Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
MeshMove.cpp
Go to the documentation of this file.
1 
2 ///////////////////////////////////////////////////////////////////////////////
3 //
4 // File MoveMeshToCriticalLayer.cpp
5 //
6 // For more information, please see: http://www.nektar.info
7 //
8 // The MIT License
9 //
10 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
11 // Department of Aeronautics, Imperial College London (UK), and Scientific
12 // Computing and Imaging Institute, University of Utah (USA).
13 //
14 // License for the specific language governing rights and limitations under
15 // Permission is hereby granted, free of charge, to any person obtaining a
16 // copy of this software and associated documentation files (the "Software"),
17 // to deal in the Software without restriction, including without limitation
18 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
19 // and/or sell copies of the Software, and to permit persons to whom the
20 // Software is furnished to do so, subject to the following conditions:
21 //
22 // The above copyright notice and this permission notice shall be included
23 // in all copies or substantial portions of the Software.
24 //
25 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
26 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
27 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
28 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
29 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
30 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
31 // DEALINGS IN THE SOFTWARE.
32 //
33 // Description: MoveMesh to critical layer
34 //
35 ///////////////////////////////////////////////////////////////////////////////
36 
37 #include <cstdio>
38 #include <cstdlib>
39 #include <iostream>
40 #include <iomanip>
41 //#include <sstream>
42 #include <MultiRegions/ExpList.h>
43 #include <MultiRegions/ExpList1D.h>
44 //#include <MultiRegions/ExpList0D.h>
45 //#include <MultiRegions/ExpList2D.h>
46 //#include <MultiRegions/ExpList3D.h>
47 //#include <MultiRegions/ExpList2DHomogeneous1D.h>
48 //#include <MultiRegions/ExpList3DHomogeneous1D.h>
49 //#include <MultiRegions/ExpList1DHomogeneous2D.h>
50 //#include <MultiRegions/ExpList3DHomogeneous2D.h>
51 #include <MultiRegions/ContField1D.h>
52 #include <MultiRegions/ContField2D.h>
53 //#include <MultiRegions/ContField3D.h>
54 //#include <MultiRegions/ContField3DHomogeneous1D.h>
55 //#include <MultiRegions/ContField3DHomogeneous2D.h>
56 #include <LocalRegions/SegExp.h>
57 #include <LocalRegions/QuadExp.h>
58 #include <LocalRegions/TriExp.h>
59 //#include </PreProcessing/NekMesh/Convert.h>
60 #include <LibUtilities/LinearAlgebra/Lapack.hpp>
61 #include <LibUtilities/BasicConst/NektarUnivTypeDefs.hpp>
62 //#include <LibUtilities/Foundations/GaussPoints.h>
63 #include <boost/lexical_cast.hpp>
64 #include <tinyxml.h>
65 
66 using namespace std;
67 using namespace Nektar;
68 
69 void OrderVertices(int nedges,SpatialDomains::MeshGraphSharedPtr graphShPt,
70  MultiRegions::ExpListSharedPtr & bndfield,
71  Array<OneD, int>& Vids, int v1,int v2 , NekDouble x_connect,
72  int & lastedge,
73  Array<OneD, NekDouble>& x, Array<OneD, NekDouble>& y);
74 void Computestreakpositions(int nvertl, MultiRegions::ExpListSharedPtr streak,
75  Array<OneD, NekDouble> xold_up, Array<OneD, NekDouble> yold_up,
76  Array<OneD, NekDouble> xold_low, Array<OneD, NekDouble> yold_low,
77  Array<OneD, NekDouble> xold_c, Array<OneD, NekDouble> yold_c,
78  Array<OneD, NekDouble> &xc, Array<OneD, NekDouble> &yc, NekDouble cr,
79  bool verts);
80 void GenerateAddPointsNewtonIt( NekDouble xi, NekDouble yi,NekDouble &xout, NekDouble &yout,
81  MultiRegions::ExpListSharedPtr function, Array<OneD, NekDouble> derfunction,
82  NekDouble cr);
83 
84 void GenerateMapEidsv1v2(MultiRegions::ExpListSharedPtr field,
85  Array<OneD, int> &V1, Array<OneD, int> &V2);
86 
87 void MappingEVids(Array<OneD, NekDouble> xoldup, Array<OneD, NekDouble> yoldup,
88  Array<OneD, NekDouble> xolddown, Array<OneD, NekDouble> yolddown,
89  Array<OneD, NekDouble> xcold, Array<OneD, NekDouble> ycold,
90  Array<OneD, int> Vids_c,
91  SpatialDomains::MeshGraphSharedPtr mesh,
92  MultiRegions::ExpListSharedPtr streak,
93  Array<OneD, int> V1, Array<OneD, int> V2,
94  int & nlays, Array<OneD, Array<OneD, int> >& Eids_lay,
95  Array<OneD, Array<OneD, int> >& Vids_lay);
96 bool checkcommonvert(Array<OneD, int> Vids_laybefore, Array<OneD, int> Vids_c, int Vid);
97 
98 void Cutrepetitions(int nedges,Array<OneD, NekDouble> inarray,
99  Array<OneD, NekDouble>& outarray);
100 
101 int DetermineclosePointxindex(NekDouble x,Array<OneD, NekDouble> xArray);
102 
103 void GenerateNeighbourArrays(int index, int neighpoints,Array<OneD, NekDouble> xArray,
104  Array<OneD, NekDouble> yArray,Array<OneD, NekDouble>& Neighbour_x,
105  Array<OneD, NekDouble>& Neighbour_y);
106 
107 NekDouble LagrangeInterpolant(NekDouble x, int npts,
108  Array<OneD,NekDouble> xpts, Array<OneD, NekDouble> funcvals);
109 
110 void EvaluateTangent(int npoints, Array<OneD, NekDouble> xcQedge,
111  Array<OneD, NekDouble> coeffsinterp,
112  Array<OneD, NekDouble> & txQedge, Array<OneD, NekDouble> & tyQedge);
113 void PolyInterp( Array<OneD, NekDouble> xpol, Array<OneD, NekDouble> ypol,
114  Array<OneD, NekDouble> & coeffsinterp,
115  Array<OneD, NekDouble> & xcout, Array<OneD, NekDouble> & ycout,
116  int edge, int npedge);
117 void PolyFit(int polyorder,int npoints,
118  Array<OneD, NekDouble> xin, Array<OneD, NekDouble> fin,
119  Array<OneD, NekDouble> & coeffsinterp,
120  Array<OneD, NekDouble> & xout, Array<OneD, NekDouble> & fout,
121  int npout);
122 void Orderfunctionx(Array<OneD, NekDouble> inarray_x,
123  Array<OneD, NekDouble> inarray_y, Array<OneD, NekDouble>& outarray_x,
124  Array<OneD, NekDouble>& outarray_y);
125 
126 void MoveLayersvertically(int nlays, int nvertl, int cntlow, int cntup,
127  Array<OneD, Array<OneD, int > > lay_Vids, Array<OneD, NekDouble> xc,
128  Array<OneD, NekDouble> yc, Array<OneD, int> Down, Array<OneD, int> Up,
129  Array<OneD, NekDouble >& xnew, Array<OneD, NekDouble>& ynew,
130  Array<OneD, Array<OneD, NekDouble > >& layers_x,
131  Array<OneD, Array<OneD, NekDouble > >& layers_y);
132 
133 void MoveLayerNfixedxpos(int nvertl, int npedge, Array<OneD, NekDouble> xcPhys,
134  Array<OneD, NekDouble> tmpx_lay, Array<OneD, NekDouble> tmpy_lay,
135  Array<OneD, int> Vids,
136  Array<OneD, NekDouble> &xlay, Array<OneD, NekDouble> &ylay,
137  Array<OneD, NekDouble> &xnew, Array<OneD, NekDouble> &ynew);
138 
139 void MoveLayerNnormpos(int nvertl, int npedge, Array<OneD, NekDouble> xcPhys,
140  Array<OneD, NekDouble> tmpx_lay, Array<OneD, NekDouble> tmpy_lay,
141  Array<OneD, int> Vids,
142  Array<OneD, NekDouble> &xlay, Array<OneD, NekDouble> &ylay,
143  Array<OneD, NekDouble> &xnew, Array<OneD, NekDouble> &ynew);
144 
145 void MoveOutsidePointsfixedxpos(int npedge, SpatialDomains::MeshGraphSharedPtr mesh,
146  Array<OneD, NekDouble> xcold,Array<OneD, NekDouble> ycold,
147  Array<OneD, NekDouble> xolddown,Array<OneD, NekDouble> yolddown,
148  Array<OneD, NekDouble> xoldup,Array<OneD, NekDouble> yoldup,
149  Array<OneD, NekDouble> ylaydown,Array<OneD, NekDouble> ylayup,
150  Array<OneD, NekDouble>& xnew,Array<OneD, NekDouble>& ynew);
151 
152 void MoveOutsidePointsNnormpos(int npedge, SpatialDomains::MeshGraphSharedPtr mesh,
153  Array<OneD, NekDouble> xcold,Array<OneD, NekDouble> ycold,
154  Array<OneD, NekDouble> xolddown,Array<OneD, NekDouble> yolddown,
155  Array<OneD, NekDouble> xoldup,Array<OneD, NekDouble> yoldup,
156  Array<OneD, NekDouble> xlaydown,Array<OneD, NekDouble> ylaydown,
157  Array<OneD, NekDouble> xlayup,Array<OneD, NekDouble> ylayup,
158  Array<OneD, NekDouble> nxPhys,Array<OneD, NekDouble> nyPhys,
159  Array<OneD, NekDouble>& xnew,Array<OneD, NekDouble>& ynew) ;
160 
161 void CheckSingularQuads( MultiRegions::ExpListSharedPtr Exp,
162  Array<OneD, int> V1, Array<OneD, int> V2,
163  Array<OneD, NekDouble>& xnew,Array<OneD, NekDouble>& ynew);
164 
165 void Replacevertices(string filename, Array<OneD, NekDouble> newx,
166  Array<OneD, NekDouble> newy,
167  Array<OneD, NekDouble> xcPhys, Array<OneD, NekDouble> ycPhys,
168  Array<OneD, int>Eids, int Npoints, string s_alp,
169  Array<OneD, Array<OneD, NekDouble> > x_lay,
170  Array<OneD, Array<OneD, NekDouble> > y_lay,
171  Array<OneD, Array<OneD, int > >lay_eids, bool curv_lay) ;
172 
173 int main(int argc, char *argv[])
174 {
175  NekDouble cr;
176  //set cr =0
177  cr=0;
178  //change argc from 6 to 5 allow the loading of cr to be optional
179  if(argc > 6 || argc < 5)
180  {
181  fprintf(stderr,
182  "Usage: ./MoveMesh meshfile fieldfile changefile alpha cr(optional)\n");
183  exit(1);
184  }
185 
186  //ATTEnTION !!! with argc=2 you impose that vSession refers to is argv[1]=meshfile!!!!!
187  LibUtilities::SessionReaderSharedPtr vSession
188  = LibUtilities::SessionReader::CreateInstance(2, argv);
189  //----------------------------------------------
190 
191  if( argc == 6 &&
192  vSession->DefinesSolverInfo("INTERFACE")
193  && vSession->GetSolverInfo("INTERFACE")=="phase" )
194  {
195  cr = boost::lexical_cast<NekDouble>(argv[argc-1]);
196  argc=5;
197  }
198 
199  // Read in mesh from input file
200  string meshfile(argv[argc-4]);
201  SpatialDomains::MeshGraphSharedPtr graphShPt = SpatialDomains::MeshGraph::Read(meshfile);
202  //----------------------------------------------
203 
204  // Also read and store the boundary conditions
205  SpatialDomains::BoundaryConditionsSharedPtr boundaryConditions;
206  boundaryConditions = MemoryManager<SpatialDomains::BoundaryConditions>
207  ::AllocateSharedPtr(vSession,graphShPt);
208  SpatialDomains::BoundaryConditions bcs(vSession, graphShPt);
209  //----------------------------------------------
210 
211 
212  //the mesh file should have 2 component: set output fields
213  //fields has to be of the SAME dimension of the mesh (that's why there is
214  //the changefile as an input)
215  //a contfield2D is needed to extract boundary conditions!!!
216 
217  // store name of the file to change
218  string changefile(argv[argc-2]);
219  //----------------------------------------------
220 
221  //store the value of alpha
222  string charalp (argv[argc-1]);
223  //NekDouble alpha = boost::lexical_cast<NekDouble>(charalp);
224  cout<<"read alpha="<<charalp<<endl;
225 
226  //---------------------------------------------
227  // Import field file.
228  string fieldfile(argv[argc-3]);
229  vector<LibUtilities::FieldDefinitionsSharedPtr> fielddef;
230  vector<vector<NekDouble> > fielddata;
231  //----------------------------------------------
232 
233  MultiRegions::ExpListSharedPtr streak;
234  streak = MemoryManager<MultiRegions::ContField2D>
235  ::AllocateSharedPtr(vSession, graphShPt, "w",true);
236 
237  LibUtilities::Import(fieldfile,fielddef,fielddata);
238 
239  for(int i=0; i<fielddata.size(); i++)
240  {
241  streak->ExtractDataToCoeffs(fielddef[i], fielddata[i], fielddef[i]->m_fields[0], streak->UpdateCoeffs());
242  }
243  streak->BwdTrans_IterPerExp(streak->GetCoeffs(), streak->UpdatePhys());
244 
245  //------------------------------------------------
246  // determine the I regions (3 region expected)
247  // hypothesys: the layes have the same number of points
248 
249  int nIregions, lastIregion;
250  const Array<OneD, SpatialDomains::BoundaryConditionShPtr> bndConditions = streak->GetBndConditions();
251  Array<OneD, int> Iregions =Array<OneD, int>(bndConditions.num_elements(),-1);
252 
253  nIregions=0;
254  int nbnd= bndConditions.num_elements();
255  for(int r=0; r<nbnd; r++)
256  {
257  if(bndConditions[r]->GetUserDefined()=="CalcBC")
258  {
259  lastIregion=r;
260  Iregions[r]=r;
261  nIregions++;
262  }
263  }
264 
265  ASSERTL0(nIregions>0,"there is any boundary region with the tag USERDEFINEDTYPE=""CalcBC"" specified");
266  cout<<"nIregions="<<nIregions<<endl;
267  //set expansion along a layers
268  Array<OneD, MultiRegions::ExpListSharedPtr> bndfieldx= streak->GetBndCondExpansions();
269  //--------------------------------------------------------
270  //determine the points in the lower and upper curve...
271  int nedges = bndfieldx[lastIregion]->GetExpSize();
272  int nvertl = nedges +1 ;
273  Array<OneD, int> Vids_low(nvertl,-10);
274  Array<OneD, NekDouble> xold_low(nvertl);
275  Array<OneD, NekDouble> yold_low(nvertl);
276  Array<OneD, NekDouble> zi(nvertl);
277 
278  //order the ids on the lower curve lastIregion starting from the id on x=0
279  NekDouble x_connect;
280  NekDouble x0,y0,z0,yt=0,zt=0;
281  int lastedge=-1;
282  int v1,v2;
283  //first point for x_connect=0(or-1.6 for the full mesh (-pi,pi) )
284  x_connect=0;
285  SpatialDomains::PointGeomSharedPtr vertex0 =
286  graphShPt->GetVertex
287  (
288  ( (bndfieldx[lastIregion]->GetExp(0)->as<LocalRegions::SegExp>())
289  ->GetGeom1D()
290  )
291  ->GetVid(0)
292  );
293  vertex0->GetCoords(x0,y0,z0);
294  if( x0 != 0.0)
295  {
296  cout<<"WARNING x0="<<x0<<endl;
297  x_connect=x0;
298  }
299  v1=0;
300  v2=1;
301  OrderVertices(nedges, graphShPt, bndfieldx[lastIregion-1],
302  Vids_low, v1, v2 , x_connect ,lastedge, xold_low,yold_low);
303  ASSERTL0(Vids_low[v2]!=-10, "Vids_low[v2] is wrong");
304  SpatialDomains::PointGeomSharedPtr vertex = graphShPt->GetVertex(Vids_low[v2]);
305 
306  //update x_connect
307  cout<<"x_conn="<<x_connect<<" yt="<<yt<<" zt="<<zt<<" vid="<<Vids_low[v2]<<endl;
308  vertex->GetCoords(x_connect,yt,zt);
309 
310  int i=2;
311  while(i<nvertl)
312  {
313  v1=i;
314  OrderVertices(nedges, graphShPt, bndfieldx[lastIregion-1],
315  Vids_low, v1, v2 , x_connect, lastedge, xold_low, yold_low );
316  SpatialDomains::PointGeomSharedPtr vertex = graphShPt->GetVertex(Vids_low[v1]);
317  //update x_connect (lastedge is updated on the OrderVertices function)
318  vertex->GetCoords(x_connect,yt,zt);
319  i++;
320  }
321 
322  //------------------------------------------------------------------------
323  //order in the same way the id of the upper curve lastIregion-1 starting from x=0:
324  Array<OneD, int> Vids_up(nvertl,-10);
325  Array<OneD,NekDouble> xold_up(nvertl);
326  Array<OneD,NekDouble> yold_up(nvertl);
327  //first point for x_connect=0 (or-1.6)
328  x_connect=0;
329  vertex0 =
330  graphShPt->GetVertex
331  (
332  ( (bndfieldx[lastIregion]->GetExp(0)->as<LocalRegions::SegExp>())
333  ->GetGeom1D()
334  )
335  ->GetVid(0)
336  );
337  vertex0->GetCoords(x0,y0,z0);
338  if( x0 != 0.0)
339  {
340  cout<<"WARNING x0="<<x0<<endl;
341  x_connect=x0;
342  }
343  lastedge=-1;
344 
345  v1=0;
346  v2=1;
347  OrderVertices(nedges, graphShPt, bndfieldx[lastIregion-2 ],
348  Vids_up, v1, v2 , x_connect ,lastedge, xold_up, yold_up);
349  SpatialDomains::PointGeomSharedPtr vertexU = graphShPt->GetVertex(Vids_up[v2]);
350  vertexU->GetCoords(x_connect,yt,zt);
351 
352  i=2;
353  while(i<nvertl)
354  {
355  v1=i;
356  OrderVertices(nedges, graphShPt, bndfieldx[lastIregion-2],
357  Vids_up, v1, v2 , x_connect, lastedge, xold_up, yold_up );
358  SpatialDomains::PointGeomSharedPtr vertex = graphShPt->GetVertex(Vids_up[v1]);
359  //cout<<"VIdup="<<Vids_up[v1]<<endl;
360  //update x_connect (lastedge is updated on the OrderVertices function)
361  vertex->GetCoords(x_connect,yt,zt);
362  i++;
363  }
364 
365  //-----------------------------------------------------------------------------------
366  //order in the same way the id of the layer curve lastIregion starting from x=0:
367  Array<OneD, int> Vids_c(nvertl,-10);
368  Array<OneD,NekDouble> xold_c(nvertl);
369  Array<OneD,NekDouble> yold_c(nvertl);
370  //first point for x_connect=0(or-1.6)
371  x_connect=0;
372  vertex0 =
373  graphShPt->GetVertex(((bndfieldx[lastIregion]->GetExp(0)
374  ->as<LocalRegions::SegExp>())->GetGeom1D())->GetVid(0));
375  vertex0->GetCoords(x0,y0,z0);
376  if( x0 != 0.0)
377  {
378  cout<<"WARNING x0="<<x0<<endl;
379  x_connect=x0;
380  }
381  lastedge=-1;
382 
383  v1=0;
384  v2=1;
385 
386  OrderVertices(nedges, graphShPt, bndfieldx[lastIregion],
387  Vids_c, v1, v2 , x_connect ,lastedge, xold_c, yold_c);
388  SpatialDomains::PointGeomSharedPtr vertexc = graphShPt->GetVertex(Vids_c[v2]);
389 
390  //update x_connect
391  vertexc->GetCoords(x_connect,yt,zt);
392 
393  i=2;
394  while(i<nvertl)
395  {
396  v1=i;
397  OrderVertices(nedges, graphShPt, bndfieldx[lastIregion],
398  Vids_c, v1, v2 , x_connect, lastedge, xold_c, yold_c );
399  SpatialDomains::PointGeomSharedPtr vertex = graphShPt->GetVertex(Vids_c[v1]);
400 //cout<<"Vids cl="<<Vids_low[v1]<<endl;
401  //update x_connect (lastedge is updated on the OrderVertices function)
402  vertex->GetCoords(x_connect,yt,zt);
403  i++;
404  }
405 
406  //calculate the distances between the layer and the upper/lower curve
407 
408  Array<OneD, NekDouble> Deltaup (nvertl, -200);
409  Array<OneD, NekDouble> Deltalow (nvertl, -200);
410 
411  for(int r=0; r<nvertl; r++)
412  {
413  //Always positive!!!
414 //cout<<"i="<<r<<" yup="<<yold_up[r]<<" yc="<<yold_c[r]<<" ylow="<<yold_low[r]<<endl;
415  Deltaup[r] = yold_up[r] - yold_c[r];
416  Deltalow[r] = yold_c[r] - yold_low[r];
417  ASSERTL0(Deltaup[r]>0, "distance between upper and layer curve is not positive");
418  ASSERTL0(Deltalow[r]>0, "distance between lower and layer curve is not positive");
419  }
420  //------------------------------------------------------------------------
421 
422 
423  //fieds to force continuity:
424  const SpatialDomains::BoundaryRegionCollection &bregions = bcs.GetBoundaryRegions();
425  MultiRegions::ContField1DSharedPtr Cont_y;
426  MultiRegions::ExpList1DSharedPtr yexp;
427 
428  yexp = MemoryManager<MultiRegions::ExpList1D>
429  ::AllocateSharedPtr(vSession,*(bregions.find(lastIregion)->second), graphShPt, true);
430  Cont_y = MemoryManager<MultiRegions::ContField1D>
431  ::AllocateSharedPtr(vSession, *yexp);
432  //--------------------------------------
433  /*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/
434  //generate additional points using the Newton iteration
435  //determine the xposition for every edge (in the middle even if it
436  // is not necessary
437  //PARAMETER which determines the number of points per edge @todo put as an input
438  int npedge;
439  if(vSession->DefinesParameter("npedge"))
440  {
441  npedge = (int)vSession->GetParameter("npedge");
442  }
443  else
444  {
445  npedge = 5;//default value
446  }
447  /*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/
448  //find the points where u=0 and determine the sign of the shift and the delta
449  int nq= streak->GetTotPoints();
450  Array<OneD, NekDouble> x(nq);
451  Array<OneD,NekDouble> y(nq);
452  streak->GetCoords(x,y);
453  Array<OneD, NekDouble> x_c(nvertl);
454  Array<OneD,NekDouble> y_c(nvertl,-200);
455  Array<OneD, NekDouble> tmp_w (nvertl, 200);
456  Array<OneD, int> Sign (nvertl,1);
457  Array<OneD, NekDouble> Delta_c(nvertl,-200);
458 
459 
460  Computestreakpositions(nvertl, streak, xold_up, yold_up,
461  xold_low, yold_low, xold_c, yold_c, x_c, y_c,cr,true);
462  // if the curve is low the old layer point, it has to shift down
463  NekDouble shift = 0;
464  for(int q=0; q<nvertl; q++)
465  {
466  if(y_c[q] < yold_c[q])
467  {
468  Sign[q] = -1;
469  }
470  //calculate delta
471  Delta_c[q] = abs(yold_c[q]-y_c[q]);
472  //check the shifting of the layer:
473  shift+= Delta_c[q];
474  cout<<x_c[q]<<" "<<y_c[q]<<endl;
475  }
476  //cout<<"shift="<<shift<<endl;
477  if(shift<0.001)
478  {
479  cout<<"Warning: the critical layer is stationary"<<endl;
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  Array<OneD, NekDouble> Addpointsx (nedges*(npedge-2), 0.0);
489  Array<OneD, NekDouble> Addpointsy (nedges*(npedge-2), 0.0);
490  //calculate the dU_dy
491  Array<OneD, NekDouble> dU(streak->GetTotPoints());
492  streak->PhysDeriv(MultiRegions::eY, streak->GetPhys(), dU);
493 
494 
495  LocalRegions::SegExpSharedPtr bndSegExp;
496 
497  int Eid,id1,id2;
498  NekDouble x1,y1,z1;
499  NekDouble x2,y2,z2;
500  SpatialDomains::PointGeomSharedPtr vertex1;
501  SpatialDomains::PointGeomSharedPtr vertex2;
502  for(int r=0; r<nedges; r++)
503  {
504 
505  bndSegExp = bndfieldx[lastIregion]->GetExp(r)
506  ->as<LocalRegions::SegExp>();
507  Eid = (bndSegExp->GetGeom1D())->GetEid();
508  id1 = (bndSegExp->GetGeom1D())->GetVid(0);
509  id2 = (bndSegExp->GetGeom1D())->GetVid(1);
510  vertex1 = graphShPt->GetVertex(id1);
511  vertex2 = graphShPt->GetVertex(id2);
512  vertex1->GetCoords(x1,y1,z1);
513  vertex2->GetCoords(x2,y2,z2);
514  //cout<<"edge="<<r<<" x1="<<x1<<" x2="<<x2<<endl;
515  //cout<<"edge="<<r<<" y1="<<y1<<" y2="<<y2<<endl;
516  cout<<"edge="<<r<<" x1="<<x1<<" y1="<<y1<<" x2="<<x2<<" y2="<<y2<<endl;
517  if(x2>x1)
518  {
519  Cpointsx[r] = x1 +(x2-x1)/2;
520  //cout<<"edge="<<r<<" x1="<<x1<<" x2="<<x2<<" Cx="<<Cpointsx[r]<<endl;
521  //cout<<"edge="<<r<<" x1="<<x1<<" y1="<<y1<<" x2="<<x2<<" y2="<<y2<<endl;
522  if( Cpointsx[r]>x2 || Cpointsx[r]< x1)
523  {
524  Cpointsx[r] = -Cpointsx[r];
525  }
526  for(int w=0; w< npedge-2; w++)
527  {
528 
529  Addpointsx[r*(npedge-2) +w] = x1 +((x2-x1)/(npedge - 1))*(w+1);
530  if( Addpointsx[r*(npedge-2) +w] > x2 || Addpointsx[r*(npedge-2) +w] < x1)
531  {
532  Addpointsx[r*(npedge-2) +w] = -Addpointsx[r*(npedge-2) +w];
533  }
534  //initial guess along the line defined by the NEW verts y_c
535  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);
536  //Addpointsy[r*(npedge-2) +w] = y1 + ((y2-y1)/(x2-x1))*(Addpointsx[r*(npedge-2) +w]-x1);
537  GenerateAddPointsNewtonIt( Addpointsx[r*(npedge-2) +w], Addpointsy[r*(npedge-2) +w],
538  Addpointsx[r*(npedge-2) +w], Addpointsy[r*(npedge-2) +w], streak, dU,cr);
539 
540  // Lay_x->UpdatePhys()[r*npedge +1 +w]= Addpointsx[r*(npedge-2) +w];
541  // Lay_y->UpdatePhys()[r*npedge +1 +w]= Addpointsy[r*(npedge-2) +w];
542  }
543  //Lay_x->UpdatePhys()[r*npedge +0] =x1;
544  //Lay_y->UpdatePhys()[r*npedge +0] =y1;
545  //Lay_x->UpdatePhys()[r*npedge +npedge-1] =x2;
546  //Lay_y->UpdatePhys()[r*npedge +npedge-1] =y2;
547 
548  }
549  else if(x1>x2)
550  {
551  Cpointsx[r] = x2+ (x1-x2)/2;
552 //cout<<"edge="<<r<<" y1="<<y1<<" y2="<<y2<<endl;
553  if( Cpointsx[r] > x1 || Cpointsx[r] < x2)
554  {
555  Cpointsx[r] = -Cpointsx[r];
556  }
557  for(int w=0; w< npedge-2; w++)
558  {
559  Addpointsx[r*(npedge-2) +w] = x2 +((x1-x2)/(npedge - 1))*(w+1);
560  if( Addpointsx[r*(npedge-2) +w] > x1 || Addpointsx[r*(npedge-2) +w] < x2)
561  {
562  Addpointsx[r*(npedge-2) +w] = -Addpointsx[r*(npedge-2) +w];
563  }
564 
565  //initial guess along the line defined by the NEW verts y_c
566  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);
567 
568  //Addpointsy[r*(npedge-2) +w] = y2 + ((y1-y2)/(x1-x2))*(Addpointsx[r*(npedge-2) +w]-x2);
569  GenerateAddPointsNewtonIt( Addpointsx[r*(npedge-2) +w], Addpointsy[r*(npedge-2) +w],
570  Addpointsx[r*(npedge-2) +w], Addpointsy[r*(npedge-2) +w], streak, dU,cr);
571  // Lay_x->UpdatePhys()[r*npedge +1]= Addpointsx[r*(npedge-2) +w];
572  // Lay_y->UpdatePhys()[r*npedge +1]= Addpointsy[r*(npedge-2) +w];
573  }
574  //Lay_x->UpdatePhys()[r*npedge +0] =x2;
575  //Lay_y->UpdatePhys()[r*npedge +0] =y2;
576  //Lay_x->UpdatePhys()[r*npedge +npedge-1] =x1;
577  //Lay_y->UpdatePhys()[r*npedge +npedge-1] =y1;
578  }
579  else
580  {
581  ASSERTL0(false, "point not generated");
582  }
583 //cout<<"calculate cpoints coords"<<endl;
584  //Cpointsy[r] = y1 + (y2-y1)/2;
585 //cout<<"central point:"<<endl;
586  //GenerateAddPointsNewtonIt( Cpointsx[r], Cpointsy[r],Cpointsx[r], Cpointsy[r],
587  // streak, dU,cr);
588  //NekDouble diff = Cpointsy[r]-Addpointsy[r*(npedge-2)];
589 //cout<<"diff="<<diff<<endl;
590  Eids[r] = Eid;
591 
592  }
593  //-------------------------------------------------------------
594 
595 
596  //fill the xPhys,yPhys array( may necessary after)
597  Array<OneD, NekDouble> xcPhys (nedges*npedge, 0.0);
598  Array<OneD, NekDouble> ycPhys (nedges*npedge, 0.0);
599 
600  for(int a=0; a<nedges; a++)
601  {
602  //v1
603  xcPhys[a*npedge+0] = x_c[a];
604  ycPhys[a*npedge+0] = y_c[a];
605  //v2
606  xcPhys[a*npedge+npedge-1] = x_c[a+1];
607  ycPhys[a*npedge+npedge-1] = y_c[a+1];
608 
609  for(int b=0; b<npedge-2; b++)
610  {
611  xcPhys[a*npedge +b+1] = Addpointsx[a*(npedge-2)+b];
612  ycPhys[a*npedge +b+1] = Addpointsy[a*(npedge-2)+b];
613  }
614  }
615 
616 cout<<"xc,yc before tanevaluate"<<endl;
617 for(int v=0; v< xcPhys.num_elements(); v++)
618 {
619 cout<<xcPhys[v]<<" "<<ycPhys[v]<<endl;
620 }
621 
622  //-------------------------------------------------
623 
624  //V1[eid],V2[eid] vertices associate with the edge Id=eid
625  Array<OneD, int> V1;
626  Array<OneD, int> V2;
627 
628  GenerateMapEidsv1v2(streak,V1,V2);
629  Array<OneD, Array<OneD, int> > lay_Eids;
630  Array<OneD, Array<OneD, int> > lay_Vids;
631  int nlays=0;
632  MappingEVids(xold_up, yold_up, xold_low, yold_low, xold_c, yold_c, Vids_c,
633  graphShPt,streak, V1, V2, nlays, lay_Eids, lay_Vids);
634 
635 
636 
637 cout<<"nlays="<<nlays<<endl;
638  Array<OneD, Array<OneD, NekDouble> > layers_y(nlays);
639  Array<OneD, Array<OneD, NekDouble> > layers_x(nlays);
640  //initialise layers_y,lay_eids
641  for(int g=0; g<nlays; g++)
642  {
643  layers_y[g]= Array<OneD, NekDouble> ( (nvertl-1)*npedge );
644  }
645 
646 
647 
648  /////////////////////////////////////////////////////////////////////////////
649  ////////////////////////////////////////////////////////////////////////////
650  //££££££££££££££££££££££££££££££££££££££££££££££££££££££££££££££££££££
651  //$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
652  //@todo set Delta0 from session file
653  NekDouble Delta0;
654  if(vSession->DefinesParameter("Delta"))
655  {
656  Delta0 = vSession->GetParameter("Delta");
657  }
658  else
659  {
660  Delta0 = 0.1;//default value
661  }
662 
663  //$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
664  //££££££££££££££££££££££££££££££££££££££££££££££££££££££££££££££££££££
665  ////////////////////////////////////////////////////////////////////////////
666  /////////////////////////////////////////////////////////////////////////////
667  //save the coords of the old vertices
668  int nVertTot = graphShPt->GetNvertices();
669  //arrays to check with the new values
670  Array<OneD, NekDouble> xold(nVertTot);
671  Array<OneD, NekDouble> yold(nVertTot);
672  //calculate the new cordinates of the vertices
673  Array<OneD, NekDouble> xnew(nVertTot);
674  Array<OneD, NekDouble> ynew(nVertTot,-20);
675  Array<OneD, int> Up(nvertl);//Vids lay Up
676  Array<OneD, int> Down(nvertl);//Vids lay Down
677  int cntup=0;
678  int cntlow=0;
679  //Vids needed only if a layers has to be moved
680  NekDouble bleft=-10;
681  NekDouble tright = 10;
682  NekDouble bright = -10;
683  NekDouble tleft = 10;
684  for(int i=0; i<nVertTot; i++)
685  {
686  bool mvpoint =false;
687  SpatialDomains::PointGeomSharedPtr vertex = graphShPt->GetVertex(i);
688  NekDouble x,y,z;
689  vertex->GetCoords(x,y,z);
690 
691  xold[i] = x;
692  yold[i] = y;
693 
694 
695  //x coord doesn't change
696  xnew[i]=x;
697  //cout<<"x="<<x<<" y="<<y<<endl;
698  //bottom, left (x=0, y<ydown)
699  if(x==0 && y< yold_low[0]
700  && y> bleft)
701  {
702  bleft = y;
703  }
704  //top, right
705  if(x== xold_c[nvertl-1] && y> yold_up[nvertl-1]
706  && y<tright)
707  {
708  tright = y;
709  }
710  //bottom, right]
711  if(x==xold_c[nvertl-1] && y<yold_low[nvertl-1]
712  && y> bright)
713  {
714  bright = y;
715  }
716  //top, left
717  if(x== 0 && y> yold_up[0]
718  && y<tleft)
719  {
720  tleft = y;
721  }
722  //find the corresponding yold_l and deltaold
723  for(int j=0; j<nvertl; j++)
724  {
725  if((xold_up[j]==x)&&(yold_up[j]==y))
726  {
727  //ratio = (1-y)*(1+y)/( (1-yold_c[j])*(1+yold_c[j]) );
728  //ynew[i] = y + Sign[j]*Delta_c[j]*ratio;
729  ynew[i] = y_c[j] +Delta0;
730  mvpoint=true;
731  Up[j] = i;
732  }
733  if((xold_low[j]==x)&&(yold_low[j]==y))
734  {
735  //ratio = (1-y)*(1+y)/( (1-yold_c[j])*(1+yold_c[j]) );
736  //ynew[i] = y + Sign[j]*Delta_c[j]*ratio;
737  ynew[i] = y_c[j] -Delta0;
738  mvpoint=true;
739  Down[j] = i;
740  }
741  if((xold_c[j]==x)&&(yold_c[j]==y))
742  {
743  ynew[i] = y_c[j];
744  mvpoint=true;
745  }
746  }
747  if(mvpoint==false)
748  {
749  //determine the closer xold_up
750  NekDouble diff=1000;
751  int qp_closer = 0;
752  for(int k=0; k<nvertl; k++)
753  {
754  if(abs(x-xold_up[k]) < diff)
755  {
756  diff = abs(x-xold_up[k]);
757  qp_closer=k;
758  }
759  }
760  if( y>yold_up[qp_closer] && y< 1)
761  {
762 
763  //ratio = (1-y)*(1+y)/( (1-yold_c[qp_closer])*(1+yold_c[qp_closer]) );
764  //ynew[i] = y + Sign[qp_closer]*Delta_c[qp_closer]*ratio;
765 
766  //ratio = (1-y)*(1+y)/( (1-yold_up[qp_closer])*(1+yold_up[qp_closer]) );
767  //distance prop to layerup
768  ynew[i] = y_c[qp_closer] +(y-yold_c[qp_closer])*
769  (1-y_c[qp_closer])/(1-yold_c[qp_closer]);
770 //cout<<"upper zone y="<<y<<" ratio="<<ratio<<endl;
771 
772 
773 
774  }
775  else if(y<yold_low[qp_closer] && y> -1)
776  {
777 
778  //ratio = (1-y)*(1+y)/( (1-yold_c[qp_closer])*(1+yold_c[qp_closer]) );
779  //ynew[i] = y + Sign[qp_closer]*Delta_c[qp_closer]*ratio;
780 
781  //ratio = (1-y)*(1+y)/( (1-yold_low[qp_closer])*(1+yold_low[qp_closer]) );
782  //distance prop to layerlow
783  ynew[i] = y_c[qp_closer] + (y-yold_c[qp_closer] )*
784  (-1-y_c[qp_closer])/(-1-yold_c[qp_closer]);
785 
786  }
787 
788  else if ( y>yold_c[qp_closer] && y < yold_up[qp_closer])
789  {
790  if(x==0){ cntup++; }
791  //ratio = (1-y)*(1+y)/( (1-yold_c[qp_closer])*(1+yold_c[qp_closer]) );
792  //ynew[i] = y + Sign[qp_closer]*Delta_c[qp_closer]*ratio;
793 
794  }
795  else if (y<yold_c[qp_closer] && y > yold_low[qp_closer])
796  {
797  if(x==0){ cntlow++; }
798  // ratio = (1-y)*(1+y)/( (1-yold_c[qp_closer])*(1+yold_c[qp_closer]) );
799  // ynew[i] = y + Sign[qp_closer]*Delta_c[qp_closer]*ratio;
800 
801  }
802  else if( y==1 || y==-1)//bcs don't move
803  {
804  ynew[i] =y;
805  //cout<<"point x="<<xnew[i]<<" y="<<y<<" closer x="<<xold_up[qp_closer]<<endl;
806  }
807 
808  //internal layers are not moved yet so...
809  if( (ynew[i]>1 || ynew[i]<-1)
810  && ( y>yold_up[qp_closer] || y<yold_low[qp_closer]) )
811  {
812  cout<<"point x="<<xnew[i]<<" y="<<y<<" closer x="<<xold_up[qp_closer]<<" ynew="<<ynew[i]<<endl;
813  ASSERTL0(false, "shifting out of range");
814  }
815 
816  }
817 
818  }
819 
820 
821 
822 
823  int nqedge = streak->GetExp(0)->GetNumPoints(0);
824  int nquad_lay = (nvertl-1)*nqedge;
825  Array<OneD, NekDouble> coeffstmp(Cont_y->GetNcoeffs(),0.0);
826  // curve the edges around the NEW critical layer (bool to turn on/off)
827  bool curv_lay=true;
828  bool move_norm=true;
829  int np_lay = (nvertl-1)*npedge;//nedges*npedge (Eq. Points!!!)
830 
831  Array<OneD, NekDouble> xcQ(nqedge*nedges,0.0);
832  Array<OneD, NekDouble> ycQ(nqedge*nedges,0.0);
833  Array<OneD, NekDouble> zcQ(nqedge*nedges,0.0);
834  Array<OneD, NekDouble> nxPhys(npedge*nedges,0.0);
835  Array<OneD, NekDouble> nyPhys(npedge*nedges,0.0);
836  Array<OneD, NekDouble> nxQ(nqedge*nedges,0.0);
837  Array<OneD, NekDouble> nyQ(nqedge*nedges,0.0);
838 
839  if( move_norm==true )
840  {
841  //np_lay = (nvertl-1)*nqedge;//nedges*nqedge (Q points!!!)
842  //extract crit lay normals (through tangents):
843  Array<OneD, NekDouble> xcPhysMOD(xcPhys.num_elements());
844  Array<OneD, NekDouble> ycPhysMOD(ycPhys.num_elements());
845  //copy(temporary the xcPhys,ycPhys into xcPhysMOD, ycPhysMOD)
846  Vmath::Vcopy(xcPhys.num_elements(),xcPhys,1,xcPhysMOD,1);
847  Vmath::Vcopy(xcPhys.num_elements(),ycPhys,1,ycPhysMOD,1);
848 
849  Array<OneD, StdRegions::StdExpansion1DSharedPtr> Edge_newcoords(2);
850 
851 cout<<"nquad per edge="<<nqedge<<endl;
852  for(int l=0; l<2; l++)
853  {
854  Edge_newcoords[l] = bndfieldx[lastIregion]->GetExp(0)
855  ->as<StdRegions::StdExpansion1D>();
856  }
857  Array<OneD, NekDouble> xnull(nqedge);
858  Array<OneD, NekDouble> ynull(nqedge);
859  Array<OneD, NekDouble> xcedgeQ(nqedge,0.0);
860  Array<OneD, NekDouble> ycedgeQ(nqedge,0.0);
861  Array<OneD, NekDouble> txedgeQ(nqedge,0.0);
862  Array<OneD, NekDouble> tyedgeQ(nqedge,0.0);
863  Array<OneD, NekDouble> normsQ(nqedge,0.0);
864 
865  Array<OneD, NekDouble> txQ(nqedge*nedges,0.0);
866  Array<OneD, NekDouble> tyQ(nqedge*nedges,0.0);
867  Array<OneD, NekDouble> tx_tyedgeQ(nqedge,0.0);
868  Array<OneD, NekDouble> nxedgeQ(nqedge,0.0);
869  Array<OneD, NekDouble> nyedgeQ(nqedge,0.0);
870  Array<OneD, const NekDouble> ntempQ(nqedge) ;
871 
872  Array<OneD, NekDouble> nxedgePhys(npedge,0.0);
873  Array<OneD, NekDouble> nyedgePhys(npedge,0.0);
874 
875 
876  Array<OneD, NekDouble> CoordsPhys(2);
877 
878 
879  bndfieldx[lastIregion]->GetCoords(xcQ, ycQ, zcQ);
880  //determine the NEW crit lay quad points values(lagrangeInterpolant):
881  //interp(from xcPhys, ycPhys).
882  Array<OneD, NekDouble>x_tmp(np_lay-(nedges-1),0.0);
883  Array<OneD, NekDouble>y_tmp(np_lay-(nedges-1),0.0);
884  Array<OneD, NekDouble>closex(4,0.0);
885  Array<OneD, NekDouble>closey(4,0.0);
886  Cutrepetitions(nedges, xcPhysMOD,x_tmp);
887  Cutrepetitions(nedges, ycPhysMOD,y_tmp);
888  for(int l=0; l< xcQ.num_elements(); l++)
889  {
890  Cutrepetitions(nedges, xcPhysMOD,x_tmp);
891  Cutrepetitions(nedges, ycPhysMOD,y_tmp);
892  int indclose = DetermineclosePointxindex( xcQ[l], x_tmp);
893 
894  //generate neighbour arrays
895  GenerateNeighbourArrays(indclose, 4,x_tmp,y_tmp,closex,closey);
896 
897  ycQ[l]= LagrangeInterpolant(
898  xcQ[l],4,closex,closey );
899 
900 
901 
902  }
903 
904 
905 
906  //force continuity
907 
908  Vmath::Vcopy(nquad_lay, ycQ,1, Cont_y->UpdatePhys(),1);
909  Array<OneD, NekDouble> coeffsy(Cont_y->GetNcoeffs(),0.0);
910 
911  Cont_y->FwdTrans_IterPerExp(Cont_y->GetPhys(), coeffsy);
912  Cont_y->BwdTrans_IterPerExp( coeffsy, Cont_y->UpdatePhys());
913  Vmath::Vcopy(nquad_lay, Cont_y->GetPhys(),1, ycQ,1);
914 
915 cout<<"xcQ, ycQ"<<endl;
916 for(int s=0; s<xcQ.num_elements(); s++)
917 {
918 cout<<xcQ[s]<<" "<<ycQ[s]<<endl;
919 }
920 //ASSERTL0(false, "dsdfs");
921  bool evaluatetan=false;
922  NekDouble incratio;
923  Array<OneD, int> edgeinterp(2);
924  int wrong=0;
925  for(int k=0; k<nedges; k++)
926  {
927  Vmath::Vcopy(nqedge, &xcQ[k*nqedge],1,&xcedgeQ[0],1);
928  Vmath::Vcopy(nqedge, &ycQ[k*nqedge],1,&ycedgeQ[0],1);
929  //calc the NEW tangent values
930  Edge_newcoords[0]->StdPhysDeriv(xcedgeQ,txedgeQ);
931  Edge_newcoords[1]->StdPhysDeriv(ycedgeQ,tyedgeQ);
932  //norms=tx*tx
933  Vmath::Vmul(nqedge,txedgeQ,1,txedgeQ,1,normsQ,1);
934  //norms=tx*tx+ty*ty
935  Vmath::Vvtvp(nqedge,tyedgeQ,1,tyedgeQ,1,normsQ,1,normsQ,1);
936 
937  Vmath::Vsqrt(nqedge, normsQ,1,normsQ,1);
938  Vmath::Sdiv(nqedge,1.0,normsQ,1,normsQ,1);
939 
940  Vmath::Vmul(nqedge,txedgeQ,1,normsQ,1,txedgeQ,1);
941  Vmath::Vmul(nqedge,tyedgeQ,1,normsQ,1,tyedgeQ,1);
942 
943  //try evaluate tangent if the incremental ratio is high
944 
945 
946 
947 
948 
949  evaluatetan=false;
950  for(int u=0; u<nqedge-1; u++)
951  {
952  incratio = (ycedgeQ[u+1]- ycedgeQ[u])/(xcedgeQ[u+1]- xcedgeQ[u]);
953 cout<<"incratio="<<incratio<<endl;
954  if(abs(incratio)> 4.0 && evaluatetan==false )
955  {
956 cout<<"wrong="<<wrong<<endl;
957  evaluatetan=true;
958  ASSERTL0(wrong<2, "number edges to change is too high!!");
959  edgeinterp[wrong]=k;
960  wrong++;
961  }
962  }
963 
964  if(evaluatetan)
965  {
966  cout<<"tan bef"<<endl;
967  for(int e=0; e< nqedge; e++)
968  {
969  cout<<xcedgeQ[e]<<" "<<ycedgeQ[e]<<" "<<txedgeQ[e]<<endl;
970  }
971 
972  //OR: fit
973  int polyorder =3;
974  Array<OneD, NekDouble> coeffsinterp(polyorder+1);
975  Array<OneD, NekDouble> yPedges(npedge,0.0);
976  Array<OneD, NekDouble> xPedges(npedge,0.0);
977  Vmath::Vcopy(npedge, &xcPhysMOD[k*npedge+0],1,&xPedges[0],1);
978  Vmath::Vcopy(npedge, &ycPhysMOD[k*npedge+0],1,&yPedges[0],1);
979 
980  PolyFit(polyorder,nqedge, xcedgeQ,ycedgeQ, coeffsinterp, xPedges,yPedges, npedge);
981  //update values
982  Vmath::Vcopy(npedge, &xPedges[0],1, &xcPhysMOD[k*npedge+0],1);
983  Vmath::Vcopy(npedge, &yPedges[0],1, &ycPhysMOD[k*npedge+0],1);
984 
985  EvaluateTangent(nqedge,xcedgeQ, coeffsinterp,txedgeQ, tyedgeQ);
986 
987 
988  }
989  //copy also the tx,ty
990  Vmath::Vcopy(nqedge, &(txedgeQ[0]), 1, &(txQ[nqedge*k]),1);
991  Vmath::Vcopy(nqedge, &(tyedgeQ[0]), 1, &(tyQ[nqedge*k]),1);
992  }
993 
994  Array<OneD, NekDouble> fz(nedges*nqedge,0.0);
995  bndfieldx[lastIregion]->PhysDeriv(MultiRegions::eX,ycQ,fz);
996  for(int w=0; w< fz.num_elements(); w++)
997  {
998  txQ[w] = cos(atan(fz[w]));
999  tyQ[w] = sin(atan(fz[w]));
1000  cout<<xcQ[w]<<" "<<ycQ[w]<<" "<<fz[w]<<endl;
1001  }
1002  //ASSERTL0(false, "bobo");
1003  //force continuity tx,ty
1004  //tx
1005  Vmath::Vcopy(nquad_lay, txQ,1, Cont_y->UpdatePhys(),1);
1006  Cont_y->FwdTrans_IterPerExp(Cont_y->GetPhys(), coeffsy);
1007  Cont_y->BwdTrans_IterPerExp( coeffsy, Cont_y->UpdatePhys());
1008  Vmath::Vcopy(nquad_lay, Cont_y->GetPhys(),1, txQ,1);
1009  //ty
1010  Vmath::Vcopy(nquad_lay, tyQ,1, Cont_y->UpdatePhys(),1);
1011  Cont_y->FwdTrans_IterPerExp(Cont_y->GetPhys(), coeffsy);
1012  Cont_y->BwdTrans_IterPerExp( coeffsy, Cont_y->UpdatePhys());
1013  Vmath::Vcopy(nquad_lay, Cont_y->GetPhys(),1, tyQ,1);
1014 
1015  //check if the tan points have the same curvature otherwise interp
1016  NekDouble inc,incbefore;
1017 
1018  //build-up the fit for the tan using the edge with
1019  //the same derivative sign (before or after)
1020 
1021  int edgebef;
1022 
1023 
1024  for(int q=0; q<2; q++)
1025  {
1026  edgebef = edgeinterp[q]-1;
1027  incbefore = (txQ[edgebef*nqedge+nqedge-1]-txQ[edgebef*nqedge])/
1028  (xcQ[edgebef*nqedge+nqedge-1]-xcQ[edgebef*nqedge]);
1029  inc = (txQ[edgeinterp[q]*nqedge+nqedge-1]-txQ[edgeinterp[q]*nqedge])/
1030  (xcQ[edgeinterp[q]*nqedge+nqedge-1]-xcQ[edgeinterp[q]*nqedge]);
1031  int npoints = 2*nqedge;
1032 
1033  Array<OneD, NekDouble> yQedges(npoints,0.0);
1034  Array<OneD, NekDouble> xQedges(npoints,0.0);
1035  Array<OneD, NekDouble> tyQedges(npoints,0.0);
1036  Array<OneD, NekDouble> txQedges(npoints,0.0);
1037  cout<<"inc="<<inc<<" incbef="<<incbefore<<endl;
1038  if( (inc/incbefore)>0. )
1039  {
1040  cout<<"before!!"<<edgebef<<endl;
1041  int polyorder =2;
1042  Array<OneD, NekDouble> coeffsinterp(polyorder+1);
1043  Vmath::Vcopy(npoints, &xcQ[edgebef*nqedge+0],1,&xQedges[0],1);
1044  Vmath::Vcopy(npoints, &ycQ[edgebef*nqedge+0],1,&yQedges[0],1);
1045  Vmath::Vcopy(npoints, &txQ[edgebef*nqedge+0],1,&txQedges[0],1);
1046  Vmath::Vcopy(npoints, &tyQ[edgebef*nqedge+0],1,&tyQedges[0],1);
1047 
1048  PolyFit(polyorder, npoints,
1049  xQedges,txQedges,
1050  coeffsinterp, xQedges,txQedges, npoints);
1051 
1052  //copy back the values:
1053  Vmath::Vcopy(npoints,&txQedges[0],1, &txQ[edgebef*nqedge+0],1);
1054 
1055  PolyFit(polyorder, npoints,
1056  xQedges,tyQedges,
1057  coeffsinterp, xQedges,tyQedges, npoints);
1058 
1059  //copy back the values:
1060  Vmath::Vcopy(npoints,&tyQedges[0],1, &tyQ[edgebef*nqedge+0],1);
1061 
1062  }
1063  else
1064  {
1065  cout<<"after!!"<<endl;
1066  int polyorder =2;
1067  Array<OneD, NekDouble> coeffsinterp(polyorder+1);
1068  Vmath::Vcopy(npoints, &xcQ[edgeinterp[q]*nqedge+0],1,&xQedges[0],1);
1069  Vmath::Vcopy(npoints, &ycQ[edgeinterp[q]*nqedge+0],1,&yQedges[0],1);
1070  Vmath::Vcopy(npoints, &txQ[edgeinterp[q]*nqedge+0],1,&txQedges[0],1);
1071  Vmath::Vcopy(npoints, &tyQ[edgeinterp[q]*nqedge+0],1,&tyQedges[0],1);
1072 
1073 
1074  PolyFit(polyorder, npoints,
1075  xQedges,txQedges,
1076  coeffsinterp, xQedges,txQedges, npoints);
1077 
1078  //copy back the values:
1079  Vmath::Vcopy(npoints,&txQedges[0],1, &txQ[edgeinterp[q]*nqedge+0],1);
1080 
1081  PolyFit(polyorder, npoints,
1082  xQedges,tyQedges,
1083  coeffsinterp, xQedges,tyQedges, npoints);
1084 
1085  //copy back the values:
1086  Vmath::Vcopy(npoints,&tyQedges[0],1, &tyQ[edgeinterp[q]*nqedge+0],1);
1087 
1088 
1089  }
1090 
1091 
1092  }
1093  //force continuity of the tangent
1094  //tyQ
1095  Vmath::Vcopy(nquad_lay, tyQ,1, Cont_y->UpdatePhys(),1);
1096  Cont_y->FwdTrans_IterPerExp(Cont_y->GetPhys(), coeffstmp);
1097  Cont_y->BwdTrans_IterPerExp( coeffstmp, Cont_y->UpdatePhys());
1098  Vmath::Vcopy(nquad_lay, Cont_y->GetPhys(),1, tyQ,1);
1099  //txQ
1100  Vmath::Vcopy(nquad_lay, txQ,1, Cont_y->UpdatePhys(),1);
1101  Cont_y->FwdTrans_IterPerExp(Cont_y->GetPhys(), coeffstmp);
1102  Cont_y->BwdTrans_IterPerExp( coeffstmp, Cont_y->UpdatePhys());
1103  Vmath::Vcopy(nquad_lay, Cont_y->GetPhys(),1, txQ,1);
1104 
1105  for(int k=0; k<nedges; k++)
1106  {
1107  //determine the normal from eqs(APART FROM SIGN):
1108  //tx^2 +ty^2= 1 = nx^2 + ny^2;
1109  //t\cdot n=0= tx*nx +ty*ny
1110  //result: nx = ( 1+(tx/ty)^2 )^(-1/2)
1111 
1112 
1113  Vmath::Vcopy(nqedge, &(txQ[nqedge*k]),1, &(txedgeQ[0]), 1);
1114  Vmath::Vcopy(nqedge, &(tyQ[nqedge*k]),1, &(tyedgeQ[0]), 1);
1115 
1116  Vmath::Vdiv(nqedge, txedgeQ,1,tyedgeQ,1,tx_tyedgeQ,1);
1117  Vmath::Vmul(nqedge, tx_tyedgeQ,1,tx_tyedgeQ,1,tx_tyedgeQ,1);
1118  Vmath::Sadd(nqedge,1.0,tx_tyedgeQ,1,nxedgeQ,1);
1119  Vmath::Vsqrt(nqedge,nxedgeQ,1,nxedgeQ,1);
1120  Vmath::Sdiv(nqedge,1.0,nxedgeQ,1,nxedgeQ,1);
1121  //normal DOWNWARDS!!! mult by -1
1122  Vmath::Smul(nqedge, -1.0,nxedgeQ,1,nxedgeQ,1);
1123  Vmath::Vcopy(nqedge, &(nxedgeQ[0]),1, &(nxQ[nqedge*k]),1);
1124  //ny = (1-nx ^2)^(1/2)
1125  Vmath::Vmul(nqedge, nxedgeQ,1,nxedgeQ,1,nyedgeQ,1);
1126  Vmath::Smul(nqedge, -1.0,nyedgeQ,1,nyedgeQ,1);
1127  Vmath::Sadd(nqedge,1.0,nyedgeQ,1,nyedgeQ,1);
1128  Vmath::Vsqrt(nqedge,nyedgeQ,1,nyedgeQ,1);
1129  //normal DOWNWARDS!!! mult by -1
1130  Vmath::Smul(nqedge, -1.0,nyedgeQ,1,nyedgeQ,1);
1131  Vmath::Vcopy(nqedge, &(nyedgeQ[0]), 1, &(nyQ[nqedge*k]),1);
1132 
1133 
1134  cout<<"edge:"<<k<<endl;
1135  cout<<"tan/normal"<<endl;
1136  for(int r=0; r<nqedge; r++)
1137  {
1138  cout<<xcQ[k*nqedge+r]<<" "<<txedgeQ[r]<<" "<<tyedgeQ[r]<<" "
1139  <<nxedgeQ[r]<<" "<<nyedgeQ[r]<<endl;
1140  }
1141  }
1142 
1143  //force continuity:
1144  //REMEMBER: the Fwd/Bwd operation get wrong with the border values!!!
1145  Vmath::Vcopy(nquad_lay, nyQ,1, Cont_y->UpdatePhys(),1);
1146 
1147  Cont_y->FwdTrans_IterPerExp(Cont_y->GetPhys(), coeffstmp);
1148  Cont_y->BwdTrans_IterPerExp( coeffstmp, Cont_y->UpdatePhys());
1149  Vmath::Vcopy(nquad_lay, Cont_y->GetPhys(),1, nyQ,1);
1150 
1151  Vmath::Zero(nquad_lay,Cont_y->UpdatePhys(),1);
1152  Vmath::Zero(Cont_y->GetNcoeffs(),Cont_y->UpdateCoeffs(),1);
1153  Vmath::Vcopy(nquad_lay, nxQ,1, Cont_y->UpdatePhys(),1);
1154  Cont_y->FwdTrans_IterPerExp(Cont_y->GetPhys(), coeffstmp);
1155  Cont_y->BwdTrans_IterPerExp( coeffstmp, Cont_y->UpdatePhys());
1156  Vmath::Vcopy(nquad_lay, Cont_y->GetPhys(),1, nxQ,1);
1157 
1158  //force the normal at interface point to be equal
1159  for(int k=0; k<nedges; k++)
1160  {
1161  if(k>0)
1162  {
1163  //nyPhys[f*npedge +0] =
1164  // (nyPhys[(f-1)*npedge+npedge-1]+nyPhys[f*npedge+0])/2.;
1165  nyQ[(k-1)*nqedge+nqedge-1]=
1166  nyQ[k*nqedge+0];
1167  //nx= (1-ny^2)^{1/2}
1168  //nxPhys[f*npedge+0]=
1169  // sqrt(1- nyPhys[f*npedge+0]*nyPhys[f*npedge+0]);
1170  nxQ[(k-1)*nqedge+nqedge-1]=
1171  nxQ[k*nqedge+0];
1172  }
1173  }
1174 
1175  Array<OneD, NekDouble>x_tmpQ(nquad_lay-(nedges-1));
1176  //Array<OneD, NekDouble>tmpnxQ(nquad_lay-(nedges-1));
1177  Array<OneD, NekDouble>tmpnyQ(nquad_lay-(nedges-1));
1178 
1179  cout<<"nx,yQbefore"<<endl;
1180  for(int u=0; u<xcQ.num_elements(); u++)
1181  {
1182  cout<<xcQ[u]<<" "<<nyQ[u]<<" "<<txQ[u]<<endl;
1183  }
1184 
1185  Cutrepetitions(nedges, xcQ,x_tmpQ);
1186  //Cutrepetitions(nedges, nxQ, tmpnxQ);
1187  Cutrepetitions(nedges, nyQ, tmpnyQ);
1188  cout<<"nx,yQ"<<endl;
1189  for(int u=0; u<x_tmpQ.num_elements(); u++)
1190  {
1191  cout<<x_tmpQ[u]<<" "<<tmpnyQ[u]<<endl;
1192  }
1193 
1194  //interpolate the values into phys points(curved points)
1195  for(int k=0; k<nedges; k++)
1196  {
1197  //cout<<"edge:"<<k<<endl;
1198  for(int a=0; a<npedge; a++)
1199  {
1200  if(a==0)//verts pos no interp necessary
1201  {
1202  nxPhys[k*npedge +a]= nxQ[k*nqedge +0];
1203  nyPhys[k*npedge +a]= nyQ[k*nqedge +0];
1204 
1205  }
1206  else if(a== npedge-1)//verts pos no interp necessary
1207  {
1208  nxPhys[k*npedge +a]= nxQ[k*nqedge +nqedge-1];
1209  nyPhys[k*npedge +a]= nyQ[k*nqedge +nqedge-1];
1210  //cout<<":last"<<nyQ[k*nqedge+a]<<endl;
1211 
1212  }
1213  else
1214  {
1215  //use lagrange interpolant to get the
1216  //normal at phys(equispaced points)
1217 
1218  //order normal functions(cut out repetitions)
1219  //QUAD POINTS
1220 
1221  int index;
1222  //determine closest index:
1223 
1224  index=
1225  DetermineclosePointxindex( Addpointsx[k*(npedge-2) +a-1], x_tmpQ);
1226 
1227  Array<OneD, NekDouble> Pxinterp(4);
1228  Array<OneD, NekDouble> Pyinterp(4);
1229 
1230  //generate neighbour arrays (y):
1231  GenerateNeighbourArrays(index, 4,x_tmpQ,tmpnyQ,Pxinterp,Pyinterp);
1232  //interp the new normal components(y)
1233  nyPhys[k*npedge +a]=
1234  LagrangeInterpolant(Addpointsx[k*(npedge-2) +a-1],4,Pxinterp,Pyinterp );
1235  /*
1236  //generate neighbour arrays (x):
1237  GenerateNeighbourArrays(index,4,x_tmpQ,tmpnxQ,Pxinterp,Pyinterp);
1238  //interp the new normal components(x)
1239  nxPhys[k*npedge +a]=
1240  LagrangeInterpolant(Addpointsx[k*(npedge-2) +a],4,Pxinterp,Pyinterp );
1241  */
1242  //nx=-(1-ny*ny){1/2} the normal is DOWNWARDS!!!
1243  nxPhys[k*npedge +a]= -sqrt(abs(1- nyPhys[k*npedge +a]*nyPhys[k*npedge +a]));
1244  /*
1245  //(put the middle points as quad points)
1246  //GaussLobattoLegendre points in the middle:
1247  nxPhys[k*npedge +a] = nxedgeQ[a];
1248  nyPhys[k*npedge +a] = nyedgeQ[a];
1249  ASSERTL0(npedge< nqedge," quad points too low");
1250  */
1251  }
1252 
1253 
1254  //force the normal at interface point to be equal
1255  if(k>0)
1256  {
1257  //nyPhys[f*npedge +0] =
1258  // (nyPhys[(f-1)*npedge+npedge-1]+nyPhys[f*npedge+0])/2.;
1259  nyPhys[(k-1)*npedge+npedge-1]=
1260  nyPhys[k*npedge+0];
1261  //nx= (1-ny^2)^{1/2}
1262  //nxPhys[f*npedge+0]=
1263  // sqrt(1- nyPhys[f*npedge+0]*nyPhys[f*npedge+0]);
1264  nxPhys[(k-1)*npedge+npedge-1]=
1265  nxPhys[k*npedge+0];
1266  }
1267  }
1268  }
1269  cout<<"xcPhys,,"<<endl;
1270  for(int s=0; s<np_lay; s++)
1271  {
1272 
1273  cout<<xcPhysMOD[s]<<" "<<ycPhysMOD[s]<<" "<<nxPhys[s]<<" "<<nyPhys[s]<<endl;
1274 
1275  }
1276 
1277  //determine the new coords of the vertices and the curve points
1278  //for each edge
1279 
1280  //int np_lay = (nvertl-1)*npedge;//nedges*npedge
1281 
1282  //NB delta=ynew-y_c DEPENDS ON the coord trnaf ynew= y+Delta_c*ratio!!!
1283  Array<OneD, NekDouble> delta(nlays);
1284  Array<OneD, NekDouble>tmpy_lay(np_lay);
1285  Array<OneD, NekDouble>tmpx_lay(np_lay);
1286  for(int m=0; m<nlays; m++)
1287  {
1288  //delta[m] = (ynew[lay_Vids[m][0]] - y_c[0])/1.0;
1289 
1290  //depends on Delta0
1291  if(m< cntlow+1)
1292  {
1293  delta[m] = -(cntlow+1-m)*Delta0/(cntlow+1);
1294  }
1295  else
1296  {
1297  delta[m] = ( m-(cntlow) )*Delta0/(cntlow+1);
1298  }
1299 
1300 
1301  layers_x[m]= Array<OneD, NekDouble>(np_lay);
1302  //cout<<"delta="<<delta[m]<<" cntlow="<<cntlow<<endl;
1303 
1304  for(int h=0; h< nvertl; h++)
1305  {
1306  //shift using the dinstance delta from the crit layer AT x=0
1307  //for each layer
1308 
1309 
1310  //cout<<m<<"Vid:"<<lay_Vids[m][h]<<" mod from y="<<ynew[lay_Vids[m][h] ]<<" to y="<<y_c[h] +delta[m]<<endl;
1311  if(move_norm==false)
1312  {
1313  ynew[lay_Vids[m][h] ]= y_c[h] +delta[m];
1314  xnew[lay_Vids[m][h] ]= x_c[h];
1315  }
1316  else
1317  {
1318  if(h==0 || h==nvertl-1 )//nx=0,ny=1 at the borders
1319  {
1320  ynew[lay_Vids[m][h] ]= y_c[h] +delta[m];
1321  xnew[lay_Vids[m][h] ]= x_c[h];
1322  }
1323  else
1324  {
1325  ynew[lay_Vids[m][h] ]= y_c[h] +delta[m]*abs(nyPhys[h*npedge+0]);
1326  xnew[lay_Vids[m][h] ]= x_c[h] +delta[m]*abs(nxPhys[h*npedge+0]);
1327  }
1328  }
1329 cout<<"Vid x="<<xnew[lay_Vids[m][h] ]<<" y="<<ynew[lay_Vids[m][h] ]<<endl;
1330  if(h< nedges
1331  //&& curv_lay==true
1332  )
1333  {
1334 cout<<"edge=="<<h<<endl;
1335  if(h>0)//check normal consistency
1336  {
1337  ASSERTL0( nyPhys[h*npedge+0]==nyPhys[(h-1)*npedge+npedge-1]," normaly wrong");
1338  ASSERTL0( nxPhys[h*npedge+0]==nxPhys[(h-1)*npedge+npedge-1]," normalx wrong");
1339 
1340  }
1341  if(move_norm==false)
1342  {
1343  //v1
1344  layers_y[m][h*npedge +0] = y_c[h] +delta[m];
1345  layers_x[m][h*npedge +0] = xnew[lay_Vids[m][h] ];
1346  //v2
1347  layers_y[m][h*npedge +npedge-1] = y_c[h+1] +delta[m];
1348  layers_x[m][h*npedge +npedge-1] = xnew[lay_Vids[m][h+1] ];
1349  //middle points (shift crit lay points by delta):
1350  for(int d=0; d< npedge-2; d++)
1351  {
1352  layers_y[m][h*npedge +d+1]= ycPhysMOD[h*npedge +d+1] +delta[m];
1353  //Addpointsy[h*(npedge-2) +d] +delta[m];
1354  layers_x[m][h*npedge +d+1]= xcPhysMOD[h*npedge +d+1];
1355  //Addpointsx[h*(npedge-2) +d];
1356  }
1357  }
1358  else
1359  {
1360  if(h==0) //nx=0,ny=1 at the borders
1361  {
1362  //v1
1363  tmpy_lay[h*npedge +0] = y_c[h] +delta[m];
1364  tmpx_lay[h*npedge +0] = xnew[lay_Vids[m][h] ];
1365  //v2
1366  tmpy_lay[h*npedge +npedge-1] =
1367  y_c[h+1] +delta[m]*abs(nyPhys[h*npedge +npedge-1]);
1368  tmpx_lay[h*npedge +npedge-1] =
1369  x_c[h+1] +delta[m]*abs(nxPhys[h*npedge +npedge-1]);
1370  }
1371  else if(h==nedges-1)//nx=0,ny=1 at the borders
1372  {
1373  //v1
1374  tmpy_lay[h*npedge +0] =
1375  y_c[h] +delta[m]*abs(nyPhys[h*npedge +0]);
1376  tmpx_lay[h*npedge +0] =
1377  x_c[h] +delta[m]*abs(nxPhys[h*npedge +0]);
1378  //v2
1379  tmpy_lay[h*npedge +npedge-1] = y_c[h+1] +delta[m];
1380  tmpx_lay[h*npedge +npedge-1] = xnew[lay_Vids[m][h+1] ];
1381  }
1382  else
1383  {
1384  //v1
1385  tmpy_lay[h*npedge +0] =
1386  y_c[h] +delta[m]*abs(nyPhys[h*npedge +0]);
1387  tmpx_lay[h*npedge +0] =
1388  x_c[h] +delta[m]*abs(nxPhys[h*npedge +0]);
1389  //v2
1390  tmpy_lay[h*npedge +npedge-1] =
1391  y_c[h+1] +delta[m]*abs(nyPhys[h*npedge +npedge-1]);
1392  tmpx_lay[h*npedge +npedge-1] =
1393  x_c[h+1] +delta[m]*abs(nxPhys[h*npedge +npedge-1]);
1394  }
1395 
1396  //middle points
1397  for(int d=0; d< npedge-2; d++)
1398  {
1399 
1400  tmpy_lay[h*npedge +d+1] = ycPhysMOD[h*npedge +d+1] +
1401  delta[m]*abs(nyPhys[h*npedge +d+1]);
1402  //Addpointsy[h*(npedge-2) +d] +
1403  // delta[m]*abs(nyPhys[h*npedge +d+1]);
1404  tmpx_lay[h*npedge +d+1]= xcPhysMOD[h*npedge +d+1] +
1405  delta[m]*abs(nxPhys[h*npedge +d+1]);
1406  //Addpointsx[h*(npedge-2) +d] +
1407  // delta[m]*abs(nxPhys[h*npedge +d+1]);
1408 
1409  //NB ycQ,xcQ refers to nqedge NOT npedge!!!
1410  //tmpy_lay[h*npedge +d+1] = ycQ[h*nqedge +d+1] +
1411  // delta[m]*abs(nyPhys[h*npedge +d+1]);
1412  //tmpx_lay[h*npedge +d+1]= xcQ[h*nqedge +d+1] +
1413  // delta[m]*abs(nxPhys[h*npedge +d+1]);
1414  //cout<<"xmoved="<<tmpx_lay[h*npedge +d+1]<<" xold="<<xcQ[h*nqedge +d+1]<<endl;
1415  //ASSERTL0(tmpx_lay[h*npedge +d+1]>0," middle point with x<0")
1416 
1417  }
1418  }
1419  }//close edges
1420  }//close verts h
1421 
1422  for(int s=0; s<np_lay; s++)
1423  {
1424  cout<<tmpx_lay[s]<<" "<<tmpy_lay[s]<<endl;
1425  }
1426  Orderfunctionx(tmpx_lay,tmpy_lay, tmpx_lay, tmpy_lay);
1427  cout<<"fisrt interp"<<endl;
1428  for(int s=0; s<np_lay; s++)
1429  {
1430  cout<<tmpx_lay[s]<<" "<<tmpy_lay[s]<<endl;
1431  }
1432 
1433  //ASSERTL0(false, "dasd");
1434  //ASSERTL0(tmpx_lay[h*npedge +0]>=0," vert 0 x<0");
1435  //ASSERTL0(tmpx_lay[h*npedge +npedge-1]>0," vert 1 x<0");
1436 
1437  //check if the x coord is 'outofbound' and calculate the
1438  //number of outofbound points
1439 
1440  //determine which boudn has been overcome:
1441  NekDouble boundleft = xcPhysMOD[0];
1442  NekDouble boundright = xcPhysMOD[np_lay-1];
1443  bool outboundleft= false;
1444  bool outboundright=false;
1445  if(tmpx_lay[1]< boundleft )
1446  {
1447  outboundleft = true;
1448  }
1449  if(tmpx_lay[np_lay-2] > boundright )
1450  {
1451  outboundright = true;
1452  }
1453 
1454  int outvert=0;
1455  int outmiddle=0;
1456  int outcount=0;
1457 
1458  Array<OneD, int> vertout(nvertl);
1459  for(int r=0; r< nedges; r++)
1460  {
1461  //check point outofboundleft
1462  if(tmpx_lay[r*npedge + npedge-1]< boundleft && outboundleft==true )//assume the neg coords start from 0
1463  {
1464  vertout[outvert]=r;
1465  outvert++;
1466 
1467  if(r<nedges-1 )
1468  {
1469  //check if after the last negvert there are neg points
1470  if( tmpx_lay[(r+1)*npedge + npedge-1]> boundleft )
1471  {
1472  for(int s=0; s<npedge-2; s++)
1473  {
1474  if(tmpx_lay[(r+1)*npedge + s+1]< boundleft)
1475  {
1476  outmiddle++;
1477  }
1478  }
1479 
1480  }
1481  }
1482  }
1483  //check point outofboundright
1484  if(tmpx_lay[r*npedge + 0]> boundright && outboundright==true )//assume the neg coords start from 0
1485  {
1486  vertout[outvert]=r;
1487  outvert++;
1488 
1489  if( r> 0)
1490  {
1491  //check if after the last outvert there are out points
1492  if( tmpx_lay[(r-1)*npedge + 0]< boundright )
1493  {
1494  for(int s=0; s<npedge-2; s++)
1495  {
1496  if(tmpx_lay[(r-1)*npedge + s+1]> boundright)
1497  {
1498  outmiddle++;
1499  }
1500  }
1501 
1502  }
1503  }
1504  }
1505  }
1506 
1507  //calc number of point to replace
1508  outcount = outvert*npedge+1+ outmiddle;
1509  //determine from(until) which index the replacement will start
1510  int replacepointsfromindex=0;
1511  for(int c=0; c<nedges; c++)
1512  {
1513  //assume at least 1 middle point per edge
1514  if(xcPhysMOD[c*npedge+npedge-1] <= tmpx_lay[c*(npedge-(npedge-2)) +2] && outboundright==true)
1515  {
1516  replacepointsfromindex = c*(npedge-(npedge-2))+2;
1517  break;
1518  }
1519 
1520  //assume at least 1 middle point per edge
1521  if(xcPhysMOD[(nedges-1 -c)*npedge+0] >= tmpx_lay[np_lay-1 -(c*(npedge-(npedge-2)) +2)] && outboundleft==true)
1522  {
1523  replacepointsfromindex = np_lay-1 -(c*(npedge-(npedge-2)) +2);
1524  break;
1525  }
1526  }
1527 
1528  //cout<<"out="<<outcount<<endl;
1529  //cout<<"replacefrom="<<replacepointsfromindex<<endl;
1530  //if xcoord is neg find the first positive xcoord
1531 
1532  if(outcount>1)
1533  {
1534  //determine x new coords:
1535  //distribute the point all over the layer
1536  int pstart,shift;
1537  NekDouble increment;
1538 
1539  if( outboundright==true)
1540  {
1541  pstart = replacepointsfromindex;
1542  shift = np_lay-outcount;
1543  increment = (xcPhysMOD[np_lay-outcount]-xcPhysMOD[pstart])/(outcount+1);
1544  outcount = outcount-1;
1545  ASSERTL0(tmpx_lay[np_lay-outcount]>xcPhysMOD[(nedges-1)*npedge+0], "no middle points in the last edge");
1546  }
1547  else
1548  {
1549  shift=1;
1550  pstart= outcount-1;
1551  increment = (xcPhysMOD[replacepointsfromindex]-xcPhysMOD[pstart])/(outcount+1);
1552  ASSERTL0(tmpx_lay[pstart]<xcPhysMOD[0*npedge +npedge-1], "no middle points in the first edge");
1553  }
1554 
1555  //interp to points between posindex and posindex-1
1556  Array<OneD, NekDouble> replace_x(outcount);
1557  Array<OneD, NekDouble> replace_y(outcount);
1558  //order normal functions(cut out repetitions)
1559  Array<OneD, NekDouble>x_tmp(np_lay-(nedges-1));
1560  Array<OneD, NekDouble>y_tmp(np_lay-(nedges-1));
1561  Array<OneD, NekDouble>tmpny(np_lay-(nedges-1));
1562  Cutrepetitions(nedges, xcPhysMOD,x_tmp);
1563  Cutrepetitions(nedges, ycPhysMOD,y_tmp);
1564  Cutrepetitions(nedges, nyPhys, tmpny);
1565  //init neigh arrays
1566  Array<OneD, NekDouble>closex(4);
1567  Array<OneD, NekDouble>closey(4);
1568  Array<OneD, NekDouble>closeny(4);
1569  NekDouble xctmp,ycinterp,nxinterp,nyinterp;
1570 
1571  for(int v=0; v<outcount;v++)
1572  {
1573  xctmp = xcPhysMOD[pstart]+(v+1)*increment;
1574 
1575  //determine closest point index:
1576  int index =
1577  DetermineclosePointxindex( xctmp, x_tmp);
1578  //cout<<" vert="<<index<<endl;
1579 
1580  //generate neighbour arrays (ny)
1581  GenerateNeighbourArrays(index, 4,x_tmp,tmpny,closex,closeny);
1582 
1583  //interp:
1584  nyinterp =
1585  LagrangeInterpolant(
1586  xctmp,4,closex,closeny );
1587 
1588  //calc nxinterp
1589  nxinterp = sqrt(abs(1-nyinterp*nyinterp));
1590 
1591  //generata neighbour arrays (yc)
1592  GenerateNeighbourArrays(index, 4,x_tmp,y_tmp,closex,closey);
1593  //interp:
1594  ycinterp = LagrangeInterpolant(xctmp,4, closex,closey);
1595  //calc new coord
1596  replace_x[v] = xctmp +delta[m]*abs(nxinterp);
1597  replace_y[v] = ycinterp +delta[m]*abs(nyinterp);
1598  tmpx_lay[ v+shift ] = replace_x[v];
1599  tmpy_lay[ v+shift ] = replace_y[v];
1600 
1601  //cout<<"xinterp="<<replace_x[v]<<" yinterp="<<replace_y[v]<<endl;
1602  }
1603  }//end outcount if
1604 
1605 
1606 
1607  /*
1608  for(int s=0; s<np_lay; s++)
1609  {
1610 
1611  cout<<tmpx_lay[s]<<" "<<tmpy_lay[s]<<endl;
1612 
1613  }
1614  if(m== 0)
1615  {
1616  //ASSERTL0(false, "ssa");
1617  }
1618  */
1619 
1620  int closepoints = 4;
1621 
1622  Array<OneD, NekDouble> Pyinterp(closepoints);
1623  Array<OneD, NekDouble> Pxinterp(closepoints);
1624 
1625  //check if any edge has less than npedge points
1626  int pointscount=0;
1627  for(int q=0; q<np_lay; q++)
1628  {
1629  for(int e=0; e<nedges; e++)
1630  {
1631  if(tmpx_lay[q]<= x_c[e+1] && tmpx_lay[q]>= x_c[e])
1632  {
1633  pointscount++;
1634  }
1635  if(q == e*npedge +npedge-1 && pointscount!=npedge )
1636  {
1637  //cout<<"edge with few points :"<<e<<endl;
1638  pointscount=0;
1639  }
1640  else if(q == e*npedge +npedge-1)
1641  {
1642  pointscount=0;
1643  }
1644  }
1645  }
1646  //----------------------------------------------------------
1647  /*
1648  cout<<"notordered"<<endl;
1649  for(int g=0; g<tmpx_lay.num_elements(); g++)
1650  {
1651  cout<<tmpx_lay[g]<<" "<<tmpy_lay[g]<<endl;
1652  }
1653  */
1654 
1655  //cout<<nedges<<"nedges"<<npedge<<" np_lay="<<np_lay<<endl;
1656  //calc lay coords
1657  //MoveLayerNfixedxpos(nvertl, npedge, xcPhysMOD, tmpx_lay, tmpy_lay,
1658  // lay_Vids[m], layers_x[m], layers_y[m],xnew,ynew);
1659  MoveLayerNnormpos(nvertl, npedge, xcPhysMOD, tmpx_lay, tmpy_lay,
1660  lay_Vids[m], layers_x[m], layers_y[m],xnew,ynew);
1661 
1662  /*
1663  //generate x,y arrays without lastedgepoint
1664  //(needed to interp correctly)
1665  Array<OneD, NekDouble>tmpx(np_lay-(nedges-1));
1666  Array<OneD, NekDouble>tmpy(np_lay-(nedges-1));
1667 
1668  Cutrepetitions(nedges, tmpx_lay, tmpx);
1669  Cutrepetitions(nedges, tmpy_lay, tmpy);
1670  //order points in x:
1671  int index;
1672  Array<OneD, NekDouble> copyarray_x(tmpx.num_elements());
1673  Array<OneD, NekDouble> copyarray_y(tmpx.num_elements());
1674  Orderfunctionx(tmpx, tmpy, tmpx, tmpy);
1675 
1676  //determine the neighbour points (-3;+3)
1677  for(int g=0; g< nvertl; g++)
1678  {
1679  //verts
1680  //cout<<"determine value for vert x="<<x_c[g]<<endl;
1681  //determine closest index:
1682 
1683  index=
1684  DetermineclosePointxindex( x_c[g], tmpx);
1685  //generate neighbour arrays:
1686  GenerateNeighbourArrays(index, closepoints,tmpx,tmpy,Pxinterp,Pyinterp);
1687  //write vert coords
1688  ynew[lay_Vids[m][g] ]= LagrangeInterpolant(x_c[g],closepoints,Pxinterp,Pyinterp );
1689  xnew[lay_Vids[m][g] ]= x_c[g];
1690 
1691 
1692  if(g<nedges)
1693  {
1694 
1695  //v1
1696  layers_y[m][g*npedge +0] = ynew[lay_Vids[m][g] ];
1697  layers_x[m][g*npedge +0] = xnew[lay_Vids[m][g] ];
1698  //v2
1699 
1700  //determine closest index:
1701  index=
1702  DetermineclosePointxindex( x_c[g+1], tmpx);
1703  //generate neighbour arrays:
1704  GenerateNeighbourArrays(index, closepoints,tmpx,tmpy,Pxinterp,Pyinterp);
1705  layers_y[m][g*npedge +npedge-1] =
1706  LagrangeInterpolant(x_c[g+1],closepoints,Pxinterp,Pyinterp );
1707  layers_x[m][g*npedge +npedge-1] = x_c[g+1];
1708 
1709  //middle points
1710  for(int r=0; r< npedge-2; r++)
1711  {
1712  //determine closest point index:
1713  index =
1714  DetermineclosePointxindex( xcPhysMOD[g*npedge +r+1], tmpx);
1715  //cout<<" vert+"<<index<<endl;
1716 
1717  ASSERTL0( index<= tmpy.num_elements()-1, " index wrong");
1718  //generate neighbour arrays Pyinterp,Pxinterp
1719  GenerateNeighbourArrays(index, closepoints,tmpx,tmpy,Pxinterp,Pyinterp);
1720 
1721  layers_y[m][g*npedge +r+1]=
1722  LagrangeInterpolant(
1723  xcPhysMOD[g*npedge +r+1],closepoints,Pxinterp,Pyinterp );
1724  //cout<<"x value="<<xcPhysMOD[g*npedge +r+1]<<endl;
1725  layers_x[m][g*npedge +r+1]= xcPhysMOD[g*npedge +r+1];
1726 
1727  }
1728 
1729 
1730  }//if edge closed g
1731  }// nvertl closed
1732  */
1733  //check if there are points out of range:
1734  //cout<<Vmath::Vmax(np_lay,layers_y[m],1)<<endl;
1735  if(curv_lay==true)
1736  {
1737  //ASSERTL0(Vmath::Vmax(np_lay,layers_y[m],1)< Vmath::Vmax(nVertTot,yold,1),"point>ymax");
1738  //ASSERTL0(Vmath::Vmin(np_lay,layers_y[m],1)> Vmath::Vmin(nVertTot,yold,1),"point<ymin");
1739  }
1740 
1741 
1742  //force Polycontinuity of the layer(3 order poly every 2 edges):
1743  int npoints = npedge;
1744  Array<OneD, NekDouble> xPedges(npoints);
1745  Array<OneD, NekDouble> yPedges(npoints);
1746  for(int f=0; f<nedges; f++)
1747  {
1748  int polyorder=2;
1749 
1750  Array<OneD, NekDouble> coeffsinterp(polyorder+1);
1751 
1752  Vmath::Vcopy(npoints, &layers_x[m][(f)*npedge+0],1,&xPedges[0],1);
1753  Vmath::Vcopy(npoints, &layers_y[m][(f)*npedge+0],1,&yPedges[0],1);
1754 
1755  PolyFit(polyorder, npoints,
1756  xPedges,yPedges,
1757  coeffsinterp, xPedges,yPedges, npoints);
1758 
1759  //copy back the values:
1760  Vmath::Vcopy(npoints,&yPedges[0],1, &layers_y[m][(f)*npedge+0],1);
1761 
1762  //verts still the same:
1763  layers_y[m][f*npedge+0]= ynew[lay_Vids[m][f]];
1764  layers_y[m][f*npedge+npedge-1]= ynew[lay_Vids[m][f+1]];
1765  }
1766 
1767  cout<<" xlay ylay lay:"<<m<<endl;
1768  for(int l=0; l<np_lay; l++)
1769  {
1770  //cout<<tmpx_lay[l]<<" "<<tmpy_lay[l]<<endl;
1771  cout<<std::setprecision(8)<<layers_x[m][l]<<" "<<layers_y[m][l]<<endl;
1772  }
1773  /*
1774  cout<<"nverts"<<endl;
1775  for(int l=0; l<nvertl; l++)
1776  {
1777  cout<<std::setprecision(8)<<xnew[lay_Vids[m][l] ]<<" "<<ynew[lay_Vids[m][l] ]<<endl;
1778  }
1779  */
1780 
1781  //ASSERTL0(false, "as");
1782 
1783  //if the layers coords are too steep use two edges verts to get an
1784  //poly interp third order
1785  /*
1786  bool polyinterp=false;
1787  for(int b=0; b< nedges; b++)
1788  {
1789  for(int u=0; u<npedge; u++)
1790  {
1791  if(
1792  abs(layers_y[m][b*npedge+u+1]- layers_y[m][b*npedge +u])/(layers_x[m][b*npedge+u+1]-layers_x[m][b*npedge+u]))>4.0 )
1793  {
1794  polyinterp=true;
1795  break;
1796  }
1797  cout<<"incratio="<<incratio<<endl;
1798  }
1799  //
1800  //Evaluatelay_tan
1801 
1802  }
1803  */
1804 
1805  cout<<"lay="<<m<<endl;
1806  ASSERTL0(Vmath::Vmin(nVertTot, yold,1)< Vmath::Vmin(np_lay,layers_y[m],1),
1807  " different layer ymin val");
1808  ASSERTL0(Vmath::Vmax(nVertTot, yold,1)> Vmath::Vmax(np_lay,layers_y[m],1),
1809  " different layer ymax val");
1810  ASSERTL0(Vmath::Vmin(nVertTot, xold,1)== Vmath::Vmin(np_lay,layers_x[m],1),
1811  " different layer xmin val");
1812  ASSERTL0(Vmath::Vmax(nVertTot, xold,1)== Vmath::Vmax(np_lay,layers_x[m],1),
1813  " different layer xmax val");
1814 
1815  }//close layers!!! m index
1816 
1817  //MoveOutsidePointsfixedxpos(npedge, graphShPt,xold_c, yold_c, xold_low, yold_low,
1818  // xold_up, yold_up, layers_y[0], layers_y[nlays-1], xnew, ynew);
1819 
1820  //lastIregion -1 = laydown
1821  //lastIregion -2 = layup
1822  MoveOutsidePointsNnormpos(npedge, graphShPt, xold_c,yold_c, xold_low,yold_low,xold_up,yold_up,
1823  layers_x[0], layers_y[0], layers_x[nlays-1], layers_y[nlays-1],nxPhys, nyPhys,xnew, ynew);
1824 
1825 
1826 
1827  /*
1828  //update vertices coords outside layers region
1829  NekDouble ratio;
1830  for(int n=0; n<nVertTot; n++)
1831  {
1832  NekDouble ratio;
1833  SpatialDomains::PointGeomSharedPtr vertex = graphShPt->GetVertex(n);
1834  NekDouble x,y,z;
1835  vertex->GetCoords(x,y,z);
1836  int qp_closer;
1837  NekDouble diff;
1838  int qp_closerup, qp_closerdown;
1839  NekDouble diffup, diffdown;
1840  //determine the closer xold_up
1841  NekDouble tmp=1000;
1842  diffup =1000;
1843  diffdown = 1000;
1844  for(int k=0; k<nvertl; k++)
1845  {
1846  if(abs(x-xold_c[k]) < tmp)
1847  {
1848  tmp = abs(x-xold_c[k]);
1849  qp_closer=k;
1850  }
1851  }
1852 
1853  //find nplay_closer
1854  int nplay_closer;
1855  if(qp_closer==0)
1856  {
1857  nplay_closer=0;//first vert
1858  }
1859  else
1860  {
1861  nplay_closer= (qp_closer-1)*npedge +npedge-1;
1862  }
1863 
1864  if( y>yold_up[qp_closer] && y<1 )//nlays-1 is layerup
1865  {
1866 
1867  // ratio = (1-layers_y[nlays-1][qp_closer])*(1-y_c[qp_closer])/
1868  // ( (1-yold_up[n])*(1-yold_c[qp_closer]) );
1869  ratio = (1-layers_y[nlays-1][nplay_closer])/
1870  ( (1-yold_up[qp_closer]) );
1871  //distance prop to layerup
1872  ynew[n] = layers_y[nlays-1][nplay_closer]
1873  + (y-yold_up[qp_closer])*ratio;
1874  xnew[n] = x;
1875 
1876  }
1877  else if( y< yold_low[qp_closer] && y>-1 )//0 is layerdown
1878  {
1879 
1880  ratio = (1+layers_y[0][nplay_closer])/
1881  ( (1+yold_low[qp_closer]) );
1882  //distance prop to layerlow
1883  ynew[n] = layers_y[0][nplay_closer]
1884  + (y-yold_low[qp_closer])*ratio;
1885  xnew[n] = x;
1886  }
1887 
1888  }
1889  */
1890 
1891  /*
1892  //update verts coords of critlay(in case EvaluateLayTnaget has been called)
1893  //it's not necessary !!!
1894  for(int e=0; e<nvertl; e++)
1895  {
1896  ynew[CritLay[e]] = y_c[e];
1897  }
1898  */
1899 
1900 
1901  }//move_norm bool
1902  else//move vertically
1903  {
1904  MoveLayersvertically(nlays, nvertl, cntlow, cntup,
1905  lay_Vids, x_c, y_c, Down, Up, xnew, ynew, layers_x, layers_y);
1906 
1907  }
1908 
1909 
1910 //check singular quads:
1911 CheckSingularQuads(streak, V1, V2,xnew, ynew);
1912 //check borders of the new mesh verts:
1913 //cout<<std::setprecision(8)<<"yoldmax="<<Vmath::Vmax(nVertTot, yold,1)<<endl;
1914 //cout<<std::setprecision(8)<<"ynewmax="<<Vmath::Vmax(nVertTot,ynew,1)<<endl;
1915 
1916 cout<<std::setprecision(8)<<"xmin="<<Vmath::Vmin(nVertTot, xnew,1)<<endl;
1917  ASSERTL0(Vmath::Vmin(nVertTot, xold,1)== Vmath::Vmin(nVertTot,xnew,1),
1918  " different xmin val");
1919  ASSERTL0(Vmath::Vmin(nVertTot, yold,1)== Vmath::Vmin(nVertTot,ynew,1),
1920  " different ymin val");
1921  ASSERTL0(Vmath::Vmax(nVertTot, xold,1)== Vmath::Vmax(nVertTot,xnew,1),
1922  " different xmax val");
1923  ASSERTL0(Vmath::Vmax(nVertTot, yold,1)== Vmath::Vmax(nVertTot,ynew,1),
1924  " different ymax val");
1925 
1926 
1927 
1928  //replace the vertices with the new ones
1929 
1930  Replacevertices(changefile, xnew , ynew, xcPhys, ycPhys, Eids, npedge, charalp, layers_x,layers_y, lay_Eids, curv_lay);
1931 
1932 
1933 }
1934 
1935 
1936 void OrderVertices (int nedges, SpatialDomains::MeshGraphSharedPtr graphShPt,
1937  MultiRegions::ExpListSharedPtr & bndfield,
1938  Array<OneD, int>& Vids, int v1,int v2, NekDouble x_connect, int & lastedge,
1939  Array<OneD, NekDouble>& x, Array<OneD, NekDouble>& y)
1940 {
1941  int nvertl = nedges+1;
1942  int edge;
1943  Array<OneD, int> Vids_temp(nvertl,-10);
1944  NekDouble x0layer = 0.0;
1945  for(int j=0; j<nedges; j++)
1946  {
1947  LocalRegions::SegExpSharedPtr bndSegExplow =
1948  bndfield->GetExp(j)->as<LocalRegions::SegExp>();
1949  edge = (bndSegExplow->GetGeom1D())->GetEid();
1950 //cout<<" edge="<<edge<<endl;
1951  for(int k=0; k<2; k++)
1952  {
1953  Vids_temp[j+k]=(bndSegExplow->GetGeom1D())->GetVid(k);
1954  SpatialDomains::PointGeomSharedPtr vertex = graphShPt->GetVertex(Vids_temp[j+k]);
1955  NekDouble x1,y1,z1;
1956  vertex->GetCoords(x1,y1,z1);
1957  if(x1==x_connect && edge!=lastedge)
1958  {
1959  //first 2 points
1960  if(x_connect==x0layer)
1961  {
1962  Vids[v1]=Vids_temp[j+k];
1963  x[v1]=x1;
1964  y[v1]=y1;
1965 
1966  if(k==0 )
1967  {
1968  Vids_temp[j+1]=(bndSegExplow->GetGeom1D())->GetVid(1);
1969  Vids[v2]=Vids_temp[j+1];
1970  SpatialDomains::PointGeomSharedPtr vertex = graphShPt->GetVertex(Vids[v2]);
1971  NekDouble x2,y2,z2;
1972  vertex->GetCoords(x2,y2,z2);
1973  x[v2]=x2;
1974  y[v2]=y2;
1975  }
1976  if(k==1)
1977  {
1978  Vids_temp[j+0]=(bndSegExplow->GetGeom1D())->GetVid(0);
1979  Vids[v2]=Vids_temp[j+0];
1980  SpatialDomains::PointGeomSharedPtr vertex = graphShPt->GetVertex(Vids[v2]);
1981  NekDouble x2,y2,z2;
1982  vertex->GetCoords(x2,y2,z2);
1983  x[v2]=x2;
1984  y[v2]=y2;
1985  }
1986  }
1987  else //other vertices
1988  {
1989  if(k==0 )
1990  {
1991  Vids_temp[j+1]=(bndSegExplow->GetGeom1D())->GetVid(1);
1992  Vids[v1]=Vids_temp[j+1];
1993  SpatialDomains::PointGeomSharedPtr vertex = graphShPt->GetVertex(Vids[v1]);
1994  NekDouble x1,y1,z1;
1995  vertex->GetCoords(x1,y1,z1);
1996  x[v1]=x1;
1997  y[v1]=y1;
1998  }
1999  if(k==1)
2000  {
2001  Vids_temp[j+0]=(bndSegExplow->GetGeom1D())->GetVid(0);
2002  Vids[v1]=Vids_temp[j+0];
2003  SpatialDomains::PointGeomSharedPtr vertex = graphShPt->GetVertex(Vids[v1]);
2004  NekDouble x1,y1,z1;
2005  vertex->GetCoords(x1,y1,z1);
2006  x[v1]=x1;
2007  y[v1]=y1;
2008  }
2009  }
2010  break;
2011  }
2012  }
2013  if(Vids[v1]!=-10)
2014  {
2015  lastedge = edge;
2016  break;
2017  }
2018  }
2019 
2020 }
2021 
2022 
2023 void Computestreakpositions(int npoints, MultiRegions::ExpListSharedPtr streak,
2024  Array<OneD, NekDouble> xold_up, Array<OneD, NekDouble> yold_up,
2025  Array<OneD, NekDouble> xold_low, Array<OneD, NekDouble> yold_low,
2026  Array<OneD, NekDouble> xold_c, Array<OneD, NekDouble> yold_c,
2027  Array<OneD, NekDouble> &xc, Array<OneD, NekDouble> &yc, NekDouble cr,
2028  bool verts)
2029 {
2030 cout<<"Computestreakpositions"<<endl;
2031  int nq = streak->GetTotPoints();
2032  Array<OneD, NekDouble> coord(2);
2033  //Array<OneD, NekDouble> stvalues(nvertl,-10);
2034  Array<OneD, NekDouble> derstreak(nq);
2035  streak->PhysDeriv(MultiRegions::eY, streak->GetPhys(), derstreak);
2036  int elmtid, offset;
2037  NekDouble U,dU;
2038  NekDouble F=1000;
2039 
2040  if(verts==true)//only for verts makes sense to init the coord values..
2041  {
2042 
2043  //start guess
2044  //yc= (yup+ydown)/2
2045  Vmath::Vadd(xc.num_elements(), yold_up,1,yold_low,1, yc,1);
2046  Vmath::Smul(xc.num_elements(), 0.5,yc,1,yc,1);
2047  Vmath::Vcopy(xc.num_elements(),xold_c,1,xc,1);
2048  }
2049  else//case of xQ,yQ
2050  {
2051  Vmath::Vcopy(xc.num_elements(), xold_c,1,xc,1);
2052  Vmath::Vcopy(xc.num_elements(), yold_c,1,yc,1);
2053  }
2054 
2055  int its;
2056  int attempt;
2057  NekDouble tol = 1e-3;
2058  NekDouble ytmp;
2059  for(int e=0; e<npoints; e++)
2060  {
2061  coord[0] =xc[e];
2062  coord[1] =yc[e];
2063 
2064  elmtid = streak->GetExpIndex(coord,0.00001);
2065  offset = streak->GetPhys_Offset(elmtid);
2066  F = 1000;
2067  its=0;
2068  attempt=0;
2069  ytmp = coord[1];
2070 //cout<<"start guess: x="<<xc[e]<<" y="<<yc[e]<<endl;
2071  while( abs(F)> 0.000000001)
2072  {
2073 
2074  elmtid = streak->GetExpIndex(coord,0.00001);
2075 
2076  //stvalues[e] = streak->GetExp(elmtid)->PhysEvaluate(coord, streak->GetPhys() +offset );
2077 //cout<<"elmtid="<<elmtid<<" x="<<coord[0]<<" y="<<coord[1]<<" stvalue="<<U<<" dU="<<dU<<endl;
2078 
2079  if( (abs(coord[1])>1 || elmtid==-1)
2080  && attempt==0 && verts==true
2081  )
2082  {
2083  //try the old crit lay position:
2084  coord[1] = yold_c[e];
2085  attempt++;
2086  }
2087  else if( (abs(coord[1])>1 || elmtid==-1) )
2088  {
2089  coord[1] = ytmp +0.01;
2090  elmtid = streak->GetExpIndex(coord,0.001);
2091  offset = streak->GetPhys_Offset(elmtid);
2092  NekDouble Utmp = streak->GetExp(elmtid)->PhysEvaluate(coord, streak->GetPhys() + offset);
2093  NekDouble dUtmp = streak->GetExp(elmtid)->PhysEvaluate(coord, derstreak + offset);
2094  coord[1] = coord[1] - (Utmp-cr)/dUtmp;
2095  if( (abs(Utmp-cr)>abs(F))||(abs(coord[1])>1) )
2096  {
2097  coord[1] = ytmp -0.01;
2098  }
2099 
2100  attempt++;
2101  }
2102  else
2103  {
2104  ASSERTL0(abs(coord[1])<= 1, " y value out of bound +/-1");
2105  }
2106  offset = streak->GetPhys_Offset(elmtid);
2107  U = streak->GetExp(elmtid)->PhysEvaluate(coord, streak->GetPhys() + offset);
2108  dU = streak->GetExp(elmtid)->PhysEvaluate(coord, derstreak + offset);
2109  coord[1] = coord[1] - (U-cr)/dU;
2110  F = U-cr;
2111  ASSERTL0( coord[0]==xc[e], " x coordinate must remain the same");
2112 
2113  its++;
2114  if(its>200 && abs(F)<0.00001)
2115  {
2116  cout<<"warning streak position obtained with precision:"<<F<<endl;
2117  break;
2118  }
2119  else if(its>1000 && abs(F)< 0.0001)
2120  {
2121  cout<<"warning streak position obtained with precision:"<<F<<endl;
2122  break;
2123  }
2124  else if(its>1000)
2125  {
2126  ASSERTL0(false, "no convergence after 1000 iterations");
2127  }
2128  }
2129  yc[e] = coord[1] - (U-cr)/dU;
2130  ASSERTL0( U<= cr + tol, "streak wrong+");
2131  ASSERTL0( U>= cr -tol, "streak wrong-");
2132  //Utilities::Zerofunction(coord[0], coord[1], xtest, ytest, streak, derstreak);
2133 cout<<"result streakvert x="<<xc[e]<<" y="<<yc[e]<<" streak="<<U<<endl;
2134  }
2135 
2136 }
2137 
2138 
2139 void GenerateAddPointsNewtonIt( NekDouble xi, NekDouble yi,NekDouble &xout, NekDouble &yout,
2140  MultiRegions::ExpListSharedPtr function, Array<OneD, NekDouble> derfunction,
2141  NekDouble cr)
2142 {
2143  int elmtid,offset;
2144  NekDouble F,U,dU;
2145  Array<OneD, NekDouble> coords(2);
2146 
2147  coords[0] = xi;
2148  coords[1] = yi;
2149  F =1000;
2150  int attempt=0;
2151  int its=0;
2152  NekDouble ytmp;
2153  ytmp = coords[1];
2154  while( abs(F)> 0.00000001)
2155  {
2156 
2157 //cout<<"generate newton it xi="<<xi<<" yi="<<yi<<endl;
2158  elmtid = function->GetExpIndex(coords, 0.01);
2159  //@to do if GetType(elmtid)==triangular WRONG!!!
2160 cout<<"gen newton xi="<<xi<<" yi="<<coords[1]<<" elmtid="<<elmtid<<" F="<<F<<endl;
2161 
2162  if( (abs(coords[1])>1 || elmtid==-1) )
2163  {
2164 
2165  coords[1] = ytmp +0.01;
2166  elmtid = function->GetExpIndex(coords,0.01);
2167  offset = function->GetPhys_Offset(elmtid);
2168  NekDouble Utmp = function->GetExp(elmtid)->PhysEvaluate(coords, function->GetPhys() + offset);
2169  NekDouble dUtmp = function->GetExp(elmtid)->PhysEvaluate(coords, derfunction + offset);
2170  coords[1] = coords[1] - (Utmp-cr)/dUtmp;
2171 cout<<"attempt:"<<coords[1]<<endl;
2172  if( (abs(Utmp-cr)>abs(F))||(abs(coords[1])>1.01) )
2173  {
2174  coords[1] = ytmp -0.01;
2175  }
2176 
2177  attempt++;
2178  }
2179  else if( abs(coords[1])<1.01 &&attempt==0)
2180  {
2181  coords[1] =1.0;
2182  attempt++;
2183  }
2184  else
2185  {
2186  ASSERTL0(abs(coords[1])<= 1.00, " y value out of bound +/-1");
2187  }
2188  offset = function->GetPhys_Offset(elmtid);
2189  U = function->GetExp(elmtid)->PhysEvaluate(coords, function->GetPhys() + offset);
2190  dU = function->GetExp(elmtid)->PhysEvaluate(coords, derfunction + offset);
2191  coords[1] = coords[1] - (U-cr)/dU;
2192 cout<<cr<<"U-cr="<<U-cr<<" tmp result y:"<<coords[1]<<" dU="<<dU<<endl;
2193  F = U-cr;
2194 
2195  its++;
2196  if(its>200 && abs(F)<0.00001)
2197  {
2198  cout<<"warning streak position obtained with precision:"<<F<<endl;
2199  break;
2200  }
2201  else if(its>1000 && abs(F)< 0.0001)
2202  {
2203  cout<<"warning streak position obtained with precision:"<<F<<endl;
2204  break;
2205  }
2206  else if(its>1000)
2207  {
2208  ASSERTL0(false, "no convergence after 1000 iterations");
2209  }
2210 
2211 
2212  ASSERTL0( coords[0]==xi, " x coordinate must remain the same");
2213  }
2214  xout = xi;
2215  yout = coords[1] - (U-cr)/dU;
2216 cout<<"NewtonIt result x="<<xout<<" y="<<coords[1]<<" U="<<U<<endl;
2217 }
2218 
2219 void GenerateMapEidsv1v2(MultiRegions::ExpListSharedPtr field,
2220  Array<OneD, int> &V1, Array<OneD, int> &V2)
2221 {
2222 
2223  const boost::shared_ptr<LocalRegions::ExpansionVector> exp2D = field->GetExp();
2224  int nel = exp2D->size();
2225  LocalRegions::QuadExpSharedPtr locQuadExp;
2226  LocalRegions::TriExpSharedPtr locTriExp;
2227  SpatialDomains::Geometry1DSharedPtr SegGeom;
2228  int id;
2229  int cnt=0;
2230  Array<OneD, int> V1tmp(4*nel, 10000);
2231  Array<OneD, int> V2tmp(4*nel, 10000);
2232  for(int i=0; i<nel; i++)
2233  {
2234  if((locQuadExp = (*exp2D)[i]->as<LocalRegions::QuadExp>()))
2235  {
2236  for(int j = 0; j < locQuadExp->GetNedges(); ++j)
2237  {
2238  SegGeom = (locQuadExp->GetGeom2D())->GetEdge(j);
2239  id = SegGeom->GetEid();
2240  if( V1tmp[id] == 10000)
2241  {
2242  V1tmp[id]= SegGeom->GetVertex(0)->GetVid();
2243  V2tmp[id]= SegGeom->GetVertex(1)->GetVid();
2244  cnt++;
2245  }
2246  }
2247  }
2248  //in the future the tri edges may be not necessary (if the nedges is known)
2249 
2250  else if((locTriExp = (*exp2D)[i]->as<LocalRegions::TriExp>()))
2251  {
2252  for(int j = 0; j < locTriExp->GetNedges(); ++j)
2253  {
2254  SegGeom = (locTriExp->GetGeom2D())->GetEdge(j);
2255  id = SegGeom->GetEid();
2256 
2257  if( V1tmp[id] == 10000)
2258  {
2259  V1tmp[id]= SegGeom->GetVertex(0)->GetVid();
2260  V2tmp[id]= SegGeom->GetVertex(1)->GetVid();
2261  cnt++;
2262  }
2263  }
2264 
2265  }
2266 
2267  }
2268 
2269  V1 = Array<OneD, int>(cnt);
2270  V2 = Array<OneD, int>(cnt);
2271  //populate the output vectors V1,V2
2272  for(int g=0; g<cnt; g++)
2273  {
2274  V1[g] = V1tmp[g];
2275  ASSERTL0(V1tmp[g]!=10000, "V1 wrong");
2276  V2[g] = V2tmp[g];
2277  ASSERTL0(V2tmp[g]!=10000, "V2 wrong");
2278  }
2279 
2280 
2281 
2282 
2283 }
2284 
2285 
2286 
2287 
2288 
2289 void MappingEVids(Array<OneD, NekDouble> xoldup, Array<OneD, NekDouble> yoldup,
2290  Array<OneD, NekDouble> xolddown, Array<OneD, NekDouble> yolddown,
2291  Array<OneD, NekDouble> xcold, Array<OneD, NekDouble> ycold,
2292  Array<OneD, int> Vids_c,
2293  SpatialDomains::MeshGraphSharedPtr mesh,
2294  MultiRegions::ExpListSharedPtr streak,
2295  Array<OneD, int> V1, Array<OneD, int> V2,
2296  int & nlays, Array<OneD, Array<OneD, int> >& Eids_lay,
2297  Array<OneD, Array<OneD, int> >& Vids_lay)
2298 {
2299 
2300  int nlay_Eids = xcold.num_elements()-1;
2301  int nlay_Vids = xcold.num_elements();
2302 
2303  int nVertsTot = mesh->GetNvertices();
2304 cout<<"nverttot="<<nVertsTot<<endl;
2305  //find the first vert of each layer
2306  //int qp_closerup,qp_closerdown;
2307  //NekDouble diffup,diffdown;
2308 cout<<"init nlays="<<nlays<<endl;
2309  //tmp first vert array (of course <100!!)
2310  Array<OneD, int> tmpVids0(100);
2311  Array<OneD, NekDouble> tmpx0(100);
2312  Array<OneD, NekDouble> tmpy0(100);
2313  Array<OneD, NekDouble> x2D(nVertsTot);
2314  Array<OneD, NekDouble> y2D(nVertsTot);
2315 cout<<"yoldup="<<yoldup[0]<<endl;
2316 cout<<"yolddown="<<yolddown[0]<<endl;
2317 
2318  for(int r=0; r< nVertsTot; r++)
2319  {
2320 
2321  SpatialDomains::PointGeomSharedPtr vertex = mesh->GetVertex(r);
2322  NekDouble x,y,z;
2323  vertex->GetCoords(x,y,z);
2324  x2D[r] = x;
2325  y2D[r] = y;
2326  if( xcold[0]==x )
2327  {
2328  //check if the vert is inside layer zone
2329  if(
2330  y<= yoldup[0] && y>= yolddown[0]
2331  && y!= ycold[0]
2332  )
2333  {
2334 //cout<<"x="<<x<<" y="<<y<<endl;
2335  tmpVids0[nlays] =r;
2336  tmpx0[nlays] = x;
2337  tmpy0[nlays] = y;
2338  nlays++;
2339  }
2340 
2341  }
2342  }
2343 cout<<"nlays="<<nlays<<endl;
2344 
2345  //reorder first verts from xlower to xhigher
2346  Array<OneD, NekDouble> tmpx(nlays);
2347  Array<OneD, NekDouble> tmpf(nlays);
2348  Array<OneD, int> tmpV(nlays);
2349  //local copy to prevent overwriting
2350  Vmath::Vcopy(nlays, tmpx0,1,tmpx,1);
2351  Vmath::Vcopy(nlays, tmpy0,1,tmpf,1);
2352  Vmath::Vcopy(nlays, tmpVids0,1,tmpV,1);
2353  NekDouble max = Vmath::Vmax(tmpf.num_elements(), tmpf,1);
2354  int index=0;
2355  for(int w=0; w< nlays; w++)
2356  {
2357  index = Vmath::Imin(tmpf.num_elements(), tmpf,1);
2358  tmpx0[w]= tmpx[index];
2359  tmpy0[w]= tmpf[index];
2360  tmpVids0[w] = tmpV[index];
2361  tmpf[index] = max+1000;
2362 
2363  }
2364 
2365 
2366 
2367 
2368 
2369  //initialise layers arrays
2370  Eids_lay = Array<OneD, Array<OneD,int> >(nlays);
2371  Vids_lay = Array<OneD, Array<OneD,int> >(nlays);
2372  for(int m=0; m<nlays; m++)
2373  {
2374  Eids_lay[m] = Array<OneD, int> (nlay_Eids);
2375  Vids_lay[m] = Array<OneD, int> (nlay_Vids);
2376  }
2377 
2378  //determine the others verts and edge for each layer
2379  NekDouble normbef = 0.0;
2380  NekDouble normtmp = 0.0;
2381  NekDouble xbef = 0.0;
2382  NekDouble ybef=0.0;
2383  NekDouble xtmp,ytmp,normnext,xnext,ynext,diff;
2384  NekDouble Ubef = 0.0, Utmp = 0.0, Unext = 0.0;
2385  Array<OneD, NekDouble> coord(2);
2386  int elmtid,offset;
2387  int nTotEdges = V1.num_elements();
2388  Array<OneD, int> edgestmp(6);
2389  for(int m=0; m<nlays; m++)
2390  {
2391  for(int g=0; g<nlay_Eids; g++)
2392  {
2393  if(g==0)
2394  {
2395  for(int h=0; h< nTotEdges; h++)
2396  {
2397 
2398  if( tmpVids0[m]== V1[h] )
2399  {
2400  SpatialDomains::PointGeomSharedPtr vertex = mesh->GetVertex(V2[h]);
2401  NekDouble x,y,z;
2402  vertex->GetCoords(x,y,z);
2403  if(x!=tmpx0[m])
2404  {
2405  Eids_lay[m][0] = h;
2406  Vids_lay[m][0] = V1[h];
2407  Vids_lay[m][1] = V2[h];
2408  SpatialDomains::PointGeomSharedPtr vertex1
2409  = mesh->GetVertex(V1[h]);
2410  NekDouble x1,y1,z1;
2411  vertex1->GetCoords(x1,y1,z1);
2412  normbef= sqrt( (y-y1)*(y-y1)+(x-x1)*(x-x1) );
2413  ybef = (y-y1);
2414  xbef = (x-x1);
2415  coord[0] = x;
2416  coord[1] = y;
2417  elmtid = streak->GetExpIndex(coord,0.00001);
2418  offset = streak->GetPhys_Offset(elmtid);
2419  Ubef = streak->GetExp(elmtid)->PhysEvaluate(coord, streak->GetPhys() + offset);
2420 
2421  }
2422 
2423  }
2424  if( tmpVids0[m]== V2[h] )
2425  {
2426  SpatialDomains::PointGeomSharedPtr vertex = mesh->GetVertex(V1[h]);
2427  NekDouble x,y,z;
2428  vertex->GetCoords(x,y,z);
2429  if(x!=tmpx0[m])
2430  {
2431  Eids_lay[m][0] = h;
2432  Vids_lay[m][0] = V2[h];
2433  Vids_lay[m][1] = V1[h];
2434  SpatialDomains::PointGeomSharedPtr vertex2
2435  = mesh->GetVertex(V2[h]);
2436  NekDouble x2=0.0,y2=0.0;
2437  normbef= sqrt( (y-y2)*(y-y2)+(x-x2)*(x-x2) );
2438  ybef = (y-y2);
2439  xbef = (x-x2);
2440  coord[0] = x;
2441  coord[1] = y;
2442  elmtid = streak->GetExpIndex(coord,0.00001);
2443  offset = streak->GetPhys_Offset(elmtid);
2444  Ubef = streak->GetExp(elmtid)->PhysEvaluate(coord, streak->GetPhys() + offset);
2445 
2446  }
2447  }
2448 
2449  }
2450 
2451 cout<<"Eid="<<Eids_lay[m][0]<<" Vids_lay0="<<Vids_lay[m][0]<<" Vidslay1="<<Vids_lay[m][1]<<endl;
2452 
2453  }
2454  else
2455  {
2456  //three verts(edges) candidates
2457  int cnt =0;
2458  for(int h=0; h< nTotEdges; h++)
2459  {
2460 //cout<<Vids_lay[m][g]<<" V1="<<V1[h]<<" V2[h]="<<V2[h]<<endl;
2461  if( (Vids_lay[m][g]==V1[h] || Vids_lay[m][g]==V2[h]) && h!= Eids_lay[m][g-1])
2462  {
2463 cout<<"edgetmp="<<h<<endl;
2464  ASSERTL0(cnt<=6, "wrong number of candidates");
2465  edgestmp[cnt]= h;
2466  cnt++;
2467  }
2468  }
2469 
2470  diff =1000;
2471  Array<OneD, NekDouble > diffarray(cnt);
2472  Array<OneD, NekDouble > diffUarray(cnt);
2473 cout<<"normbef="<<normbef<<endl;
2474 cout<<"Ubefcc="<<Ubef<<endl;
2475  //choose the right candidate
2476  for(int e=0; e< cnt; e++)
2477  {
2478  SpatialDomains::PointGeomSharedPtr vertex1 = mesh->GetVertex(V1[edgestmp[e]]);
2479  NekDouble x1,y1,z1;
2480  vertex1->GetCoords(x1,y1,z1);
2481  SpatialDomains::PointGeomSharedPtr vertex2 = mesh->GetVertex(V2[edgestmp[e]]);
2482  NekDouble x2,y2,z2;
2483  vertex2->GetCoords(x2,y2,z2);
2484 
2485  normtmp= sqrt( (y2-y1)*(y2-y1)+(x2-x1)*(x2-x1) );
2486 
2487 cout<<"edgetmp1="<<edgestmp[e]<<endl;
2488 cout<<"V1 x1="<<x1<<" y1="<<y1<<endl;
2489 cout<<"V2 x2="<<x2<<" y2="<<y2<<endl;
2490  if( Vids_lay[m][g]==V1[edgestmp[e]] )
2491  {
2492 
2493 
2494  ytmp = (y2-y1);
2495  xtmp = (x2-x1);
2496  coord[0] = x2;
2497  coord[1] = y2;
2498  elmtid = streak->GetExpIndex(coord,0.00001);
2499  offset = streak->GetPhys_Offset(elmtid);
2500  Utmp = streak->GetExp(elmtid)->PhysEvaluate(coord, streak->GetPhys() + offset);
2501  diffarray[e] = abs((xtmp*xbef+ytmp*ybef)/(normtmp*normbef)-1);
2502  diffUarray[e] = abs(Ubef-Utmp);
2503 cout<<" normtmp="<<normtmp<<endl;
2504 cout<<" Utmpcc="<<Utmp<<endl;
2505 cout<<xtmp<<" ytmp="<<ytmp<<" diff="<<abs(((xtmp*xbef+ytmp*ybef)/(normtmp*normbef))-1)<<endl;
2506  if(
2507  abs( (xtmp*xbef+ytmp*ybef)/(normtmp*normbef)-1)<diff
2508  && y2<= yoldup[g+1] && y2>= yolddown[g+1]
2509  && y1<= yoldup[g] && y1>= yolddown[g]
2510  )
2511  {
2512 
2513  Eids_lay[m][g] = edgestmp[e];
2514  Vids_lay[m][g+1] = V2[edgestmp[e]];
2515  diff = abs((xtmp*xbef+ytmp*ybef)/(normtmp*normbef)-1);
2516  normnext =normtmp;
2517  ynext = ytmp;
2518  xnext = xtmp;
2519  Unext = Utmp;
2520  }
2521  }
2522  else if( Vids_lay[m][g]==V2[edgestmp[e]] )
2523  {
2524 
2525 
2526  ytmp = (y1-y2);
2527  xtmp = (x1-x2);
2528  coord[0] = x1;
2529  coord[1] = y1;
2530  elmtid = streak->GetExpIndex(coord,0.00001);
2531  offset = streak->GetPhys_Offset(elmtid);
2532  Utmp = streak->GetExp(elmtid)->PhysEvaluate(coord, streak->GetPhys() + offset);
2533  diffarray[e] = abs((xtmp*xbef+ytmp*ybef)/(normtmp*normbef)-1);
2534  diffUarray[e] = abs(Ubef-Utmp);
2535 cout<<" normtmp="<<normtmp<<endl;
2536 cout<<" Utmpcc="<<Utmp<<endl;
2537 cout<<xtmp<<" ytmp="<<ytmp<<" diff="<<abs(((xtmp*xbef+ytmp*ybef)/(normtmp*normbef))-1)<<endl;
2538  if(
2539  abs((xtmp*xbef+ytmp*ybef)/(normtmp*normbef)-1)<diff
2540  && y2<= yoldup[g] && y2>= yolddown[g]
2541  && y1<= yoldup[g+1] && y1>= yolddown[g+1]
2542  )
2543  {
2544  Eids_lay[m][g] = edgestmp[e];
2545  Vids_lay[m][g+1] = V1[edgestmp[e]];
2546  diff = abs((xtmp*xbef+ytmp*ybef)/(normtmp*normbef)-1);
2547  normnext =normtmp;
2548  ynext = ytmp;
2549  xnext = xtmp;
2550  Unext = Utmp;
2551  }
2552 
2553  }
2554  else
2555  {
2556  ASSERTL0(false, "eid not found");
2557  }
2558 
2559 
2560 
2561  }
2562 cout<<"Eid before check="<<Eids_lay[m][g]<<endl;
2563 for(int q=0; q<cnt; q++)
2564 {
2565 cout<<q<<" diff"<<diffarray[q]<<endl;
2566 }
2567  //check if the eid has a vert in common with another layer
2568  bool check =false;
2569  if(m>0 && m< nlays)
2570  {
2571  check = checkcommonvert(Vids_lay[m-1],Vids_c,Vids_lay[m][g+1]);
2572  }
2573  if(check == true)
2574  {
2575 cout<<"COMMON VERT"<<endl;
2576  int eid = Vmath::Imin(cnt, diffarray,1);
2577  diffarray[eid]=1000;
2578  eid = Vmath::Imin(cnt,diffarray,1);
2579 
2580 
2581  SpatialDomains::PointGeomSharedPtr vertex1 = mesh->GetVertex(V1[edgestmp[eid]]);
2582  NekDouble x1,y1,z1;
2583  vertex1->GetCoords(x1,y1,z1);
2584  SpatialDomains::PointGeomSharedPtr vertex2 = mesh->GetVertex(V2[edgestmp[eid]]);
2585  NekDouble x2,y2,z2;
2586  vertex2->GetCoords(x2,y2,z2);
2587 
2588  normtmp= sqrt( (y2-y1)*(y2-y1)+(x2-x1)*(x2-x1) );
2589 
2590  Eids_lay[m][g] = edgestmp[eid];
2591  if(Vids_lay[m][g] == V1[edgestmp[eid]])
2592  {
2593  coord[0] = x2;
2594  coord[1] = y2;
2595  elmtid = streak->GetExpIndex(coord,0.00001);
2596  offset = streak->GetPhys_Offset(elmtid);
2597  Utmp = streak->GetExp(elmtid)->PhysEvaluate(coord, streak->GetPhys() + offset);
2598  Vids_lay[m][g+1] = V2[edgestmp[eid]];
2599  normnext =normtmp;
2600  ynext = (y2-y1);
2601  xnext = (x2-x1);;
2602  Unext = Utmp;
2603 
2604 
2605  }
2606  if(Vids_lay[m][g] == V2[edgestmp[eid]])
2607  {
2608  coord[0] = x1;
2609  coord[1] = y1;
2610  elmtid = streak->GetExpIndex(coord,0.00001);
2611  offset = streak->GetPhys_Offset(elmtid);
2612  Utmp = streak->GetExp(elmtid)->PhysEvaluate(coord, streak->GetPhys() + offset);
2613  Vids_lay[m][g+1] = V1[edgestmp[eid]];
2614  normnext =normtmp;
2615  ynext = (y1-y2);
2616  xnext = (x1-x2);;
2617  Unext = Utmp;
2618 
2619 
2620  }
2621 
2622  }
2623 
2624 cout<<m<<"edge aft:"<<Eids_lay[m][g]<<" Vid="<<Vids_lay[m][g+1]<<endl;
2625  normbef = normnext;
2626  ybef = ynext;
2627  xbef = xnext;
2628  Ubef = Unext;
2629 
2630 cout<<"endelse"<<normtmp<<endl;
2631  } //end else
2632 
2633  }//end g it
2634 
2635 
2636 
2637 
2638 
2639 
2640  } //end m it
2641 
2642 for(int w=0; w< nlays; w++)
2643 {
2644 for(int f=0; f< nlay_Eids; f++)
2645 {
2646 cout<<"check="<<w<<" Eid:"<<Eids_lay[w][f]<<endl;
2647 }
2648 }
2649 
2650 }
2651 
2652 
2653 bool checkcommonvert(Array<OneD, int> Vids_laybefore, Array<OneD, int> Vids_c, int Vid)
2654 {
2655  bool check=false;
2656  for(int u=0; u< Vids_laybefore.num_elements(); u++)
2657  {
2658  if( Vids_laybefore[u]==Vid || Vids_c[u]==Vid)
2659  {
2660  check =true;
2661  }
2662 cout<<Vid<<" Vert test="<<Vids_laybefore[u]<<endl;
2663  }
2664  return check;
2665 
2666 }
2667 
2668 
2669 void Cutrepetitions(int nedges,Array<OneD, NekDouble> inarray,
2670  Array<OneD, NekDouble>& outarray)
2671 {
2672 
2673  //determine npedge:
2674  int np_lay = inarray.num_elements();
2675  ASSERTL0(inarray.num_elements()%nedges==0," something on number npedge");
2676  //cut out the repetitions:
2677 
2678  int cnt=0;
2679  for(int w=0; w< np_lay; w++)
2680  {
2681 
2682  if(w< np_lay-1)
2683  {
2684  if(inarray[w] ==inarray[w+1])
2685  {
2686  continue;
2687  }
2688  }
2689  outarray[cnt]= inarray[w];
2690  cnt++;
2691  }
2692 
2693 
2694  ASSERTL0( cnt== np_lay-(nedges-1), "wrong cut");
2695 
2696 }
2697 
2698 int DetermineclosePointxindex(NekDouble x,Array<OneD, NekDouble> xArray)
2699 {
2700  int npts = xArray.num_elements();
2701  Array<OneD, NekDouble> xcopy(npts);
2702  Vmath::Vcopy(npts,xArray,1,xcopy,1);
2703  //subtract xpoint and abs
2704 
2705  Vmath::Sadd(npts, -x, xcopy,1, xcopy,1);
2706  Vmath::Vabs(npts, xcopy,1,xcopy,1);
2707 
2708  int index = Vmath::Imin(npts, xcopy,1);
2709  if(xArray[index]> x)//assume always x[index]< x(HYPHOTHESIS)
2710  {
2711  index = index-1;
2712  }
2713  return index;
2714 }
2715 
2716 void GenerateNeighbourArrays(int index,int neighpoints, Array<OneD, NekDouble> xArray,
2717  Array<OneD, NekDouble> yArray,Array<OneD, NekDouble>& Neighbour_x,
2718  Array<OneD, NekDouble>& Neighbour_y)
2719 {
2720  ASSERTL0( neighpoints%2==0,"number of neighbour points should be even");
2721  int leftpoints = (neighpoints/2)-1;//-1 because xArray[index]< x
2722  int rightpoints = neighpoints/2;
2723  int diff ;
2724  int start;
2725 //cout<<"index="<<index<<" left="<<leftpoints<<" right="<<rightpoints<<endl;
2726  if(index-leftpoints<0)
2727  {
2728 //cout<"case0"<<endl;
2729  diff = index-leftpoints;
2730  start= 0;
2731  Vmath::Vcopy(neighpoints, &yArray[0],1,&Neighbour_y[0],1);
2732  Vmath::Vcopy(neighpoints, &xArray[0],1,&Neighbour_x[0],1);
2733  }
2734  else if( (yArray.num_elements()-1)-index < rightpoints)
2735  {
2736 //cout"case1 closest="<<xArray[index]<<endl;
2737  int rpoints = (yArray.num_elements()-1)-index;//
2738  diff = rightpoints-rpoints;
2739 //cout<"start index="<<index-leftpoints-diff<<endl;
2740  start = index-leftpoints-diff;
2741  Vmath::Vcopy(neighpoints, &yArray[start],1,&Neighbour_y[0],1);
2742  Vmath::Vcopy(neighpoints, &xArray[start],1,&Neighbour_x[0],1);
2743  }
2744  else
2745  {
2746 //cout<<"caseaa"<<endl;
2747  start = index-leftpoints;
2748  Vmath::Vcopy(neighpoints, &yArray[start],1,&Neighbour_y[0],1);
2749  Vmath::Vcopy(neighpoints, &xArray[start],1,&Neighbour_x[0],1);
2750  }
2751 /*
2752 for(int t= start; t<start+neighpoints; t++)
2753 {
2754 cout<<"Px="<<xArray[t]<<" "<<yArray[t]<<endl;
2755 }
2756 */
2757  //check if there is any repetition
2758  for(int f=1; f< neighpoints; f++)
2759  {
2760  ASSERTL0(Neighbour_x[f]!=Neighbour_x[f-1]," repetition on NeighbourArrays");
2761  }
2762 
2763 }
2764 
2765 NekDouble LagrangeInterpolant(NekDouble x, int npts,
2766  Array<OneD,NekDouble> xpts, Array<OneD, NekDouble> funcvals)
2767 {
2768  NekDouble sum = 0.0;
2769  NekDouble LagrangePoly;
2770 //cout<<"lagrange"<<endl;
2771  for(int pt=0;pt<npts;++pt)
2772  {
2773  NekDouble h=1.0;
2774 
2775  for(int j=0;j<pt; ++j)
2776  {
2777  h = h * (x - xpts[j])/(xpts[pt]-xpts[j]);
2778  }
2779 
2780  for(int k=pt+1;k<npts;++k)
2781  {
2782  h = h * (x - xpts[k])/(xpts[pt]-xpts[k]);
2783  }
2784  LagrangePoly=h;
2785 
2786  sum += funcvals[pt]*LagrangePoly;
2787  }
2788 //cout<<"result :"<<sum<<endl;
2789  return sum;
2790 }
2791 
2792 
2793 void EvaluateTangent(int npoints, Array<OneD, NekDouble> xcQedge,
2794  Array<OneD, NekDouble> coeffsinterp,
2795  Array<OneD, NekDouble> & txQedge, Array<OneD, NekDouble> & tyQedge)
2796 {
2797  Array<OneD, NekDouble> yprime(npoints,0.0);
2798  int np_pol= coeffsinterp.num_elements();
2799 cout<<"evaluatetan with "<<np_pol<<endl;
2800 
2801  //calc derivative poly (cut last entry on b array that is the const)
2802  int derorder;
2803  Array<OneD, NekDouble> yinterp(npoints,0.0);
2804  int polorder;
2805  for(int q=0; q< npoints; q++)
2806  {
2807  polorder=np_pol-1;
2808  derorder=np_pol-2;
2809  yprime[q] =0;
2810  for(int d=0; d< np_pol-1; d++)
2811  {
2812  yprime[q] += (derorder +1)*coeffsinterp[d]*std::pow(xcQedge[q],derorder);
2813  derorder--;
2814  }
2815  //test
2816  for(int a=0; a< np_pol; a++)
2817  {
2818  yinterp[q] += coeffsinterp[a]*std::pow(xcQedge[q],polorder);
2819 //cout<<"coeff*x^n="<<b[a]*std::pow(xcQedge[q],polorder)<<" sum="<<yinterp[q]<<endl;
2820  polorder--;
2821  }
2822 
2823  }
2824 
2825  //transf yprime into tx,ty:
2826  for(int n=0; n< npoints; n++)
2827  {
2828  //ABS???!!
2829  txQedge[n]=0;
2830  txQedge[n] = cos((atan((yprime[n]))));
2831  tyQedge[n] = sin((atan((yprime[n]))));
2832 cout<<xcQedge[n]<<" "<<yinterp[n]<<" "<<yprime[n]<<" "<<txQedge[n]<<" "<<tyQedge[n]<<endl;
2833  }
2834 }
2835 
2836 void PolyInterp( Array<OneD, NekDouble> xpol, Array<OneD, NekDouble> ypol,
2837  Array<OneD, NekDouble> & coeffsinterp,
2838  Array<OneD, NekDouble> & xcout, Array<OneD, NekDouble> & ycout,
2839  int edge, int npedge)
2840 {
2841  int np_pol = xpol.num_elements();
2842  int N = np_pol;
2843  Array<OneD, NekDouble> A (N*N,1.0);
2844  Array<OneD, NekDouble> b (N);
2845  int row=0;
2846  //fill column by column
2847  for(int e=0; e<N; e++)
2848  {
2849  row=0;
2850  for(int w=0; w < N; w++)
2851  {
2852  A[N*e+row] = std::pow( xpol[w], N-1-e);
2853  row++;
2854  }
2855  }
2856  row=0;
2857  for(int r= 0; r< np_pol; r++)
2858  {
2859  b[row] = ypol[r];
2860 //cout<<"b="<<b[row]<<" y_c="<<y_c[r]<<endl;
2861  row++;
2862  }
2863 
2864 //cout<<"A elements="<<A.num_elements()<<endl;
2865  Array<OneD, int> ipivot (N);
2866  int info =0;
2867  //Lapack::Dgesv( N, 1, A.get(), N, ipivot.get(), b.get(), N, info);
2868  Lapack::Dgetrf( N, N, A.get(), N, ipivot.get(), info);
2869  if( info < 0 )
2870  {
2871  std::string message = "ERROR: The " + boost::lexical_cast<std::string>(-info) +
2872  "th parameter had an illegal parameter for dgetrf";
2873  ASSERTL0(false, message.c_str());
2874  }
2875  else if( info > 0 )
2876  {
2877  std::string message = "ERROR: Element u_" + boost::lexical_cast<std::string>(info) +
2878  boost::lexical_cast<std::string>(info) + " is 0 from dgetrf";
2879  ASSERTL0(false, message.c_str());
2880  }
2881 
2882  // N means no transponse (direct matrix)
2883  int ncolumns_b =1;
2884  Lapack::Dgetrs( 'N', N, ncolumns_b , A.get() , N, ipivot.get(), b.get(), N, info);
2885  if( info < 0 )
2886  {
2887  std::string message = "ERROR: The " + boost::lexical_cast<std::string>(-info) +
2888  "th parameter had an illegal parameter for dgetrf";
2889  ASSERTL0(false, message.c_str());
2890  }
2891  else if( info > 0 )
2892  {
2893  std::string message = "ERROR: Element u_" + boost::lexical_cast<std::string>(info) +
2894  boost::lexical_cast<std::string>(info) + " is 0 from dgetrf";
2895  ASSERTL0(false, message.c_str());
2896  }
2897 /*
2898 for(int h=0; h<np_pol; h++)
2899 {
2900 cout<<"coeff:"<<b[h]<<endl;
2901 }
2902 */
2903 
2904 
2905 
2906  //ovewrite the ycPhysMOD:
2907  int polorder;
2908  for(int c=0; c< npedge; c++)
2909  {
2910  polorder=np_pol-1;
2911  //put ycPhysMOD to zero
2912  ycout[edge*(npedge)+c+1]=0;
2913  for(int d=0; d< np_pol; d++)
2914  {
2915  ycout[edge*(npedge)+c+1] += b[d]
2916  *std::pow(xcout[edge*(npedge)+c+1],polorder);
2917  polorder--;
2918  }
2919  }
2920 
2921  //write coeffs
2922  Vmath::Vcopy(np_pol, b,1,coeffsinterp,1);
2923 
2924 }
2925 
2926 
2927 void PolyFit(int polyorder,int npoints,
2928  Array<OneD, NekDouble> xin, Array<OneD, NekDouble> fin,
2929  Array<OneD, NekDouble> & coeffsinterp,
2930  Array<OneD, NekDouble> & xout, Array<OneD, NekDouble> & fout,
2931  int npout)
2932 {
2933 
2934  int N = polyorder+1;
2935  Array<OneD, NekDouble> A (N*N,0.0);
2936  Array<OneD, NekDouble> b (N,0.0);
2937 cout<<npoints<<endl;
2938 for(int u=0; u<npoints; u++)
2939 {
2940 cout<<"c="<<xin[u]<<" "<<
2941 fin[u]<<endl;
2942 }
2943  //fill column by column
2944  //e counts cols
2945  for(int e=0; e<N; e++)
2946  {
2947 
2948  for(int row=0; row<N; row++)
2949  {
2950  for(int w=0; w < npoints; w++)
2951  {
2952  A[N*e+row] += std::pow( xin[w], e+row);
2953  }
2954  }
2955  }
2956 
2957  for(int row= 0; row< N; row++)
2958  {
2959  for(int h=0; h< npoints; h++)
2960  {
2961  b[row] += fin[h]*std::pow(xin[h],row);
2962 //cout<<b="<<b[row]<<" y_c="<<y_c[r]<<endl;
2963  }
2964 
2965  }
2966 
2967 
2968 
2969 
2970 
2971 //cout<<"A elements="<<A.num_elements()<<endl;
2972  Array<OneD, int> ipivot (N);
2973  int info =0;
2974  //Lapack::Dgesv( N, 1, A.get(), N, ipivot.get(), b.get(), N, info);
2975  Lapack::Dgetrf( N, N, A.get(), N, ipivot.get(), info);
2976 
2977  if( info < 0 )
2978  {
2979  std::string message = "ERROR: The " + boost::lexical_cast<std::string>(-info) +
2980  "th parameter had an illegal parameter for dgetrf";
2981  ASSERTL0(false, message.c_str());
2982  }
2983  else if( info > 0 )
2984  {
2985  std::string message = "ERROR: Element u_" + boost::lexical_cast<std::string>(info) +
2986  boost::lexical_cast<std::string>(info) + " is 0 from dgetrf";
2987  ASSERTL0(false, message.c_str());
2988  }
2989  // N means no transponse (direct matrix)
2990  int ncolumns_b =1;
2991  Lapack::Dgetrs( 'N', N, ncolumns_b , A.get() , N, ipivot.get(), b.get(), N, info);
2992  if( info < 0 )
2993  {
2994  std::string message = "ERROR: The " + boost::lexical_cast<std::string>(-info) +
2995  "th parameter had an illegal parameter for dgetrf";
2996  ASSERTL0(false, message.c_str());
2997  }
2998  else if( info > 0 )
2999  {
3000  std::string message = "ERROR: Element u_" + boost::lexical_cast<std::string>(info) +
3001  boost::lexical_cast<std::string>(info) + " is 0 from dgetrf";
3002  ASSERTL0(false, message.c_str());
3003  }
3004 
3005  //the lower coeff the lower is the grad
3006  //reverse:
3007  Array<OneD, NekDouble> tmp(N);
3008  Vmath::Vcopy(N,b,1,tmp,1);
3009  int cnt =N;
3010  for(int j=0; j<N; j++)
3011  {
3012  b[j]= tmp[cnt-1];
3013  cnt--;
3014  }
3015 
3016 for(int h=0; h<N; h++)
3017 {
3018 cout<<"coeff:"<<b[h]<<endl;
3019 }
3020 
3021  //ovewrite the ycPhysMOD:
3022  int polorder;
3023  for(int c=0; c< npout; c++)
3024  {
3025  polorder=polyorder;
3026  //put ycPhysMOD to zero
3027  fout[c]=0;
3028  for(int d=0; d< N; d++)
3029  {
3030 
3031  fout[c] += b[d]
3032  *std::pow(xout[c],polorder);
3033  polorder--;
3034 
3035  }
3036 
3037  }
3038 
3039  //write coeffs
3040  Vmath::Vcopy(N, b,1,coeffsinterp,1);
3041 
3042 }
3043 
3044 
3045 void Orderfunctionx(Array<OneD, NekDouble> inarray_x,
3046  Array<OneD, NekDouble> inarray_y, Array<OneD, NekDouble>& outarray_x,
3047  Array<OneD, NekDouble>& outarray_y)
3048 {
3049  Array<OneD, NekDouble>tmpx(inarray_x.num_elements());
3050  Array<OneD, NekDouble>tmpy(inarray_x.num_elements());
3051  //local copy to prevent overwriting
3052  Vmath::Vcopy(inarray_x.num_elements() , inarray_x,1,tmpx,1);
3053  Vmath::Vcopy(inarray_x.num_elements() , inarray_y,1,tmpy,1);
3054 
3055  //order function with respect to x
3056  int index;
3057  NekDouble max = Vmath::Vmax(tmpx.num_elements(), tmpx,1);
3058  for(int w=0; w<tmpx.num_elements(); w++)
3059  {
3060  index = Vmath::Imin(tmpx.num_elements(), tmpx,1);
3061  outarray_x[w]= tmpx[index];
3062  outarray_y[w]= tmpy[index];
3063  if(w< tmpx.num_elements()-1)//case of repetitions
3064  {
3065  if(tmpx[index] == tmpx[index+1])
3066  {
3067  outarray_x[w+1]= tmpx[index+1];
3068  outarray_y[w+1]= tmpy[index+1];
3069  tmpx[index+1] = max+1000;
3070  w++;
3071  }
3072  }
3073 /*
3074  if(w>0)//case of repetitions
3075  {
3076  if(inarray_x[index] == tmpx[index-1])
3077  {
3078  outarray_x[w+1]= tmpx[index-1];
3079  outarray_y[w+1]= tmpy[index-1];
3080  tmpx[index-1] = max+1000;
3081  w++;
3082  }
3083  }
3084 */
3085  tmpx[index] = max+1000;
3086 
3087  }
3088 }
3089 
3090 void MoveLayersvertically(int nlays, int nvertl, int cntlow, int cntup,
3091  Array<OneD, Array<OneD, int > > lay_Vids, Array<OneD, NekDouble> xc,
3092  Array<OneD, NekDouble> yc, Array<OneD, int> Down, Array<OneD, int> Up,
3093  Array<OneD, NekDouble >& xnew,Array<OneD, NekDouble>& ynew,
3094  Array<OneD, Array<OneD, NekDouble > >& layers_x,
3095  Array<OneD, Array<OneD, NekDouble > >& layers_y)
3096 {
3097  int np_lay = layers_y[0].num_elements();
3098  // 0<h<nlays-1 to fill only the 'internal' layers(no up,low);
3099  for(int h=1; h<nlays-1; h++)
3100  {
3101  layers_x[h]= Array<OneD, NekDouble>(np_lay);
3102  for(int s=0; s<nvertl; s++)
3103  {
3104  //check if ynew is still empty
3105  ASSERTL0(ynew[ lay_Vids[h][s] ]==-20, "ynew layers not empty");
3106  if(h<cntlow+1)//layers under the crit lay
3107  {
3108  //y= ylow+delta
3109  ynew[ lay_Vids[h][s] ] = ynew[Down[s]]+ h*abs(ynew[Down[s]] - yc[s])/(cntlow+1);
3110  //put the layer vertical
3111  xnew[lay_Vids[h][s] ] = xc[s];
3112 //cout<<"ynew="<<ynew[ lay_Vids[h][s] ]<<" ydown="<<ynew[Down[s]]<<
3113 //" delta="<<abs(ynew[Down[s]] - y_c[s])/(cntlow+1)<<endl;
3114  //until now layers_y=yold
3115  layers_y[h][s] = ynew[ lay_Vids[h][s] ];
3116  layers_x[h][s] = xnew[ lay_Vids[h][s] ];
3117  }
3118  else
3119  {
3120  //y = yc+delta
3121  ynew[ lay_Vids[h][s] ] = yc[s] + (h-cntlow)*abs(ynew[Up[s]] - yc[s])/(cntup+1);
3122  //put the layer vertical
3123  xnew[lay_Vids[h][s] ] = xc[s];
3124  //until now layers_y=yold
3125  layers_y[h][s] = ynew[ lay_Vids[h][s] ];
3126  layers_x[h][s] = xnew[ lay_Vids[h][s] ];
3127  }
3128  }
3129  }
3130 
3131 }
3132 
3133 
3134 void MoveLayerNfixedxpos(int nvertl, int npedge, Array<OneD, NekDouble> xcPhys,
3135  Array<OneD, NekDouble> tmpx_lay, Array<OneD, NekDouble> tmpy_lay,
3136  Array<OneD, int> Vids,
3137  Array<OneD, NekDouble> &xlay, Array<OneD, NekDouble> &ylay,
3138  Array<OneD, NekDouble> &xnew, Array<OneD, NekDouble> &ynew)
3139 {
3140  int np_lay = xcPhys.num_elements();
3141  int nedges = nvertl-1;
3142  Array<OneD, NekDouble>tmpx(np_lay-(nedges-1));
3143  Array<OneD, NekDouble>tmpy(np_lay-(nedges-1));
3144  Cutrepetitions(nedges, tmpx_lay, tmpx);
3145  Cutrepetitions(nedges, tmpy_lay, tmpy);
3146  //order points in x:
3147  int index;
3148  int closepoints = 4;
3149  Array<OneD, NekDouble>Pxinterp(closepoints);
3150  Array<OneD, NekDouble>Pyinterp(closepoints);
3151  Orderfunctionx(tmpx, tmpy, tmpx, tmpy);
3152  //determine the neighbour points (-3;+3)
3153  for(int g=0; g< nedges; g++)
3154  {
3155  //write vert coords
3156  //v1
3157  index=
3158  DetermineclosePointxindex( xcPhys[g*npedge+0], tmpx);
3159  //generate neighbour arrays:
3160  GenerateNeighbourArrays(index, closepoints,tmpx,tmpy,Pxinterp,Pyinterp);
3161  ynew[Vids[g] ]= LagrangeInterpolant(xcPhys[g*npedge+0],closepoints,Pxinterp,Pyinterp );
3162  xnew[Vids[g] ]= xcPhys[g*npedge+0];
3163  ylay[g*npedge +0] = ynew[ Vids[g] ];
3164  xlay[g*npedge +0] = xnew[ Vids[g] ];
3165 
3166  //v2
3167  //determine closest index:
3168  index=
3169  DetermineclosePointxindex( xcPhys[g*npedge +npedge-1], tmpx);
3170  //generate neighbour arrays:
3171  GenerateNeighbourArrays(index, closepoints,tmpx,tmpy,Pxinterp,Pyinterp);
3172  ynew[Vids[g+1] ]= LagrangeInterpolant(xcPhys[g*npedge +npedge-1],closepoints,Pxinterp,Pyinterp );
3173  xnew[Vids[g+1] ]= xcPhys[g*npedge +npedge-1];
3174  ylay[g*npedge +npedge-1] = ynew[Vids[g+1] ];
3175  xlay[g*npedge +npedge-1] = xnew[Vids[g+1] ];
3176 
3177 
3178 
3179  //middle points
3180  for(int r=0; r< npedge-2; r++)
3181  {
3182 
3183  //determine closest point index:
3184  index =
3185  DetermineclosePointxindex( xcPhys[g*npedge +r+1], tmpx);
3186 //cout<<" vert+"<<index<<endl;
3187 
3188  ASSERTL0( index<= tmpy.num_elements()-1, " index wrong");
3189  //generate neighbour arrays Pyinterp,Pxinterp
3190  GenerateNeighbourArrays(index, closepoints,tmpx,tmpy,Pxinterp,Pyinterp);
3191 
3192  ylay[g*npedge +r+1]=
3193  LagrangeInterpolant(
3194  xcPhys[g*npedge +r+1],closepoints,Pxinterp,Pyinterp );
3195 //cout<<"x value="<<xcPhysMOD[g*npedge +r+1]<<endl;
3196  xlay[g*npedge +r+1]= xcPhys[g*npedge +r+1];
3197 
3198 
3199 
3200 
3201 /*
3202 for(int t=0; t<6; t++)
3203 {
3204 cout<<"Px="<<Pxinterp[t]<<" "<<Pyinterp[t]<<endl;
3205 }
3206 */
3207  }
3208 
3209  }//edge closed g
3210 
3211 }
3212 
3213 void MoveLayerNnormpos(int nvertl, int npedge, Array<OneD, NekDouble> xcPhys,
3214  Array<OneD, NekDouble> tmpx_lay, Array<OneD, NekDouble> tmpy_lay,
3215  Array<OneD, int> Vids,
3216  Array<OneD, NekDouble> &xlay, Array<OneD, NekDouble> &ylay,
3217  Array<OneD, NekDouble> &xnew, Array<OneD, NekDouble> &ynew)
3218 {
3219  int np_lay = xcPhys.num_elements();
3220  int nedges = nvertl-1;
3221  NekDouble x0,x1, xtmp;
3222  Array<OneD, NekDouble>tmpx(np_lay-(nedges-1));
3223  Array<OneD, NekDouble>tmpy(np_lay-(nedges-1));
3224  Cutrepetitions(nedges, tmpx_lay, tmpx);
3225  Cutrepetitions(nedges, tmpy_lay, tmpy);
3226  //order points in x:
3227  int index;
3228  int closepoints = 4;
3229  Array<OneD, NekDouble>Pxinterp(closepoints);
3230  Array<OneD, NekDouble>Pyinterp(closepoints);
3231  Orderfunctionx(tmpx, tmpy, tmpx, tmpy);
3232  //determine the neighbour points (-3;+3)
3233  for(int g=0; g< nedges; g++)
3234  {
3235  //write vert coords
3236  //v1
3237  ynew[Vids[g] ]= tmpy_lay[g*npedge+0];
3238  xnew[Vids[g] ]= tmpx_lay[g*npedge+0];
3239 
3240 
3241  ylay[g*npedge +0] = ynew[ Vids[g] ];
3242  xlay[g*npedge +0] = xnew[ Vids[g] ];
3243 
3244  //v2
3245  ynew[Vids[g+1] ]= tmpy_lay[g*npedge+npedge-1];
3246  xnew[Vids[g+1] ]= tmpx_lay[g*npedge+npedge-1];
3247  ylay[g*npedge +npedge-1] = ynew[Vids[g+1] ];
3248  xlay[g*npedge +npedge-1] = xnew[Vids[g+1] ];
3249 
3250 
3251 
3252  //middle points
3253  for(int r=0; r< npedge-2; r++)
3254  {
3255  x0 = xlay[g*npedge +0];
3256  x1 = xlay[g*npedge +npedge-1];
3257  xtmp = x0 + r*(x1-x0)/(npedge-1);
3258  //determine closest point index:
3259  index =
3260  DetermineclosePointxindex( xtmp, tmpx);
3261 //cout<<" vert+"<<index<<endl;
3262 
3263  ASSERTL0( index<= tmpy.num_elements()-1, " index wrong");
3264  //generate neighbour arrays Pyinterp,Pxinterp
3265  GenerateNeighbourArrays(index, closepoints,tmpx,tmpy,Pxinterp,Pyinterp);
3266 
3267  ylay[g*npedge +r+1]=
3268  LagrangeInterpolant(
3269  xtmp,closepoints,Pxinterp,Pyinterp );
3270 //cout<<"x value="<<xcPhysMOD[g*npedge +r+1]<<endl;
3271  xlay[g*npedge +r+1]= xtmp;
3272 
3273 
3274  }
3275 
3276  }//edge closed g
3277 }
3278 
3279 void MoveOutsidePointsfixedxpos(int npedge, SpatialDomains::MeshGraphSharedPtr mesh,
3280  Array<OneD, NekDouble> xcold,Array<OneD, NekDouble> ycold,
3281  Array<OneD, NekDouble> xolddown,Array<OneD, NekDouble> yolddown,
3282  Array<OneD, NekDouble> xoldup,Array<OneD, NekDouble> yoldup,
3283  Array<OneD, NekDouble> ylaydown,Array<OneD, NekDouble> ylayup,
3284  Array<OneD, NekDouble>& xnew,Array<OneD, NekDouble>& ynew)
3285 {
3286  //update vertices coords outside layers region
3287  int nvertl = ycold.num_elements();
3288  int nVertTot = mesh->GetNvertices();
3289  for(int n=0; n<nVertTot; n++)
3290  {
3291  NekDouble ratio;
3292  SpatialDomains::PointGeomSharedPtr vertex = mesh->GetVertex(n);
3293  NekDouble x,y,z;
3294  vertex->GetCoords(x,y,z);
3295  int qp_closer;
3296  //determine the closer xold_up
3297  NekDouble tmp=1000;
3298 
3299  for(int k=0; k<nvertl; k++)
3300  {
3301  if(abs(x-xcold[k]) < tmp)
3302  {
3303  tmp = abs(x-xcold[k]);
3304  qp_closer=k;
3305  }
3306 
3307  }
3308  //find nplay_closer
3309  int nplay_closer;
3310  if(qp_closer==0)
3311  {
3312  nplay_closer=0;//first vert
3313  }
3314  else
3315  {
3316  nplay_closer= (qp_closer-1)*npedge +npedge-1;
3317  }
3318 
3319 
3320  if( y>yoldup[qp_closer] && y<1 )//nlays-1 is layerup
3321  {
3322 
3323 // ratio = (1-layers_y[nlays-1][qp_closer])*(1-y_c[qp_closer])/
3324 // ( (1-yold_up[n])*(1-yold_c[qp_closer]) );
3325  ratio = (1-ylayup[nplay_closer])/
3326  ( (1-yoldup[qp_closer]) );
3327  //distance prop to layerup
3328  ynew[n] = ylayup[nplay_closer]
3329  + (y-yoldup[qp_closer])*ratio;
3330  xnew[n] = x;
3331 
3332  }
3333  else if( y< yolddown[qp_closer] && y>-1 )//0 is layerdown
3334  {
3335 
3336  ratio = (1+ylaydown[nplay_closer])/
3337  ( (1+yolddown[qp_closer]) );
3338  //distance prop to layerlow
3339  ynew[n] = ylaydown[nplay_closer]
3340  + (y-yolddown[qp_closer])*ratio;
3341  xnew[n] = x;
3342  }
3343 
3344  }
3345 }
3346 
3347 void MoveOutsidePointsNnormpos(int npedge, SpatialDomains::MeshGraphSharedPtr mesh,
3348  Array<OneD, NekDouble> xcold,Array<OneD, NekDouble> ycold,
3349  Array<OneD, NekDouble> xolddown,Array<OneD, NekDouble> yolddown,
3350  Array<OneD, NekDouble> xoldup,Array<OneD, NekDouble> yoldup,
3351  Array<OneD, NekDouble> xlaydown,Array<OneD, NekDouble> ylaydown,
3352  Array<OneD, NekDouble> xlayup,Array<OneD, NekDouble> ylayup,
3353  Array<OneD, NekDouble> nxPhys,Array<OneD, NekDouble> nyPhys,
3354  Array<OneD, NekDouble>& xnew,Array<OneD, NekDouble>& ynew)
3355 {
3356 /*
3357  int nq1D =bndfieldup->GetTotPoints();
3358  Array<OneD, NekDouble> xlayoldup(nq1D);
3359  Array<OneD, NekDouble> xlayolddown(nq1D);
3360  Array<OneD, NekDouble> ylayoldup(nq1D);
3361  Array<OneD, NekDouble> ylayolddown(nq1D);
3362  Array<OneD, NekDouble> zlayoldup(nq1D);
3363  Array<OneD, NekDouble> zlayolddown(nq1D);
3364  bndfielddown->GetCoords( xlayolddown, ylayolddown,zlayolddown);
3365  bndfieldup->GetCoords( xlayoldup, ylayoldup,zlayoldup);
3366 
3367  NekDouble xmax = Vmath::Vmax(nq1D, xlayoldup,1);
3368  NekDouble xmin = Vmath::Vmin(nq1D, xlayoldup,1);
3369 */
3370  //determine the new verts up/down pos:
3371  int nvertl = xoldup.num_elements();
3372  int nedges = nvertl-1;
3373  Array<OneD, NekDouble> xnew_down(nvertl);
3374  Array<OneD, NekDouble> ynew_down(nvertl);
3375  Array<OneD, NekDouble> xnew_up(nvertl);
3376  Array<OneD, NekDouble> ynew_up(nvertl);
3377  Array<OneD, NekDouble> nxvert(nvertl);
3378  Array<OneD, NekDouble> nyvert(nvertl);
3379  Array<OneD, NekDouble> norm(nvertl);
3380  Array<OneD, NekDouble> tmp(nvertl);
3381  for(int a=0; a< nedges;a++)
3382  {
3383  if(a==0)
3384  {
3385  //v1
3386  xnew_down[a] = xlaydown[a*npedge+0];
3387  ynew_down[a] = ylaydown[a*npedge+0];
3388  xnew_up[a] = xlayup[a*npedge+0];
3389  ynew_up[a] = ylayup[a*npedge+0];
3390  nxvert[a] = nxPhys[a*npedge+0];
3391  nyvert[a] = nyPhys[a*npedge+0];
3392  //v2
3393  xnew_down[a+1] = xlaydown[a*npedge+npedge-1];
3394  ynew_down[a+1] = ylaydown[a*npedge+npedge-1];
3395  xnew_up[a+1] = xlayup[a*npedge+npedge-1];
3396  ynew_up[a+1] = ylayup[a*npedge+npedge-1];
3397  nxvert[a+1] = nxPhys[a*npedge+npedge-1];
3398  nyvert[a+1] = nyPhys[a*npedge+npedge-1];
3399  }
3400  else
3401  {
3402  //v2
3403  xnew_down[a+1] = xlaydown[a*npedge+npedge-1];
3404  ynew_down[a+1] = ylaydown[a*npedge+npedge-1];
3405  xnew_up[a+1] = xlayup[a*npedge+npedge-1];
3406  ynew_up[a+1] = ylayup[a*npedge+npedge-1];
3407  nxvert[a+1] = nxPhys[a*npedge+npedge-1];
3408  nyvert[a+1] = nyPhys[a*npedge+npedge-1];
3409  }
3410 
3411  }
3412 
3413  NekDouble xmax = Vmath::Vmax(nvertl, xoldup,1);
3414  NekDouble xmin = Vmath::Vmin(nvertl, xoldup,1);
3415 
3416  //update vertices coords outside layers region
3417  NekDouble ratiox;
3418  //NekDouble ratioy;
3419 
3420  int nVertTot = mesh->GetNvertices();
3421  for(int n=0; n<nVertTot; n++)
3422  {
3423  NekDouble ratio;
3424  SpatialDomains::PointGeomSharedPtr vertex = mesh->GetVertex(n);
3425  NekDouble x,y,z;
3426  vertex->GetCoords(x,y,z);
3427  int qp_closeroldup = 0, qp_closerolddown = 0;
3428  NekDouble diffup, diffdown;
3429  //determine the closer xold_up,down
3430  diffdown =1000;
3431  diffup = 1000;
3432 
3433 
3434 
3435 
3436  for(int k=0; k<nvertl; k++)
3437  {
3438  if(abs(x-xolddown[k]) < diffdown)
3439  {
3440  diffdown = abs(x-xolddown[k]);
3441  qp_closerolddown=k;
3442  }
3443  if(abs(x-xoldup[k]) < diffup)
3444  {
3445  diffup = abs(x-xoldup[k]);
3446  qp_closeroldup=k;
3447  }
3448 
3449  }
3450 
3451  //find nplay_closer
3452  diffdown =1000;
3453  diffup = 1000;
3454 
3455  int qp_closerup = 0, qp_closerdown = 0;
3456 
3457  for(int f=0; f< nvertl; f++)
3458  {
3459  if(abs(x-xnew_down[f]) < diffdown)
3460  {
3461  diffdown = abs(x-xnew_down[f]);
3462  qp_closerdown=f;
3463  }
3464  if(abs(x-xnew_up[f]) < diffup)
3465  {
3466  diffup = abs(x-xnew_up[f]);
3467  qp_closerup=f;
3468  }
3469  }
3470 
3471 // int qp_closernormoldup;
3472  Vmath::Sadd(nvertl, -x,xoldup,1, tmp,1);
3473  Vmath::Vmul(nvertl, tmp,1,tmp,1,tmp,1);
3474  Vmath::Sadd(nvertl,-y,yoldup,1,norm,1);
3475  Vmath::Vmul(nvertl, norm,1,norm,1,norm,1);
3476  Vmath::Vadd(nvertl, norm,1,tmp,1,norm,1);
3477 // qp_closernormoldup = Vmath::Imin(nvertl, norm,1);
3478 
3479  Vmath::Zero(nvertl, norm,1);
3480  Vmath::Zero(nvertl, tmp,1);
3481 
3482 // int qp_closernormolddown;
3483  Vmath::Sadd(nvertl, -x,xolddown,1, tmp,1);
3484  Vmath::Vmul(nvertl, tmp,1,tmp,1,tmp,1);
3485  Vmath::Sadd(nvertl,-y,yolddown,1,norm,1);
3486  Vmath::Vmul(nvertl, norm,1,norm,1,norm,1);
3487  Vmath::Vadd(nvertl, norm,1,tmp,1,norm,1);
3488 // qp_closernormolddown = Vmath::Imin(nvertl, norm,1);
3489 
3490  Vmath::Zero(nvertl, norm,1);
3491  Vmath::Zero(nvertl, tmp,1);
3492 
3493  int qp_closernormup;
3494  Vmath::Sadd(nvertl, -x,xnew_up,1, tmp,1);
3495  Vmath::Vmul(nvertl, tmp,1,tmp,1,tmp,1);
3496  Vmath::Sadd(nvertl,-y,ynew_up,1,norm,1);
3497  Vmath::Vmul(nvertl, norm,1,norm,1,norm,1);
3498  Vmath::Vadd(nvertl, norm,1,tmp,1,norm,1);
3499  qp_closernormup = Vmath::Imin(nvertl, norm,1);
3500 
3501  Vmath::Zero(nvertl, norm,1);
3502  Vmath::Zero(nvertl, tmp,1);
3503 
3504  int qp_closernormdown;
3505  Vmath::Sadd(nvertl, -x,xnew_down,1, tmp,1);
3506  Vmath::Vmul(nvertl, tmp,1,tmp,1,tmp,1);
3507  Vmath::Sadd(nvertl,-y,ynew_down,1,norm,1);
3508  Vmath::Vmul(nvertl, norm,1,norm,1,norm,1);
3509  Vmath::Vadd(nvertl, norm,1,tmp,1,norm,1);
3510  qp_closernormdown = Vmath::Imin(nvertl, norm,1);
3511 
3512 
3513 
3514 
3515  if( y>yoldup[qp_closeroldup] && y<1 )
3516  {
3517 
3518 // ratio = (1-layers_y[nlays-1][qp_closer])*(1-y_c[qp_closer])/
3519 // ( (1-yold_up[n])*(1-yold_c[qp_closer]) );
3520  ratio = (1-ynew_up[qp_closerup])/
3521  ( (1-yoldup[qp_closeroldup]) );
3522 
3523 // ratioy = (1-ynew_up[qp_closernormup])/
3524 // ( (1-yoldup[qp_closernormoldup]) );
3525  //distance prop to layerup
3526  ynew[n] = ynew_up[qp_closerup]
3527  + (y-yoldup[qp_closeroldup])*ratio;
3528  //ynew[n] = y +abs(nyvert[qp_closernormup])*(ynew_up[qp_closeroldup]-yoldup[qp_closeroldup])*ratioy;
3529 
3530  //ynew[n] = y + 0.3*(ynew_up[qp_closerup]-yoldup[qp_closerup]);
3531  //xnew[n] = x + abs(nxvert[qp_closeroldup])*(xnew_up[qp_closeroldup]-xoldup[qp_closeroldup]);
3532 
3533  if(x> (xmax-xmin)/2. && x< xmax)
3534  {
3535  ratiox = (xmax-xnew_up[qp_closernormup])/
3536  (xmax-xoldup[qp_closernormup]) ;
3537  if( (xmax-xoldup[qp_closernormup])==0)
3538  {
3539  ratiox = 1.0;
3540  }
3541 
3542  //xnew[n] = xnew_up[qp_closerup]
3543  // + (x-xoldup[qp_closerup])*ratiox;
3544  xnew[n] = x + abs(nxvert[qp_closernormup])*(xnew_up[qp_closeroldup]-xoldup[qp_closeroldup])*ratiox;
3545  ASSERTL0(x>xmin," x value <xmin up second half");
3546  ASSERTL0(x<xmax," x value >xmax up second half");
3547  }
3548  else if( x> xmin && x<= (xmax-xmin)/2.)
3549  {
3550 //cout<<"up close normold="<<qp_closernormoldup<<" closenorm="<<qp_closernormup<<endl;
3551  ratiox = (xnew_up[qp_closernormup]-xmin)/
3552  ( (xoldup[qp_closernormup]-xmin) );
3553  if( (xoldup[qp_closernormup]-xmin)==0)
3554  {
3555  ratiox = 1.0;
3556  }
3557  //xnew[n] = xnew_up[qp_closerup]
3558  // + (x-xoldup[qp_closeroldup])*ratiox;
3559  xnew[n] = x + abs(nxvert[qp_closernormup])*(xnew_up[qp_closeroldup]-xoldup[qp_closeroldup])*ratiox;
3560 //cout<<"up xold="<<x<<" xnew="<<xnew[n]<<endl;
3561  ASSERTL0(x>xmin," x value <xmin up first half");
3562  ASSERTL0(x<xmax," x value >xmax up first half");
3563  }
3564 
3565 
3566  }
3567  else if( y< yolddown[qp_closerolddown] && y>-1 )
3568  {
3569 
3570  ratio = (1+ynew_down[qp_closerdown])/
3571  ( (1+yolddown[qp_closerolddown]) );
3572 
3573 // ratioy = (1-ynew_down[qp_closernormdown])/
3574 // ( (1-yolddown[qp_closernormolddown]) );
3575 
3576  //distance prop to layerlow
3577  ynew[n] = ynew_down[qp_closerdown]
3578  + (y-yolddown[qp_closerolddown])*ratio;
3579  //ynew[n] = y +abs(nyvert[qp_closernormdown])*
3580  // (ynew_down[qp_closerolddown]-yolddown[qp_closerolddown])*ratioy;
3581  //ynew[n] = y + 0.3*(ynew_down[qp_closerdown]-yolddown[qp_closerdown]);
3582  //xnew[n] = x + abs(nxvert[qp_closerolddown])*(xnew_down[qp_closerolddown]-xolddown[qp_closerolddown]);
3583 /*
3584 if(n==74)
3585 {
3586 cout<<qp_closerolddown<<" nplaydown="<<qp_closerdown<<endl;
3587 cout<<"xolddown="<<xolddown[qp_closerolddown]<<" xnewdown="<<xnew_down[qp_closerdown]<<endl;
3588 cout<<"xold+"<<x<<" xnew+"<<xnew[n]<<endl;
3589 }
3590 */
3591 
3592  if(x> (xmax-xmin)/2. && x <xmax)
3593  {
3594  ratiox = (xmax-xnew_down[qp_closernormdown])/
3595  ( (xmax-xolddown[qp_closernormdown]) );
3596  if( (xmax-xolddown[qp_closernormdown])==0)
3597  {
3598  ratiox = 1.0;
3599  }
3600  //xnew[n] = xnew_down[qp_closerdown]
3601  // + (x-xolddown[qp_closerolddown])*ratiox;
3602  xnew[n] = x +
3603  abs(nxvert[qp_closernormdown])*(xnew_down[qp_closerolddown]-xolddown[qp_closerolddown])*ratiox;
3604  ASSERTL0(x>xmin," x value <xmin down second half");
3605  ASSERTL0(x<xmax," x value >xmax down second half");
3606  }
3607  else if( x>xmin && x<= (xmax-xmin)/2.)
3608  {
3609  ratiox = (xnew_down[qp_closernormdown]-xmin)/
3610  ( (xolddown[qp_closernormdown]-xmin) );
3611  if( (xolddown[qp_closernormdown]-xmin)==0)
3612  {
3613  ratiox = 1.0;
3614  }
3615  //xnew[n] = xnew_down[qp_closerdown]
3616  // + (x-xolddown[qp_closerolddown])*ratiox;
3617  xnew[n] = x +
3618  abs(nxvert[qp_closernormdown])*(xnew_down[qp_closerolddown]-xolddown[qp_closerolddown])*ratiox;
3619  ASSERTL0(x>xmin," x value <xmin down first half");
3620  ASSERTL0(x<xmax," x value >xmax down first half");
3621  }
3622 
3623  }
3624 
3625 cout<<"xold"<<x<<" xnew="<<xnew[n]<<endl;
3626  ASSERTL0(xnew[n] >= xmin, "newx < xmin");
3627  ASSERTL0(xnew[n]<= xmax, "newx > xmax");
3628 
3629  }// verts closed
3630 }
3631 
3632 void CheckSingularQuads( MultiRegions::ExpListSharedPtr Exp,
3633  Array<OneD, int> V1, Array<OneD, int> V2,
3634  Array<OneD, NekDouble>& xnew,Array<OneD, NekDouble>& ynew)
3635 {
3636  const boost::shared_ptr<LocalRegions::ExpansionVector> exp2D = Exp->GetExp();
3637  int nel = exp2D->size();
3638  LocalRegions::QuadExpSharedPtr locQuadExp;
3639  LocalRegions::TriExpSharedPtr locTriExp;
3640  SpatialDomains::Geometry1DSharedPtr SegGeom;
3641  int idbef, idnext;
3642  NekDouble xV1, yV1, xV2,yV2;
3643  NekDouble slopebef,slopenext,slopenew;
3644  Array<OneD, int> locEids(4);
3645  for(int i=0; i<nel; i++)
3646  {
3647  if((locQuadExp = (*exp2D)[i]->as<LocalRegions::QuadExp>()))
3648  {
3649  SegGeom = (locQuadExp->GetGeom2D())->GetEdge(0);
3650  idbef = SegGeom->GetEid();
3651  if(xnew[ V1[idbef] ]<= xnew[ V2[idbef] ])
3652  {
3653  xV1 = xnew[ V1[idbef] ];
3654  yV1 = ynew[ V1[idbef] ];
3655  xV2 = xnew[ V2[idbef] ];
3656  yV2 = ynew[ V2[idbef] ];
3657  slopebef = (yV2 -yV1)/(xV2 -xV1);
3658  }
3659  else
3660  {
3661  xV1 = xnew[ V2[idbef] ];
3662  yV1 = ynew[ V2[idbef] ];
3663  xV2 = xnew[ V1[idbef] ];
3664  yV2 = ynew[ V1[idbef] ];
3665  slopebef = (yV2 -yV1)/(xV2 -xV1);
3666  }
3667 //cout<<"00 V1 x="<<xnew[ V1[idbef] ]<<" y="<<ynew[ V1[idbef] ]<<endl;
3668 //cout<<"00 V2 x="<<xnew[ V2[idbef] ]<<" y="<<ynew[ V2[idbef] ]<<endl;
3669  for(int j = 1; j < locQuadExp->GetNedges(); ++j)
3670  {
3671  SegGeom = (locQuadExp->GetGeom2D())->GetEdge(j);
3672  idnext = SegGeom->GetEid();
3673 //cout<<"id="<<idnext<<" locid="<<j<<endl;
3674 //cout<<" V1 x="<<xnew[ V1[idnext] ]<<" y="<<ynew[ V1[idnext] ]<<endl;
3675 //cout<<" V2 x="<<xnew[ V2[idnext] ]<<" y="<<ynew[ V2[idnext] ]<<endl;
3676  if(xV1 == xnew[ V1[idnext] ] && yV1 == ynew[ V1[idnext] ] )
3677  {
3678  xV1 = xnew[ V1[idnext] ];
3679  yV1 = ynew[ V1[idnext] ];
3680  xV2 = xnew[ V2[idnext] ];
3681  yV2 = ynew[ V2[idnext] ];
3682  slopenext = (yV2 -yV1)/(xV2 -xV1);
3683 if(i==23)
3684 {
3685 cout<<"case1 x0="<<xV1<<" x1="<<xV2<<endl;
3686 cout<<idnext<<" 11slope bef ="<<slopebef<<" slopenext="<<slopenext<<endl;
3687 }
3688  //compare with slope before
3689  if( slopebef/slopenext>0.84 && slopebef/slopenext <1.18)
3690  {
3691  xnew[ V1[idnext] ] = xnew[ V1[idnext] ] -0.01;
3692  slopenew = (yV2-yV1)/(xV2- xnew[ V1[idnext] ]);
3693 
3694  if( abs(slopebef-slopenew) < abs(slopebef-slopenext) )
3695  {
3696  xnew[ V1[idnext] ] = xnew[ V1[idnext] ] +0.02;
3697  slopenew = (yV2-yV1)/(xV2- xnew[ V1[idnext] ]);
3698  }
3699  slopenext = slopenew;
3700 cout<<"slopenew="<<slopenew<<endl;
3701 cout<<"moved x="<<xnew[ V1[idnext] ]<<endl;
3702  }
3703  }
3704  else if(xV2 == xnew[ V2[idnext] ] && yV2 == ynew[ V2[idnext] ] )
3705  {
3706  xV1 = xnew[ V2[idnext] ];
3707  yV1 = ynew[ V2[idnext] ];
3708  xV2 = xnew[ V1[idnext] ];
3709  yV2 = ynew[ V1[idnext] ];
3710  slopenext = (yV2 -yV1)/(xV2 -xV1);
3711 if(i==23)
3712 {
3713 cout<<"case2 x0="<<xV1<<" x1="<<xV2<<endl;
3714 cout<<idnext<<" 22slope bef ="<<slopebef<<" slopenext="<<slopenext<<endl;
3715 }
3716  //compare with slope before
3717  if( slopebef/slopenext>0.84 && slopebef/slopenext <1.18)
3718  {
3719  xnew[ V2[idnext] ] = xnew[ V2[idnext] ] -0.01;
3720  slopenew = (yV2-yV1)/(xV2- xnew[ V2[idnext] ]);
3721 
3722  if( abs(slopebef-slopenew) < abs(slopebef-slopenext) )
3723  {
3724  xnew[ V2[idnext] ] = xnew[ V2[idnext] ] +0.02;
3725  slopenew = (yV2-yV1)/(xV2- xnew[ V2[idnext] ]);
3726 
3727  }
3728  slopenext = slopenew;
3729 cout<<"slopenew="<<slopenew<<endl;
3730 cout<<"moved x="<<xnew[ V2[idnext] ]<<endl;
3731  }
3732  }
3733  else if(xV1 == xnew[ V2[idnext] ] && yV1 == ynew[ V2[idnext] ] )
3734  {
3735  xV1 = xnew[ V2[idnext] ];
3736  yV1 = ynew[ V2[idnext] ];
3737  xV2 = xnew[ V1[idnext] ];
3738  yV2 = ynew[ V1[idnext] ];
3739  slopenext = (yV2 -yV1)/(xV2 -xV1);
3740 if(i==23)
3741 {
3742 cout<<"case3 x0="<<xV1<<" x1="<<xV2<<endl;
3743 cout<<idnext<<" 22slope bef ="<<slopebef<<" slopenext="<<slopenext<<endl;
3744 }
3745  //compare with slope before
3746  if( slopebef/slopenext>0.84 && slopebef/slopenext <1.18)
3747  {
3748  xnew[ V2[idnext] ] = xnew[ V2[idnext] ] -0.01;
3749  slopenew = (yV2-yV1)/(xV2- xnew[ V2[idnext] ]);
3750 
3751  if( abs(slopebef-slopenew) < abs(slopebef-slopenext) )
3752  {
3753  xnew[ V2[idnext] ] = xnew[ V2[idnext] ] +0.02;
3754  slopenew = (yV2-yV1)/(xV2- xnew[ V2[idnext] ]);
3755  }
3756  slopenext = slopenew;
3757 cout<<"slopenew="<<slopenew<<endl;
3758 cout<<"moved x="<<xnew[ V2[idnext] ]<<endl;
3759  }
3760 
3761  }
3762  else if(xV2 == xnew[ V1[idnext] ] && yV2 == ynew[ V1[idnext] ] )
3763  {
3764  xV1 = xnew[ V1[idnext] ];
3765  yV1 = ynew[ V1[idnext] ];
3766  xV2 = xnew[ V2[idnext] ];
3767  yV2 = ynew[ V2[idnext] ];
3768  slopenext = (yV2 -yV1)/(xV2 -xV1);
3769 if(i==23)
3770 {
3771 cout<<"case4 x0="<<xV1<<" x1="<<xV2<<endl;
3772 cout<<idnext<<" 22slope bef ="<<slopebef<<" slopenext="<<slopenext<<endl;
3773 }
3774  //compare with slope before
3775  if( slopebef/slopenext>0.84 && slopebef/slopenext <1.18)
3776  {
3777  xnew[ V1[idnext] ] = xnew[ V1[idnext] ] -0.01;
3778  slopenew = (yV2-yV1)/(xV2- xnew[ V1[idnext] ]);
3779 
3780  if( abs(slopebef-slopenew) < abs(slopebef-slopenext) )
3781  {
3782  xnew[ V1[idnext] ] = xnew[ V1[idnext] ] +0.02;
3783  slopenew = (yV2-yV1)/(xV2- xnew[ V1[idnext] ]);
3784  }
3785  slopenext = slopenew;
3786 cout<<"slopenew="<<slopenew<<endl;
3787 cout<<"moved x="<<xnew[ V1[idnext] ]<<endl;
3788  }
3789 
3790  }
3791  else
3792  {
3793  ASSERTL0(false, "edge not connected");
3794  }
3795  slopebef = slopenext;
3796 
3797 
3798 
3799  }
3800  }
3801  }
3802 }
3803 
3804 void Replacevertices(string filename, Array<OneD, NekDouble> newx,
3805  Array<OneD, NekDouble> newy,
3806  Array<OneD, NekDouble> xcPhys, Array<OneD, NekDouble> ycPhys,
3807  Array<OneD, int>Eids, int Npoints, string s_alp,
3808  Array<OneD, Array<OneD, NekDouble> > x_lay,
3809  Array<OneD, Array<OneD, NekDouble> > y_lay,
3810  Array<OneD, Array<OneD, int > >lay_eids, bool curv_lay)
3811 {
3812  //load existing file
3813  string newfile;
3814  TiXmlDocument doc(filename);
3815  //load xscale parameter (if exists)
3816  TiXmlElement* master = doc.FirstChildElement("NEKTAR");
3817  TiXmlElement* mesh = master->FirstChildElement("GEOMETRY");
3818  TiXmlElement* element = mesh->FirstChildElement("VERTEX");
3819  NekDouble xscale = 1.0;
3820  LibUtilities::AnalyticExpressionEvaluator expEvaluator;
3821  const char *xscal = element->Attribute("XSCALE");
3822  if(xscal)
3823  {
3824  std::string xscalstr = xscal;
3825  int expr_id = expEvaluator.DefineFunction("",xscalstr);
3826  xscale = expEvaluator.Evaluate(expr_id);
3827  }
3828 
3829 
3830  // Save a new XML file.
3831  newfile = filename.substr(0, filename.find_last_of("."))+"_moved.xml";
3832  doc.SaveFile( newfile );
3833 
3834  //write the new vertices
3835  TiXmlDocument docnew(newfile);
3836  bool loadOkaynew = docnew.LoadFile();
3837 
3838  std::string errstr = "Unable to load file: ";
3839  errstr += newfile;
3840  ASSERTL0(loadOkaynew, errstr.c_str());
3841 
3842  TiXmlHandle docHandlenew(&docnew);
3843  TiXmlElement* meshnew = NULL;
3844  TiXmlElement* masternew = NULL;
3845  TiXmlElement* condnew = NULL;
3846  TiXmlElement* Parsnew = NULL;
3847  TiXmlElement* parnew = NULL;
3848 
3849  // Master tag within which all data is contained.
3850 
3851 
3852  masternew = docnew.FirstChildElement("NEKTAR");
3853  ASSERTL0(masternew, "Unable to find NEKTAR tag in file.");
3854 
3855  //set the alpha value
3856  string alphastring;
3857  condnew = masternew->FirstChildElement("CONDITIONS");
3858  Parsnew = condnew->FirstChildElement("PARAMETERS");
3859 cout<<"alpha="<<s_alp<<endl;
3860  parnew = Parsnew->FirstChildElement("P");
3861  while(parnew)
3862  {
3863  TiXmlNode *node = parnew->FirstChild();
3864  if (node)
3865  {
3866  // Format is "paramName = value"
3867  std::string line = node->ToText()->Value();
3868  std::string lhs;
3869  std::string rhs;
3870  /// Pull out lhs and rhs and eliminate any spaces.
3871  int beg = line.find_first_not_of(" ");
3872  int end = line.find_first_of("=");
3873  // Check for no parameter name
3874  if (beg == end) throw 1;
3875  // Check for no parameter value
3876  if (end != line.find_last_of("=")) throw 1;
3877  // Check for no equals sign
3878  if (end == std::string::npos) throw 1;
3879  lhs = line.substr(line.find_first_not_of(" "), end-beg);
3880  lhs = lhs.substr(0, lhs.find_last_not_of(" ")+1);
3881 
3882  //rhs = line.substr(line.find_last_of("=")+1);
3883  //rhs = rhs.substr(rhs.find_first_not_of(" "));
3884  //rhs = rhs.substr(0, rhs.find_last_not_of(" ")+1);
3885 
3886  boost::to_upper(lhs);
3887  if(lhs == "ALPHA")
3888  {
3889  alphastring = "Alpha = "+ s_alp;
3890  parnew->RemoveChild(node);
3891  parnew->LinkEndChild(new TiXmlText(alphastring) );
3892  }
3893  }
3894 
3895  parnew = parnew->NextSiblingElement("P");
3896  }
3897 
3898 
3899  // Find the Mesh tag and same the dim and space attributes
3900  meshnew = masternew->FirstChildElement("GEOMETRY");
3901 
3902  ASSERTL0(meshnew, "Unable to find GEOMETRY tag in file.");
3903  // Now read the vertices
3904  TiXmlElement* elementnew = meshnew->FirstChildElement("VERTEX");
3905  ASSERTL0(elementnew, "Unable to find mesh VERTEX tag in file.");
3906  //set xscale 1!!
3907  if(xscale!=1.0)
3908  {
3909  elementnew->SetAttribute("XSCALE",1.0);
3910  }
3911  TiXmlElement *vertexnew = elementnew->FirstChildElement("V");
3912 
3913 
3914 
3915  int indx;
3916  int err, numPts;
3917  int nextVertexNumber = -1;
3918 
3919  while (vertexnew)
3920  {
3921  nextVertexNumber++;
3922  //delete the old one
3923  TiXmlAttribute *vertexAttr = vertexnew->FirstAttribute();
3924  std::string attrName(vertexAttr->Name());
3925  ASSERTL0(attrName == "ID", (std::string("Unknown attribute name: ") + attrName).c_str());
3926 
3927  err = vertexAttr->QueryIntValue(&indx);
3928  ASSERTL0(err == TIXML_SUCCESS, "Unable to read attribute ID.");
3929  ASSERTL0(indx == nextVertexNumber, "Vertex IDs must begin with zero and be sequential.");
3930 
3931  std::string vertexBodyStr;
3932  // Now read body of vertex
3933  TiXmlNode *vertexBody = vertexnew->FirstChild();
3934  // Accumulate all non-comment body data.
3935  if (vertexBody->Type() == TiXmlNode::TINYXML_TEXT)
3936  {
3937  vertexBodyStr += vertexBody->ToText()->Value();
3938  vertexBodyStr += " ";
3939  }
3940  ASSERTL0(!vertexBodyStr.empty(), "Vertex definitions must contain vertex data.");
3941  //remove the old coordinates
3942  vertexnew->RemoveChild(vertexBody);
3943  //write the new one
3944 //cout<<"writing.. v:"<<nextVertexNumber<<endl;
3945  stringstream s;
3946  //we need at least 5 digits (setprecision 5) to get the streak position with
3947  // precision 10^-10
3948  s << std::scientific << std::setprecision(8) << newx[nextVertexNumber] << " "
3949  << newy[nextVertexNumber] << " " << 0.0;
3950  vertexnew->LinkEndChild(new TiXmlText(s.str()));
3951  //TiXmlNode *newvertexBody = vertexnew->FirstChild();
3952  //string newvertexbodystr= newvertexBody->SetValue(s.str());
3953  //vertexnew->ReplaceChild(vertexBody,new TiXmlText(newvertexbodystr));
3954 
3955  vertexnew = vertexnew->NextSiblingElement("V");
3956  }
3957 
3958 
3959 
3960  //read the curved tag
3961  TiXmlElement* curvednew = meshnew->FirstChildElement("CURVED");
3962  ASSERTL0(curvednew, "Unable to find mesh CURVED tag in file.");
3963  TiXmlElement *edgenew = curvednew->FirstChildElement("E");
3964  int cnt =-1;
3965  //ID is different from index...
3966  std::string charindex;
3967  int eid;
3968  int index;
3969  int indexeid;
3970  int neids_lay = lay_eids[0].num_elements();
3971  //if edgenew belongs to the crit lay replace it, else delete it.
3972  while (edgenew)
3973  {
3974  indexeid =-1;
3975  cnt++;
3976  //get the index...
3977  TiXmlAttribute *edgeAttr = edgenew->FirstAttribute();
3978  std::string attrName(edgeAttr->Name());
3979  charindex = edgeAttr->Value();
3980  std::istringstream iss(charindex);
3981  iss >> std::dec >> index;
3982  //get the eid
3983  edgenew->QueryIntAttribute("EDGEID", &eid);
3984 //cout<<"eid="<<eid<<" neid="<<Eids.num_elements()<<endl;
3985  //find the corresponding index curve point
3986  for(int u=0; u<Eids.num_elements(); u++)
3987  {
3988 //cout<<"Eids="<<Eids[u]<<" eid="<<eid<<endl;
3989  if(Eids[u]==eid)
3990  {
3991  indexeid = u;
3992  }
3993  }
3994  if(indexeid==-1)
3995  {
3996  curvednew->RemoveChild(edgenew);
3997  //ASSERTL0(false, "edge to update not found");
3998  }
3999  else
4000  {
4001 
4002  std::string edgeBodyStr;
4003  //read the body of the edge
4004  TiXmlNode *edgeBody = edgenew->FirstChild();
4005  if(edgeBody->Type() == TiXmlNode::TINYXML_TEXT)
4006  {
4007  edgeBodyStr += edgeBody->ToText()->Value();
4008  edgeBodyStr += " ";
4009  }
4010  ASSERTL0(!edgeBodyStr.empty(), "Edge definitions must contain edge data");
4011  //remove the old coordinates
4012  edgenew->RemoveChild(edgeBody);
4013  //write the new points coordinates
4014  //we need at least 5 digits (setprecision 5) to get the streak position with
4015  // precision 10^-10
4016 
4017  //Determine the number of points
4018  err = edgenew->QueryIntAttribute("NUMPOINTS", &numPts);
4019  ASSERTL0(err == TIXML_SUCCESS, "Unable to read curve attribute NUMPOINTS.");
4020 
4021  stringstream st;
4022  edgenew->SetAttribute("NUMPOINTS", Npoints);
4023  for(int u=0; u< Npoints; u++)
4024  {
4025  st << std::scientific <<
4026  std::setprecision(8) <<xcPhys[cnt*Npoints+u]
4027  << " " << ycPhys[cnt*Npoints+u] << " " << 0.000<<" ";
4028  }
4029 
4030  edgenew->LinkEndChild(new TiXmlText(st.str()));
4031 
4032 
4033 /*
4034  st << std::scientific << std::setprecision(8) << x_crit[v1] << " "
4035  << y_crit[v1] << " " << 0.000<<" ";
4036  for(int a=0; a< Npoints-2; a++)
4037  {
4038  st << std::scientific << std::setprecision(8) <<
4039  " "<<Pcurvx[indexeid*(Npoints-2) +a]<<" "<<Pcurvy[indexeid*(Npoints-2) +a]
4040  <<" "<<0.000<<" ";
4041 
4042  }
4043  st << std::scientific << std::setprecision(8) <<
4044  " "<<x_crit[v2]<<" "<< y_crit[v2] <<" "<< 0.000;
4045  edgenew->LinkEndChild(new TiXmlText(st.str()));
4046 
4047 */
4048 
4049 
4050  }
4051 
4052  edgenew = edgenew->NextSiblingElement("E");
4053 
4054  }
4055 
4056  //write also the others layers curve points
4057  if(curv_lay == true)
4058  {
4059 cout<<"write other curved edges"<<endl;
4060  TiXmlElement * curved = meshnew->FirstChildElement("CURVED");
4061  int idcnt = 300;
4062  int nlays = lay_eids.num_elements();
4063 
4064  //TiXmlComment * comment = new TiXmlComment();
4065  //comment->SetValue(" new edges ");
4066  //curved->LinkEndChild(comment);
4067  for (int g=0; g< nlays; ++g)
4068  {
4069  for(int p=0; p< neids_lay; p++)
4070  {
4071  stringstream st;
4072  TiXmlElement * e = new TiXmlElement( "E" );
4073  e->SetAttribute("ID", idcnt++);
4074  e->SetAttribute("EDGEID", lay_eids[g][p]);
4075  e->SetAttribute("NUMPOINTS", Npoints);
4076  e->SetAttribute("TYPE", "PolyEvenlySpaced");
4077  for(int c=0; c< Npoints; c++)
4078  {
4079  st << std::scientific << std::setprecision(8) <<x_lay[g][p*Npoints +c]
4080  << " " << y_lay[g][p*Npoints +c] << " " << 0.000<<" ";
4081 
4082  }
4083 
4084  TiXmlText * t0 = new TiXmlText(st.str());
4085  e->LinkEndChild(t0);
4086  curved->LinkEndChild(e);
4087  }
4088 
4089  }
4090  }
4091 
4092 
4093 
4094  docnew.SaveFile( newfile );
4095 
4096  cout<<"new file: "<<newfile<<endl;
4097 }
PolyFit
void PolyFit(int polyorder, int npoints, Array< OneD, NekDouble > xin, Array< OneD, NekDouble > fin, Array< OneD, NekDouble > &coeffsinterp, Array< OneD, NekDouble > &xout, Array< OneD, NekDouble > &fout, int npout)
Definition: MeshMove.cpp:2927
checkcommonvert
bool checkcommonvert(Array< OneD, int > Vids_laybefore, Array< OneD, int > Vids_c, int Vid)
Definition: MeshMove.cpp:2653
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:188
Nektar::MultiRegions::ContField1DSharedPtr
boost::shared_ptr< ContField1D > ContField1DSharedPtr
Definition: ContField1D.h:237
OrderVertices
void OrderVertices(int nedges, SpatialDomains::MeshGraphSharedPtr graphShPt, MultiRegions::ExpListSharedPtr &bndfield, Array< OneD, int > &Vids, int v1, int v2, NekDouble x_connect, int &lastedge, Array< OneD, NekDouble > &x, Array< OneD, NekDouble > &y)
Definition: MeshMove.cpp:1936
GenerateMapEidsv1v2
void GenerateMapEidsv1v2(MultiRegions::ExpListSharedPtr field, Array< OneD, int > &V1, Array< OneD, int > &V2)
Definition: MeshMove.cpp:2219
Vmath::Vsqrt
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
Definition: Vmath.cpp:394
MoveOutsidePointsNnormpos
void MoveOutsidePointsNnormpos(int npedge, SpatialDomains::MeshGraphSharedPtr mesh, Array< OneD, NekDouble > xcold, Array< OneD, NekDouble > ycold, Array< OneD, NekDouble > xolddown, Array< OneD, NekDouble > yolddown, Array< OneD, NekDouble > xoldup, Array< OneD, NekDouble > yoldup, Array< OneD, NekDouble > xlaydown, Array< OneD, NekDouble > ylaydown, Array< OneD, NekDouble > xlayup, Array< OneD, NekDouble > ylayup, Array< OneD, NekDouble > nxPhys, Array< OneD, NekDouble > nyPhys, Array< OneD, NekDouble > &xnew, Array< OneD, NekDouble > &ynew)
Definition: MeshMove.cpp:3347
GenerateAddPointsNewtonIt
void GenerateAddPointsNewtonIt(NekDouble xi, NekDouble yi, NekDouble &xout, NekDouble &yout, MultiRegions::ExpListSharedPtr function, Array< OneD, NekDouble > derfunction, NekDouble cr)
Definition: MeshMove.cpp:2139
Vmath::Vmax
T Vmax(int n, const T *x, const int incx)
Return the maximum element in x – called vmax to avoid conflict with max.
Definition: Vmath.cpp:765
Nektar::MemoryManager
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
Definition: NekMemoryManager.hpp:102
Vmath::Vmin
T Vmin(int n, const T *x, const int incx)
Return the minimum element in x - called vmin to avoid conflict with min.
Definition: Vmath.cpp:857
ContField2D.h
MoveLayerNnormpos
void MoveLayerNnormpos(int nvertl, int npedge, Array< OneD, NekDouble > xcPhys, Array< OneD, NekDouble > tmpx_lay, Array< OneD, NekDouble > tmpy_lay, Array< OneD, int > Vids, Array< OneD, NekDouble > &xlay, Array< OneD, NekDouble > &ylay, Array< OneD, NekDouble > &xnew, Array< OneD, NekDouble > &ynew)
Definition: MeshMove.cpp:3213
Vmath::Imin
int Imin(int n, const T *x, const int incx)
Return the index of the minimum element in x.
Definition: Vmath.cpp:833
Vmath::Vvtvp
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.cpp:428
std
STL namespace.
Vmath::Sdiv
void Sdiv(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha/y.
Definition: Vmath.cpp:257
Vmath::Vdiv
void Vdiv(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x/y.
Definition: Vmath.cpp:227
Nektar::LibUtilities::SessionReaderSharedPtr
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
Vmath::Vabs
void Vabs(int n, const T *x, const int incx, T *y, const int incy)
vabs: y = |x|
Definition: Vmath.cpp:410
MappingEVids
void MappingEVids(Array< OneD, NekDouble > xoldup, Array< OneD, NekDouble > yoldup, Array< OneD, NekDouble > xolddown, Array< OneD, NekDouble > yolddown, Array< OneD, NekDouble > xcold, Array< OneD, NekDouble > ycold, Array< OneD, int > Vids_c, SpatialDomains::MeshGraphSharedPtr mesh, MultiRegions::ExpListSharedPtr streak, Array< OneD, int > V1, Array< OneD, int > V2, int &nlays, Array< OneD, Array< OneD, int > > &Eids_lay, Array< OneD, Array< OneD, int > > &Vids_lay)
Definition: MeshMove.cpp:2289
Nektar::StdRegions::StdExpansion::GetCoords
void GetCoords(Array< OneD, NekDouble > &coords_1, Array< OneD, NekDouble > &coords_2=NullNekDouble1DArray, Array< OneD, NekDouble > &coords_3=NullNekDouble1DArray)
this function returns the physical coordinates of the quadrature points of the expansion ...
Definition: StdExpansion.h:676
Nektar::LibUtilities::Import
void Import(const std::string &infilename, std::vector< FieldDefinitionsSharedPtr > &fielddefs, std::vector< std::vector< NekDouble > > &fielddata, FieldMetaDataMap &fieldinfomap, const Array< OneD, int > ElementiDs)
Imports an FLD file.
Definition: FieldIO.cpp:115
Cutrepetitions
void Cutrepetitions(int nedges, Array< OneD, NekDouble > inarray, Array< OneD, NekDouble > &outarray)
Definition: MeshMove.cpp:2669
ExpList1D.h
NektarUnivTypeDefs.hpp
main
int main(int argc, char *argv[])
Definition: MeshMove.cpp:173
Nektar::LocalRegions::SegExpSharedPtr
boost::shared_ptr< SegExp > SegExpSharedPtr
Definition: SegExp.h:266
TriExp.h
Vmath::Smul
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:199
MoveOutsidePointsfixedxpos
void MoveOutsidePointsfixedxpos(int npedge, SpatialDomains::MeshGraphSharedPtr mesh, Array< OneD, NekDouble > xcold, Array< OneD, NekDouble > ycold, Array< OneD, NekDouble > xolddown, Array< OneD, NekDouble > yolddown, Array< OneD, NekDouble > xoldup, Array< OneD, NekDouble > yoldup, Array< OneD, NekDouble > ylaydown, Array< OneD, NekDouble > ylayup, Array< OneD, NekDouble > &xnew, Array< OneD, NekDouble > &ynew)
Definition: MeshMove.cpp:3279
ExpList.h
Nektar::MultiRegions::ExpList1DSharedPtr
boost::shared_ptr< ExpList1D > ExpList1DSharedPtr
Shared pointer to an ExpList1D object.
Definition: ExpList1D.h:50
Lapack.hpp
LagrangeInterpolant
NekDouble LagrangeInterpolant(NekDouble x, int npts, Array< OneD, NekDouble > xpts, Array< OneD, NekDouble > funcvals)
Definition: MeshMove.cpp:2765
CheckSingularQuads
void CheckSingularQuads(MultiRegions::ExpListSharedPtr Exp, Array< OneD, int > V1, Array< OneD, int > V2, Array< OneD, NekDouble > &xnew, Array< OneD, NekDouble > &ynew)
Definition: MeshMove.cpp:3632
EvaluateTangent
void EvaluateTangent(int npoints, Array< OneD, NekDouble > xcQedge, Array< OneD, NekDouble > coeffsinterp, Array< OneD, NekDouble > &txQedge, Array< OneD, NekDouble > &tyQedge)
Definition: MeshMove.cpp:2793
Nektar::SpatialDomains::BoundaryRegionCollection
std::map< int, BoundaryRegionShPtr > BoundaryRegionCollection
Definition: Conditions.h:206
Nektar::MultiRegions::ExpListSharedPtr
boost::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
Definition: LocTraceToTraceMap.h:56
npts
static std::string npts
Definition: InputFld.cpp:43
Computestreakpositions
void Computestreakpositions(int nvertl, MultiRegions::ExpListSharedPtr streak, Array< OneD, NekDouble > xold_up, Array< OneD, NekDouble > yold_up, Array< OneD, NekDouble > xold_low, Array< OneD, NekDouble > yold_low, Array< OneD, NekDouble > xold_c, Array< OneD, NekDouble > yold_c, Array< OneD, NekDouble > &xc, Array< OneD, NekDouble > &yc, NekDouble cr, bool verts)
Definition: MeshMove.cpp:2023
Nektar::LocalRegions::QuadExpSharedPtr
boost::shared_ptr< QuadExp > QuadExpSharedPtr
Definition: QuadExp.h:293
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:44
SegExp.h
QuadExp.h
Vertex
Represents a vertex in the mesh.
Definition: makeEllipsoidNektar.c:22
Vmath::Sadd
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add vector y = alpha + x.
Definition: Vmath.cpp:301
MoveLayerNfixedxpos
void MoveLayerNfixedxpos(int nvertl, int npedge, Array< OneD, NekDouble > xcPhys, Array< OneD, NekDouble > tmpx_lay, Array< OneD, NekDouble > tmpy_lay, Array< OneD, int > Vids, Array< OneD, NekDouble > &xlay, Array< OneD, NekDouble > &ylay, Array< OneD, NekDouble > &xnew, Array< OneD, NekDouble > &ynew)
Definition: MeshMove.cpp:3134
Orderfunctionx
void Orderfunctionx(Array< OneD, NekDouble > inarray_x, Array< OneD, NekDouble > inarray_y, Array< OneD, NekDouble > &outarray_x, Array< OneD, NekDouble > &outarray_y)
Definition: MeshMove.cpp:3045
Nektar::SpatialDomains::BoundaryConditionsSharedPtr
boost::shared_ptr< BoundaryConditions > BoundaryConditionsSharedPtr
Definition: Conditions.h:271
Nektar::SpatialDomains::Geometry1DSharedPtr
boost::shared_ptr< Geometry1D > Geometry1DSharedPtr
Definition: Geometry1D.h:48
ContField1D.h
Nektar::StdRegions::find
InputIterator find(InputIterator first, InputIterator last, InputIterator startingpoint, const EqualityComparable &value)
Definition: StdRegions.hpp:315
DetermineclosePointxindex
int DetermineclosePointxindex(NekDouble x, Array< OneD, NekDouble > xArray)
Definition: MeshMove.cpp:2698
Nektar::SpatialDomains::BoundaryConditions::GetBoundaryRegions
const BoundaryRegionCollection & GetBoundaryRegions(void) const
Definition: Conditions.h:227
MoveLayersvertically
void MoveLayersvertically(int nlays, int nvertl, int cntlow, int cntup, Array< OneD, Array< OneD, int > > lay_Vids, Array< OneD, NekDouble > xc, Array< OneD, NekDouble > yc, Array< OneD, int > Down, Array< OneD, int > Up, Array< OneD, NekDouble > &xnew, Array< OneD, NekDouble > &ynew, Array< OneD, Array< OneD, NekDouble > > &layers_x, Array< OneD, Array< OneD, NekDouble > > &layers_y)
Definition: MeshMove.cpp:3090
PolyInterp
void PolyInterp(Array< OneD, NekDouble > xpol, Array< OneD, NekDouble > ypol, Array< OneD, NekDouble > &coeffsinterp, Array< OneD, NekDouble > &xcout, Array< OneD, NekDouble > &ycout, int edge, int npedge)
Definition: MeshMove.cpp:2836
Vmath::Zero
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:359
Nektar::SpatialDomains::MeshGraphSharedPtr
boost::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:442
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047
Replacevertices
void Replacevertices(string filename, Array< OneD, NekDouble > newx, Array< OneD, NekDouble > newy, Array< OneD, NekDouble > xcPhys, Array< OneD, NekDouble > ycPhys, Array< OneD, int >Eids, int Npoints, string s_alp, Array< OneD, Array< OneD, NekDouble > > x_lay, Array< OneD, Array< OneD, NekDouble > > y_lay, Array< OneD, Array< OneD, int > >lay_eids, bool curv_lay)
Definition: MeshMove.cpp:3804
Nektar::SpatialDomains::PointGeomSharedPtr
boost::shared_ptr< PointGeom > PointGeomSharedPtr
Definition: Geometry.h:60
Vmath::Vadd
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:285
Vmath::Vmul
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:169
GenerateNeighbourArrays
void GenerateNeighbourArrays(int index, int neighpoints, Array< OneD, NekDouble > xArray, Array< OneD, NekDouble > yArray, Array< OneD, NekDouble > &Neighbour_x, Array< OneD, NekDouble > &Neighbour_y)
Definition: MeshMove.cpp:2716
Nektar::LocalRegions::TriExpSharedPtr
boost::shared_ptr< TriExp > TriExpSharedPtr
Definition: TriExp.h:291
Generated on Fri Dec 2 2016 13:10:49 for Nektar++ by   doxygen 1.8.8