Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
MoveMesh.cpp
Go to the documentation of this file.
1 #include <cstdio>
2 #include <cstdlib>
3 #include <iostream>
4 #include <iomanip>
5 #include <MultiRegions/ExpList.h>
6 #include <MultiRegions/ContField1D.h>
7 #include <MultiRegions/ContField2D.h>
8 #include <MultiRegions/ContField3D.h>
9 #include <LocalRegions/SegExp.h>
10 #include <LocalRegions/TriExp.h>
11 #include <LocalRegions/QuadExp.h>
12 #include <LibUtilities/LinearAlgebra/Lapack.hpp>
13 #include <LibUtilities/BasicConst/NektarUnivTypeDefs.hpp>
14 #include <boost/lexical_cast.hpp>
15 #include <tinyxml.h>
16 
17 
18 using namespace Nektar;
19 
20 int main(int argc, char *argv[])
21 {
22  void SetFields(SpatialDomains::MeshGraphSharedPtr &mesh,
23  SpatialDomains::BoundaryConditionsSharedPtr &boundaryConditions,
24  LibUtilities::SessionReaderSharedPtr &session,
25  Array<OneD,MultiRegions::ExpListSharedPtr> &Exp,int nvariables);
26  void OrderVertices(int nedges,SpatialDomains::MeshGraphSharedPtr graphShPt,
27  MultiRegions::ExpListSharedPtr & bndfield,
28  Array<OneD, int>& Vids, int v1,int v2 , NekDouble x_connect,
29  int & lastedge,
30  Array<OneD, NekDouble>& x, Array<OneD, NekDouble>& y);
31  void Computestreakpositions(int nvertl, MultiRegions::ExpListSharedPtr &streak,
32  Array<OneD, NekDouble> x, Array<OneD, NekDouble> y,
33  Array<OneD, NekDouble> &xold_up, Array<OneD, NekDouble> &yold_up,
34  Array<OneD, NekDouble> &xold_low, Array<OneD, NekDouble> &yold_low,
35  Array<OneD, NekDouble> &xold_c, Array<OneD, NekDouble> &yold_c,
36  Array<OneD, NekDouble> &xc, Array<OneD, NekDouble> &yc, NekDouble cr,
37  bool verts);
38  void GenerateAddPointsNewtonIt( NekDouble &xi, NekDouble &yi,NekDouble &x0, NekDouble &y0,
39  MultiRegions::ExpListSharedPtr &function, Array<OneD, NekDouble> &derfunction,
40  NekDouble cr);
41  void GenerateCurveSp(Array<OneD, NekDouble> &x_c, Array<OneD, NekDouble> &y_c,
42  Array<OneD, NekDouble> &x_lay, Array<OneD, NekDouble> &y_lay);
43  void GenerateCurve(int npoints, int npused, Array<OneD, NekDouble> &x_c,
44  Array<OneD, NekDouble> &y_c, Array<OneD, NekDouble> &curve);
45  void GenerateAddPoints(int region, SpatialDomains::MeshGraphSharedPtr &mesh,int np, int npused,
46  Array<OneD, NekDouble> &curve, MultiRegions::ExpListSharedPtr & bndfield,
47  Array<OneD, NekDouble>& outx, Array<OneD, NekDouble>& outy,
48  Array<OneD, int>&Eids);
49  void MapEdgeVertices(MultiRegions::ExpListSharedPtr field, Array<OneD, int> &V1,
50  Array<OneD, int> &V2);
51  void Findlay_eids(Array<OneD, Array<OneD, NekDouble> >& layers_y,
52  Array<OneD, int> V1, Array<OneD, int> V2,
53  Array<OneD, Array<OneD, int > >& lay_Vids,
54  Array<OneD, Array<OneD, int > >& lay_eids);
55  void MoveLayersvertically(int nlays, int nvertl, int cntlow, int cntup,
56  Array<OneD, Array<OneD, int > > lay_Vids, Array<OneD, NekDouble> xc,
57  Array<OneD, NekDouble> yc, Array<OneD, int> Down, Array<OneD, int> Up,
58  Array<OneD, NekDouble >& xnew, Array<OneD, NekDouble>& ynew,
59  Array<OneD, Array<OneD, NekDouble > >& layers_x,
60  Array<OneD, Array<OneD, NekDouble > >& layers_y);
61  void Cutrepetitions(int nedges,Array<OneD, NekDouble> inarray,
62  Array<OneD, NekDouble>& outarray);
63  void Orderfunctionx(Array<OneD, NekDouble> inarray_x,
64  Array<OneD, NekDouble> inarray_y, Array<OneD, NekDouble>& outarray_x,
65  Array<OneD, NekDouble>& outarray_y);
66  int DetermineclosePointxindex(NekDouble x,Array<OneD, NekDouble> xArray);
67  void GenerateNeighbourArrays(int index, int neighpoints,Array<OneD, NekDouble> xArray,
68  Array<OneD, NekDouble> yArray,Array<OneD, NekDouble>& Neighbour_x,
69  Array<OneD, NekDouble>& Neighbour_y);
70  NekDouble LagrangeInterpolant(NekDouble x, int npts,
71  Array<OneD,NekDouble> xpts, Array<OneD, NekDouble> funcvals);
72  void ChangeLayerspos(Array<OneD, NekDouble> & ynew,
73  Array<OneD, NekDouble> yc,
74  Array<OneD, NekDouble> Addpointsy,
75  Array<OneD, Array<OneD, NekDouble > >& layers_y,
76  Array<OneD, Array<OneD, int> >lay_Vids,
77  NekDouble delt, NekDouble delt_opp,string movelay, int npedge);
78 
79  void Replacevertices(string filename, Array<OneD, NekDouble> newx,
80  Array<OneD, NekDouble> newy,
81  Array<OneD, NekDouble> x_crit, Array<OneD, NekDouble> y_crit,
82  Array<OneD, NekDouble> & Pcurvx,
83  Array<OneD, NekDouble> & Pcurvy,
84  Array<OneD, int>&Eids, int Npoints, string s_alp,
85  Array<OneD, Array<OneD, NekDouble> > x_lay,
86  Array<OneD, Array<OneD, NekDouble> > y_lay,
87  Array<OneD, Array<OneD, int > >lay_eids, bool curv_lay);
88 
89 
90 
91  int i,j;
92 
93 //ATTEnTION !!! with argc=2 you impose that vSession refers to is argv[1]=meshfile!!!!!
94  LibUtilities::SessionReaderSharedPtr vSession
95  = LibUtilities::SessionReader::CreateInstance(2, argv);
96  //----------------------------------------------
97  NekDouble cr;
98  //set cr =0
99  cr=0;
100  //change argc from 6 to 5 allow the loading of cr to be optional
101  if(argc > 6 || argc <5 )
102  {
103  fprintf(stderr,
104  "Usage: ./MoveMesh meshfile fieldfile changefile alpha cr(optional)\n");
105  exit(1);
106  }
107  else if( argc == 6 &&
108  vSession->DefinesSolverInfo("INTERFACE")
109  && vSession->GetSolverInfo("INTERFACE")=="phase" )
110  {
111  cr = boost::lexical_cast<double>(argv[argc-1]);
112  argc=5;
113  }
114 
115  // Read in mesh from input file
116  string meshfile(argv[argc-4]);
117  SpatialDomains::MeshGraphSharedPtr graphShPt = SpatialDomains::MeshGraph::Read(vSession);//meshfile);
118  //----------------------------------------------
119 
120  // Also read and store the boundary conditions
121  SpatialDomains::BoundaryConditionsSharedPtr boundaryConditions;
122  boundaryConditions = MemoryManager<SpatialDomains::BoundaryConditions>
123  ::AllocateSharedPtr(vSession,graphShPt);
124  SpatialDomains::BoundaryConditions bcs(vSession, graphShPt);
125  //----------------------------------------------
126 
127  // Define Expansion
128  Array<OneD, MultiRegions::ExpListSharedPtr> fields;
129  int nfields;
130 cout<<"o2k"<<endl;
131  //the mesh file should have 2 component: set output fields
132  //fields has to be of the SAME dimension of the mesh (that's why there is
133  //the changefile as an input)
134  //a contfield2D is needed to extract boundary conditions!!!
135  nfields=3;
136  SetFields(graphShPt,boundaryConditions,vSession,fields,nfields);
137  //---------------------------------------------------------------
138 cout<<"ok"<<endl;
139  // store name of the file to change
140  string changefile(argv[argc-2]);
141  //----------------------------------------------
142 
143  //store the value of alpha
144  string charalp (argv[argc-1]);
145  //NekDouble alpha = boost::lexical_cast<double>(charalp);
146 cout<<"read alpha="<<charalp<<endl;
147  // Import field file.
148  string fieldfile(argv[argc-3]);
149  vector<LibUtilities::FieldDefinitionsSharedPtr> fielddef;
150  vector<vector<NekDouble> > fielddata;
151  LibUtilities::Import(fieldfile,fielddef,fielddata);
152  //----------------------------------------------
153 
154 //cout<<"dim st="<<fielddef[0]->m_fields.size()<<endl;
155  //fill a vector with the streak sol
156 /*
157  Array<OneD, MultiRegions::ExpListSharedPtr> streak;
158  SetFields(graphShPt, boundaryConditions, vSession, streak, 1);
159  for(int i = 0; i < fielddata.size(); ++i)
160  {
161  streak[0]->ExtractDataToCoeffs(fielddef[i],fielddata[i],fielddef[i]->m_fields[0]);
162  }
163  streak[0]->BwdTrans(streak[0]->GetCoeffs(), streak[0]->UpdatePhys());
164 */
165  MultiRegions::ExpListSharedPtr streak;
166  streak = MemoryManager<MultiRegions::ContField2D>
167  ::AllocateSharedPtr(vSession, graphShPt, "w",true);
168 
169  for(int i=0; i<fielddata.size(); i++)
170  {
171  streak->ExtractDataToCoeffs(fielddef[i], fielddata[i], fielddef[i]->m_fields[0], streak->UpdateCoeffs());
172  }
173  streak->BwdTrans_IterPerExp(streak->GetCoeffs(), streak->UpdatePhys());
174 
175  //------------------------------------------------
176 /*
177 static int cnt=0;
178 int nquad = streak->GetTotPoints();
179 Array<OneD, NekDouble> xs(nquad);
180 Array<OneD, NekDouble> ys(nquad);
181 streak->GetCoords(xs,ys);
182 if(
183 //abs(streak->GetPhys()[j])< 0.01
184 x[j]==1.6
185 )
186 {
187 cnt++;
188 cout<<"cnt="<<cnt<<" x="<<xs[j]<<" y="<<ys[j]<<" U="<<streak->GetPhys()[j]<<endl;
189 }
190 */
191  //----------------------------------------------
192 /*
193  // Copy data from file:fill fields with the fielddata
194  for(j = 0; j < nfields; ++j)
195  {
196  for(int i = 0; i < fielddata.size(); ++i)
197  {
198  fields[j]->ExtractDataToCoeffs(fielddef[i],fielddata[i],fielddef[i]->m_fields[j]);
199  }
200  fields[j]->BwdTrans(fields[j]->GetCoeffs(),fields[j]->UpdatePhys());
201  }
202  //----------------------------------------------
203 */
204  // determine the I regions (3 region expected)
205  //hypothesys: the layes have the same number of points
206 
207  int nIregions, lastIregion;
208  const Array<OneD, SpatialDomains::BoundaryConditionShPtr> bndConditions = fields[0]->GetBndConditions();
209  Array<OneD, int> Iregions =Array<OneD, int>(bndConditions.num_elements(),-1);
210 
211  nIregions=0;
212  int nbnd= bndConditions.num_elements();
213  for(int r=0; r<nbnd; r++)
214  {
215  if(bndConditions[r]->GetUserDefined()==SpatialDomains::eCalcBC)
216  {
217  lastIregion=r;
218  Iregions[r]=r;
219  nIregions++;
220  }
221  }
222 
223  ASSERTL0(nIregions>0,"there is any boundary region with the tag USERDEFINEDTYPE=""CalcBC"" specified");
224 cout<<"nIregions="<<nIregions<<endl;
225  //set expansion along a layers
226  Array<OneD, MultiRegions::ExpListSharedPtr> bndfieldx= fields[0]->GetBndCondExpansions();
227  //--------------------------------------------------------
228 
229  //determine the points in the lower and upper curve...
230  int nq1D= bndfieldx[lastIregion]->GetTotPoints();
231  int nedges = bndfieldx[lastIregion]->GetExpSize();
232  int nvertl = nedges +1 ;
233  Array<OneD, int> Vids_low(nvertl,-10);
234  Array<OneD, NekDouble> xold_low(nvertl);
235  Array<OneD, NekDouble> yold_low(nvertl);
236  Array<OneD, NekDouble> zi(nvertl);
237 
238  //order the ids on the lower curve lastIregion starting from the id on x=0
239  NekDouble x_connect;
240  NekDouble x0,y0,z0,xt,yt=0.0,zt=0.0;
241  int lastedge=-1;
242  int v1,v2;
243  //first point for x_connect=0(or-1.6 for the full mesh (-pi,pi) )
244  x_connect=0;
245  SpatialDomains::PointGeomSharedPtr vertex0 =
246  graphShPt->GetVertex
247  (
248  ( (bndfieldx[lastIregion]->GetExp(0)->as<LocalRegions::SegExp>())
249  ->GetGeom1D()
250  )
251  ->GetVid(0)
252  );
253  vertex0->GetCoords(x0,y0,z0);
254  if( x0 != 0.0)
255  {
256 cout<<"WARNING x0="<<x0<<endl;
257  x_connect=x0;
258  }
259  v1=0;
260  v2=1;
261  OrderVertices(nedges, graphShPt, bndfieldx[lastIregion-1],
262  Vids_low, v1, v2 , x_connect ,lastedge, xold_low,yold_low);
263  ASSERTL0(Vids_low[v2]!=-10, "Vids_low[v2] is wrong");
264  SpatialDomains::PointGeomSharedPtr vertex = graphShPt->GetVertex(Vids_low[v2]);
265 
266  //update x_connect
267 cout<<"x_conn="<<x_connect<<" yt="<<yt<<" zt="<<zt<<" vid="<<Vids_low[v2]<<endl;
268  vertex->GetCoords(x_connect,yt,zt);
269 
270  i=2;
271  while(i<nvertl)
272  {
273  v1=i;
274  OrderVertices(nedges, graphShPt, bndfieldx[lastIregion-1],
275  Vids_low, v1, v2 , x_connect, lastedge, xold_low, yold_low );
276  SpatialDomains::PointGeomSharedPtr vertex = graphShPt->GetVertex(Vids_low[v1]);
277 //cout<<"Vids low="<<Vids_low[v1]<<endl;
278  //update x_connect (lastedge is updated on the OrderVertices function)
279  vertex->GetCoords(x_connect,yt,zt);
280  i++;
281  }
282 
283  //------------------------------------------------------------------------
284 
285  //order in the same way the id of the upper curve lastIregion-1 starting from x=0:
286  Array<OneD, int> Vids_up(nvertl,-10);
287  Array<OneD,NekDouble> xold_up(nvertl);
288  Array<OneD,NekDouble> yold_up(nvertl);
289  //first point for x_connect=0 (or-1.6)
290  x_connect=0;
291  vertex0 =
292  graphShPt->GetVertex
293  (
294  ( (bndfieldx[lastIregion]->GetExp(0)->as<LocalRegions::SegExp>())
295  ->GetGeom1D()
296  )
297  ->GetVid(0)
298  );
299  vertex0->GetCoords(x0,y0,z0);
300  if( x0 != 0.0)
301  {
302 cout<<"WARNING x0="<<x0<<endl;
303  x_connect=x0;
304  }
305  lastedge=-1;
306 
307  v1=0;
308  v2=1;
309  OrderVertices(nedges, graphShPt, bndfieldx[lastIregion-2 ],
310  Vids_up, v1, v2 , x_connect ,lastedge, xold_up, yold_up);
311  SpatialDomains::PointGeomSharedPtr vertexU = graphShPt->GetVertex(Vids_up[v2]);
312 //cout<<"VIdup="<<Vids_up[v2]<<endl;
313  //update x_connect
314  vertexU->GetCoords(x_connect,yt,zt);
315 
316  i=2;
317  while(i<nvertl)
318  {
319  v1=i;
320  OrderVertices(nedges, graphShPt, bndfieldx[lastIregion-2],
321  Vids_up, v1, v2 , x_connect, lastedge, xold_up, yold_up );
322  SpatialDomains::PointGeomSharedPtr vertex = graphShPt->GetVertex(Vids_up[v1]);
323 //cout<<"VIdup="<<Vids_up[v1]<<endl;
324  //update x_connect (lastedge is updated on the OrderVertices function)
325  vertex->GetCoords(x_connect,yt,zt);
326  i++;
327  }
328 
329  //-----------------------------------------------------------------------------------
330 
331 
332  //order in the same way the id of the layer curve lastIregion starting from x=0:
333  Array<OneD, int> Vids_c(nvertl,-10);
334  Array<OneD,NekDouble> xold_c(nvertl);
335  Array<OneD,NekDouble> yold_c(nvertl);
336  //first point for x_connect=0(or-1.6)
337  x_connect=0;
338  vertex0 =
339  graphShPt->GetVertex
340  (
341  ( (bndfieldx[lastIregion]->GetExp(0)->as<LocalRegions::SegExp>())
342  ->GetGeom1D()
343  )
344  ->GetVid(0)
345  );
346  vertex0->GetCoords(x0,y0,z0);
347  if( x0 != 0.0)
348  {
349 cout<<"WARNING x0="<<x0<<endl;
350  x_connect=x0;
351  }
352  lastedge=-1;
353 
354  v1=0;
355  v2=1;
356 
357  OrderVertices(nedges, graphShPt, bndfieldx[lastIregion],
358  Vids_c, v1, v2 , x_connect ,lastedge, xold_c, yold_c);
359  SpatialDomains::PointGeomSharedPtr vertexc = graphShPt->GetVertex(Vids_c[v2]);
360 
361  //update x_connect
362  vertexc->GetCoords(x_connect,yt,zt);
363 
364  i=2;
365  while(i<nvertl)
366  {
367  v1=i;
368  OrderVertices(nedges, graphShPt, bndfieldx[lastIregion],
369  Vids_c, v1, v2 , x_connect, lastedge, xold_c, yold_c );
370  SpatialDomains::PointGeomSharedPtr vertex = graphShPt->GetVertex(Vids_c[v1]);
371 //cout<<"Vids cl="<<Vids_low[v1]<<endl;
372  //update x_connect (lastedge is updated on the OrderVertices function)
373  vertex->GetCoords(x_connect,yt,zt);
374  i++;
375  }
376 
377  //calculate the distances between the layer and the upper/lower curve
378 
379  Array<OneD, NekDouble> Deltaup (nvertl, -200);
380  Array<OneD, NekDouble> Deltalow (nvertl, -200);
381 
382  for(int r=0; r<nvertl; r++)
383  {
384  //Always positive!!!
385 //cout<<"i="<<r<<" yup="<<yold_up[r]<<" yc="<<yold_c[r]<<" ylow="<<yold_low[r]<<endl;
386  Deltaup[r] = yold_up[r] - yold_c[r];
387  Deltalow[r] = yold_c[r] - yold_low[r];
388  ASSERTL0(Deltaup[r]>0, "distance between upper and layer curve is not positive");
389  ASSERTL0(Deltalow[r]>0, "distance between lower and layer curve is not positive");
390  }
391  //------------------------------------------------------------------------
392 
393  //-----------------------------------------------------------------------------------
394 
395 
396  //fieds to force continuity:
397 /*
398  MultiRegions::ContField1DSharedPtr Lay_x;
399  MultiRegions::ContField1DSharedPtr Lay_y;
400  //initialie fields
401  const SpatialDomains::BoundaryRegionCollection &bregions = bcs.GetBoundaryRegions();
402  MultiRegions::ExpList1DSharedPtr xtmp;
403  MultiRegions::ExpList1DSharedPtr ytmp;
404  xtmp = MemoryManager<MultiRegions::ExpList1D>
405  ::AllocateSharedPtr(*(bregions[lastIregion]), graphShPt, true);
406  ytmp = MemoryManager<MultiRegions::ExpList1D>
407  ::AllocateSharedPtr(*(bregions[lastIregion]), graphShPt, true);
408 
409  Lay_x = MemoryManager<MultiRegions::ContField1D>
410  ::AllocateSharedPtr(vSession, *xtmp);
411 
412  Lay_y = MemoryManager<MultiRegions::ContField1D>
413  ::AllocateSharedPtr(vSession, *ytmp);
414 */
415  //--------------------------------------
416  /*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/
417  //generate additional points using the Newton iteration
418  //determine the xposition for every edge (in the middle even if it
419  // is not necessary
420  //PARAMETER which determines the number of points per edge @todo put as an input
421  int npedge;
422  if(vSession->DefinesParameter("npedge"))
423  {
424  npedge = (int)vSession->GetParameter("npedge");
425  }
426  else
427  {
428  npedge = 5;//default value
429  }
430  /*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/
431  //find the points where u=0 and determine the sign of the shift and the delta
432  int nq= streak->GetTotPoints();
433  Array<OneD, NekDouble> x(nq);
434  Array<OneD,NekDouble> y(nq);
435  fields[0]->GetCoords(x,y);
436  Array<OneD, NekDouble> x_c(nvertl);
437  Array<OneD,NekDouble> y_c(nvertl,-200);
438  Array<OneD, NekDouble> tmp_w (nvertl, 200);
439  Array<OneD, int> Sign (nvertl,1);
440  Array<OneD, NekDouble> Delta_c(nvertl,-200);
441 
442 
443  Computestreakpositions(nvertl, streak, x,y,xold_up, yold_up,
444  xold_low, yold_low, xold_c, yold_c, x_c, y_c,cr,true);
445  // if the curve is low the old layer point, it has to shift down
446  NekDouble shift;
447  for(int q=0; q<nvertl; q++)
448  {
449  if(y_c[q] < yold_c[q])
450  {
451  Sign[q] = -1;
452  }
453  //calculate delta
454  Delta_c[q] = abs(yold_c[q]-y_c[q]);
455  //fill layer arrays:
456 /*
457  if(q==0 || q == nvertl-1)
458  {
459  Lay_x->UpdatePhys()[q]=x_c[q];
460  Lay_y->UpdatePhys()[q]=y_c[q];
461  }
462  else if(q<nvertl && q!=nvertl-1)
463  {
464  Lay_x->UpdatePhys()[ q*(npedge-1) ]=x_c[q];
465  Lay_y->UpdatePhys()[ q*(npedge-1) ]=y_c[q];
466  Lay_x->UpdatePhys()[ q*(npedge-1) +1 ]=x_c[q];
467  Lay_y->UpdatePhys()[ q*(npedge-1) +1 ]=y_c[q];
468  }
469 */
470  //check the shifting of the layer:
471  shift+= Delta_c[q];
472 cout<<x_c[q]<<" "<<y_c[q]<<endl;
473  }
474 //cout<<"shift="<<shift<<endl;
475  if(shift<0.001)
476  {
477  cout<<"Warning: the critical layer is stationary"<<endl;
478  }
479  //------------------------------------------------------------------
480 
481 
482 
483 
484  //additional points arrays
485  Array<OneD, NekDouble> Cpointsx (nedges);
486  Array<OneD, NekDouble> Cpointsy (nedges, 0.0);
487  Array<OneD, int> Eids (nedges);
488 
489 /*
490  int ncurves=1;
491 
492  int ntotcurvecoeffs=0;
493  //array storing ALL the curvecoeffs
494  //hypothesis: totcoeffs< 2*nvertl; ncurves< nedges;
495  Array<OneD, NekDouble> totcurvecoeffs (2*nvertl,0.0);
496  Array<OneD, int> ncurvepoints (nedges);
497 
498  //calculate the approximated curve coeffs
499  //since the number of degree of freedom is too large
500  // we split the layer in as many curves as many changing of the rate
501  // \Delta( (\Delta y /\Delta x) ) (second derivative)
502  //for each change of rate, calculate one curve and the additional point for
503  // the curved tag of the session file
504 
505  //change of curve second derivative sign quantities
506  NekDouble der,tmp;
507  tmp = (y_c[1]-y_c[0])/(x_c[1]-x_c[0]) ;
508  der = (y_c[2]-y_c[1])/(x_c[2]-x_c[1]) ;
509  int sign2der=1;
510  int Npoints=0;
511  int Npused=0;
512  int cnt=0;
513  if(der < tmp)
514  {
515  sign2der =-1;
516  }
517 //cout<<"sign="<<sign2der<<" tmp="<<tmp<<" der="<<der<<endl;
518  for(int r=2; r< nedges; r++)
519  {
520  tmp = (y_c[r+1]-y_c[r])/(x_c[r+1]-x_c[r]);
521 
522  if(sign2der==1)
523  {
524  if(tmp < der)
525  {
526  cnt++;
527  //Array<OneD, NekDouble> curvecoeffs (nvertl);
528 //cout<<"change 1"<<endl;
529  Npoints = r-(ncurves-1)-Npoints;
530  ncurvepoints[ncurves-1] = Npoints;
531  ncurves++;
532  sign2der = -1*sign2der;
533 
534 //cout<<"Npoints="<<Npoints<<" r="<<r<<endl;
535  if(Npused==0)
536  {
537  ntotcurvecoeffs += Npoints;
538  GenerateCurve(Npoints, Npused, x_c, y_c, totcurvecoeffs);
539  GenerateAddPoints(lastIregion, graphShPt, Npoints, Npused, totcurvecoeffs,
540  bndfieldx[lastIregion], Cpointsx, Cpointsy, Eids);
541 
542  }
543  else
544  {
545  ntotcurvecoeffs += Npoints+1;
546  GenerateCurve(Npoints+1, Npused-1, x_c, y_c, totcurvecoeffs);
547  GenerateAddPoints(lastIregion, graphShPt, Npoints+1, Npused-1, totcurvecoeffs,
548  bndfieldx[lastIregion], Cpointsx, Cpointsy, Eids);
549  }
550  Npused +=Npoints;
551  }
552  }
553  else
554  {
555  if(tmp > der)
556  {
557  cnt++;
558 //cout<<"change -1"<<endl;
559  Npoints = r-(ncurves-1)-Npoints;
560  ncurvepoints[ncurves-1] = Npoints;
561  ncurves++;
562  sign2der = -1*sign2der;
563 //cout<<"npoints="<<Npoints<<" r="<<r<<endl;
564  if(Npused==0)
565  {
566  ntotcurvecoeffs += Npoints;
567  GenerateCurve(Npoints, Npused, x_c, y_c, totcurvecoeffs);
568  GenerateAddPoints(lastIregion, graphShPt, Npoints, Npused, totcurvecoeffs,
569  bndfieldx[lastIregion], Cpointsx, Cpointsy, Eids);
570  }
571  else
572  {
573  ntotcurvecoeffs += Npoints+1;
574  GenerateCurve(Npoints+1, Npused-1, x_c, y_c, totcurvecoeffs);
575  GenerateAddPoints(lastIregion, graphShPt, Npoints+1, Npused-1, totcurvecoeffs,
576  bndfieldx[lastIregion], Cpointsx, Cpointsy, Eids);
577  }
578  Npused +=Npoints;
579  }
580  }
581 //cout<<"sign="<<sign2der<<" tmp="<<tmp<<" der="<<der<<endl;
582 
583  der =tmp;
584  }
585 
586  // generate last curve adding a point to fill the gap
587  //between the previus calculated curve
588  int Nplast = nvertl - Npused;
589  if(Npused==0)
590  {
591  ntotcurvecoeffs += Nplast;
592  ncurvepoints[ncurves-1] = Nplast;
593  GenerateCurve(Npoints, Npused, x_c, y_c, totcurvecoeffs);
594  GenerateAddPoints(lastIregion, graphShPt, Npoints, Npused, totcurvecoeffs,
595  bndfieldx[lastIregion], Cpointsx, Cpointsy, Eids);
596  }
597  else
598  {
599  ntotcurvecoeffs += Nplast+1;
600 //cout<<"last curve npoints="<<Nplast<<" Npused="<<Npused<<endl;
601 //cout<<"ntotcurvecoeffs="<<ntotcurvecoeffs<<endl;
602  GenerateCurve(Nplast+1, Npused-1, x_c, y_c, totcurvecoeffs);
603  GenerateAddPoints(lastIregion, graphShPt, Nplast+1, Npused-1, totcurvecoeffs,
604  bndfieldx[lastIregion], Cpointsx, Cpointsy, Eids);
605  }
606 //cout<<" final curves coeffs values"<<endl;
607 //cout<<"ntotcurvecoeffs="<<ntotcurvecoeffs<<endl;
608  ASSERTL0( ntotcurvecoeffs <= 2*nvertl, "numer of curves coeffs exceeds 2*nvertices");
609  for(int u=0; u <2*nvertl; u++)
610  {
611 //cout<<"coeffs = "<<totcurvecoeffs[u]<<endl;
612  }
613 
614 */
615 
616  //generate additional points using the Newton iteration
617  //determine the xposition for every edge (in the middle even if it
618  // is not necessary
619  //PARAMETER which determines the number of points @todo put as an input
620  //int npedge=5;
621  //additional points arrays
622 /*
623  //NB Double array LEAK!!!!
624  Array<OneD, Array<OneD, NekDouble> >Addpointsx;
625  Array<OneD, Array<OneD, NekDouble> >Addpointsy;
626  Addpointsx = Array<OneD, Array<OneD, NekDouble> > (nedges);
627  Addpointsy = Array<OneD, Array<OneD, NekDouble> > (nedges);
628  for(int q=0; q< npedge-2; q++)
629  {
630  // npedge-2 num points per edge
631  Addpointsx[q] = Array<OneD, NekDouble> (npedge-2, 0.0);
632  Addpointsy[q] = Array<OneD, NekDouble> (npedge-2, 0.0);
633  }
634 
635 */
636 
637  //Array<OneD, NekDouble> Ycoords (nedges*npedge-2, 0.0);
638 
639 
640  Array<OneD, NekDouble> Addpointsx (nedges*(npedge-2), 0.0);
641  Array<OneD, NekDouble> Addpointsy (nedges*(npedge-2), 0.0);
642  //calculate the dU_dy
643  Array<OneD, NekDouble> dU(streak->GetTotPoints());
644  streak->PhysDeriv(MultiRegions::eY, streak->GetPhys(), dU);
645 
646 
647  LocalRegions::SegExpSharedPtr bndSegExp;
648 
649  int Eid,id1,id2;
650  NekDouble x1,y1,z1;
651  NekDouble x2,y2,z2;
652  SpatialDomains::PointGeomSharedPtr vertex1;
653  SpatialDomains::PointGeomSharedPtr vertex2;
654  for(int r=0; r<nedges; r++)
655  {
656 
657  bndSegExp = bndfieldx[lastIregion]->GetExp(r)
658  ->as<LocalRegions::SegExp>();
659  Eid = (bndSegExp->GetGeom1D())->GetEid();
660  id1 = (bndSegExp->GetGeom1D())->GetVid(0);
661  id2 = (bndSegExp->GetGeom1D())->GetVid(1);
662  vertex1 = graphShPt->GetVertex(id1);
663  vertex2 = graphShPt->GetVertex(id2);
664  vertex1->GetCoords(x1,y1,z1);
665  vertex2->GetCoords(x2,y2,z2);
666 //cout<<"edge="<<r<<" x1="<<x1<<" x2="<<x2<<endl;
667 //cout<<"edge="<<r<<" y1="<<y1<<" y2="<<y2<<endl;
668 cout<<"edge="<<r<<" x1="<<x1<<" y1="<<y1<<" x2="<<x2<<" y2="<<y2<<endl;
669  if(x2>x1)
670  {
671  Cpointsx[r] = x1 +(x2-x1)/2;
672 //cout<<"edge="<<r<<" x1="<<x1<<" x2="<<x2<<" Cx="<<Cpointsx[r]<<endl;
673 //cout<<"edge="<<r<<" x1="<<x1<<" y1="<<y1<<" x2="<<x2<<" y2="<<y2<<endl;
674  if( Cpointsx[r]>x2 || Cpointsx[r]< x1)
675  {
676  Cpointsx[r] = -Cpointsx[r];
677  }
678  for(int w=0; w< npedge-2; w++)
679  {
680 
681  Addpointsx[r*(npedge-2) +w] = x1 +((x2-x1)/(npedge - 1))*(w+1);
682  if( Addpointsx[r*(npedge-2) +w] > x2 || Addpointsx[r*(npedge-2) +w] < x1)
683  {
684  Addpointsx[r*(npedge-2) +w] = -Addpointsx[r*(npedge-2) +w];
685  }
686  //initial guess along the line defined by the NEW verts y_c
687  Addpointsy[r*(npedge-2) +w] = y_c[r] + ((y_c[r+1]-y_c[r])/(x_c[r+1]-x_c[r]))*(Addpointsx[r*(npedge-2) +w]-x1);
688  //Addpointsy[r*(npedge-2) +w] = y1 + ((y2-y1)/(x2-x1))*(Addpointsx[r*(npedge-2) +w]-x1);
689  GenerateAddPointsNewtonIt( Addpointsx[r*(npedge-2) +w], Addpointsy[r*(npedge-2) +w],
690  Addpointsx[r*(npedge-2) +w], Addpointsy[r*(npedge-2) +w], streak, dU,cr);
691 
692  // Lay_x->UpdatePhys()[r*npedge +1 +w]= Addpointsx[r*(npedge-2) +w];
693  // Lay_y->UpdatePhys()[r*npedge +1 +w]= Addpointsy[r*(npedge-2) +w];
694  }
695  //Lay_x->UpdatePhys()[r*npedge +0] =x1;
696  //Lay_y->UpdatePhys()[r*npedge +0] =y1;
697  //Lay_x->UpdatePhys()[r*npedge +npedge-1] =x2;
698  //Lay_y->UpdatePhys()[r*npedge +npedge-1] =y2;
699 
700  }
701  else if(x1>x2)
702  {
703  Cpointsx[r] = x2+ (x1-x2)/2;
704 //cout<<"edge="<<r<<" y1="<<y1<<" y2="<<y2<<endl;
705  if( Cpointsx[r] > x1 || Cpointsx[r] < x2)
706  {
707  Cpointsx[r] = -Cpointsx[r];
708  }
709  for(int w=0; w< npedge-2; w++)
710  {
711  Addpointsx[r*(npedge-2) +w] = x2 +((x1-x2)/(npedge - 1))*(w+1);
712  if( Addpointsx[r*(npedge-2) +w] > x1 || Addpointsx[r*(npedge-2) +w] < x2)
713  {
714  Addpointsx[r*(npedge-2) +w] = -Addpointsx[r*(npedge-2) +w];
715  }
716 
717  //initial guess along the line defined by the NEW verts y_c
718  Addpointsy[r*(npedge-2) +w] = y_c[r+1] + ((y_c[r]-y_c[r+1])/(x_c[r]-x_c[r+1]))*(Addpointsx[r*(npedge-2) +w]-x2);
719 
720  //Addpointsy[r*(npedge-2) +w] = y2 + ((y1-y2)/(x1-x2))*(Addpointsx[r*(npedge-2) +w]-x2);
721  GenerateAddPointsNewtonIt( Addpointsx[r*(npedge-2) +w], Addpointsy[r*(npedge-2) +w],
722  Addpointsx[r*(npedge-2) +w], Addpointsy[r*(npedge-2) +w], streak, dU,cr);
723  // Lay_x->UpdatePhys()[r*npedge +1]= Addpointsx[r*(npedge-2) +w];
724  // Lay_y->UpdatePhys()[r*npedge +1]= Addpointsy[r*(npedge-2) +w];
725  }
726  //Lay_x->UpdatePhys()[r*npedge +0] =x2;
727  //Lay_y->UpdatePhys()[r*npedge +0] =y2;
728  //Lay_x->UpdatePhys()[r*npedge +npedge-1] =x1;
729  //Lay_y->UpdatePhys()[r*npedge +npedge-1] =y1;
730  }
731  else
732  {
733  ASSERTL0(false, "point not generated");
734  }
735 //cout<<"calculate cpoints coords"<<endl;
736  //Cpointsy[r] = y1 + (y2-y1)/2;
737 //cout<<"central point:"<<endl;
738  //GenerateAddPointsNewtonIt( Cpointsx[r], Cpointsy[r],Cpointsx[r], Cpointsy[r],
739  // streak, dU,cr);
740  //NekDouble diff = Cpointsy[r]-Addpointsy[r*(npedge-2)];
741 //cout<<"diff="<<diff<<endl;
742  Eids[r] = Eid;
743 
744  }
745  //-------------------------------------------------------------
746 
747 
748 
749  //fill the xPhys,yPhys array( may necessary after)
750  Array<OneD, NekDouble> xcPhys (nedges*npedge, 0.0);
751  Array<OneD, NekDouble> ycPhys (nedges*npedge, 0.0);
752 
753  for(int a=0; a<nedges; a++)
754  {
755  //v1
756  xcPhys[a*npedge+0] = x_c[a];
757  ycPhys[a*npedge+0] = y_c[a];
758  //v2
759  xcPhys[a*npedge+npedge-1] = x_c[a+1];
760  ycPhys[a*npedge+npedge-1] = y_c[a+1];
761 
762  for(int b=0; b<npedge-2; b++)
763  {
764  xcPhys[a*npedge +b+1] = Addpointsx[a*(npedge-2)+b];
765  ycPhys[a*npedge +b+1] = Addpointsy[a*(npedge-2)+b];
766  }
767  }
768 
769  //force the continuity of the layer
770 
771 /*
772 cout<<"force"<<endl;
773  int Nregcoeffs = Lay_x->GetNcoeffs();
774  Array<OneD, NekDouble> coeffs (Nregcoeffs);
775  //Lay_x->FwdTrans(Lay_x->GetPhys(), coeffs);
776  //Lay_x->BwdTrans(coeffs, Lay_x->UpdatePhys());
777  //Array<OneD, NekDouble> dery(npedge*nedges);
778  //Lay_y->PhysDeriv(MultiRegions::eX, Lay_y->GetPhys(), dery);
779 
780 
781  Vmath::Zero(Nregcoeffs,coeffs,1);
782 
783  //Lay_y->FwdTrans(Lay_y->GetPhys(), coeffs);
784  //Lay_y->BwdTrans(coeffs, Lay_y->UpdatePhys());
785  for(int r=0; r<nedges; r++)
786  {
787 
788  for(int w=0; w< npedge-2; w++)
789  {
790  // Addpointsx[r*(npedge-2) +w] = Lay_x->GetPhys()[r*npedge +1 +w];
791  // Addpointsy[r*(npedge-2) +w] = Lay_y->GetPhys()[r*npedge +1 +w];
792  }
793  }
794 */
795  Array<OneD, NekDouble> curve(nvertl);
796 cout<<"gen curve"<<endl;
797  //GenerateCurve(nvertl,0, x_c, y_c, curve);
798 
799 
800  //-------------------------------------------------
801 
802 
803 
804 /*
805  //generate the closest curve to the critical layer positions taking
806  // npedge-2 points per edge
807  Array<OneD, NekDouble> Crlay_pointsx(nedges*(npedge-2)+nedges+1);
808  Array<OneD, NekDouble> Crlay_pointsy(nedges*(npedge-2)+nedges+1);
809  int cnt=0;
810  int cnt1=0;
811  int cnt2=0;
812  //HYPOTHESIS: Cpoints ALREADY generated
813 //cout<<"Cr x--y"<<endl;
814  for(int u=0; u<Crlay_pointsx.num_elements(); u++)
815  {
816  if( u== cnt)// u pari
817  {
818  Crlay_pointsx[u] = x_c[cnt1];
819  Crlay_pointsy[u] = y_c[cnt1];
820  cnt =cnt +npedge-1;
821  cnt1++;
822  }
823  else//u dispari
824  {
825  //Crlay_pointsx[u] = Cpointsx[cnt1];
826  //Crlay_pointsy[u] = Cpointsy[cnt1];
827  Crlay_pointsx[u] = Addpointsx[cnt2];
828  Crlay_pointsy[u] = Addpointsy[cnt2];
829  cnt2++;
830  }
831 
832 //cout<<u<<" "<<Crlay_pointsx[u]<<" "<<Crlay_pointsy[u]<<endl;
833  }
834 */
835 //cout<<"num px="<<Crlay_pointsx.num_elements()<<endl;
836 /*
837  for(int r=0; r< Crlay_pointsx.num_elements(); r++)
838  {
839 cout<<r<<" "<<Crlay_pointsx[r]<<" "<<Crlay_pointsy[r]
840 <<" "<<sqrt((Crlay_pointsx[r+1]- Crlay_pointsx[r])
841 *(Crlay_pointsx[r+1]- Crlay_pointsx[r]) +
842 (Crlay_pointsy[r+1]- Crlay_pointsy[r])
843 *(Crlay_pointsy[r+1]- Crlay_pointsy[r])) <<endl;
844  }
845 
846 */
847 
848 
849  Array<OneD, NekDouble> curve_coeffs (nedges*(npedge-2)+nedges+1);
850  //generate the curve
851  //NB:: for a high number of points (>20) it does not work for all
852  // points
853 //cout<<"CALL generate curve"<<nedges*(npedge-2)+nedges+1<<endl;
854  //GenerateCurve(nedges*(npedge-2)+nedges+1, 0, Crlay_pointsx, Crlay_pointsy, curve_coeffs);
855  //calculate the coords of the additional points as y=poly(x)
856  //HYPOTHESES:
857  //x coords ALREADY CALCULATED
858  //npedge is the number of points per edge
859  //Eid[s] ALREADY GENERATED
860 /*
861  int polorder;
862  double xl;
863  //put zero Addpointsy!!!!!
864  Vmath::Zero(nedges*(npedge-2), Addpointsy,1);
865  for(int e=0; e< nedges*(npedge-2); e++)
866  {
867  xl = Addpointsx[e];
868  polorder = curve_coeffs.num_elements()-1;
869 
870  for(int g= 0; g< curve_coeffs.num_elements(); g++)
871  {
872  Addpointsy[e] += curve_coeffs[g]*(std::pow( xl , polorder));
873 //cout<<g<<" coeff="<<curve_coeffs[g]<<" x^exp="<<(std::pow( xl , polorder))<<" exp="<<polorder<<endl;
874  ASSERTL0(polorder >=0, " polynomial with one negative exponent");
875  polorder--;
876  }
877 cout<<xl<<" "<<Addpointsy[e]<<" "<<Crlay_pointsx[e]<<" "
878 <<Crlay_pointsy[e]<<" "
879 <<curve_coeffs[e]<<endl;
880 
881  }
882 cout<<"LAST coeff==y(x==0)="<<curve_coeffs[nedges*(npedge-2)+nedges] <<endl;
883 */
884  //------------------------------------------------------------
885 
886  // determine the number of layers
887  int nlays=0;
888 
889 
890  int cnt=0;
891  //find coords and Vids of the layers vertices
892  int nVertTot = graphShPt->GetNvertices();
893  //arrays to check with the new values
894  Array<OneD, NekDouble> xold(nVertTot);
895  Array<OneD, NekDouble> yold(nVertTot);
896  cnt=0;
897  for(int n=0; n<nVertTot; n++)
898  {
899  SpatialDomains::PointGeomSharedPtr vertex = graphShPt->GetVertex(n);
900  NekDouble x,y,z;
901  vertex->GetCoords(x,y,z);
902  xold[n] =x;
903  yold[n] =y;
904  //tollerance 5e-2 to get the points
905  if(x<= xold_c[0] +5e-2 && x >= xold_c[0] -5e-2
906  && y<= yold_up[0] && y>=yold_low[0]
907  && y!= yold_c[0])
908  {
909  cnt++;
910 //cout<<"vert index="<<m<<" vertID="<<n<<endl;
911  }
912  }
913 
914  nlays =cnt;
915 cout<<"nlays="<<nlays<<endl;
916  Array<OneD, Array<OneD, NekDouble> > layers_y(nlays);
917  Array<OneD, Array<OneD, NekDouble> > layers_x(nlays);
918  Array<OneD, Array<OneD, int > >lay_Vids(nlays);
919  Array<OneD, Array<OneD, int > >lay_eids(nlays);
920  //initialise layers_y,lay_eids
921  for(int g=0; g<nlays; g++)
922  {
923  layers_y[g]= Array<OneD, NekDouble> ( (nvertl-1)*npedge );
924  lay_Vids[g]= Array<OneD, int> (nvertl);
925  lay_eids[g]= Array<OneD, int> (nvertl-1);
926  }
927 
928  //------------------------------------------------------------
929 
930  //calculate the new cordinates of the vertices
931  Array<OneD, NekDouble> xnew(nVertTot);
932  Array<OneD, NekDouble> ynew(nVertTot,-20);
933  Array<OneD, int> Up(nvertl);//Vids lay Up
934  Array<OneD, int> Down(nvertl);//Vids lay Down
935  /////////////////////////////////////////////////////////////////////////////
936  ////////////////////////////////////////////////////////////////////////////
937  //££££££££££££££££££££££££££££££££££££££££££££££££££££££££££££££££££££
938  //$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
939  //@todo set Delta0 from session file
940  NekDouble Delta0;
941  if(vSession->DefinesParameter("Delta"))
942  {
943  Delta0 = vSession->GetParameter("Delta");
944  }
945  else
946  {
947  Delta0 = 0.1;//default value
948  }
949 
950  //$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
951  //££££££££££££££££££££££££££££££££££££££££££££££££££££££££££££££££££££
952  ////////////////////////////////////////////////////////////////////////////
953  /////////////////////////////////////////////////////////////////////////////
954  int cntup=0;
955  int cntlow=0;
956  Array<OneD, int> Acntlay(nvertl,0);
957  //Vids needed only if a layers has to be moved
958  NekDouble bleft=-10;
959  NekDouble tright = 10;
960  NekDouble bright = -10;
961  NekDouble tleft = 10;
962  cnt=0;
963  int bottomleft, topright,bottomright,topleft;
964  for(int i=0; i<nVertTot; i++)
965  {
966  bool mvpoint =false;
967  SpatialDomains::PointGeomSharedPtr vertex = graphShPt->GetVertex(i);
968  NekDouble x,y,z;
969  vertex->GetCoords(x,y,z);
970  //x coord doesn't change
971  xnew[i]=x;
972 //cout<<"x="<<x<<" y="<<y<<endl;
973  //bottom, left (x=0, y<ydown)
974  if(x==0 && y< yold_low[0]
975  && y> bleft)
976  {
977  bleft = y;
978  bottomleft = i;
979  }
980  //top, right
981  if(x== xold_c[nvertl-1] && y> yold_up[nvertl-1]
982  && y<tright)
983  {
984  tright = y;
985  topright = i;
986  }
987  //bottom, right]
988  if(x==xold_c[nvertl-1] && y<yold_low[nvertl-1]
989  && y> bright)
990  {
991  bright = y;
992  bottomright = i;
993  }
994  //top, left
995  if(x== 0 && y> yold_up[0]
996  && y<tleft)
997  {
998  tleft = y;
999  topleft = i;
1000  }
1001  //find the corresponding yold_l and deltaold
1002  for(int j=0; j<nvertl; j++)
1003  {
1004  if((xold_up[j]==x)&&(yold_up[j]==y))
1005  {
1006  //ratio = (1-y)*(1+y)/( (1-yold_c[j])*(1+yold_c[j]) );
1007  //ynew[i] = y + Sign[j]*Delta_c[j]*ratio;
1008  ynew[i] = y_c[j] +Delta0;
1009  mvpoint=true;
1010  Up[j] = i;
1011  }
1012  if((xold_low[j]==x)&&(yold_low[j]==y))
1013  {
1014  //ratio = (1-y)*(1+y)/( (1-yold_c[j])*(1+yold_c[j]) );
1015  //ynew[i] = y + Sign[j]*Delta_c[j]*ratio;
1016  ynew[i] = y_c[j] -Delta0;
1017  mvpoint=true;
1018  Down[j] = i;
1019  }
1020  if((xold_c[j]==x)&&(yold_c[j]==y))
1021  {
1022  ynew[i] = y_c[j];
1023  mvpoint=true;
1024  }
1025  }
1026  if(mvpoint==false)
1027  {
1028  //determine the closer xold_up
1029  NekDouble diff=1000;
1030  int qp_closer;
1031  for(int k=0; k<nvertl; k++)
1032  {
1033  if(abs(x-xold_up[k]) < diff)
1034  {
1035  diff = abs(x-xold_up[k]);
1036  qp_closer=k;
1037  }
1038  }
1039  if( y>yold_up[qp_closer] && y< 1)
1040  {
1041 
1042  //ratio = (1-y)*(1+y)/( (1-yold_c[qp_closer])*(1+yold_c[qp_closer]) );
1043  //ynew[i] = y + Sign[qp_closer]*Delta_c[qp_closer]*ratio;
1044 
1045  //ratio = (1-y)*(1+y)/( (1-yold_up[qp_closer])*(1+yold_up[qp_closer]) );
1046  //distance prop to layerup
1047  ynew[i] = y_c[qp_closer] +(y-yold_c[qp_closer])*
1048  (1-y_c[qp_closer])/(1-yold_c[qp_closer]);
1049 //cout<<"upper zone y="<<y<<" ratio="<<ratio<<endl;
1050 
1051 
1052 
1053  }
1054  else if(y<yold_low[qp_closer] && y> -1)
1055  {
1056 
1057  //ratio = (1-y)*(1+y)/( (1-yold_c[qp_closer])*(1+yold_c[qp_closer]) );
1058  //ynew[i] = y + Sign[qp_closer]*Delta_c[qp_closer]*ratio;
1059 
1060  //ratio = (1-y)*(1+y)/( (1-yold_low[qp_closer])*(1+yold_low[qp_closer]) );
1061  //distance prop to layerlow
1062  ynew[i] = y_c[qp_closer] + (y-yold_c[qp_closer] )*
1063  (-1-y_c[qp_closer])/(-1-yold_c[qp_closer]);
1064 
1065  }
1066 
1067  else if ( y>yold_c[qp_closer] && y < yold_up[qp_closer])
1068  {
1069  if(x==0){ cntup++; }
1070  //ratio = (1-y)*(1+y)/( (1-yold_c[qp_closer])*(1+yold_c[qp_closer]) );
1071  //ynew[i] = y + Sign[qp_closer]*Delta_c[qp_closer]*ratio;
1072 
1073  }
1074  else if (y<yold_c[qp_closer] && y > yold_low[qp_closer])
1075  {
1076  if(x==0){ cntlow++; }
1077  // ratio = (1-y)*(1+y)/( (1-yold_c[qp_closer])*(1+yold_c[qp_closer]) );
1078  // ynew[i] = y + Sign[qp_closer]*Delta_c[qp_closer]*ratio;
1079 
1080  }
1081  else if( y==1 || y==-1)//bcs don't move
1082  {
1083  ynew[i] =y;
1084  //cout<<"point x="<<xnew[i]<<" y="<<y<<" closer x="<<xold_up[qp_closer]<<endl;
1085  }
1086 
1087  //internal layers are not moved yet so...
1088  if( (ynew[i]>1 || ynew[i]<-1)
1089  && ( y>yold_up[qp_closer] || y<yold_low[qp_closer]) )
1090  {
1091  cout<<"point x="<<xnew[i]<<" y="<<y<<" closer x="<<xold_up[qp_closer]<<" ynew="<<ynew[i]<<endl;
1092  ASSERTL0(false, "shifting out of range");
1093  }
1094 
1095  }
1096 //cout<<i<<" "<<xnew[i]<<" "<<ynew[i]<<endl;
1097  //count number of layers nlays(again) initialise vectors
1098  //(crit lay not incl but uplay and down lay included)
1099  //nb tollerance 10-5 to get the points
1100  for(int m=0; m<nvertl; m++)
1101  {
1102 
1103  if( x>= x_c[m] -5e-2 && x<= x_c[m] +5e-2
1104  && y<= yold_up[m] && y>=yold_low[m]
1105  && y!= yold_c[m])
1106  {
1107 //cout<<"LAY x"<<x<<" y="<<y<<endl;
1108 //cout<<"yup="<<yold_up[4]<<" ylow="<<yold_low[4]<<endl;
1109 
1110  //need to fill layers_y to order it(through Findlay_eids)
1111  layers_y[Acntlay[m]][m]= y;
1112  lay_Vids[Acntlay[m]][m]= i;
1113 
1114  Acntlay[m]++;
1115  }
1116 
1117  }
1118 
1119  }
1120 cout<<"cntlow="<<cntlow<<endl;
1121  ASSERTL0(Acntlay[4]==nlays, "something wrong with the number of layers");
1122  //order layers coords:
1123  //V1[eid],V2[eid] vertices associate with the edge Id=eid
1124  Array<OneD, int> V1;
1125  Array<OneD, int> V2;
1126  MapEdgeVertices(streak,V1,V2);
1127  //find eids which belong to each layer(essentially order the
1128  // lay_Vids[n][m] on m layers)
1129  Findlay_eids(layers_y, V1, V2, lay_Vids, lay_eids);
1130  //determine the new coords of the vertices and the curve points
1131  //for each edge
1132 
1133 
1134  //------------------------------------------------------------------
1135 
1136 
1137 
1138 
1139 
1140 /*
1141 //possible interface to alglib!!!!!!!!
1142  Array<OneD, NekDouble> xalg(5, 1.0); //= "[-1.0,-0.5,0.0,+0.5,+1.0]";
1143  Array<OneD, NekDouble> yalg(5, 1.0); //= "[+1.0,0.25,0.0,0.25,+1.0]";
1144  Array<OneD, NekDouble> calg(1,3.0);
1145  double t = 0.25;
1146  double v;
1147  double diffstep = 0.001;
1148  lsfitstate state;
1149  alglib::Lsfitcreatef(xalg,yalg,calg,diffstep,state);
1150 
1151  //alglib::spline1dinterpolant s;
1152 
1153 */
1154 
1155 
1156 
1157 
1158 
1159 
1160 
1161 
1162 
1163 
1164 
1165 
1166 
1167 
1168 
1169 
1170  // curve the edges around the NEW critical layer (bool to turn on/off)
1171  bool curv_lay=true;
1172  bool move_norm=true;
1173  int np_lay = (nvertl-1)*npedge;//nedges*npedge (Eq. Points!!!)
1174  int nqedge = streak->GetExp(0)->GetNumPoints(0);
1175  Array<OneD, NekDouble> xcQ(nqedge*nedges,0.0);
1176  Array<OneD, NekDouble> ycQ(nqedge*nedges,0.0);
1177  Array<OneD, NekDouble> zcQ(nqedge*nedges,0.0);
1178  Array<OneD, NekDouble> nxPhys(npedge*nedges,0.0);
1179  Array<OneD, NekDouble> nyPhys(npedge*nedges,0.0);
1180  Array<OneD, NekDouble> nxQ(nqedge*nedges,0.0);
1181  Array<OneD, NekDouble> nyQ(nqedge*nedges,0.0);
1182  if( move_norm==true)
1183  {
1184  //np_lay = (nvertl-1)*nqedge;//nedges*nqedge (Q points!!!)
1185  //extract crit lay normals (through tangents):
1186 
1187  Array<OneD, StdRegions::StdExpansion1DSharedPtr> Edge_newcoords(2);
1188 
1189 cout<<"nquad per edge="<<nqedge<<endl;
1190  for(int l=0; l<2; l++)
1191  {
1192  Edge_newcoords[l] = bndfieldx[lastIregion]->GetExp(0)
1193  ->as<StdRegions::StdExpansion1D>();
1194  }
1195  Array<OneD, NekDouble> xnull(nqedge);
1196  Array<OneD, NekDouble> ynull(nqedge);
1197  Array<OneD, NekDouble> xcedgeQ(nqedge,0.0);
1198  Array<OneD, NekDouble> ycedgeQ(nqedge,0.0);
1199  Array<OneD, NekDouble> txedgeQ(nqedge,0.0);
1200  Array<OneD, NekDouble> tyedgeQ(nqedge,0.0);
1201  Array<OneD, NekDouble> normsQ(nqedge,0.0);
1202 
1203  Array<OneD, NekDouble> txQ(nqedge*nedges,0.0);
1204  Array<OneD, NekDouble> tyQ(nqedge*nedges,0.0);
1205  Array<OneD, NekDouble> tx_tyedgeQ(nqedge,0.0);
1206  Array<OneD, NekDouble> nxedgeQ(nqedge,0.0);
1207  Array<OneD, NekDouble> nyedgeQ(nqedge,0.0);
1208  Array<OneD, const NekDouble> ntempQ(nqedge) ;
1209 
1210  Array<OneD, NekDouble> nxedgePhys(npedge,0.0);
1211  Array<OneD, NekDouble> nyedgePhys(npedge,0.0);
1212 
1213 
1214  Array<OneD, NekDouble> CoordsPhys(2);
1215 
1216 
1217  bndfieldx[lastIregion]->GetCoords(xcQ, ycQ, zcQ);
1218  //determine the NEW crit lay quad points values(lagrangeInterpolant):
1219 /*
1220  for(int w=0; w< nqedge*nedges; w++)
1221  {
1222 
1223  }
1224 */
1225  //determine the NEW crit lay quad points values(Newtonit):
1226 cout<<"xcQ, ycQ"<<endl;
1227  Computestreakpositions(nqedge*nedges, streak, xcQ, ycQ,xnull,ynull,xnull,xnull,
1228  xnull,ynull, xcQ,ycQ, cr,false);
1229 /*
1230 for(int s=0; s<xcQ.num_elements(); s++)
1231 {
1232 cout<<xcQ[s]<<" "<<ycQ[s]<<endl;
1233 }
1234 ASSERTL0(false, "dsdfs");
1235 */
1236 
1237  for(int k=0; k<nedges; k++)
1238  {
1239  Vmath::Vcopy(nqedge, &xcQ[k*nqedge],1,&xcedgeQ[0],1);
1240  Vmath::Vcopy(nqedge, &ycQ[k*nqedge],1,&ycedgeQ[0],1);
1241  //calc the NEW tangent values
1242  Edge_newcoords[0]->StdPhysDeriv(xcedgeQ,txedgeQ);
1243  Edge_newcoords[1]->StdPhysDeriv(ycedgeQ,tyedgeQ);
1244  //norms=tx*tx
1245  Vmath::Vmul(nqedge,txedgeQ,1,txedgeQ,1,normsQ,1);
1246  //norms=tx*tx+ty*ty
1247  Vmath::Vvtvp(nqedge,tyedgeQ,1,tyedgeQ,1,normsQ,1,normsQ,1);
1248 
1249  Vmath::Vsqrt(nqedge, normsQ,1,normsQ,1);
1250  Vmath::Sdiv(nqedge,1.0,normsQ,1,normsQ,1);
1251 
1252  Vmath::Vmul(nqedge,txedgeQ,1,normsQ,1,txedgeQ,1);
1253  Vmath::Vmul(nqedge,tyedgeQ,1,normsQ,1,tyedgeQ,1);
1254 
1255 
1256  //determine the normal from eqs(APART FROM SIGN):
1257  //tx^2 +ty^2= 1 = nx^2 + ny^2;
1258  //t\cdot n=0= tx*nx +ty*ny
1259  //result: nx = ( 1+(tx/ty)^2 )^(-1/2)
1260  Vmath::Vdiv(nqedge, txedgeQ,1,tyedgeQ,1,tx_tyedgeQ,1);
1261  Vmath::Vmul(nqedge, tx_tyedgeQ,1,tx_tyedgeQ,1,tx_tyedgeQ,1);
1262  Vmath::Sadd(nqedge,1.0,tx_tyedgeQ,1,nxedgeQ,1);
1263  Vmath::Vsqrt(nqedge,nxedgeQ,1,nxedgeQ,1);
1264  Vmath::Sdiv(nqedge,1.0,nxedgeQ,1,nxedgeQ,1);
1265  //normal DOWNWARDS!!! mult by -1
1266  Vmath::Smul(nqedge, -1.0,nxedgeQ,1,nxedgeQ,1);
1267  Vmath::Vcopy(nqedge, &(nxedgeQ[0]),1, &(nxQ[nqedge*k]),1);
1268  //ny = (1-nx ^2)^(1/2)
1269  Vmath::Vmul(nqedge, nxedgeQ,1,nxedgeQ,1,nyedgeQ,1);
1270  Vmath::Smul(nqedge, -1.0,nyedgeQ,1,nyedgeQ,1);
1271  Vmath::Sadd(nqedge,1.0,nyedgeQ,1,nyedgeQ,1);
1272  Vmath::Vsqrt(nqedge,nyedgeQ,1,nyedgeQ,1);
1273  //normal DOWNWARDS!!! mult by -1
1274  Vmath::Smul(nqedge, -1.0,nyedgeQ,1,nyedgeQ,1);
1275  Vmath::Vcopy(nqedge, &(nyedgeQ[0]), 1, &(nyQ[nqedge*k]),1);
1276 /*
1277 cout<<"edge:"<<k<<endl;
1278 cout<<"tan/normal"<<endl;
1279 for(int r=0; r<nqedge; r++)
1280 {
1281 cout<<xcQ[k*nqedge+r]<<" "<<txedgeQ[r]<<" "<<tyedgeQ[r]<<" "
1282 <<nxedgeQ[r]<<" "<<nyedgeQ[r]<<endl;
1283 }
1284 */
1285  //force the normal at interface point to be equal
1286  if(k>0)
1287  {
1288  //nyPhys[f*npedge +0] =
1289  // (nyPhys[(f-1)*npedge+npedge-1]+nyPhys[f*npedge+0])/2.;
1290  nyQ[(k-1)*nqedge+nqedge-1]=
1291  nyQ[k*nqedge+0];
1292  //nx= (1-ny^2)^{1/2}
1293  //nxPhys[f*npedge+0]=
1294  // sqrt(1- nyPhys[f*npedge+0]*nyPhys[f*npedge+0]);
1295  nxQ[(k-1)*nqedge+nqedge-1]=
1296  nxQ[k*nqedge+0];
1297  }
1298  }
1299 
1300  int nquad_lay = (nvertl-1)*nqedge;
1301  Array<OneD, NekDouble>x_tmpQ(nquad_lay-(nedges-1));
1302  //Array<OneD, NekDouble>tmpnxQ(nquad_lay-(nedges-1));
1303  Array<OneD, NekDouble>tmpnyQ(nquad_lay-(nedges-1));
1304 
1305  Cutrepetitions(nedges, xcQ,x_tmpQ);
1306  //Cutrepetitions(nedges, nxQ, tmpnxQ);
1307  Cutrepetitions(nedges, nyQ, tmpnyQ);
1308 /*
1309 for(int u=0; u<x_tmpQ.num_elements(); u++)
1310 {
1311 cout<<x_tmpQ[u]<<" "<<tmpnyQ[u]<<endl;
1312 }
1313 */
1314 
1315 
1316  //interpolate the values into phys points(curved points)
1317  for(int k=0; k<nedges; k++)
1318  {
1319 //cout<<"edge:"<<k<<endl;
1320  for(int a=0; a<npedge; a++)
1321  {
1322  if(a==0)//verts pos no interp necessary
1323  {
1324  nxPhys[k*npedge +a]= nxQ[k*nqedge +0];
1325  nyPhys[k*npedge +a]= nyQ[k*nqedge +0];
1326 
1327  }
1328  else if(a== npedge-1)//verts pos no interp necessary
1329  {
1330  nxPhys[k*npedge +a]= nxQ[k*nqedge +nqedge-1];
1331  nyPhys[k*npedge +a]= nyQ[k*nqedge +nqedge-1];
1332 //cout<<":last"<<nyQ[k*nqedge+a]<<endl;
1333 
1334  }
1335 
1336  else
1337  {
1338  //use lagrange interpolant to get the
1339  //normal at phys(equispaced points)
1340 
1341  //order normal functions(cut out repetitions)
1342  //QUAD POINTS
1343 
1344 
1345  int index;
1346  //determine closest index:
1347 
1348  index=
1349  DetermineclosePointxindex( Addpointsx[k*(npedge-2) +a-1], x_tmpQ);
1350 
1351  Array<OneD, NekDouble> Pxinterp(4);
1352  Array<OneD, NekDouble> Pyinterp(4);
1353 
1354  //generate neighbour arrays (y):
1355  GenerateNeighbourArrays(index, 4,x_tmpQ,tmpnyQ,Pxinterp,Pyinterp);
1356  //interp the new normal components(y)
1357  nyPhys[k*npedge +a]=
1358  LagrangeInterpolant(Addpointsx[k*(npedge-2) +a-1],4,Pxinterp,Pyinterp );
1359 /*
1360  //generate neighbour arrays (x):
1361  GenerateNeighbourArrays(index,4,x_tmpQ,tmpnxQ,Pxinterp,Pyinterp);
1362  //interp the new normal components(x)
1363  nxPhys[k*npedge +a]=
1364  LagrangeInterpolant(Addpointsx[k*(npedge-2) +a],4,Pxinterp,Pyinterp );
1365 */
1366  //nx=-(1-ny*ny){1/2} the normal is DOWNWARDS!!!
1367  nxPhys[k*npedge +a]= -sqrt(abs(1- nyPhys[k*npedge +a]*nyPhys[k*npedge +a]));
1368 /*
1369  //(put the middle points as quad points)
1370  //GaussLobattoLegendre points in the middle:
1371  nxPhys[k*npedge +a] = nxedgeQ[a];
1372  nyPhys[k*npedge +a] = nyedgeQ[a];
1373  ASSERTL0(npedge< nqedge," quad points too low");
1374 */
1375  }
1376 
1377 
1378  //force the normal at interface point to be equal
1379  if(k>0)
1380  {
1381  //nyPhys[f*npedge +0] =
1382  // (nyPhys[(f-1)*npedge+npedge-1]+nyPhys[f*npedge+0])/2.;
1383  nyPhys[(k-1)*npedge+npedge-1]=
1384  nyPhys[k*npedge+0];
1385  //nx= (1-ny^2)^{1/2}
1386  //nxPhys[f*npedge+0]=
1387  // sqrt(1- nyPhys[f*npedge+0]*nyPhys[f*npedge+0]);
1388  nxPhys[(k-1)*npedge+npedge-1]=
1389  nxPhys[k*npedge+0];
1390  }
1391 
1392 
1393 
1394  }
1395 
1396 
1397 
1398 
1399 
1400  }
1401 /*
1402 for(int s=0; s<np_lay; s++)
1403 {
1404 
1405 cout<<xcPhys[s]<<" "<<nxPhys[s]<<" "<<nyPhys[s]<<endl;
1406 
1407 }
1408 cout<<"xcPhys,,"<<endl;
1409 */
1410 
1411 
1412  //determine the new coords of the vertices and the curve points
1413  //for each edge
1414 
1415  //int np_lay = (nvertl-1)*npedge;//nedges*npedge
1416 
1417  //NB delta=ynew-y_c DEPENDS ON the coord trnaf ynew= y+Delta_c*ratio!!!
1418  Array<OneD, NekDouble> delta(nlays);
1419  Array<OneD, NekDouble>tmpy_lay(np_lay);
1420  Array<OneD, NekDouble>tmpx_lay(np_lay);
1421  for(int m=0; m<nlays; m++)
1422  {
1423  //delta[m] = (ynew[lay_Vids[m][0]] - y_c[0])/1.0;
1424 
1425  //depends on Delta0
1426  if(m< cntlow+1)
1427  {
1428  delta[m] = -(cntlow+1-m)*Delta0/(cntlow+1);
1429  }
1430  else
1431  {
1432  delta[m] = ( m-(cntlow) )*Delta0/(cntlow+1);
1433  }
1434 
1435 
1436  layers_x[m]= Array<OneD, NekDouble>(np_lay);
1437 //cout<<"delta="<<delta[m]<<" cntlow="<<cntlow<<endl;
1438 
1439  for(int h=0; h< nvertl; h++)
1440  {
1441  //shift using the dinstance delta from the crit layer AT x=0
1442  //for each layer
1443 
1444 
1445 //cout<<m<<"Vid:"<<lay_Vids[m][h]<<" mod from y="<<ynew[lay_Vids[m][h] ]<<" to y="<<y_c[h] +delta[m]<<endl;
1446  if(move_norm==false)
1447  {
1448  ynew[lay_Vids[m][h] ]= y_c[h] +delta[m];
1449  xnew[lay_Vids[m][h] ]= x_c[h];
1450  }
1451  else
1452  {
1453  if(h==0 || h==nvertl-1 )//nx=0,ny=1 at the borders
1454  {
1455  ynew[lay_Vids[m][h] ]= y_c[h] +delta[m];
1456  xnew[lay_Vids[m][h] ]= x_c[h];
1457  }
1458  else
1459  {
1460  ynew[lay_Vids[m][h] ]= y_c[h] +delta[m]*abs(nyPhys[h*npedge+0]);
1461  xnew[lay_Vids[m][h] ]= x_c[h] +delta[m]*abs(nxPhys[h*npedge+0]);
1462  }
1463  }
1464 //cout<<"Vid x="<<xnew[lay_Vids[m][h] ]<<" y="<<ynew[lay_Vids[m][h] ]<<endl;
1465  if(h< nedges && curv_lay==true)
1466  {
1467 //cout<<"edge=="<<h<<endl;
1468  if(h>0)//check normal consistency
1469  {
1470  ASSERTL0( nyPhys[h*npedge+0]==nyPhys[(h-1)*npedge+npedge-1]," normaly wrong");
1471  ASSERTL0( nxPhys[h*npedge+0]==nxPhys[(h-1)*npedge+npedge-1]," normalx wrong");
1472 
1473  }
1474  if(move_norm==false)
1475  {
1476  //v1
1477  layers_y[m][h*npedge +0] = y_c[h] +delta[m];
1478  layers_x[m][h*npedge +0] = xnew[lay_Vids[m][h] ];
1479  //v2
1480  layers_y[m][h*npedge +npedge-1] = y_c[h+1] +delta[m];
1481  layers_x[m][h*npedge +npedge-1] = xnew[lay_Vids[m][h+1] ];
1482  //middle points (shift crit lay points by delta):
1483  for(int d=0; d< npedge-2; d++)
1484  {
1485  layers_y[m][h*npedge +d+1]= Addpointsy[h*(npedge-2) +d] +delta[m];
1486  layers_x[m][h*npedge +d+1]= Addpointsx[h*(npedge-2) +d];
1487  }
1488  }
1489  else
1490  {
1491  if(h==0) //nx=0,ny=1 at the borders
1492  {
1493  //v1
1494  tmpy_lay[h*npedge +0] = y_c[h] +delta[m];
1495  tmpx_lay[h*npedge +0] = xnew[lay_Vids[m][h] ];
1496  //v2
1497  tmpy_lay[h*npedge +npedge-1] =
1498  y_c[h+1] +delta[m]*abs(nyPhys[h*npedge +npedge-1]);
1499  tmpx_lay[h*npedge +npedge-1] =
1500  x_c[h+1] +delta[m]*abs(nxPhys[h*npedge +npedge-1]);
1501  }
1502  else if(h==nedges-1)//nx=0,ny=1 at the borders
1503  {
1504  //v1
1505  tmpy_lay[h*npedge +0] =
1506  y_c[h] +delta[m]*abs(nyPhys[h*npedge +0]);
1507  tmpx_lay[h*npedge +0] =
1508  x_c[h] +delta[m]*abs(nxPhys[h*npedge +0]);
1509  //v2
1510  tmpy_lay[h*npedge +npedge-1] = y_c[h+1] +delta[m];
1511  tmpx_lay[h*npedge +npedge-1] = xnew[lay_Vids[m][h+1] ];
1512  }
1513  else
1514  {
1515  //v1
1516  tmpy_lay[h*npedge +0] =
1517  y_c[h] +delta[m]*abs(nyPhys[h*npedge +0]);
1518  tmpx_lay[h*npedge +0] =
1519  x_c[h] +delta[m]*abs(nxPhys[h*npedge +0]);
1520  //v2
1521  tmpy_lay[h*npedge +npedge-1] =
1522  y_c[h+1] +delta[m]*abs(nyPhys[h*npedge +npedge-1]);
1523  tmpx_lay[h*npedge +npedge-1] =
1524  x_c[h+1] +delta[m]*abs(nxPhys[h*npedge +npedge-1]);
1525  }
1526 
1527  //middle points
1528  for(int d=0; d< npedge-2; d++)
1529  {
1530 
1531  tmpy_lay[h*npedge +d+1] = Addpointsy[h*(npedge-2) +d] +
1532  delta[m]*abs(nyPhys[h*npedge +d+1]);
1533  tmpx_lay[h*npedge +d+1]= Addpointsx[h*(npedge-2) +d] +
1534  delta[m]*abs(nxPhys[h*npedge +d+1]);
1535 
1536  //NB ycQ,xcQ refers to nqedge NOT npedge!!!
1537  //tmpy_lay[h*npedge +d+1] = ycQ[h*nqedge +d+1] +
1538  // delta[m]*abs(nyPhys[h*npedge +d+1]);
1539  //tmpx_lay[h*npedge +d+1]= xcQ[h*nqedge +d+1] +
1540  // delta[m]*abs(nxPhys[h*npedge +d+1]);
1541 //cout<<"xmoved="<<tmpx_lay[h*npedge +d+1]<<" xold="<<xcQ[h*nqedge +d+1]<<endl;
1542  //ASSERTL0(tmpx_lay[h*npedge +d+1]>0," middle point with x<0")
1543 
1544  }
1545 
1546  }
1547 
1548 
1549 
1550  }//close edges
1551  }//close verts h
1552 /*
1553 for(int s=0; s<np_lay; s++)
1554 {
1555 
1556 cout<<tmpx_lay[s]<<" "<<tmpy_lay[s]<<endl;
1557 
1558 }
1559 //ASSERTL0(false, "dasd");
1560 */
1561 
1562 
1563  //ASSERTL0(tmpx_lay[h*npedge +0]>=0," vert 0 x<0");
1564  //ASSERTL0(tmpx_lay[h*npedge +npedge-1]>0," vert 1 x<0");
1565 
1566 
1567  //Orderfunctionx(tmpx_lay,tmpy_lay, Array<OneD, NekDouble>& outarray_x,
1568  // Array<OneD, NekDouble>& outarray_y);
1569 
1570 
1571 
1572  //check if the x coord is 'outofbound' and calculate the
1573  //number of outofbound points
1574 
1575  //determine which boudn has been overcome:
1576  NekDouble boundleft = xcPhys[0];
1577  NekDouble boundright = xcPhys[np_lay-1];
1578  bool outboundleft= false;
1579  bool outboundright=false;
1580  if(tmpx_lay[1]< boundleft )
1581  {
1582  outboundleft = true;
1583  }
1584  if(tmpx_lay[np_lay-2] > boundright )
1585  {
1586  outboundright = true;
1587  }
1588 
1589 
1590 
1591  int outvert=0;
1592  int outmiddle=0;
1593  int outcount=0;
1594 
1595  Array<OneD, int> vertout(nvertl);
1596  for(int r=0; r< nedges; r++)
1597  {
1598  //check point outofboundleft
1599  if(tmpx_lay[r*npedge + npedge-1]< boundleft && outboundleft==true )//assume the neg coords start from 0
1600  {
1601  vertout[outvert]=r;
1602  outvert++;
1603 
1604  if(r<nedges-1 )
1605  {
1606  //check if after the last negvert there are neg points
1607  if( tmpx_lay[(r+1)*npedge + npedge-1]> boundleft )
1608  {
1609 
1610  for(int s=0; s<npedge-2; s++)
1611  {
1612  if(tmpx_lay[(r+1)*npedge + s+1]< boundleft)
1613  {
1614  outmiddle++;
1615  }
1616  }
1617 
1618  }
1619  }
1620  }
1621  //check point outofboundright
1622  if(tmpx_lay[r*npedge + 0]> boundright && outboundright==true )//assume the neg coords start from 0
1623  {
1624  vertout[outvert]=r;
1625  outvert++;
1626 
1627  if( r> 0)
1628  {
1629  //check if after the last outvert there are out points
1630  if( tmpx_lay[(r-1)*npedge + 0]< boundright )
1631  {
1632 
1633  for(int s=0; s<npedge-2; s++)
1634  {
1635  if(tmpx_lay[(r-1)*npedge + s+1]> boundright)
1636  {
1637  outmiddle++;
1638  }
1639  }
1640 
1641  }
1642  }
1643  }
1644  }
1645  //calc number of point to replace
1646  outcount = outvert*npedge+1+ outmiddle;
1647  //determine from(until) which index the replacement will start
1648  int replacepointsfromindex=0;
1649  for(int c=0; c<nedges; c++)
1650  {
1651  //assume at least 1 middle point per edge
1652  if(xcPhys[c*npedge+npedge-1] <= tmpx_lay[c*(npedge-(npedge-2)) +2] && outboundright==true)
1653  {
1654  replacepointsfromindex = c*(npedge-(npedge-2))+2;
1655  break;
1656  }
1657 
1658 
1659  //assume at least 1 middle point per edge
1660  if(xcPhys[(nedges-1 -c)*npedge+0] >= tmpx_lay[np_lay-1 -(c*(npedge-(npedge-2)) +2)]
1661  && outboundleft==true)
1662  {
1663  replacepointsfromindex = np_lay-1 -(c*(npedge-(npedge-2)) +2);
1664  break;
1665  }
1666 
1667 
1668  }
1669 
1670 
1671 
1672 
1673 //cout<<"out="<<outcount<<endl;
1674 //cout<<"replacefrom="<<replacepointsfromindex<<endl;
1675 
1676 
1677  //if xcoord is neg find the first positive xcoord
1678 
1679 
1680  if(outcount>1)
1681  {
1682  //determine x new coords:
1683  //distribute the point all over the layer
1684  int pstart,shift;
1685  NekDouble increment;
1686 
1687  if( outboundright==true)
1688  {
1689  pstart = replacepointsfromindex;
1690  shift = np_lay-outcount;
1691  increment = (xcPhys[np_lay-outcount]-xcPhys[pstart])/(outcount+1);
1692  outcount = outcount-1;
1693  ASSERTL0(tmpx_lay[np_lay-outcount]>xcPhys[(nedges-1)*npedge+0], "no middle points in the last edge");
1694  }
1695  else
1696  {
1697  shift=1;
1698  pstart= outcount-1;
1699  increment = (xcPhys[replacepointsfromindex]-xcPhys[pstart])/(outcount+1);
1700  ASSERTL0(tmpx_lay[pstart]<xcPhys[0*npedge +npedge-1], "no middle points in the first edge");
1701  }
1702 
1703  //interp to points between posindex and posindex-1
1704  Array<OneD, NekDouble> replace_x(outcount);
1705  Array<OneD, NekDouble> replace_y(outcount);
1706  //order normal functions(cut out repetitions)
1707  Array<OneD, NekDouble>x_tmp(np_lay-(nedges-1));
1708  Array<OneD, NekDouble>y_tmp(np_lay-(nedges-1));
1709  Array<OneD, NekDouble>tmpny(np_lay-(nedges-1));
1710  Cutrepetitions(nedges, xcPhys,x_tmp);
1711  Cutrepetitions(nedges, ycPhys,y_tmp);
1712  Cutrepetitions(nedges, nyPhys, tmpny);
1713  //init neigh arrays
1714  Array<OneD, NekDouble>closex(4);
1715  Array<OneD, NekDouble>closey(4);
1716  Array<OneD, NekDouble>closeny(4);
1717  NekDouble xctmp,ycinterp,nxinterp,nyinterp;
1718 
1719 
1720 
1721  for(int v=0; v<outcount;v++)
1722  {
1723  xctmp = xcPhys[pstart]+(v+1)*increment;
1724 
1725 
1726 
1727  //determine closest point index:
1728  int index =
1729  DetermineclosePointxindex( xctmp, x_tmp);
1730 //cout<<" vert="<<index<<endl;
1731 
1732 
1733  //generate neighbour arrays (ny)
1734  GenerateNeighbourArrays(index, 4,x_tmp,tmpny,closex,closeny);
1735 
1736 
1737  //interp:
1738  nyinterp =
1739  LagrangeInterpolant(
1740  xctmp,4,closex,closeny );
1741 
1742  //calc nxinterp
1743  nxinterp = sqrt(abs(1-nyinterp*nyinterp));
1744 
1745  //generata neighbour arrays (yc)
1746  GenerateNeighbourArrays(index, 4,x_tmp,y_tmp,closex,closey);
1747  //interp:
1748  ycinterp = LagrangeInterpolant(xctmp,4, closex,closey);
1749  //calc new coord
1750  replace_x[v] = xctmp +delta[m]*abs(nxinterp);
1751  replace_y[v] = ycinterp +delta[m]*abs(nyinterp);
1752  tmpx_lay[ v+shift ] = replace_x[v];
1753  tmpy_lay[ v+shift ] = replace_y[v];
1754 
1755 //cout<<"xinterp="<<replace_x[v]<<" yinterp="<<replace_y[v]<<endl;
1756 
1757 
1758 
1759 
1760  }
1761  }
1762 
1763 
1764 
1765 /*
1766 for(int s=0; s<np_lay; s++)
1767 {
1768 
1769 cout<<tmpx_lay[s]<<" "<<tmpy_lay[s]<<endl;
1770 
1771 }
1772 if(m== 0)
1773 {
1774 //ASSERTL0(false, "ssa");
1775 }
1776 */
1777 
1778 
1779 
1780 
1781 
1782 
1783 
1784  int closepoints = 4;
1785 
1786  Array<OneD, NekDouble> Pyinterp(closepoints);
1787  Array<OneD, NekDouble> Pxinterp(closepoints);
1788 
1789  //check if any edge has less than npedge points
1790  int pointscount=0;
1791  for(int q=0; q<np_lay; q++)
1792  {
1793  for(int e=0; e<nedges; e++)
1794  {
1795  if(tmpx_lay[q]<= x_c[e+1] && tmpx_lay[q]>= x_c[e])
1796  {
1797  pointscount++;
1798  }
1799  if(q == e*npedge +npedge-1 && pointscount!=npedge )
1800  {
1801 //cout<<"edge with few points :"<<e<<endl;
1802  pointscount=0;
1803  }
1804  else if(q == e*npedge +npedge-1)
1805  {
1806  pointscount=0;
1807  }
1808  }
1809  }
1810  //----------------------------------------------------------
1811 /*
1812 cout<<"notordered"<<endl;
1813 for(int g=0; g<tmpx_lay.num_elements(); g++)
1814 {
1815 cout<<tmpx_lay[g]<<" "<<tmpy_lay[g]<<endl;
1816 }
1817 */
1818 
1819 
1820 cout<<nedges<<"nedges"<<npedge<<" np_lay="<<np_lay<<endl;
1821  //generate x,y arrays without lastedgepoint
1822  //(needed to interp correctly)
1823  Array<OneD, NekDouble>tmpx(np_lay-(nedges-1));
1824  Array<OneD, NekDouble>tmpy(np_lay-(nedges-1));
1825 /*
1826  for(int b=0; b<nedges; b++)
1827  {
1828  Vmath::Vcopy( npedge-1, &tmpy_lay[b*(npedge)],1,&tmpy[b*(npedge-1)],1);
1829  Vmath::Vcopy( npedge-1, &tmpx_lay[b*(npedge)],1,&tmpx[b*(npedge-1)],1);
1830  if(b== nedges-1)
1831  {
1832  tmpy[b*(npedge-1)+npedge-1] = tmpy_lay[b*npedge+npedge-1];
1833  tmpx[b*(npedge-1)+npedge-1] = tmpx_lay[b*npedge+npedge-1];
1834  }
1835  }
1836 
1837  //check if there is any repetition
1838  for(int f=1; f< tmpx.num_elements(); f++)
1839  {
1840  string error= "repetition on arrays interp:"+boost::lexical_cast<string>(m);
1841  ASSERTL0(tmpx[f]!=tmpx[f-1],error);
1842  }
1843 */
1844 
1845  Cutrepetitions(nedges, tmpx_lay, tmpx);
1846  Cutrepetitions(nedges, tmpy_lay, tmpy);
1847 
1848 
1849  //order points in x:
1850  int index;
1851  Array<OneD, NekDouble> copyarray_x(tmpx.num_elements());
1852  Array<OneD, NekDouble> copyarray_y(tmpx.num_elements());
1853 /*
1854  //NB: it may be possible that some edges have
1855  //more/less points than npedge
1856  //@todo: write an orderfunction
1857 
1858  Vmath::Vcopy(tmpx.num_elements(), tmpx,1, copyarray_x,1);
1859  Vmath::Vcopy(tmpx.num_elements(), tmpy,1, copyarray_y,1);
1860  NekDouble max = Vmath::Vmax(tmpx.num_elements(), tmpx,1);
1861  for(int w=0; w<tmpx.num_elements(); w++)
1862  {
1863  index = Vmath::Imin(tmpx.num_elements(), copyarray_x,1);
1864  tmpx[w]= copyarray_x[index];
1865  tmpy[w]= copyarray_y[index];
1866  copyarray_x[index] = max+1000;
1867  }
1868 */
1869 
1870 
1871  Orderfunctionx(tmpx, tmpy, tmpx, tmpy);
1872 /*
1873 cout<<"ordered"<<endl;
1874 for(int g=0; g<tmpx.num_elements(); g++)
1875 {
1876 cout<<tmpx[g]<<" "<<tmpy[g]<<endl;
1877 }
1878 */
1879 
1880 
1881 
1882  //determine the neighbour points (-3;+3)
1883  for(int g=0; g< nvertl; g++)
1884  {
1885  //verts
1886 //cout<<"determine value for vert x="<<x_c[g]<<endl;
1887  //determine closest index:
1888 
1889  index=
1890  DetermineclosePointxindex( x_c[g], tmpx);
1891  //generate neighbour arrays:
1892  GenerateNeighbourArrays(index, closepoints,tmpx,tmpy,Pxinterp,Pyinterp);
1893  //write vert coords
1894  ynew[lay_Vids[m][g] ]= LagrangeInterpolant(x_c[g],closepoints,Pxinterp,Pyinterp );
1895  xnew[lay_Vids[m][g] ]= x_c[g];
1896 
1897 
1898  if(g<nedges)
1899  {
1900 
1901 
1902 
1903  //v1
1904  layers_y[m][g*npedge +0] = ynew[lay_Vids[m][g] ];
1905  layers_x[m][g*npedge +0] = xnew[lay_Vids[m][g] ];
1906  //v2
1907 
1908  //determine closest index:
1909  index=
1910  DetermineclosePointxindex( x_c[g+1], tmpx);
1911  //generate neighbour arrays:
1912  GenerateNeighbourArrays(index, closepoints,tmpx,tmpy,Pxinterp,Pyinterp);
1913  layers_y[m][g*npedge +npedge-1] =
1914  LagrangeInterpolant(x_c[g+1],closepoints,Pxinterp,Pyinterp );
1915  layers_x[m][g*npedge +npedge-1] = x_c[g+1];
1916 
1917 
1918 
1919  //middle points
1920  for(int r=0; r< npedge-2; r++)
1921  {
1922  //cout<<"determine value for x="<<Addpointsx[g*(npedge-2) +r]<<endl;
1923  //determine closest point index:
1924  index =
1925  DetermineclosePointxindex( Addpointsx[g*(npedge-2) +r], tmpx);
1926  //cout<<" vert+"<<index<<endl;
1927 
1928  ASSERTL0( index<= tmpy.num_elements()-1, " index wrong");
1929  //generate neighbour arrays Pyinterp,Pxinterp
1930  GenerateNeighbourArrays(index, closepoints,tmpx,tmpy,Pxinterp,Pyinterp);
1931 
1932  layers_y[m][g*npedge +r+1]=
1933  LagrangeInterpolant(
1934  Addpointsx[g*(npedge-2) +r],closepoints,Pxinterp,Pyinterp );
1935  //cout<<"x value="<<Addpointsx[g*(npedge-2) +r]<<endl;
1936  layers_x[m][g*npedge +r+1]= Addpointsx[g*(npedge-2) +r];
1937 
1938  /*
1939  for(int t=0; t<6; t++)
1940  {
1941  cout<<"Px="<<Pxinterp[t]<<" "<<Pyinterp[t]<<endl;
1942  }
1943  */
1944  }
1945 
1946  }//if edge closed g
1947  }
1948  //check if there are points out of range:
1949  ASSERTL0(Vmath::Vmax(np_lay,layers_y[m],1)< Vmath::Vmax(nVertTot,yold,1),"point>ymax");
1950  ASSERTL0(Vmath::Vmin(np_lay,layers_y[m],1)> Vmath::Vmin(nVertTot,yold,1),"point<ymin");
1951 
1952 
1953 /*
1954 cout<<" xlay ylay"<<endl;
1955 for(int l=0; l<np_lay; l++)
1956 {
1957 //cout<<tmpx_lay[l]<<" "<<tmpy_lay[l]<<endl;
1958 cout<<std::setprecision(8)<<layers_x[m][l]<<" "<<layers_y[m][l]<<endl;
1959 }
1960 
1961 cout<<"nverts"<<endl;
1962 for(int l=0; l<nvertl; l++)
1963 {
1964 cout<<std::setprecision(8)<<xnew[lay_Vids[m][l] ]<<" "<<ynew[lay_Vids[m][l] ]<<endl;
1965 }
1966 */
1967 
1968 //ASSERTL0(false, "as");
1969 
1970  }//close layers!!! m index
1971 
1972 
1973 
1974  //update vertices coords outside layers region
1975  NekDouble ratio;
1976  for(int n=0; n<nVertTot; n++)
1977  {
1978  NekDouble ratio;
1979  SpatialDomains::PointGeomSharedPtr vertex = graphShPt->GetVertex(n);
1980  NekDouble x,y,z;
1981  vertex->GetCoords(x,y,z);
1982  int qp_closer;
1983  //determine the closer xold_up
1984  NekDouble tmp=1000;
1985 
1986  for(int k=0; k<nvertl; k++)
1987  {
1988  if(abs(x-xold_c[k]) < tmp)
1989  {
1990  tmp = abs(x-xold_c[k]);
1991  qp_closer=k;
1992  }
1993 
1994  }
1995  //find nplay_closer
1996  int nplay_closer;
1997  if(qp_closer==0)
1998  {
1999  nplay_closer=0;//first vert
2000  }
2001  else
2002  {
2003  nplay_closer= (qp_closer-1)*npedge +npedge-1;
2004  }
2005 
2006 
2007  if( y>yold_up[qp_closer] && y<1 )//nlays-1 is layerup
2008  {
2009 
2010 // ratio = (1-layers_y[nlays-1][qp_closer])*(1-y_c[qp_closer])/
2011 // ( (1-yold_up[n])*(1-yold_c[qp_closer]) );
2012  ratio = (1-layers_y[nlays-1][nplay_closer])/
2013  ( (1-yold_up[qp_closer]) );
2014  //distance prop to layerup
2015  ynew[n] = layers_y[nlays-1][nplay_closer]
2016  + (y-yold_up[qp_closer])*ratio;
2017  xnew[n] = x;
2018 
2019  }
2020  else if( y< yold_low[qp_closer] && y>-1 )//0 is layerdown
2021  {
2022 
2023  ratio = (1+layers_y[0][nplay_closer])/
2024  ( (1+yold_low[qp_closer]) );
2025  //distance prop to layerlow
2026  ynew[n] = layers_y[0][nplay_closer]
2027  + (y-yold_low[qp_closer])*ratio;
2028  xnew[n] = x;
2029  }
2030 
2031  }
2032 
2033  }//move_norm bool
2034  else//move vertically
2035  {
2036  MoveLayersvertically(nlays, nvertl, cntlow, cntup,
2037  lay_Vids, x_c, y_c, Down, Up, xnew, ynew, layers_x, layers_y);
2038 
2039  }
2040 
2041 
2042  //check borders of the new mesh verts:
2043 //cout<<std::setprecision(8)<<"yoldmax="<<Vmath::Vmax(nVertTot, yold,1)<<endl;
2044 //cout<<std::setprecision(8)<<"ynewmax="<<Vmath::Vmax(nVertTot,ynew,1)<<endl;
2045  ASSERTL0(Vmath::Vmin(nVertTot, xold,1)== Vmath::Vmin(nVertTot,xnew,1),
2046  " different xmin val");
2047  ASSERTL0(Vmath::Vmin(nVertTot, yold,1)== Vmath::Vmin(nVertTot,ynew,1),
2048  " different ymin val");
2049  ASSERTL0(Vmath::Vmax(nVertTot, xold,1)== Vmath::Vmax(nVertTot,xnew,1),
2050  " different xmax val");
2051  ASSERTL0(Vmath::Vmax(nVertTot, yold,1)== Vmath::Vmax(nVertTot,ynew,1),
2052  " different ymax val");
2053 
2054 
2055 
2056  //check x,y coords
2057 /*
2058  NekDouble tmpdiff;
2059  for(int f=0; f< nlays; f++)
2060  {
2061 //cout<<"delta lay="<<f<<" is:"<<delta[f]<<endl;
2062  for(int u=0; u< nvertl; u++)
2063  {
2064  if(u==0)
2065  {
2066 //cout<<"vert x="<<xnew[lay_Vids[f][u] ]<<" curv x="<<layers_x[f][0]<<endl;
2067  ASSERTL0(xnew[lay_Vids[f][u] ]==layers_x[f][0],"wro");
2068 //cout<<"vert y="<<ynew[lay_Vids[f][u] ]<<" curv y="<<layers_y[f][0]<<endl;
2069  ASSERTL0(ynew[lay_Vids[f][u] ]==layers_y[f][0],"wroy");
2070  }
2071  else if(u== nvertl-1)
2072  {
2073 //cout<<"vert x="<<xnew[lay_Vids[f][u] ]<<" curv x="<<layers_x[f][(nedges-1)*npedge +npedge-1]<<endl;
2074  ASSERTL0(
2075  xnew[lay_Vids[f][u] ]==layers_x[f][(nedges-1)*npedge +npedge-1],"wroLAST");
2076 
2077 
2078 //cout<<"vert y="<<ynew[lay_Vids[f][u] ]<<" curv y="<<layers_y[f][(nedges-1)*npedge +npedge-1]<<endl;
2079  ASSERTL0(ynew[lay_Vids[f][u] ]==layers_y[f][(nedges-1)*npedge +npedge-1],"wroyLAST");
2080  }
2081  else
2082  {
2083 //cout<<"vert x="<<xnew[lay_Vids[f][u] ]<<" curv x1="<<layers_x[f][u*npedge +0]<<" curv x2="
2084 //<<layers_x[f][(u-1)*npedge +npedge-1]<<endl;
2085  ASSERTL0(xnew[lay_Vids[f][u] ]== layers_x[f][u*npedge +0],"wrongg");
2086  ASSERTL0(layers_x[f][u*npedge +0]== layers_x[f][(u-1)*npedge +npedge-1],"wrong1");
2087 
2088 //cout<<"vert y="<<ynew[lay_Vids[f][u] ]<<" curv y1="<<layers_y[f][u*npedge +0]<<" curv y2="
2089 //<<layers_y[f][(u-1)*npedge +npedge-1]<<endl;
2090  ASSERTL0(ynew[lay_Vids[f][u] ]== layers_y[f][u*npedge +0],"wronggy");
2091  ASSERTL0(layers_y[f][u*npedge +0]== layers_y[f][(u-1)*npedge +npedge-1],"wrong1y");
2092  }
2093  }
2094  }
2095 */
2096 
2097 
2098  //------------------------------------------------------------
2099 
2100  //check if the layers need to be shifted
2101 /*
2102  string movelay;
2103  int Vid_di;
2104  NekDouble delt, delt_opp;
2105  NekDouble tmpleft=10;
2106  NekDouble tmpright = -10;
2107 
2108  //check top left (x=0, y>yup)
2109  //last layer nlays-1
2110  if(
2111  ( abs(ynew[topleft]-ynew[ lay_Vids[nlays-1][0] ])<0.01)||
2112  ( ynew[topleft]< ynew[lay_Vids[nlays-1][0] ])
2113  )
2114  {
2115 
2116  movelay = "layup";
2117  //calculate the new delta based on the space at the border:
2118  delt = abs(y_c[0]-ynew[topleft])/(nlays/2);
2119  delt_opp = abs(y_c[0]-ynew[bottomleft])/(nlays/2);
2120 cout<<"move layerup to down:"<<endl;
2121  }
2122  //check bottom right (x=pi/2, y<ydown)
2123  //first lay: 0
2124  if(
2125  ( abs(ynew[bottomright]-ynew[ lay_Vids[0][(nvertl-1)] ])<0.01 ) ||
2126  ( ynew[bottomright]> ynew[ lay_Vids[0][(nvertl-1)] ] )
2127  )
2128  {
2129 
2130  movelay= "laydown";
2131  //calculate the new delta based on the space at the border:
2132  delt = abs(y_c[nvertl-1]-ynew[bottomright])/(nlays/2);
2133  delt_opp = abs(y_c[nvertl-1]-ynew[topright])/(nlays/2);
2134 cout<<"move layerdown to up:"<<endl;
2135 
2136  }
2137 
2138 
2139  //ChangeLayerspos(ynew, y_c, Addpointsy, layers_y, lay_Vids, delt, delt_opp, movelay, npedge);
2140 */
2141  //--------------------------------------------------------------
2142 
2143 
2144 
2145  //replace the vertices with the new ones
2146 
2147  Replacevertices(changefile, xnew , ynew, x_c, y_c, Addpointsx, Addpointsy, Eids, npedge, charalp, layers_x,layers_y, lay_eids, curv_lay);
2148 
2149 }
2150 
2151  // Define Expansion
2152  void SetFields(SpatialDomains::MeshGraphSharedPtr &mesh,
2153  SpatialDomains::BoundaryConditionsSharedPtr &boundaryConditions,
2154  LibUtilities::SessionReaderSharedPtr &session,
2155  Array<OneD,MultiRegions::ExpListSharedPtr> &Exp,int nvariables)
2156  {
2157 
2158  // Setting parameteres for homogenous problems
2159  MultiRegions::GlobalSysSolnType solnType;
2160  NekDouble LhomX; ///< physical length in X direction (if homogeneous)
2161  NekDouble LhomY; ///< physical length in Y direction (if homogeneous)
2162  NekDouble LhomZ; ///< physical length in Z direction (if homogeneous)
2163 
2164  bool DeclareCoeffPhysArrays = true;
2165  int npointsX; ///< number of points in X direction (if homogeneous)
2166  int npointsY; ///< number of points in Y direction (if homogeneous)
2167  int npointsZ; ///< number of points in Z direction (if homogeneous)
2168  int HomoDirec = 0;
2169  bool useFFT = false;
2170  ///Parameter for homogeneous expansions
2171  enum HomogeneousType
2172  {
2173  eHomogeneous1D,
2174  eHomogeneous2D,
2175  eHomogeneous3D,
2176  eNotHomogeneous
2177  };
2178 
2179  enum HomogeneousType HomogeneousType = eNotHomogeneous;
2180 
2181  if(session->DefinesSolverInfo("HOMOGENEOUS"))
2182  {
2183  std::string HomoStr = session->GetSolverInfo("HOMOGENEOUS");
2184  //m_spacedim = 3;
2185 
2186  if((HomoStr == "HOMOGENEOUS1D")||(HomoStr == "Homogeneous1D")||
2187  (HomoStr == "1D")||(HomoStr == "Homo1D"))
2188  {
2189  HomogeneousType = eHomogeneous1D;
2190  npointsZ = session->GetParameter("HomModesZ");
2191  LhomZ = session->GetParameter("LZ");
2192  HomoDirec = 1;
2193  }
2194 
2195  if((HomoStr == "HOMOGENEOUS2D")||(HomoStr == "Homogeneous2D")||
2196  (HomoStr == "2D")||(HomoStr == "Homo2D"))
2197  {
2198  HomogeneousType = eHomogeneous2D;
2199  npointsY = session->GetParameter("HomModesY");
2200  LhomY = session->GetParameter("LY");
2201  npointsZ = session->GetParameter("HomModesZ");
2202  LhomZ = session->GetParameter("LZ");
2203  HomoDirec = 2;
2204  }
2205 
2206  if((HomoStr == "HOMOGENEOUS3D")||(HomoStr == "Homogeneous3D")||
2207  (HomoStr == "3D")||(HomoStr == "Homo3D"))
2208  {
2209  HomogeneousType = eHomogeneous3D;
2210  npointsX = session->GetParameter("HomModesX");
2211  LhomX = session->GetParameter("LX");
2212  npointsY = session->GetParameter("HomModesY");
2213  LhomY = session->GetParameter("LY");
2214  npointsZ = session->GetParameter("HomModesZ");
2215  LhomZ = session->GetParameter("LZ");
2216  HomoDirec = 3;
2217  }
2218 
2219  if(session->DefinesSolverInfo("USEFFT"))
2220  {
2221  useFFT = true;
2222  }
2223  }
2224 
2225  int i;
2226  int expdim = mesh->GetMeshDimension();
2227  Exp= Array<OneD, MultiRegions::ExpListSharedPtr>(nvariables);
2228  // Continuous Galerkin projection
2229 
2230 
2231  switch(expdim)
2232  {
2233  case 1:
2234  {
2235 
2236 /*
2237  if(HomogeneousType == eHomogeneous2D)
2238  {
2239  const LibUtilities::PointsKey PkeyY(npointsY,LibUtilities::eFourierEvenlySpaced);
2240  const LibUtilities::BasisKey BkeyY(LibUtilities::eFourier,npointsY,PkeyY);
2241  const LibUtilities::PointsKey PkeyZ(npointsZ,LibUtilities::eFourierEvenlySpaced);
2242  const LibUtilities::BasisKey BkeyZ(LibUtilities::eFourier,npointsZ,PkeyZ);
2243 
2244  for(i = 0 ; i < Exp.num_elements(); i++)
2245  {
2246  Exp[i] = MemoryManager<MultiRegions::ContField3DHomogeneous2D>
2247  ::AllocateSharedPtr(session,BkeyY,BkeyZ,LhomY,LhomZ,useFFT,mesh,session->GetVariable(i));
2248  }
2249  }
2250  else
2251  {
2252 */
2253  for(i = 0 ; i < Exp.num_elements(); i++)
2254  {
2255  Exp[i] = MemoryManager<MultiRegions::ContField1D>
2256  ::AllocateSharedPtr(session,mesh,session->GetVariable(i));
2257  }
2258  //}
2259 
2260  break;
2261  }
2262  case 2:
2263  {
2264 /*
2265  if(HomogeneousType == eHomogeneous1D)
2266  {
2267  const LibUtilities::PointsKey PkeyZ(npointsZ,LibUtilities::eFourierEvenlySpaced);
2268  const LibUtilities::BasisKey BkeyZ(LibUtilities::eFourier,npointsZ,PkeyZ);
2269 
2270  for(i = 0 ; i < Exp.num_elements(); i++)
2271  {
2272  Exp[i] = MemoryManager<MultiRegions::ContField3DHomogeneous1D>
2273  ::AllocateSharedPtr(session,BkeyZ,LhomZ,useFFT,mesh,session->GetVariable(i));
2274  }
2275  }
2276  else
2277  {
2278 */
2279 
2280  i = 0;
2281  MultiRegions::ContField2DSharedPtr firstfield;
2282 
2283  firstfield = MemoryManager<MultiRegions::ContField2D>
2284  ::AllocateSharedPtr(session,mesh,session->GetVariable(i),DeclareCoeffPhysArrays);
2285 
2286  Exp[0] = firstfield;
2287  for(i = 1 ; i < Exp.num_elements(); i++)
2288  {
2289  Exp[i] = MemoryManager<MultiRegions::ContField2D>
2290  ::AllocateSharedPtr(*firstfield,mesh,session->GetVariable(i),DeclareCoeffPhysArrays);
2291  }
2292  // }
2293 
2294  break;
2295  }
2296  case 3:
2297  {
2298 /*
2299  if(HomogeneousType == eHomogeneous3D)
2300  {
2301  ASSERTL0(false,"3D fully periodic problems not implemented yet");
2302  }
2303  else
2304  {
2305 */
2306  i = 0;
2307  MultiRegions::ContField3DSharedPtr firstfield =
2308  MemoryManager<MultiRegions::ContField3D>
2309  ::AllocateSharedPtr(session,mesh,session->GetVariable(i));
2310 
2311  Exp[0] = firstfield;
2312  for(i = 1 ; i < Exp.num_elements(); i++)
2313  {
2314  Exp[i] = MemoryManager<MultiRegions::ContField3D>
2315  ::AllocateSharedPtr(*firstfield,mesh,session->GetVariable(i));
2316  }
2317  //}
2318  break;
2319  }
2320  default:
2321  ASSERTL0(false,"Expansion dimension not recognised");
2322  break;
2323  }
2324  }
2325 
2326  void OrderVertices (int nedges, SpatialDomains::MeshGraphSharedPtr graphShPt,
2327  MultiRegions::ExpListSharedPtr & bndfield,
2328  Array<OneD, int>& Vids, int v1,int v2, NekDouble x_connect, int & lastedge,
2329  Array<OneD, NekDouble>& x, Array<OneD, NekDouble>& y)
2330  {
2331  int nvertl = nedges+1;
2332  int edge;
2333  Array<OneD, int> Vids_temp(nvertl,-10);
2334  NekDouble x0layer = 0.0;
2335  for(int j=0; j<nedges; j++)
2336  {
2337  LocalRegions::SegExpSharedPtr bndSegExplow =
2338  bndfield->GetExp(j)->as<LocalRegions::SegExp>() ;
2339  edge = (bndSegExplow->GetGeom1D())->GetEid();
2340 //cout<<" edge="<<edge<<endl;
2341  for(int k=0; k<2; k++)
2342  {
2343  Vids_temp[j+k]=(bndSegExplow->GetGeom1D())->GetVid(k);
2344  SpatialDomains::PointGeomSharedPtr vertex = graphShPt->GetVertex(Vids_temp[j+k]);
2345  NekDouble x1,y1,z1;
2346  vertex->GetCoords(x1,y1,z1);
2347  if(x1==x_connect && edge!=lastedge)
2348  {
2349  //first 2 points
2350  if(x_connect==x0layer)
2351  {
2352  Vids[v1]=Vids_temp[j+k];
2353  x[v1]=x1;
2354  y[v1]=y1;
2355 
2356  if(k==0 )
2357  {
2358  Vids_temp[j+1]=(bndSegExplow->GetGeom1D())->GetVid(1);
2359  Vids[v2]=Vids_temp[j+1];
2360  SpatialDomains::PointGeomSharedPtr vertex = graphShPt->GetVertex(Vids[v2]);
2361  NekDouble x2,y2,z2;
2362  vertex->GetCoords(x2,y2,z2);
2363  x[v2]=x2;
2364  y[v2]=y2;
2365  }
2366  if(k==1)
2367  {
2368  Vids_temp[j+0]=(bndSegExplow->GetGeom1D())->GetVid(0);
2369  Vids[v2]=Vids_temp[j+0];
2370  SpatialDomains::PointGeomSharedPtr vertex = graphShPt->GetVertex(Vids[v2]);
2371  NekDouble x2,y2,z2;
2372  vertex->GetCoords(x2,y2,z2);
2373  x[v2]=x2;
2374  y[v2]=y2;
2375  }
2376 
2377  }
2378  else //other vertices
2379  {
2380  if(k==0 )
2381  {
2382  Vids_temp[j+1]=(bndSegExplow->GetGeom1D())->GetVid(1);
2383  Vids[v1]=Vids_temp[j+1];
2384  SpatialDomains::PointGeomSharedPtr vertex = graphShPt->GetVertex(Vids[v1]);
2385  NekDouble x1,y1,z1;
2386  vertex->GetCoords(x1,y1,z1);
2387  x[v1]=x1;
2388  y[v1]=y1;
2389  }
2390  if(k==1)
2391  {
2392  Vids_temp[j+0]=(bndSegExplow->GetGeom1D())->GetVid(0);
2393  Vids[v1]=Vids_temp[j+0];
2394  SpatialDomains::PointGeomSharedPtr vertex = graphShPt->GetVertex(Vids[v1]);
2395  NekDouble x1,y1,z1;
2396  vertex->GetCoords(x1,y1,z1);
2397  x[v1]=x1;
2398  y[v1]=y1;
2399  }
2400  }
2401  break;
2402  }
2403  }
2404  if(Vids[v1]!=-10)
2405  {
2406  lastedge = edge;
2407  break;
2408  }
2409  }
2410 
2411  }
2412 
2413  void Computestreakpositions(int npoints, MultiRegions::ExpListSharedPtr &streak,
2414  Array<OneD, NekDouble> x, Array<OneD, NekDouble> y,
2415  Array<OneD, NekDouble> &xold_up, Array<OneD, NekDouble> &yold_up,
2416  Array<OneD, NekDouble> &xold_low, Array<OneD, NekDouble> &yold_low,
2417  Array<OneD, NekDouble> &xold_c, Array<OneD, NekDouble> &yold_c,
2418  Array<OneD, NekDouble> &xc, Array<OneD, NekDouble> &yc, NekDouble cr,
2419  bool verts)
2420  {
2421 cout<<"Computestreakpositions"<<endl;
2422  int nq = streak->GetTotPoints();
2423  Array<OneD, NekDouble> coord(2);
2424  //Array<OneD, NekDouble> stvalues(nvertl,-10);
2425  Array<OneD, NekDouble> derstreak(nq);
2426  streak->PhysDeriv(MultiRegions::eY, streak->GetPhys(), derstreak);
2427  int elmtid, offset;
2428  NekDouble U,dU;
2429  NekDouble F=1000;
2430 
2431  if(verts==true)//only for verts makes sense to init the coord values..
2432  {
2433 /*
2434 //we can't rely on the layer down/up position to calc the streak position !!!(criti lay
2435 //EXTERNAL to the layer zone)
2436  for(int q=0; q<npoints; q++)
2437  {
2438  NekDouble streaktmp =200;
2439  NekDouble streaktmppos=200;
2440  NekDouble streaktmpneg=-200;
2441  int ipos,ineg, it =0;
2442  NekDouble weightpos, weightneg;
2443 
2444 cout<<"nq="<<nq<<" xold_down="<<xold_low[q]<<" xold_up="<<xold_up[q]<<" xold_c="<<xold_c[q]<<endl;
2445 
2446  ASSERTL0(xold_low[q]==xold_up[q], "layer region non valid");
2447 cout<<q<<" xup="<<xold_up[q]<<" yup="<<yold_up[q]<<" ydown="<<yold_low[q]<<endl;
2448  //algorithm to determine the closest points to the streak positions
2449  //using the weighted mean
2450  for(int j=0; j<nq; j++)
2451  {
2452  if(x[j]==xold_c[q] && y[j]<= yold_up[q] && y[j]>= yold_low[q])
2453  {
2454  if(streak->GetPhys()[j]< streaktmppos && streak->GetPhys()[j]>cr)
2455  {
2456  streaktmppos = streak->GetPhys()[j];
2457  ipos =j;
2458  }
2459  static int cnt=0;
2460 
2461  if( abs(streak->GetPhys()[j]) < -streaktmpneg && streak->GetPhys()[j]<cr)
2462  {
2463 //cout<<"test streakneg="<<streak->GetPhys()[j]<<endl;
2464  streaktmpneg = streak->GetPhys()[j];
2465  ineg =j;
2466  }
2467 
2468 //cout<<" x="<<x[j]<<" y="<<y[j]<<" streak="<<streak->GetPhys()[j]<<endl;
2469  }
2470  }
2471 cout<<"ipos="<<ipos<<" ineg="<<ineg<<endl;
2472 //cout<<"closer streak points ypos="<<y[ipos]<<" yneg="<<y[ineg]<<endl;
2473  //determine the streak y position as the result of the weighted mean
2474  if(streaktmppos< 200 && streaktmpneg>-200)
2475  {
2476  xc[q]= x[ipos];
2477  weightpos = 1/(streaktmppos*streaktmppos);
2478  weightneg = 1/(streaktmpneg*streaktmpneg);
2479  yc[q]= ( (y[ipos]*weightpos) + (y[ineg]*weightneg) )/(weightpos+weightneg);
2480  }
2481  else if(streaktmppos< 200)
2482  {
2483  xc[q]= x[ipos];
2484  yc[q]= y[ipos];
2485  }
2486  else if(streaktmpneg> -200)
2487  {
2488  xc[q]= x[ineg];
2489  yc[q]= y[ineg];
2490  }
2491  else
2492  {
2493  ASSERTL0(false, "streak not found");
2494  }
2495 
2496 cout<<" streak x="<<xc[q]<<" y="<<yc[q]<<endl;
2497  }
2498 */
2499  //start guess
2500  //yc= (yup+ydown)/2
2501  Vmath::Vadd(xc.num_elements(), yold_up,1,yold_low,1, yc,1);
2502  Vmath::Smul(xc.num_elements(), 0.5,yc,1,yc,1);
2503  Vmath::Vcopy(xc.num_elements(),xold_c,1,xc,1);
2504  }
2505  else//case of xQ,yQ
2506  {
2507  Vmath::Vcopy(xc.num_elements(), x,1,xc,1);
2508  Vmath::Vcopy(xc.num_elements(), y,1,yc,1);
2509  }
2510 
2511  int its;
2512  int attempt;
2513  NekDouble tol = 1e-3;
2514  NekDouble ytmp;
2515  for(int e=0; e<npoints; e++)
2516  {
2517  coord[0] =xc[e];
2518  coord[1] =yc[e];
2519 
2520  elmtid = streak->GetExpIndex(coord,0.00001);
2521  offset = streak->GetPhys_Offset(elmtid);
2522  F = 1000;
2523  its=0;
2524  attempt=0;
2525 //cout<<"start guess: x="<<xc[e]<<" y="<<yc[e]<<endl;
2526  while( abs(F)> 0.000000001)
2527  {
2528  ytmp = coord[1];
2529  elmtid = streak->GetExpIndex(coord,0.00001);
2530  offset = streak->GetPhys_Offset(elmtid);
2531  U = streak->GetExp(elmtid)->PhysEvaluate(coord, streak->GetPhys() + offset);
2532  dU = streak->GetExp(elmtid)->PhysEvaluate(coord, derstreak + offset);
2533  coord[1] = coord[1] - (U-cr)/dU;
2534  F = U-cr;
2535  ASSERTL0( coord[0]==xc[e], " x coordinate must remain the same");
2536  //stvalues[e] = streak->GetExp(elmtid)->PhysEvaluate(coord, streak->GetPhys() +offset );
2537 //cout<<"elmtid="<<elmtid<<" x="<<coord[0]<<" y="<<coord[1]<<" stvalue="<<U<<" dU="<<dU<<endl;
2538  if(abs(coord[1])>1
2539  && attempt==0 && verts==true
2540  )
2541  {
2542 
2543  //try the old crit lay position:
2544  coord[1] = yold_c[e];
2545  attempt++;
2546  }
2547  else if(abs(coord[1])>1 )
2548  {
2549  coord[1] = ytmp +0.01;
2550  elmtid = streak->GetExpIndex(coord,0.00001);
2551  offset = streak->GetPhys_Offset(elmtid);
2552  NekDouble Utmp = streak->GetExp(elmtid)->PhysEvaluate(coord, streak->GetPhys() + offset);
2553  NekDouble dUtmp = streak->GetExp(elmtid)->PhysEvaluate(coord, derstreak + offset);
2554  coord[1] = coord[1] - (Utmp-cr)/dUtmp;
2555 
2556  if( (abs(Utmp-cr)>abs(F))||(abs(coord[1])>1) )
2557  {
2558  coord[1] = ytmp -0.01;
2559  }
2560 
2561  attempt++;
2562  }
2563  else
2564  {
2565  ASSERTL0(abs(coord[1])<= 1, " y value out of bound +/-1");
2566  }
2567 
2568  its++;
2569  if(its>1000 && F< 0.0001)
2570  {
2571  cout<<"warning streak position obtained with precision:"<<F<<endl;
2572  break;
2573  }
2574  else if(its>1000)
2575  {
2576  ASSERTL0(false, "no convergence after 1000 iterations");
2577  }
2578  }
2579  yc[e] = coord[1] - (U-cr)/dU;
2580 
2581  ASSERTL0( U<= cr + tol, "streak wrong+");
2582  ASSERTL0( U>= cr -tol, "streak wrong-");
2583  //Utilities::Zerofunction(coord[0], coord[1], xtest, ytest, streak, derstreak);
2584 cout<<"result streakvert x="<<xc[e]<<" y="<<yc[e]<<" streak="<<U<<endl;
2585  }
2586 
2587 
2588  }
2589 
2590 
2591  void GenerateAddPointsNewtonIt( NekDouble &xi, NekDouble &yi,NekDouble &x0, NekDouble &y0,
2592  MultiRegions::ExpListSharedPtr &function, Array<OneD, NekDouble> &derfunction,
2593  NekDouble cr)
2594  {
2595  int elmtid,offset;
2596  NekDouble F,U,dU;
2597  Array<OneD, NekDouble> coords(2);
2598 
2599  coords[0] = xi;
2600  coords[1] = yi;
2601  F =1000;
2602  int attempt=0;
2603  int its=0;
2604  NekDouble ytmp;
2605  while( abs(F)> 0.00000001)
2606  {
2607  ytmp = coords[1];
2608 //cout<<"generate newton it xi="<<xi<<" yi="<<yi<<endl;
2609  elmtid = function->GetExpIndex(coords, 0.00001);
2610  //@to do if GetType(elmtid)==triangular WRONG!!!
2611 //cout<<"gen newton xi="<<xi<<" yi="<<coords[1]<<" elmtid="<<elmtid<<" F="<<F<<endl;
2612  offset = function->GetPhys_Offset(elmtid);
2613 
2614  U = function->GetExp(elmtid)->PhysEvaluate(coords, function->GetPhys() + offset);
2615  dU = function->GetExp(elmtid)->PhysEvaluate(coords, derfunction + offset);
2616  coords[1] = coords[1] - (U-cr)/dU;
2617 //cout<<cr<<"U-cr="<<U-cr<<" tmp result y:"<<coords[1]<<" dU="<<dU<<endl;
2618  F = U-cr;
2619 
2620  if( abs(coords[1])>1
2621  //&& attempt==0
2622  )
2623  {
2624 
2625  coords[1] = ytmp +0.01;
2626  elmtid = function->GetExpIndex(coords,0.00001);
2627  offset = function->GetPhys_Offset(elmtid);
2628  NekDouble Utmp = function->GetExp(elmtid)->PhysEvaluate(coords, function->GetPhys() + offset);
2629  NekDouble dUtmp = function->GetExp(elmtid)->PhysEvaluate(coords, derfunction + offset);
2630  coords[1] = coords[1] - (Utmp-cr)/dUtmp;
2631 cout<<"attempt:"<<coords[1]<<endl;
2632  if( (abs(Utmp-cr)>abs(F))||(abs(coords[1])>1) )
2633  {
2634  coords[1] = ytmp -0.01;
2635  }
2636 
2637  attempt++;
2638  }
2639  else
2640  {
2641  ASSERTL0(abs(coords[1])<= 1, " y value out of bound +/-1");
2642  }
2643 
2644  its++;
2645  if(its>1000 && F< 0.0001)
2646  {
2647  cout<<"warning streak position obtained with precision:"<<F<<endl;
2648  break;
2649  }
2650  else if(its>1000)
2651  {
2652  ASSERTL0(false, "no convergence after 1000 iterations");
2653  }
2654 
2655 
2656  ASSERTL0( coords[0]==xi, " x coordinate must remain the same");
2657  }
2658  x0 = xi;
2659  y0 = coords[1] - (U-cr)/dU;
2660 cout<<"NewtonIt result x="<<x0<<" y="<<coords[1]<<" U="<<U<<endl;
2661  }
2662 
2663  void GenerateCurve(int npoints, int npused, Array<OneD, NekDouble> &x_c,
2664  Array<OneD, NekDouble> &y_c, Array<OneD, NekDouble> &curve)
2665  {
2666  int N= npoints;
2667  int totpoints = npoints + npused;
2668 //cout<<"totppoints="<<totpoints<<endl;
2669  Array<OneD, NekDouble> A (N*N,1.0);
2670  Array<OneD, NekDouble> b (N);
2671  int row=0;
2672  //fill column by column
2673  for(int e=0; e<N; e++)
2674  {
2675  row=0;
2676  for(int w=npused; w < totpoints; w++)
2677  {
2678  A[N*e+row] = std::pow( x_c[w], N-1-e);
2679  row++;
2680  }
2681  }
2682  row=0;
2683  for(int r= npused; r< totpoints; r++)
2684  {
2685  b[row] = y_c[r];
2686 //cout<<"b="<<b[row]<<" y_c="<<y_c[r]<<endl;
2687  row++;
2688  }
2689 
2690 //cout<<"A elements="<<A.num_elements()<<endl;
2691  Array<OneD, int> ipivot (N);
2692  int info =0;
2693  //Lapack::Dgesv( N, 1, A.get(), N, ipivot.get(), b.get(), N, info);
2694  Lapack::Dgetrf( N, N, A.get(), N, ipivot.get(), info);
2695 
2696  if( info < 0 )
2697  {
2698  std::string message = "ERROR: The " + boost::lexical_cast<std::string>(-info) +
2699  "th parameter had an illegal parameter for dgetrf";
2700  ASSERTL0(false, message.c_str());
2701  }
2702  else if( info > 0 )
2703  {
2704  std::string message = "ERROR: Element u_" + boost::lexical_cast<std::string>(info) +
2705  boost::lexical_cast<std::string>(info) + " is 0 from dgetrf";
2706  ASSERTL0(false, message.c_str());
2707  }
2708 
2709  // N means no transponse (direct matrix)
2710  int ncolumns_b =1;
2711  Lapack::Dgetrs( 'N', N, ncolumns_b , A.get() , N, ipivot.get(), b.get(), N, info);
2712  if( info < 0 )
2713  {
2714  std::string message = "ERROR: The " + boost::lexical_cast<std::string>(-info) +
2715  "th parameter had an illegal parameter for dgetrf";
2716  ASSERTL0(false, message.c_str());
2717  }
2718  else if( info > 0 )
2719  {
2720  std::string message = "ERROR: Element u_" + boost::lexical_cast<std::string>(info) +
2721  boost::lexical_cast<std::string>(info) + " is 0 from dgetrf";
2722  ASSERTL0(false, message.c_str());
2723  }
2724 //cout<<"coeffs a,b,c,d..."<<endl;
2725  //fill the coeffs array
2726  if(npused==0)
2727  {
2728  Vmath::Vcopy(totpoints, &(b[0]), 1, &(curve[npused]), 1);
2729 //cout<<b[1]<<" ipiv="<<ipivot[1]<<endl;
2730  }
2731  else
2732  {
2733  row=0;
2734  for(int a= npused+1; a< totpoints+1; a++)
2735  {
2736  curve[a] = b[row];
2737  row++;
2738 //cout<<curve[a]<<" ipiv="<<ipivot[row-1]<<endl;
2739  }
2740  }
2741 
2742 
2743 
2744  }
2745 
2746 
2747  void GenerateCurveSp(Array<OneD, NekDouble> &x_c, Array<OneD, NekDouble> &y_c,
2748  Array<OneD, NekDouble> &x_lay, Array<OneD, NekDouble> &y_lay)
2749  {
2750  Array<OneD, NekDouble> k(x_c.num_elements());
2751  Array<OneD, NekDouble> b(x_c.num_elements());
2752  NekDouble dx,dx1,dy,dy1;
2753  //determine the values of k
2754  //remark fix the derivative=0 at the border of the domain: k_0=k_n=0
2755  for(int s=0; s< k.num_elements(); s++)
2756  {
2757  dx = x_c[s] - x_c[s-1];
2758  dy = y_c[s] - y_c[s-1];
2759  dx1 = x_c[s+1] - x_c[s];
2760  dy1 = y_c[s+1] - y_c[s];
2761  if(s==0 || s==k.num_elements()-1)
2762  {
2763  k[s]=0;
2764 
2765  }
2766  else
2767  {
2768  b[s] = 3*( dy/(dx*dx) +dy1/(dx1*dx1) );
2769  }
2770  }
2771  }
2772 
2773 
2774 
2775 
2776  void GenerateAddPoints(int region, SpatialDomains::MeshGraphSharedPtr &mesh,int np, int npused,
2777  Array<OneD, NekDouble> &curve, MultiRegions::ExpListSharedPtr & bndfield,
2778  Array<OneD, NekDouble>& outx, Array<OneD, NekDouble>& outy, Array<OneD, int>&Eids)
2779  {
2780  LocalRegions::SegExpSharedPtr bndSegExp;
2781 
2782  int Eid,id1,id2;
2783  NekDouble x1,y1,z1;
2784  NekDouble x2,y2,z2;
2785  SpatialDomains::PointGeomSharedPtr vertex1;
2786  SpatialDomains::PointGeomSharedPtr vertex2;
2787 
2788  int firstedge=0;
2789  int firstcoeff=0;
2790  if(npused!=0)
2791  {
2792  firstedge = npused;
2793  firstcoeff = firstedge+1;
2794 
2795  }
2796  int lastedge = firstedge + (np-1);
2797  int lastcoeff = firstcoeff + np;
2798 //cout<<"firstedge="<<firstedge<<" lastedge="<<lastedge<<" firstcoeff="<<firstcoeff
2799 // <<" lastcoeff="<<lastcoeff<<endl;
2800  int polorder;
2801  for(int s= firstedge; s< lastedge; s++)
2802  {
2803  bndSegExp = bndfield->GetExp(s)->as<LocalRegions::SegExp>();
2804  Eid = (bndSegExp->GetGeom1D())->GetEid();
2805  id1 = (bndSegExp->GetGeom1D())->GetVid(0);
2806  id2 = (bndSegExp->GetGeom1D())->GetVid(1);
2807  vertex1 = mesh->GetVertex(id1);
2808  vertex2 = mesh->GetVertex(id2);
2809  vertex1->GetCoords(x1,y1,z1);
2810  vertex2->GetCoords(x2,y2,z2);
2811 //cout<<"x1="<<x1<<" x2="<<x2<<endl;
2812  if(x2>x1)
2813  {
2814  outx[s] = x1 +(x2-x1)/2;
2815  if( outx[s]>x2 || outx[s]< x1)
2816  {
2817  outx[s] = -outx[s];
2818  }
2819  }
2820  else if(x1>x2)
2821  {
2822  outx[s] = x2+ (x1-x2)/2;
2823  if( outx[s] > x1 || outx[s] <x2)
2824  {
2825  outx[s] = -outx[s];
2826  }
2827  }
2828  else
2829  {
2830  ASSERTL0(false, "point not generated");
2831  }
2832  polorder = np-1;
2833  for(int g= firstcoeff; g< lastcoeff; g++)
2834  {
2835  outy[s] += curve[g]*(std::pow(outx[s], polorder));
2836 //cout<<"coeff*x^i="<<outy[s]<<" coeff="<<curve[g]<<" exp="<<polorder<<endl;
2837  ASSERTL0(polorder >=0, " polynomial with one negative exponent");
2838  polorder--;
2839  }
2840  Eids[s] = Eid;
2841 
2842 //cout<<"eid="<<Eids[s]<<" xadd="<<outx[s]<<" yadd="<<outy[s]<<endl;
2843  }
2844 
2845  }
2846 
2847  void MapEdgeVertices(MultiRegions::ExpListSharedPtr field, Array<OneD, int> &V1,
2848  Array<OneD, int> &V2)
2849  {
2850  const boost::shared_ptr<LocalRegions::ExpansionVector> exp2D = field->GetExp();
2851  int nel = exp2D->size();
2852  LocalRegions::QuadExpSharedPtr locQuadExp;
2853  LocalRegions::TriExpSharedPtr locTriExp;
2854  SpatialDomains::Geometry1DSharedPtr SegGeom;
2855  int id;
2856  int cnt=0;
2857  Array<OneD, int> V1tmp(4*nel, 10000);
2858  Array<OneD, int> V2tmp(4*nel, 10000);
2859  for(int i=0; i<nel; i++)
2860  {
2861  if((locQuadExp = (*exp2D)[i]->as<LocalRegions::QuadExp>()))
2862  {
2863  for(int j = 0; j < locQuadExp->GetNedges(); ++j)
2864  {
2865  SegGeom = (locQuadExp->GetGeom2D())->GetEdge(j);
2866  id = SegGeom->GetEid();
2867  if( V1tmp[id] == 10000)
2868  {
2869  V1tmp[id]= SegGeom->GetVertex(0)->GetVid();
2870  V2tmp[id]= SegGeom->GetVertex(1)->GetVid();
2871  cnt++;
2872  }
2873  }
2874  }
2875  //in the future the tri edges may be not necessary (if the nedges is known)
2876 
2877  else if((locTriExp = (*exp2D)[i]->as<LocalRegions::TriExp>()))
2878  {
2879  for(int j = 0; j < locTriExp->GetNedges(); ++j)
2880  {
2881  SegGeom = (locTriExp->GetGeom2D())->GetEdge(j);
2882  id = SegGeom->GetEid();
2883 
2884  if( V1tmp[id] == 10000)
2885  {
2886  V1tmp[id]= SegGeom->GetVertex(0)->GetVid();
2887  V2tmp[id]= SegGeom->GetVertex(1)->GetVid();
2888  cnt++;
2889  }
2890  }
2891 
2892  }
2893 
2894  }
2895 
2896  V1 = Array<OneD, int>(cnt);
2897  V2 = Array<OneD, int>(cnt);
2898  //populate the output vectors V1,V2
2899  for(int g=0; g<cnt; g++)
2900  {
2901  V1[g] = V1tmp[g];
2902  ASSERTL0(V1tmp[g]!=10000, "V1 wrong");
2903  V2[g] = V2tmp[g];
2904  ASSERTL0(V2tmp[g]!=10000, "V2 wrong");
2905  }
2906 
2907  }
2908 
2909  void Findlay_eids(Array<OneD, Array<OneD, NekDouble> >& layers_y,
2910  Array<OneD, int> V1, Array<OneD, int> V2,
2911  Array<OneD, Array<OneD, int > >& lay_Vids,
2912  Array<OneD, Array<OneD, int > >& lay_eids)
2913  {
2914  int neids = V1.num_elements();
2915  int nlays = layers_y.num_elements();
2916  int nvertl = lay_Vids[0].num_elements();
2917  Array<OneD, NekDouble > tmpy(nlays,-1000);
2918  Array<OneD, int > tmpVids(nlays,-1000);
2919  Array<OneD, int > tmpIndex(nlays,-1000);
2920 
2921  //sort layers_y low to high
2922  for(int a=0; a<nvertl; a++)
2923  {
2924  tmpy = Array<OneD, NekDouble> (nlays,-1000);
2925  //copy into array:(not possible use Vmath with matrices)
2926  for(int t=0; t<nlays; t++)
2927  {
2928  tmpy[t] = layers_y[t][a];
2929  tmpVids[t] = lay_Vids[t][a];
2930  }
2931 
2932 
2933  for(int b=0; b<nlays; b++)
2934  {
2935  tmpIndex[b] = Vmath::Imin(nlays,tmpy,1);
2936  layers_y[b][a] = tmpy[tmpIndex[b]];
2937  lay_Vids[b][a] = tmpVids[tmpIndex[b]];
2938  //tmpy[b] = Vmath::Vmin(nlays,&(layers_y[0][a]),1);
2939  ASSERTL0(layers_y[b][a]==Vmath::Vmin(nlays,tmpy,1),
2940  "sort failed ");
2941  tmpy[tmpIndex[b]] = 1000;
2942 
2943  }
2944 
2945  }
2946 
2947  for(int a=0; a< nvertl-1; a++)
2948  {
2949  for(int b=0; b<nlays; b++)
2950  {
2951  for(int p=0; p< neids; p++)
2952  {
2953  if(
2954  (lay_Vids[b][a] == V1[p] && lay_Vids[b][a+1]==V2[p])||
2955  (lay_Vids[b][a] == V2[p] && lay_Vids[b][a+1]==V1[p])
2956  )
2957  {
2958  lay_eids[b][a] = p;
2959  break;
2960  }
2961  }
2962  }
2963  }
2964 
2965  }
2966 
2967  void MoveLayersvertically(int nlays, int nvertl, int cntlow, int cntup,
2968  Array<OneD, Array<OneD, int > > lay_Vids, Array<OneD, NekDouble> xc,
2969  Array<OneD, NekDouble> yc, Array<OneD, int> Down, Array<OneD, int> Up,
2970  Array<OneD, NekDouble >& xnew,Array<OneD, NekDouble>& ynew,
2971  Array<OneD, Array<OneD, NekDouble > >& layers_x,
2972  Array<OneD, Array<OneD, NekDouble > >& layers_y)
2973  {
2974  int np_lay = layers_y[0].num_elements();
2975  // 0<h<nlays-1 to fill only the 'internal' layers(no up,low);
2976  for(int h=1; h<nlays-1; h++)
2977  {
2978  layers_x[h]= Array<OneD, NekDouble>(np_lay);
2979  for(int s=0; s<nvertl; s++)
2980  {
2981  //check if ynew is still empty
2982  ASSERTL0(ynew[ lay_Vids[h][s] ]==-20, "ynew layers not empty");
2983  if(h<cntlow+1)//layers under the crit lay
2984  {
2985  //y= ylow+delta
2986  ynew[ lay_Vids[h][s] ] = ynew[Down[s]]+ h*abs(ynew[Down[s]] - yc[s])/(cntlow+1);
2987  //put the layer vertical
2988  xnew[lay_Vids[h][s] ] = xc[s];
2989 //cout<<"ynew="<<ynew[ lay_Vids[h][s] ]<<" ydown="<<ynew[Down[s]]<<
2990 //" delta="<<abs(ynew[Down[s]] - y_c[s])/(cntlow+1)<<endl;
2991  //until now layers_y=yold
2992  layers_y[h][s] = ynew[ lay_Vids[h][s] ];
2993  layers_x[h][s] = xnew[ lay_Vids[h][s] ];
2994  }
2995  else
2996  {
2997  //y = yc+delta
2998  ynew[ lay_Vids[h][s] ] = yc[s] + (h-cntlow)*abs(ynew[Up[s]] - yc[s])/(cntup+1);
2999  //put the layer vertical
3000  xnew[lay_Vids[h][s] ] = xc[s];
3001  //until now layers_y=yold
3002  layers_y[h][s] = ynew[ lay_Vids[h][s] ];
3003  layers_x[h][s] = xnew[ lay_Vids[h][s] ];
3004  }
3005  }
3006  }
3007 
3008  }
3009  void Cutrepetitions(int nedges,Array<OneD, NekDouble> inarray,
3010  Array<OneD, NekDouble>& outarray)
3011  {
3012 
3013  //determine npedge:
3014  int np_lay = inarray.num_elements();
3015 
3016  int npedge = np_lay/nedges;
3017  ASSERTL0(inarray.num_elements()%nedges==0," something on number npedge");
3018  //cut out the repetitions:
3019 
3020  cout<<nedges<<"nedges"<<npedge<<" nplay="<<np_lay<<endl;
3021  //generate x,y arrays without lastedgepoint
3022  //(needed to interp correctly)
3023 
3024  for(int b=0; b<nedges; b++)
3025  {
3026  Vmath::Vcopy( npedge-1, &inarray[b*(npedge)],1,&outarray[b*(npedge-1)],1);
3027  if(b== nedges-1)
3028  {
3029  outarray[b*(npedge-1)+npedge-1] = inarray[b*npedge+npedge-1];
3030 
3031  }
3032  }
3033 
3034 
3035 
3036  }
3037 
3038  void Orderfunctionx(Array<OneD, NekDouble> inarray_x,
3039  Array<OneD, NekDouble> inarray_y, Array<OneD, NekDouble>& outarray_x,
3040  Array<OneD, NekDouble>& outarray_y)
3041  {
3042 
3043  Array<OneD, NekDouble>tmpx(inarray_x.num_elements());
3044  //local copy to prevent overwriting
3045  Vmath::Vcopy(inarray_x.num_elements() , inarray_x,1,tmpx,1);
3046 
3047  //order function with respect to x
3048  int index;
3049 
3050  NekDouble max = Vmath::Vmax(tmpx.num_elements(), tmpx,1);
3051  for(int w=0; w<tmpx.num_elements(); w++)
3052  {
3053  index = Vmath::Imin(tmpx.num_elements(), tmpx,1);
3054  outarray_x[w]= tmpx[index];
3055  outarray_y[w]= inarray_y[index];
3056  tmpx[index] = max+1000;
3057  }
3058 
3059 
3060  }
3061 
3062  int DetermineclosePointxindex(NekDouble x,Array<OneD, NekDouble> xArray)
3063  {
3064  int npts = xArray.num_elements();
3065  Array<OneD, NekDouble> xcopy(npts);
3066  Vmath::Vcopy(npts,xArray,1,xcopy,1);
3067  //subtract xpoint and abs
3068 
3069  Vmath::Sadd(npts, -x, xcopy,1, xcopy,1);
3070  Vmath::Vabs(npts, xcopy,1,xcopy,1);
3071 
3072  int index = Vmath::Imin(npts, xcopy,1);
3073  if(xArray[index]> x)//assume always x[index]< x(HYPHOTHESIS)
3074  {
3075  index = index-1;
3076  }
3077 
3078  return index;
3079  }
3080 
3081  void GenerateNeighbourArrays(int index,int neighpoints, Array<OneD, NekDouble> xArray,
3082  Array<OneD, NekDouble> yArray,Array<OneD, NekDouble>& Neighbour_x,
3083  Array<OneD, NekDouble>& Neighbour_y)
3084  {
3085  ASSERTL0( neighpoints%2==0,"number of neighbour points should be even");
3086  int leftpoints = (neighpoints/2)-1;//-1 because xArray[index]< x
3087  int rightpoints = neighpoints/2;
3088  int diff ;
3089  int start;
3090 //cout<<"index="<<index<<" left="<<leftpoints<<" right="<<rightpoints<<endl;
3091  if(index-leftpoints<0)
3092  {
3093 //cout<<"case0"<<endl;
3094  diff = index-leftpoints;
3095  start= 0;
3096  Vmath::Vcopy(neighpoints, &yArray[0],1,&Neighbour_y[0],1);
3097  Vmath::Vcopy(neighpoints, &xArray[0],1,&Neighbour_x[0],1);
3098  }
3099  else if( (yArray.num_elements()-1)-index < rightpoints)
3100  {
3101 //cout<<"case1 closest="<<xArray[index]<<endl;
3102  int rpoints = (yArray.num_elements()-1)-index;//
3103  diff = rightpoints-rpoints;
3104 //cout<<"start index="<<index-leftpoints-diff<<endl;
3105  start = index-leftpoints-diff;
3106  Vmath::Vcopy(neighpoints, &yArray[start],1,&Neighbour_y[0],1);
3107  Vmath::Vcopy(neighpoints, &xArray[start],1,&Neighbour_x[0],1);
3108  }
3109  else
3110  {
3111 //cout<<"caseaa"<<endl;
3112  start = index-leftpoints;
3113  Vmath::Vcopy(neighpoints, &yArray[start],1,&Neighbour_y[0],1);
3114  Vmath::Vcopy(neighpoints, &xArray[start],1,&Neighbour_x[0],1);
3115  }
3116 /*
3117 for(int t= start; t<start+neighpoints; t++)
3118 {
3119 cout<<"Px="<<xArray[t]<<" "<<yArray[t]<<endl;
3120 }
3121 */
3122  //check if there is any repetition
3123  for(int f=1; f< neighpoints; f++)
3124  {
3125  ASSERTL0(Neighbour_x[f]!=Neighbour_x[f-1]," repetition on NeighbourArrays");
3126  }
3127 
3128 
3129  }
3130 
3131  NekDouble LagrangeInterpolant(NekDouble x, int npts,
3132  Array<OneD,NekDouble> xpts, Array<OneD, NekDouble> funcvals)
3133  {
3134  NekDouble sum = 0.0;
3135  NekDouble LagrangePoly;
3136 //cout<<"lagrange"<<endl;
3137  for(int pt=0;pt<npts;++pt)
3138  {
3139  NekDouble h=1.0;
3140 
3141  for(int j=0;j<pt; ++j)
3142  {
3143  h = h * (x - xpts[j])/(xpts[pt]-xpts[j]);
3144  }
3145 
3146  for(int k=pt+1;k<npts;++k)
3147  {
3148  h = h * (x - xpts[k])/(xpts[pt]-xpts[k]);
3149  }
3150  LagrangePoly=h;
3151 
3152  sum += funcvals[pt]*LagrangePoly;
3153  }
3154 //cout<<"result :"<<sum<<endl;
3155  return sum;
3156  }
3157 
3158  void ChangeLayerspos(Array<OneD, NekDouble> & ynew,
3159  Array<OneD, NekDouble> yc,
3160  Array<OneD, NekDouble> Addpointsy,
3161  Array<OneD, Array<OneD, NekDouble > >& layers_y,
3162  Array<OneD, Array<OneD, int> >lay_Vids,
3163  NekDouble delt, NekDouble delt_opp,string movelay, int npedge)
3164  {
3165  int nvertl = yc.num_elements();
3166  int nedges = nvertl-1;
3167  int nlays = layers_y.num_elements();
3168  int np_lay = layers_y[0].num_elements();
3169  Array<OneD, Array<OneD, NekDouble> > tmpy(layers_y.num_elements());
3170  Array<OneD, Array<OneD, NekDouble> > tmpynew(layers_y.num_elements());
3171  int cntup=1;
3172  NekDouble delta;
3173  //cp layers_y into tmpy and move the first one
3174 cout<<"movelay"<<movelay<<endl;
3175  if(movelay == "laydown")
3176  {
3177  delta = abs(layers_y[nlays-1][(nvertl-2)*npedge+npedge-1]-yc[nvertl-1]);
3178 
3179  for(int g=0; g< nlays; g++)
3180  {
3181  tmpy[g] = Array<OneD, NekDouble> (np_lay);
3182  tmpynew[g] = Array<OneD, NekDouble>(nvertl);
3183  for(int h=0; h< nvertl; h++)
3184  {
3185 
3186  //move layers_up
3187  if(g==0 )
3188  {
3189  tmpynew[g][h]= yc[h] + delta + delta/(nlays/2);
3190 
3191  if(h<nedges)
3192  {
3193 
3194  //v1
3195  tmpy[g][h*npedge +0] = yc[h] + delta+ delta/(nlays/2);
3196  //v2
3197  tmpy[g][h*npedge +npedge-1] = yc[h+1] + delta+delta/(nlays/2);
3198  //middle points (shift crit lay points by delta):
3199  for(int d=0; d< npedge-2; d++)
3200  {
3201  tmpy[g][h*npedge +d+1]=
3202  Addpointsy[h*(npedge-2) +d] +delta + delta/(nlays/2);
3203 
3204  }
3205  }
3206  }
3207  else
3208  {
3209  tmpynew[g][h]= ynew[ lay_Vids[g][h] ];
3210  if(h<nedges)
3211  {
3212  //v1
3213  tmpy[g][h*npedge +0] = layers_y[g][h*npedge +0];
3214  //v2
3215  tmpy[g][h*npedge +npedge-1] = layers_y[g][h*npedge +npedge-1];
3216  //middle points (shift crit lay points by delta):
3217  for(int d=0; d< npedge-2; d++)
3218  {
3219  tmpy[g][h*npedge +d+1]= layers_y[g][h*npedge +d+1];
3220  }
3221  }
3222  }
3223 
3224  }
3225  }
3226 
3227  //copy back in a different order
3228  int cnt=-1;
3229  int vcnt=-1;
3230  for(int g=0; g< nlays; g++)
3231  {
3232 cout<<"move y coords"<<g<<endl;
3233  for(int h=0; h<np_lay; h++)
3234  {
3235  //first becomes last
3236  if(g==0)
3237  {
3238  layers_y[nlays-1][h]= tmpy[0][h];
3239  }
3240  else
3241  {
3242  layers_y[cnt][h] = tmpy[g][h];
3243  }
3244  }
3245  cnt++;
3246 
3247  for(int v=0; v<nvertl; v++)
3248  {
3249  if(g==0)
3250  {
3251  ynew[lay_Vids[nlays-1][v]] = tmpynew[0][v];
3252  }
3253  else
3254  {
3255  ynew[lay_Vids[vcnt][v] ] = tmpynew[g][v];
3256  }
3257  }
3258  vcnt++;
3259 
3260 
3261 
3262  }
3263  }
3264  if( movelay =="layup")
3265  {
3266  ASSERTL0(false, "not implemented layup");
3267  }
3268 
3269 
3270 
3271 
3272  }
3273 
3274  void Replacevertices(string filename, Array<OneD, NekDouble> newx,
3275  Array<OneD, NekDouble> newy,
3276  Array<OneD, NekDouble> x_crit, Array<OneD, NekDouble> y_crit,
3277  Array<OneD, NekDouble> & Pcurvx,
3278  Array<OneD, NekDouble> & Pcurvy,
3279  Array<OneD, int>&Eids, int Npoints, string s_alp,
3280  Array<OneD, Array<OneD, NekDouble> > x_lay,
3281  Array<OneD, Array<OneD, NekDouble> > y_lay,
3282  Array<OneD, Array<OneD, int > >lay_eids, bool curv_lay)
3283  {
3284  //load existing file
3285  string newfile;
3286  TiXmlDocument doc(filename);
3287  bool loadOkay = doc.LoadFile();
3288  //load xscale parameter (if exists)
3289  TiXmlElement* master = doc.FirstChildElement("NEKTAR");
3290  TiXmlElement* mesh = master->FirstChildElement("GEOMETRY");
3291  TiXmlElement* element = mesh->FirstChildElement("VERTEX");
3292  NekDouble xscale = 1.0;
3293  LibUtilities::AnalyticExpressionEvaluator expEvaluator;
3294  const char *xscal = element->Attribute("XSCALE");
3295  if(xscal)
3296  {
3297  std::string xscalstr = xscal;
3298  int expr_id = expEvaluator.DefineFunction("",xscalstr);
3299  xscale = expEvaluator.Evaluate(expr_id);
3300  }
3301 
3302 
3303  // Save a new XML file.
3304  newfile = filename.substr(0, filename.find_last_of("."))+"_moved.xml";
3305 
3306  doc.SaveFile( newfile );
3307 
3308  //write the new vertices
3309  TiXmlDocument docnew(newfile);
3310  bool loadOkaynew = docnew.LoadFile();
3311 
3312  std::string errstr = "Unable to load file: ";
3313  errstr += newfile;
3314  ASSERTL0(loadOkaynew, errstr.c_str());
3315 
3316  TiXmlHandle docHandlenew(&docnew);
3317  TiXmlNode* nodenew = NULL;
3318  TiXmlElement* meshnew = NULL;
3319  TiXmlElement* masternew = NULL;
3320  TiXmlElement* condnew = NULL;
3321  TiXmlElement* Parsnew = NULL;
3322  TiXmlElement* parnew = NULL;
3323 
3324  // Master tag within which all data is contained.
3325 
3326 
3327  masternew = docnew.FirstChildElement("NEKTAR");
3328  ASSERTL0(masternew, "Unable to find NEKTAR tag in file.");
3329 
3330  //set the alpha value
3331  string alphastring;
3332  condnew = masternew->FirstChildElement("CONDITIONS");
3333  Parsnew = condnew->FirstChildElement("PARAMETERS");
3334 cout<<"alpha="<<s_alp<<endl;
3335  parnew = Parsnew->FirstChildElement("P");
3336  while(parnew)
3337  {
3338  TiXmlNode *node = parnew->FirstChild();
3339  if (node)
3340  {
3341  // Format is "paramName = value"
3342  std::string line = node->ToText()->Value();
3343  std::string lhs;
3344  std::string rhs;
3345  /// Pull out lhs and rhs and eliminate any spaces.
3346  int beg = line.find_first_not_of(" ");
3347  int end = line.find_first_of("=");
3348  // Check for no parameter name
3349  if (beg == end) throw 1;
3350  // Check for no parameter value
3351  if (end != line.find_last_of("=")) throw 1;
3352  // Check for no equals sign
3353  if (end == std::string::npos) throw 1;
3354 
3355  lhs = line.substr(line.find_first_not_of(" "), end-beg);
3356  lhs = lhs.substr(0, lhs.find_last_not_of(" ")+1);
3357 
3358  //rhs = line.substr(line.find_last_of("=")+1);
3359  //rhs = rhs.substr(rhs.find_first_not_of(" "));
3360  //rhs = rhs.substr(0, rhs.find_last_not_of(" ")+1);
3361 
3362  boost::to_upper(lhs);
3363  if(lhs == "ALPHA")
3364  {
3365  alphastring = "Alpha = "+ s_alp;
3366  parnew->RemoveChild(node);
3367  parnew->LinkEndChild(new TiXmlText(alphastring) );
3368  }
3369  }
3370 
3371  parnew = parnew->NextSiblingElement("P");
3372  }
3373 
3374 
3375  // Find the Mesh tag and same the dim and space attributes
3376  meshnew = masternew->FirstChildElement("GEOMETRY");
3377 
3378  ASSERTL0(meshnew, "Unable to find GEOMETRY tag in file.");
3379  TiXmlAttribute *attrnew = meshnew->FirstAttribute();
3380  // Now read the vertices
3381  TiXmlElement* elementnew = meshnew->FirstChildElement("VERTEX");
3382  ASSERTL0(elementnew, "Unable to find mesh VERTEX tag in file.");
3383  //set xscale 1!!
3384  if(xscale!=1.0)
3385  {
3386  elementnew->SetAttribute("XSCALE",1.0);
3387  }
3388  TiXmlElement *vertexnew = elementnew->FirstChildElement("V");
3389 
3390 
3391 
3392  int indx;
3393  int err, numPts;
3394  int nextVertexNumber = -1;
3395 
3396  while (vertexnew)
3397  {
3398  nextVertexNumber++;
3399  //delete the old one
3400  TiXmlAttribute *vertexAttr = vertexnew->FirstAttribute();
3401  std::string attrName(vertexAttr->Name());
3402  ASSERTL0(attrName == "ID", (std::string("Unknown attribute name: ") + attrName).c_str());
3403 
3404  err = vertexAttr->QueryIntValue(&indx);
3405  ASSERTL0(err == TIXML_SUCCESS, "Unable to read attribute ID.");
3406  ASSERTL0(indx == nextVertexNumber, "Vertex IDs must begin with zero and be sequential.");
3407 
3408  std::string vertexBodyStr;
3409  // Now read body of vertex
3410  TiXmlNode *vertexBody = vertexnew->FirstChild();
3411  // Accumulate all non-comment body data.
3412  if (vertexBody->Type() == TiXmlNode::TEXT)
3413  {
3414  vertexBodyStr += vertexBody->ToText()->Value();
3415  vertexBodyStr += " ";
3416  }
3417  ASSERTL0(!vertexBodyStr.empty(), "Vertex definitions must contain vertex data.");
3418  //remove the old coordinates
3419  vertexnew->RemoveChild(vertexBody);
3420  //write the new one
3421 //cout<<"writing.. v:"<<nextVertexNumber<<endl;
3422  stringstream s;
3423  //we need at least 5 digits (setprecision 5) to get the streak position with
3424  // precision 10^-10
3425  s << std::scientific << std::setprecision(8) << newx[nextVertexNumber] << " "
3426  << newy[nextVertexNumber] << " " << 0.0;
3427  vertexnew->LinkEndChild(new TiXmlText(s.str()));
3428  //TiXmlNode *newvertexBody = vertexnew->FirstChild();
3429  //string newvertexbodystr= newvertexBody->SetValue(s.str());
3430  //vertexnew->ReplaceChild(vertexBody,new TiXmlText(newvertexbodystr));
3431 
3432  vertexnew = vertexnew->NextSiblingElement("V");
3433  }
3434 
3435 
3436 
3437  //read the curved tag
3438  TiXmlElement* curvednew = meshnew->FirstChildElement("CURVED");
3439  ASSERTL0(curvednew, "Unable to find mesh CURVED tag in file.");
3440  TiXmlElement *edgenew = curvednew->FirstChildElement("E");
3441  int cnt =-1;
3442  //ID is different from index...
3443  std::string charindex;
3444  int eid;
3445  int index;
3446  int indexeid, v1,v2;
3447  NekDouble x1,x2,y1,y2;
3448  //if edgenew belongs to the crit lay replace it, else delete it.
3449  while (edgenew)
3450  {
3451  indexeid =-1;
3452  cnt++;
3453  //get the index...
3454  TiXmlAttribute *edgeAttr = edgenew->FirstAttribute();
3455  std::string attrName(edgeAttr->Name());
3456  charindex = edgeAttr->Value();
3457  std::istringstream iss(charindex);
3458  iss >> std::dec >> index;
3459  //get the eid
3460  edgenew->QueryIntAttribute("EDGEID", &eid);
3461 //cout<<"eid="<<eid<<" neid="<<Eids.num_elements()<<endl;
3462  //find the corresponding index curve point
3463  for(int u=0; u<Eids.num_elements(); u++)
3464  {
3465 //cout<<"Eids="<<Eids[u]<<" eid="<<eid<<endl;
3466  if(Eids[u]==eid)
3467  {
3468  indexeid = u;
3469  }
3470  }
3471  if(indexeid==-1)
3472  {
3473  curvednew->RemoveChild(edgenew);
3474  //ASSERTL0(false, "edge to update not found");
3475  }
3476  else
3477  {
3478 
3479 
3480  std::string edgeBodyStr;
3481  //read the body of the edge
3482  TiXmlNode *edgeBody = edgenew->FirstChild();
3483  if(edgeBody->Type() == TiXmlNode::TEXT)
3484  {
3485  edgeBodyStr += edgeBody->ToText()->Value();
3486  edgeBodyStr += " ";
3487  }
3488  ASSERTL0(!edgeBodyStr.empty(), "Edge definitions must contain edge data");
3489  //remove the old coordinates
3490  edgenew->RemoveChild(edgeBody);
3491  //write the new points coordinates
3492  if(x_crit[cnt]< x_crit[cnt+1])
3493  {
3494  //x1 = x_crit[cnt];
3495  //x2 = x_crit[cnt+1];
3496  //y1 = y_crit[cnt];
3497  //y2 = y_crit[cnt+1];
3498  v1 = cnt;
3499  v2 = cnt+1;
3500  }
3501  else
3502  {
3503  //x1 = x_crit[cnt+1];
3504  //x2 = x_crit[cnt];
3505  //y1 = y_crit[cnt+1];
3506  //y2 = y_crit[cnt];
3507  v1 = cnt+1;
3508  v2 = cnt;
3509  cout<<"Warning: the edges can have the vertices in the reverse order";
3510  }
3511  //we need at least 5 digits (setprecision 5) to get the streak position with
3512  // precision 10^-10
3513 
3514  //Determine the number of points
3515  err = edgenew->QueryIntAttribute("NUMPOINTS", &numPts);
3516  ASSERTL0(err == TIXML_SUCCESS, "Unable to read curve attribute NUMPOINTS.");
3517  edgenew->SetAttribute("NUMPOINTS", Npoints);
3518  stringstream st;
3519  st << std::scientific << std::setprecision(8) << x_crit[v1] << " "
3520  << y_crit[v1] << " " << 0.000<<" ";
3521 //cout<<x_crit[v1]<<" "<<y_crit[v1]<<endl;
3522  for(int a=0; a< Npoints-2; a++)
3523  {
3524  st << std::scientific << std::setprecision(8) <<
3525  " "<<Pcurvx[indexeid*(Npoints-2) +a]<<" "<<Pcurvy[indexeid*(Npoints-2) +a]
3526  <<" "<<0.000<<" ";
3527 //cout<<Pcurvx[indexeid*(Npoints-2) +a]<<" "<<Pcurvy[indexeid*(Npoints-2)+a]<<endl;
3528  }
3529  st << std::scientific << std::setprecision(8) <<
3530  " "<<x_crit[v2]<<" "<< y_crit[v2] <<" "<< 0.000;
3531 //cout<<x_crit[v2]<<" "<<y_crit[v2]<<endl;
3532  edgenew->LinkEndChild(new TiXmlText(st.str()));
3533 
3534 //cout<<st.str()<<endl;
3535 
3536  //if(x_crit[cnt] < x_crit[cnt+1])
3537  //{
3538  // cout<<"Warning: the edges can have the vertices in the reverse order";
3539  //}
3540  }
3541 
3542  edgenew = edgenew->NextSiblingElement("E");
3543 
3544  }
3545 
3546  //write also the others layers curve points
3547  if(curv_lay == true)
3548  {
3549 cout<<"write other curved edges"<<endl;
3550  TiXmlElement * curved = meshnew->FirstChildElement("CURVED");
3551  int idcnt = 300;
3552  int nlays = lay_eids.num_elements();
3553  int neids_lay = lay_eids[0].num_elements();
3554  //TiXmlComment * comment = new TiXmlComment();
3555  //comment->SetValue(" new edges ");
3556  //curved->LinkEndChild(comment);
3557 
3558  for (int g=0; g< nlays; ++g)
3559  {
3560  for(int p=0; p< neids_lay; p++)
3561  {
3562  stringstream st;
3563  TiXmlElement * e = new TiXmlElement( "E" );
3564  e->SetAttribute("ID", idcnt++);
3565  e->SetAttribute("EDGEID", lay_eids[g][p]);
3566  e->SetAttribute("NUMPOINTS", Npoints);
3567  e->SetAttribute("TYPE", "PolyEvenlySpaced");
3568  for(int c=0; c< Npoints; c++)
3569  {
3570  st << std::scientific << std::setprecision(8) <<x_lay[g][p*Npoints +c]
3571  << " " << y_lay[g][p*Npoints +c] << " " << 0.000<<" ";
3572 
3573  }
3574 
3575  TiXmlText * t0 = new TiXmlText(st.str());
3576  e->LinkEndChild(t0);
3577  curved->LinkEndChild(e);
3578  }
3579 
3580  }
3581  }
3582 
3583 
3584 
3585  docnew.SaveFile( newfile );
3586 
3587  cout<<"new file: "<<newfile<<endl;
3588 
3589  }
3590 
Generated on Sun Mar 15 2015 00:11:33 for Nektar++ by   doxygen 1.8.1.2