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