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