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