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