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