Nektar++
MeshMove.cpp
Go to the documentation of this file.
1 
2 ///////////////////////////////////////////////////////////////////////////////
3 //
4 // File: MeshMove.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 
1892  // MoveOutsidePointsfixedxpos(npedge, graphShPt,xold_c, yold_c,
1893  // xold_low, yold_low, xold_up, yold_up, layers_y[0],
1894  // layers_y[nlays-1],
1895  // xnew, ynew);
1896  // lastIregion -1 = laydown
1897  // lastIregion -2 = layup
1898  MoveOutsidePointsNnormpos(
1899  npedge, graphShPt, xold_c, yold_c, xold_low, yold_low, xold_up,
1900  yold_up, layers_x[0], layers_y[0], layers_x[nlays - 1],
1901  layers_y[nlays - 1], nxPhys, nyPhys, xnew, ynew);
1902 
1903  /*
1904  //update vertices coords outside layers region
1905  NekDouble ratio;
1906  for(int n=0; n<nVertTot; n++)
1907  {
1908  NekDouble ratio;
1909  SpatialDomains::PointGeomSharedPtr vertex = graphShPt->GetVertex(n);
1910  NekDouble x,y,z;
1911  vertex->GetCoords(x,y,z);
1912  int qp_closer;
1913  NekDouble diff;
1914  int qp_closerup, qp_closerdown;
1915  NekDouble diffup, diffdown;
1916  //determine the closer xold_up
1917  NekDouble tmp=1000;
1918  diffup =1000;
1919  diffdown = 1000;
1920  for(int k=0; k<nvertl; k++)
1921  {
1922  if(abs(x-xold_c[k]) < tmp)
1923  {
1924  tmp = abs(x-xold_c[k]);
1925  qp_closer=k;
1926  }
1927  }
1928 
1929  //find nplay_closer
1930  int nplay_closer;
1931  if(qp_closer==0)
1932  {
1933  nplay_closer=0;//first vert
1934  }
1935  else
1936  {
1937  nplay_closer= (qp_closer-1)*npedge +npedge-1;
1938  }
1939 
1940  if( y>yold_up[qp_closer] && y<1 )//nlays-1 is layerup
1941  {
1942 
1943  // ratio =
1944  (1-layers_y[nlays-1][qp_closer])*(1-y_c[qp_closer])/
1945  // ( (1-yold_up[n])*(1-yold_c[qp_closer]) );
1946  ratio = (1-layers_y[nlays-1][nplay_closer])/
1947  ( (1-yold_up[qp_closer]) );
1948  //distance prop to layerup
1949  ynew[n] = layers_y[nlays-1][nplay_closer]
1950  + (y-yold_up[qp_closer])*ratio;
1951  xnew[n] = x;
1952 
1953  }
1954  else if( y< yold_low[qp_closer] && y>-1 )//0 is layerdown
1955  {
1956 
1957  ratio = (1+layers_y[0][nplay_closer])/
1958  ( (1+yold_low[qp_closer]) );
1959  //distance prop to layerlow
1960  ynew[n] = layers_y[0][nplay_closer]
1961  + (y-yold_low[qp_closer])*ratio;
1962  xnew[n] = x;
1963  }
1964 
1965  }
1966  */
1967 
1968  /*
1969  //update verts coords of critlay(in case EvaluateLayTnaget has been
1970  called)
1971  //it's not necessary !!!
1972  for(int e=0; e<nvertl; e++)
1973  {
1974  ynew[CritLay[e]] = y_c[e];
1975  }
1976  */
1977 
1978  } // move_norm bool
1979  else // move vertically
1980  {
1981  MoveLayersvertically(nlays, nvertl, cntlow, cntup, lay_Vids, x_c, y_c,
1982  Down, Up, xnew, ynew, layers_x, layers_y);
1983  }
1984 
1985  // check singular quads:
1986  CheckSingularQuads(streak, V1, V2, xnew, ynew);
1987  // check borders of the new mesh verts:
1988  // cout<<std::setprecision(8)<<"yoldmax="<<Vmath::Vmax(nVertTot,
1989  // yold,1)<<endl;
1990  // cout<<std::setprecision(8)<<"ynewmax="<<Vmath::Vmax(nVertTot,ynew,1)<<endl;
1991 
1992  cout << std::setprecision(8) << "xmin=" << Vmath::Vmin(nVertTot, xnew, 1)
1993  << endl;
1994  ASSERTL0(Vmath::Vmin(nVertTot, xold, 1) == Vmath::Vmin(nVertTot, xnew, 1),
1995  " different xmin val");
1996  ASSERTL0(Vmath::Vmin(nVertTot, yold, 1) == Vmath::Vmin(nVertTot, ynew, 1),
1997  " different ymin val");
1998  ASSERTL0(Vmath::Vmax(nVertTot, xold, 1) == Vmath::Vmax(nVertTot, xnew, 1),
1999  " different xmax val");
2000  ASSERTL0(Vmath::Vmax(nVertTot, yold, 1) == Vmath::Vmax(nVertTot, ynew, 1),
2001  " different ymax val");
2002 
2003  // replace the vertices with the new ones
2004 
2005  Replacevertices(changefile, xnew, ynew, xcPhys, ycPhys, Eids, npedge,
2006  charalp, layers_x, layers_y, lay_Eids, curv_lay);
2007 }
2008 
2009 void OrderVertices(int nedges, SpatialDomains::MeshGraphSharedPtr graphShPt,
2010  MultiRegions::ExpListSharedPtr &bndfield,
2011  Array<OneD, int> &Vids, int v1, int v2, NekDouble x_connect,
2012  int &lastedge, Array<OneD, NekDouble> &x,
2013  Array<OneD, NekDouble> &y)
2014 {
2015  int nvertl = nedges + 1;
2016  int edge;
2017  Array<OneD, int> Vids_temp(nvertl, -10);
2018  NekDouble x0layer = 0.0;
2019  for (int j = 0; j < nedges; j++)
2020  {
2021  LocalRegions::SegExpSharedPtr bndSegExplow =
2022  bndfield->GetExp(j)->as<LocalRegions::SegExp>();
2023  edge = (bndSegExplow->GetGeom1D())->GetGlobalID();
2024  // cout<<" edge="<<edge<<endl;
2025  for (int k = 0; k < 2; k++)
2026  {
2027  Vids_temp[j + k] = (bndSegExplow->GetGeom1D())->GetVid(k);
2028  SpatialDomains::PointGeomSharedPtr vertex =
2029  graphShPt->GetVertex(Vids_temp[j + k]);
2030  NekDouble x1, y1, z1;
2031  vertex->GetCoords(x1, y1, z1);
2032  if (x1 == x_connect && edge != lastedge)
2033  {
2034  // first 2 points
2035  if (x_connect == x0layer)
2036  {
2037  Vids[v1] = Vids_temp[j + k];
2038  x[v1] = x1;
2039  y[v1] = y1;
2040 
2041  if (k == 0)
2042  {
2043  Vids_temp[j + 1] =
2044  (bndSegExplow->GetGeom1D())->GetVid(1);
2045  Vids[v2] = Vids_temp[j + 1];
2046  SpatialDomains::PointGeomSharedPtr vertex =
2047  graphShPt->GetVertex(Vids[v2]);
2048  NekDouble x2, y2, z2;
2049  vertex->GetCoords(x2, y2, z2);
2050  x[v2] = x2;
2051  y[v2] = y2;
2052  }
2053  if (k == 1)
2054  {
2055  Vids_temp[j + 0] =
2056  (bndSegExplow->GetGeom1D())->GetVid(0);
2057  Vids[v2] = Vids_temp[j + 0];
2058  SpatialDomains::PointGeomSharedPtr vertex =
2059  graphShPt->GetVertex(Vids[v2]);
2060  NekDouble x2, y2, z2;
2061  vertex->GetCoords(x2, y2, z2);
2062  x[v2] = x2;
2063  y[v2] = y2;
2064  }
2065  }
2066  else // other vertices
2067  {
2068  if (k == 0)
2069  {
2070  Vids_temp[j + 1] =
2071  (bndSegExplow->GetGeom1D())->GetVid(1);
2072  Vids[v1] = Vids_temp[j + 1];
2073  SpatialDomains::PointGeomSharedPtr vertex =
2074  graphShPt->GetVertex(Vids[v1]);
2075  NekDouble x1, y1, z1;
2076  vertex->GetCoords(x1, y1, z1);
2077  x[v1] = x1;
2078  y[v1] = y1;
2079  }
2080  if (k == 1)
2081  {
2082  Vids_temp[j + 0] =
2083  (bndSegExplow->GetGeom1D())->GetVid(0);
2084  Vids[v1] = Vids_temp[j + 0];
2085  SpatialDomains::PointGeomSharedPtr vertex =
2086  graphShPt->GetVertex(Vids[v1]);
2087  NekDouble x1, y1, z1;
2088  vertex->GetCoords(x1, y1, z1);
2089  x[v1] = x1;
2090  y[v1] = y1;
2091  }
2092  }
2093  break;
2094  }
2095  }
2096  if (Vids[v1] != -10)
2097  {
2098  lastedge = edge;
2099  break;
2100  }
2101  }
2102 }
2103 
2104 void Computestreakpositions(
2105  int npoints, MultiRegions::ExpListSharedPtr streak,
2106  Array<OneD, NekDouble> xold_up, Array<OneD, NekDouble> yold_up,
2107  Array<OneD, NekDouble> xold_low, Array<OneD, NekDouble> yold_low,
2108  Array<OneD, NekDouble> xold_c, Array<OneD, NekDouble> yold_c,
2109  Array<OneD, NekDouble> &xc, Array<OneD, NekDouble> &yc, NekDouble cr,
2110  bool verts)
2111 {
2112  boost::ignore_unused(xold_up, xold_low);
2113 
2114  cout << "Computestreakpositions" << endl;
2115  int nq = streak->GetTotPoints();
2116  Array<OneD, NekDouble> coord(2);
2117  // Array<OneD, NekDouble> stvalues(nvertl,-10);
2118  Array<OneD, NekDouble> derstreak(nq);
2119  streak->PhysDeriv(MultiRegions::eY, streak->GetPhys(), derstreak);
2120  int elmtid, offset;
2121  NekDouble U = 0.0, dU = 0.0;
2122  NekDouble F = 1000;
2123 
2124  if (verts == true) // only for verts makes sense to init the coord values..
2125  {
2126 
2127  // start guess
2128  // yc= (yup+ydown)/2
2129  Vmath::Vadd(xc.size(), yold_up, 1, yold_low, 1, yc, 1);
2130  Vmath::Smul(xc.size(), 0.5, yc, 1, yc, 1);
2131  Vmath::Vcopy(xc.size(), xold_c, 1, xc, 1);
2132  }
2133  else // case of xQ,yQ
2134  {
2135  Vmath::Vcopy(xc.size(), xold_c, 1, xc, 1);
2136  Vmath::Vcopy(xc.size(), yold_c, 1, yc, 1);
2137  }
2138 
2139  int its;
2140  int attempt;
2141  NekDouble tol = 1e-3;
2142  NekDouble ytmp;
2143  for (int e = 0; e < npoints; e++)
2144  {
2145  coord[0] = xc[e];
2146  coord[1] = yc[e];
2147 
2148  elmtid = streak->GetExpIndex(coord, 0.00001);
2149  offset = streak->GetPhys_Offset(elmtid);
2150  F = 1000;
2151  its = 0;
2152  attempt = 0;
2153  ytmp = coord[1];
2154  // cout<<"start guess: x="<<xc[e]<<" y="<<yc[e]<<endl;
2155  while (abs(F) > 0.000000001)
2156  {
2157 
2158  elmtid = streak->GetExpIndex(coord, 0.00001);
2159 
2160  // stvalues[e] = streak->GetExp(elmtid)->PhysEvaluate(coord,
2161  // streak->GetPhys() +offset );
2162  // cout<<"elmtid="<<elmtid<<" x="<<coord[0]<<" y="<<coord[1]<<"
2163  // stvalue="<<U<<" dU="<<dU<<endl;
2164 
2165  if ((abs(coord[1]) > 1 || elmtid == -1) && attempt == 0 &&
2166  verts == true)
2167  {
2168  // try the old crit lay position:
2169  coord[1] = yold_c[e];
2170  attempt++;
2171  }
2172  else if ((abs(coord[1]) > 1 || elmtid == -1))
2173  {
2174  coord[1] = ytmp + 0.01;
2175  elmtid = streak->GetExpIndex(coord, 0.001);
2176  offset = streak->GetPhys_Offset(elmtid);
2177  NekDouble Utmp = streak->GetExp(elmtid)->PhysEvaluate(
2178  coord, streak->GetPhys() + offset);
2179  NekDouble dUtmp = streak->GetExp(elmtid)->PhysEvaluate(
2180  coord, derstreak + offset);
2181  coord[1] = coord[1] - (Utmp - cr) / dUtmp;
2182  if ((abs(Utmp - cr) > abs(F)) || (abs(coord[1]) > 1))
2183  {
2184  coord[1] = ytmp - 0.01;
2185  }
2186 
2187  attempt++;
2188  }
2189  else
2190  {
2191  ASSERTL0(abs(coord[1]) <= 1, " y value out of bound +/-1");
2192  }
2193  offset = streak->GetPhys_Offset(elmtid);
2194  U = streak->GetExp(elmtid)->PhysEvaluate(coord, streak->GetPhys() +
2195  offset);
2196  dU =
2197  streak->GetExp(elmtid)->PhysEvaluate(coord, derstreak + offset);
2198  coord[1] = coord[1] - (U - cr) / dU;
2199  F = U - cr;
2200  ASSERTL0(coord[0] == xc[e], " x coordinate must remain the same");
2201 
2202  its++;
2203  if (its > 200 && abs(F) < 0.00001)
2204  {
2205  cout << "warning streak position obtained with precision:" << F
2206  << endl;
2207  break;
2208  }
2209  else if (its > 1000 && abs(F) < 0.0001)
2210  {
2211  cout << "warning streak position obtained with precision:" << F
2212  << endl;
2213  break;
2214  }
2215  else if (its > 1000)
2216  {
2217  ASSERTL0(false, "no convergence after 1000 iterations");
2218  }
2219  }
2220  yc[e] = coord[1] - (U - cr) / dU;
2221  ASSERTL0(U <= cr + tol, "streak wrong+");
2222  ASSERTL0(U >= cr - tol, "streak wrong-");
2223  // Utilities::Zerofunction(coord[0], coord[1], xtest, ytest, streak,
2224  // derstreak);
2225  cout << "result streakvert x=" << xc[e] << " y=" << yc[e]
2226  << " streak=" << U << endl;
2227  }
2228 }
2229 
2230 void GenerateAddPointsNewtonIt(NekDouble xi, NekDouble yi, NekDouble &xout,
2231  NekDouble &yout,
2232  MultiRegions::ExpListSharedPtr function,
2233  Array<OneD, NekDouble> derfunction, NekDouble cr)
2234 {
2235  int elmtid, offset;
2236  NekDouble F = 1000;
2237  NekDouble U = 0.0;
2238  NekDouble dU = 0.0;
2239 
2240  Array<OneD, NekDouble> coords(2);
2241  coords[0] = xi;
2242  coords[1] = yi;
2243 
2244  int attempt = 0;
2245  int its = 0;
2246  NekDouble ytmp = coords[1];
2247  while (abs(F) > 0.00000001)
2248  {
2249 
2250  // cout<<"generate newton it xi="<<xi<<" yi="<<yi<<endl;
2251  elmtid = function->GetExpIndex(coords, 0.01);
2252  //@to do if GetType(elmtid)==triangular WRONG!!!
2253  cout << "gen newton xi=" << xi << " yi=" << coords[1]
2254  << " elmtid=" << elmtid << " F=" << F << endl;
2255 
2256  if ((abs(coords[1]) > 1 || elmtid == -1))
2257  {
2258 
2259  coords[1] = ytmp + 0.01;
2260  elmtid = function->GetExpIndex(coords, 0.01);
2261  offset = function->GetPhys_Offset(elmtid);
2262  NekDouble Utmp = function->GetExp(elmtid)->PhysEvaluate(
2263  coords, function->GetPhys() + offset);
2264  NekDouble dUtmp = function->GetExp(elmtid)->PhysEvaluate(
2265  coords, derfunction + offset);
2266  coords[1] = coords[1] - (Utmp - cr) / dUtmp;
2267  cout << "attempt:" << coords[1] << endl;
2268  if ((abs(Utmp - cr) > abs(F)) || (abs(coords[1]) > 1.01))
2269  {
2270  coords[1] = ytmp - 0.01;
2271  }
2272 
2273  attempt++;
2274  }
2275  else if (abs(coords[1]) < 1.01 && attempt == 0)
2276  {
2277  coords[1] = 1.0;
2278  attempt++;
2279  }
2280  else
2281  {
2282  ASSERTL0(abs(coords[1]) <= 1.00, " y value out of bound +/-1");
2283  }
2284  offset = function->GetPhys_Offset(elmtid);
2285  U = function->GetExp(elmtid)->PhysEvaluate(coords, function->GetPhys() +
2286  offset);
2287  dU = function->GetExp(elmtid)->PhysEvaluate(coords,
2288  derfunction + offset);
2289  coords[1] = coords[1] - (U - cr) / dU;
2290  cout << cr << "U-cr=" << U - cr << " tmp result y:" << coords[1]
2291  << " dU=" << dU << endl;
2292  F = U - cr;
2293 
2294  its++;
2295  if (its > 200 && abs(F) < 0.00001)
2296  {
2297  cout << "warning streak position obtained with precision:" << F
2298  << endl;
2299  break;
2300  }
2301  else if (its > 1000 && abs(F) < 0.0001)
2302  {
2303  cout << "warning streak position obtained with precision:" << F
2304  << endl;
2305  break;
2306  }
2307  else if (its > 1000)
2308  {
2309  ASSERTL0(false, "no convergence after 1000 iterations");
2310  }
2311 
2312  ASSERTL0(coords[0] == xi, " x coordinate must remain the same");
2313  }
2314  xout = xi;
2315  yout = coords[1] - (U - cr) / dU;
2316  cout << "NewtonIt result x=" << xout << " y=" << coords[1] << " U=" << U
2317  << endl;
2318 }
2319 
2320 void GenerateMapEidsv1v2(MultiRegions::ExpListSharedPtr field,
2321  Array<OneD, int> &V1, Array<OneD, int> &V2)
2322 {
2323 
2324  const std::shared_ptr<LocalRegions::ExpansionVector> exp2D =
2325  field->GetExp();
2326  int nel = exp2D->size();
2327  LocalRegions::QuadExpSharedPtr locQuadExp;
2328  LocalRegions::TriExpSharedPtr locTriExp;
2329  SpatialDomains::Geometry1DSharedPtr SegGeom;
2330  int id;
2331  int cnt = 0;
2332  Array<OneD, int> V1tmp(4 * nel, 10000);
2333  Array<OneD, int> V2tmp(4 * nel, 10000);
2334  for (int i = 0; i < nel; i++)
2335  {
2336  if ((locQuadExp = (*exp2D)[i]->as<LocalRegions::QuadExp>()))
2337  {
2338  for (int j = 0; j < locQuadExp->GetNtraces(); ++j)
2339  {
2340  SegGeom = (locQuadExp->GetGeom2D())->GetEdge(j);
2341  id = SegGeom->GetGlobalID();
2342 
2343  if (V1tmp[id] == 10000)
2344  {
2345  V1tmp[id] = SegGeom->GetVertex(0)->GetVid();
2346  V2tmp[id] = SegGeom->GetVertex(1)->GetVid();
2347  cnt++;
2348  }
2349  }
2350  }
2351  // in the future the tri edges may be not necessary (if the nedges is
2352  // known)
2353 
2354  else if ((locTriExp = (*exp2D)[i]->as<LocalRegions::TriExp>()))
2355  {
2356  for (int j = 0; j < locTriExp->GetNtraces(); ++j)
2357  {
2358  SegGeom = (locTriExp->GetGeom2D())->GetEdge(j);
2359  id = SegGeom->GetGlobalID();
2360 
2361  if (V1tmp[id] == 10000)
2362  {
2363  V1tmp[id] = SegGeom->GetVertex(0)->GetVid();
2364  V2tmp[id] = SegGeom->GetVertex(1)->GetVid();
2365  cnt++;
2366  }
2367  }
2368  }
2369  }
2370 
2371  V1 = Array<OneD, int>(cnt);
2372  V2 = Array<OneD, int>(cnt);
2373  // populate the output vectors V1,V2
2374  for (int g = 0; g < cnt; g++)
2375  {
2376  V1[g] = V1tmp[g];
2377  ASSERTL0(V1tmp[g] != 10000, "V1 wrong");
2378  V2[g] = V2tmp[g];
2379  ASSERTL0(V2tmp[g] != 10000, "V2 wrong");
2380  }
2381 }
2382 
2383 void MappingEVids(Array<OneD, NekDouble> xoldup, Array<OneD, NekDouble> yoldup,
2384  Array<OneD, NekDouble> xolddown,
2385  Array<OneD, NekDouble> yolddown, Array<OneD, NekDouble> xcold,
2386  Array<OneD, NekDouble> ycold, Array<OneD, int> Vids_c,
2387  SpatialDomains::MeshGraphSharedPtr mesh,
2388  MultiRegions::ExpListSharedPtr streak, Array<OneD, int> V1,
2389  Array<OneD, int> V2, int &nlays,
2390  Array<OneD, Array<OneD, int>> &Eids_lay,
2391  Array<OneD, Array<OneD, int>> &Vids_lay)
2392 {
2393  boost::ignore_unused(xoldup, xolddown);
2394 
2395  int nlay_Eids = xcold.size() - 1;
2396  int nlay_Vids = xcold.size();
2397 
2398  int nVertsTot = mesh->GetNvertices();
2399  cout << "nverttot=" << nVertsTot << endl;
2400  // find the first vert of each layer
2401  // int qp_closerup,qp_closerdown;
2402  // NekDouble diffup,diffdown;
2403  cout << "init nlays=" << nlays << endl;
2404  // tmp first vert array (of course <100!!)
2405  Array<OneD, int> tmpVids0(100);
2406  Array<OneD, NekDouble> tmpx0(100);
2407  Array<OneD, NekDouble> tmpy0(100);
2408  Array<OneD, NekDouble> x2D(nVertsTot);
2409  Array<OneD, NekDouble> y2D(nVertsTot);
2410  cout << "yoldup=" << yoldup[0] << endl;
2411  cout << "yolddown=" << yolddown[0] << endl;
2412 
2413  for (int r = 0; r < nVertsTot; r++)
2414  {
2415 
2416  SpatialDomains::PointGeomSharedPtr vertex = mesh->GetVertex(r);
2417  NekDouble x, y, z;
2418  vertex->GetCoords(x, y, z);
2419  x2D[r] = x;
2420  y2D[r] = y;
2421  if (xcold[0] == x)
2422  {
2423  // check if the vert is inside layer zone
2424  if (y <= yoldup[0] && y >= yolddown[0] && y != ycold[0])
2425  {
2426  // cout<<"x="<<x<<" y="<<y<<endl;
2427  tmpVids0[nlays] = r;
2428  tmpx0[nlays] = x;
2429  tmpy0[nlays] = y;
2430  nlays++;
2431  }
2432  }
2433  }
2434  cout << "nlays=" << nlays << endl;
2435 
2436  // reorder first verts from xlower to xhigher
2437  Array<OneD, NekDouble> tmpx(nlays);
2438  Array<OneD, NekDouble> tmpf(nlays);
2439  Array<OneD, int> tmpV(nlays);
2440  // local copy to prevent overwriting
2441  Vmath::Vcopy(nlays, tmpx0, 1, tmpx, 1);
2442  Vmath::Vcopy(nlays, tmpy0, 1, tmpf, 1);
2443  Vmath::Vcopy(nlays, tmpVids0, 1, tmpV, 1);
2444  NekDouble max = Vmath::Vmax(tmpf.size(), tmpf, 1);
2445  int index = 0;
2446  for (int w = 0; w < nlays; w++)
2447  {
2448  index = Vmath::Imin(tmpf.size(), tmpf, 1);
2449  tmpx0[w] = tmpx[index];
2450  tmpy0[w] = tmpf[index];
2451  tmpVids0[w] = tmpV[index];
2452  tmpf[index] = max + 1000;
2453  }
2454 
2455  // initialise layers arrays
2456  Eids_lay = Array<OneD, Array<OneD, int>>(nlays);
2457  Vids_lay = Array<OneD, Array<OneD, int>>(nlays);
2458  for (int m = 0; m < nlays; m++)
2459  {
2460  Eids_lay[m] = Array<OneD, int>(nlay_Eids);
2461  Vids_lay[m] = Array<OneD, int>(nlay_Vids);
2462  }
2463 
2464  // determine the others verts and edge for each layer
2465  NekDouble normbef = 0.0;
2466  NekDouble normtmp = 0.0;
2467  NekDouble xbef = 0.0;
2468  NekDouble ybef = 0.0;
2469  NekDouble xtmp, ytmp, normnext = 0.0, xnext = 0.0, ynext = 0.0, diff;
2470  NekDouble Ubef = 0.0, Utmp = 0.0, Unext = 0.0;
2471  Array<OneD, NekDouble> coord(2);
2472  int elmtid, offset;
2473  int nTotEdges = V1.size();
2474  Array<OneD, int> edgestmp(6);
2475  for (int m = 0; m < nlays; m++)
2476  {
2477  for (int g = 0; g < nlay_Eids; g++)
2478  {
2479  if (g == 0)
2480  {
2481  for (int h = 0; h < nTotEdges; h++)
2482  {
2483 
2484  if (tmpVids0[m] == V1[h])
2485  {
2486  SpatialDomains::PointGeomSharedPtr vertex =
2487  mesh->GetVertex(V2[h]);
2488  NekDouble x, y, z;
2489  vertex->GetCoords(x, y, z);
2490  if (x != tmpx0[m])
2491  {
2492  Eids_lay[m][0] = h;
2493  Vids_lay[m][0] = V1[h];
2494  Vids_lay[m][1] = V2[h];
2495  SpatialDomains::PointGeomSharedPtr vertex1 =
2496  mesh->GetVertex(V1[h]);
2497  NekDouble x1, y1, z1;
2498  vertex1->GetCoords(x1, y1, z1);
2499  normbef =
2500  sqrt((y - y1) * (y - y1) + (x - x1) * (x - x1));
2501  ybef = (y - y1);
2502  xbef = (x - x1);
2503  coord[0] = x;
2504  coord[1] = y;
2505  elmtid = streak->GetExpIndex(coord, 0.00001);
2506  offset = streak->GetPhys_Offset(elmtid);
2507  Ubef = streak->GetExp(elmtid)->PhysEvaluate(
2508  coord, streak->GetPhys() + offset);
2509  }
2510  }
2511  if (tmpVids0[m] == V2[h])
2512  {
2513  SpatialDomains::PointGeomSharedPtr vertex =
2514  mesh->GetVertex(V1[h]);
2515  NekDouble x, y, z;
2516  vertex->GetCoords(x, y, z);
2517  if (x != tmpx0[m])
2518  {
2519  Eids_lay[m][0] = h;
2520  Vids_lay[m][0] = V2[h];
2521  Vids_lay[m][1] = V1[h];
2522  SpatialDomains::PointGeomSharedPtr vertex2 =
2523  mesh->GetVertex(V2[h]);
2524  NekDouble x2 = 0.0, y2 = 0.0;
2525  normbef =
2526  sqrt((y - y2) * (y - y2) + (x - x2) * (x - x2));
2527  ybef = (y - y2);
2528  xbef = (x - x2);
2529  coord[0] = x;
2530  coord[1] = y;
2531  elmtid = streak->GetExpIndex(coord, 0.00001);
2532  offset = streak->GetPhys_Offset(elmtid);
2533  Ubef = streak->GetExp(elmtid)->PhysEvaluate(
2534  coord, streak->GetPhys() + offset);
2535  }
2536  }
2537  }
2538 
2539  cout << "Eid=" << Eids_lay[m][0]
2540  << " Vids_lay0=" << Vids_lay[m][0]
2541  << " Vidslay1=" << Vids_lay[m][1] << endl;
2542  }
2543  else
2544  {
2545  // three verts(edges) candidates
2546  int cnt = 0;
2547  for (int h = 0; h < nTotEdges; h++)
2548  {
2549  // cout<<Vids_lay[m][g]<<" V1="<<V1[h]<<"
2550  // V2[h]="<<V2[h]<<endl;
2551  if ((Vids_lay[m][g] == V1[h] || Vids_lay[m][g] == V2[h]) &&
2552  h != Eids_lay[m][g - 1])
2553  {
2554  cout << "edgetmp=" << h << endl;
2555  ASSERTL0(cnt <= 6, "wrong number of candidates");
2556  edgestmp[cnt] = h;
2557  cnt++;
2558  }
2559  }
2560 
2561  diff = 1000;
2562  Array<OneD, NekDouble> diffarray(cnt);
2563  Array<OneD, NekDouble> diffUarray(cnt);
2564  cout << "normbef=" << normbef << endl;
2565  cout << "Ubefcc=" << Ubef << endl;
2566  // choose the right candidate
2567  for (int e = 0; e < cnt; e++)
2568  {
2569  SpatialDomains::PointGeomSharedPtr vertex1 =
2570  mesh->GetVertex(V1[edgestmp[e]]);
2571  NekDouble x1, y1, z1;
2572  vertex1->GetCoords(x1, y1, z1);
2573  SpatialDomains::PointGeomSharedPtr vertex2 =
2574  mesh->GetVertex(V2[edgestmp[e]]);
2575  NekDouble x2, y2, z2;
2576  vertex2->GetCoords(x2, y2, z2);
2577 
2578  normtmp =
2579  sqrt((y2 - y1) * (y2 - y1) + (x2 - x1) * (x2 - x1));
2580 
2581  cout << "edgetmp1=" << edgestmp[e] << endl;
2582  cout << "V1 x1=" << x1 << " y1=" << y1 << endl;
2583  cout << "V2 x2=" << x2 << " y2=" << y2 << endl;
2584  if (Vids_lay[m][g] == V1[edgestmp[e]])
2585  {
2586 
2587  ytmp = (y2 - y1);
2588  xtmp = (x2 - x1);
2589  coord[0] = x2;
2590  coord[1] = y2;
2591  elmtid = streak->GetExpIndex(coord, 0.00001);
2592  offset = streak->GetPhys_Offset(elmtid);
2593  Utmp = streak->GetExp(elmtid)->PhysEvaluate(
2594  coord, streak->GetPhys() + offset);
2595  diffarray[e] = abs((xtmp * xbef + ytmp * ybef) /
2596  (normtmp * normbef) -
2597  1);
2598  diffUarray[e] = abs(Ubef - Utmp);
2599  cout << " normtmp=" << normtmp << endl;
2600  cout << " Utmpcc=" << Utmp << endl;
2601  cout << xtmp << " ytmp=" << ytmp << " diff="
2602  << abs(((xtmp * xbef + ytmp * ybef) /
2603  (normtmp * normbef)) -
2604  1)
2605  << endl;
2606  if (abs((xtmp * xbef + ytmp * ybef) /
2607  (normtmp * normbef) -
2608  1) < diff &&
2609  y2 <= yoldup[g + 1] && y2 >= yolddown[g + 1] &&
2610  y1 <= yoldup[g] && y1 >= yolddown[g])
2611  {
2612 
2613  Eids_lay[m][g] = edgestmp[e];
2614  Vids_lay[m][g + 1] = V2[edgestmp[e]];
2615  diff = abs((xtmp * xbef + ytmp * ybef) /
2616  (normtmp * normbef) -
2617  1);
2618  normnext = normtmp;
2619  ynext = ytmp;
2620  xnext = xtmp;
2621  Unext = Utmp;
2622  }
2623  }
2624  else if (Vids_lay[m][g] == V2[edgestmp[e]])
2625  {
2626 
2627  ytmp = (y1 - y2);
2628  xtmp = (x1 - x2);
2629  coord[0] = x1;
2630  coord[1] = y1;
2631  elmtid = streak->GetExpIndex(coord, 0.00001);
2632  offset = streak->GetPhys_Offset(elmtid);
2633  Utmp = streak->GetExp(elmtid)->PhysEvaluate(
2634  coord, streak->GetPhys() + offset);
2635  diffarray[e] = abs((xtmp * xbef + ytmp * ybef) /
2636  (normtmp * normbef) -
2637  1);
2638  diffUarray[e] = abs(Ubef - Utmp);
2639  cout << " normtmp=" << normtmp << endl;
2640  cout << " Utmpcc=" << Utmp << endl;
2641  cout << xtmp << " ytmp=" << ytmp << " diff="
2642  << abs(((xtmp * xbef + ytmp * ybef) /
2643  (normtmp * normbef)) -
2644  1)
2645  << endl;
2646  if (abs((xtmp * xbef + ytmp * ybef) /
2647  (normtmp * normbef) -
2648  1) < diff &&
2649  y2 <= yoldup[g] && y2 >= yolddown[g] &&
2650  y1 <= yoldup[g + 1] && y1 >= yolddown[g + 1])
2651  {
2652  Eids_lay[m][g] = edgestmp[e];
2653  Vids_lay[m][g + 1] = V1[edgestmp[e]];
2654  diff = abs((xtmp * xbef + ytmp * ybef) /
2655  (normtmp * normbef) -
2656  1);
2657  normnext = normtmp;
2658  ynext = ytmp;
2659  xnext = xtmp;
2660  Unext = Utmp;
2661  }
2662  }
2663  else
2664  {
2665  ASSERTL0(false, "eid not found");
2666  }
2667  }
2668  cout << "Eid before check=" << Eids_lay[m][g] << endl;
2669  for (int q = 0; q < cnt; q++)
2670  {
2671  cout << q << " diff" << diffarray[q] << endl;
2672  }
2673  // check if the eid has a vert in common with another layer
2674  bool check = false;
2675  if (m > 0 && m < nlays)
2676  {
2677  check = checkcommonvert(Vids_lay[m - 1], Vids_c,
2678  Vids_lay[m][g + 1]);
2679  }
2680  if (check == true)
2681  {
2682  cout << "COMMON VERT" << endl;
2683  int eid = Vmath::Imin(cnt, diffarray, 1);
2684  diffarray[eid] = 1000;
2685  eid = Vmath::Imin(cnt, diffarray, 1);
2686 
2687  SpatialDomains::PointGeomSharedPtr vertex1 =
2688  mesh->GetVertex(V1[edgestmp[eid]]);
2689  NekDouble x1, y1, z1;
2690  vertex1->GetCoords(x1, y1, z1);
2691  SpatialDomains::PointGeomSharedPtr vertex2 =
2692  mesh->GetVertex(V2[edgestmp[eid]]);
2693  NekDouble x2, y2, z2;
2694  vertex2->GetCoords(x2, y2, z2);
2695 
2696  normtmp =
2697  sqrt((y2 - y1) * (y2 - y1) + (x2 - x1) * (x2 - x1));
2698 
2699  Eids_lay[m][g] = edgestmp[eid];
2700  if (Vids_lay[m][g] == V1[edgestmp[eid]])
2701  {
2702  coord[0] = x2;
2703  coord[1] = y2;
2704  elmtid = streak->GetExpIndex(coord, 0.00001);
2705  offset = streak->GetPhys_Offset(elmtid);
2706  Utmp = streak->GetExp(elmtid)->PhysEvaluate(
2707  coord, streak->GetPhys() + offset);
2708  Vids_lay[m][g + 1] = V2[edgestmp[eid]];
2709  normnext = normtmp;
2710  ynext = (y2 - y1);
2711  xnext = (x2 - x1);
2712  ;
2713  Unext = Utmp;
2714  }
2715  if (Vids_lay[m][g] == V2[edgestmp[eid]])
2716  {
2717  coord[0] = x1;
2718  coord[1] = y1;
2719  elmtid = streak->GetExpIndex(coord, 0.00001);
2720  offset = streak->GetPhys_Offset(elmtid);
2721  Utmp = streak->GetExp(elmtid)->PhysEvaluate(
2722  coord, streak->GetPhys() + offset);
2723  Vids_lay[m][g + 1] = V1[edgestmp[eid]];
2724  normnext = normtmp;
2725  ynext = (y1 - y2);
2726  xnext = (x1 - x2);
2727  ;
2728  Unext = Utmp;
2729  }
2730  }
2731 
2732  cout << m << "edge aft:" << Eids_lay[m][g]
2733  << " Vid=" << Vids_lay[m][g + 1] << endl;
2734  normbef = normnext;
2735  ybef = ynext;
2736  xbef = xnext;
2737  Ubef = Unext;
2738 
2739  cout << "endelse" << normtmp << endl;
2740  } // end else
2741 
2742  } // end g it
2743 
2744  } // end m it
2745 
2746  for (int w = 0; w < nlays; w++)
2747  {
2748  for (int f = 0; f < nlay_Eids; f++)
2749  {
2750  cout << "check=" << w << " Eid:" << Eids_lay[w][f] << endl;
2751  }
2752  }
2753 }
2754 
2755 bool checkcommonvert(Array<OneD, int> Vids_laybefore, Array<OneD, int> Vids_c,
2756  int Vid)
2757 {
2758  bool check = false;
2759  for (int u = 0; u < Vids_laybefore.size(); u++)
2760  {
2761  if (Vids_laybefore[u] == Vid || Vids_c[u] == Vid)
2762  {
2763  check = true;
2764  }
2765  cout << Vid << " Vert test=" << Vids_laybefore[u] << endl;
2766  }
2767  return check;
2768 }
2769 
2770 void Cutrepetitions(int nedges, Array<OneD, NekDouble> inarray,
2771  Array<OneD, NekDouble> &outarray)
2772 {
2773 
2774  // determine npedge:
2775  int np_lay = inarray.size();
2776  ASSERTL0(inarray.size() % nedges == 0, " something on number npedge");
2777  // cut out the repetitions:
2778 
2779  int cnt = 0;
2780  for (int w = 0; w < np_lay; w++)
2781  {
2782 
2783  if (w < np_lay - 1)
2784  {
2785  if (inarray[w] == inarray[w + 1])
2786  {
2787  continue;
2788  }
2789  }
2790  outarray[cnt] = inarray[w];
2791  cnt++;
2792  }
2793 
2794  ASSERTL0(cnt == np_lay - (nedges - 1), "wrong cut");
2795 }
2796 
2797 int DetermineclosePointxindex(NekDouble x, Array<OneD, NekDouble> xArray)
2798 {
2799  int npts = xArray.size();
2800  Array<OneD, NekDouble> xcopy(npts);
2801  Vmath::Vcopy(npts, xArray, 1, xcopy, 1);
2802  // subtract xpoint and abs
2803 
2804  Vmath::Sadd(npts, -x, xcopy, 1, xcopy, 1);
2805  Vmath::Vabs(npts, xcopy, 1, xcopy, 1);
2806 
2807  int index = Vmath::Imin(npts, xcopy, 1);
2808  if (xArray[index] > x) // assume always x[index]< x(HYPHOTHESIS)
2809  {
2810  index = index - 1;
2811  }
2812  return index;
2813 }
2814 
2815 void GenerateNeighbourArrays(int index, int neighpoints,
2816  Array<OneD, NekDouble> xArray,
2817  Array<OneD, NekDouble> yArray,
2818  Array<OneD, NekDouble> &Neighbour_x,
2819  Array<OneD, NekDouble> &Neighbour_y)
2820 {
2821  ASSERTL0(neighpoints % 2 == 0, "number of neighbour points should be even");
2822  int leftpoints = (neighpoints / 2) - 1; //-1 because xArray[index]< x
2823  int rightpoints = neighpoints / 2;
2824  int diff;
2825  int start;
2826  // cout<<"index="<<index<<" left="<<leftpoints<<"
2827  // right="<<rightpoints<<endl;
2828  if (index - leftpoints < 0)
2829  {
2830  // cout<"case0"<<endl;
2831  diff = index - leftpoints;
2832  start = 0;
2833  Vmath::Vcopy(neighpoints, &yArray[0], 1, &Neighbour_y[0], 1);
2834  Vmath::Vcopy(neighpoints, &xArray[0], 1, &Neighbour_x[0], 1);
2835  }
2836  else if ((yArray.size() - 1) - index < rightpoints)
2837  {
2838  // cout"case1 closest="<<xArray[index]<<endl;
2839  int rpoints = (yArray.size() - 1) - index; //
2840  diff = rightpoints - rpoints;
2841  // cout<"start index="<<index-leftpoints-diff<<endl;
2842  start = index - leftpoints - diff;
2843  Vmath::Vcopy(neighpoints, &yArray[start], 1, &Neighbour_y[0], 1);
2844  Vmath::Vcopy(neighpoints, &xArray[start], 1, &Neighbour_x[0], 1);
2845  }
2846  else
2847  {
2848  // cout<<"caseaa"<<endl;
2849  start = index - leftpoints;
2850  Vmath::Vcopy(neighpoints, &yArray[start], 1, &Neighbour_y[0], 1);
2851  Vmath::Vcopy(neighpoints, &xArray[start], 1, &Neighbour_x[0], 1);
2852  }
2853  /*
2854  for(int t= start; t<start+neighpoints; t++)
2855  {
2856  cout<<"Px="<<xArray[t]<<" "<<yArray[t]<<endl;
2857  }
2858  */
2859  // check if there is any repetition
2860  for (int f = 1; f < neighpoints; f++)
2861  {
2862  ASSERTL0(Neighbour_x[f] != Neighbour_x[f - 1],
2863  " repetition on NeighbourArrays");
2864  }
2865 }
2866 
2867 NekDouble LagrangeInterpolant(NekDouble x, int npts,
2868  Array<OneD, NekDouble> xpts,
2869  Array<OneD, NekDouble> funcvals)
2870 {
2871  NekDouble sum = 0.0;
2872  NekDouble LagrangePoly;
2873  // cout<<"lagrange"<<endl;
2874  for (int pt = 0; pt < npts; ++pt)
2875  {
2876  NekDouble h = 1.0;
2877 
2878  for (int j = 0; j < pt; ++j)
2879  {
2880  h = h * (x - xpts[j]) / (xpts[pt] - xpts[j]);
2881  }
2882 
2883  for (int k = pt + 1; k < npts; ++k)
2884  {
2885  h = h * (x - xpts[k]) / (xpts[pt] - xpts[k]);
2886  }
2887  LagrangePoly = h;
2888 
2889  sum += funcvals[pt] * LagrangePoly;
2890  }
2891  // cout<<"result :"<<sum<<endl;
2892  return sum;
2893 }
2894 
2895 void EvaluateTangent(int npoints, Array<OneD, NekDouble> xcQedge,
2896  Array<OneD, NekDouble> coeffsinterp,
2897  Array<OneD, NekDouble> &txQedge,
2898  Array<OneD, NekDouble> &tyQedge)
2899 {
2900  Array<OneD, NekDouble> yprime(npoints, 0.0);
2901  int np_pol = coeffsinterp.size();
2902  cout << "evaluatetan with " << np_pol << endl;
2903 
2904  // calc derivative poly (cut last entry on b array that is the const)
2905  int derorder;
2906  Array<OneD, NekDouble> yinterp(npoints, 0.0);
2907  int polorder;
2908  for (int q = 0; q < npoints; q++)
2909  {
2910  polorder = np_pol - 1;
2911  derorder = np_pol - 2;
2912  yprime[q] = 0;
2913  for (int d = 0; d < np_pol - 1; d++)
2914  {
2915  yprime[q] += (derorder + 1) * coeffsinterp[d] *
2916  std::pow(xcQedge[q], derorder);
2917  derorder--;
2918  }
2919  // test
2920  for (int a = 0; a < np_pol; a++)
2921  {
2922  yinterp[q] += coeffsinterp[a] * std::pow(xcQedge[q], polorder);
2923  // cout<<"coeff*x^n="<<b[a]*std::pow(xcQedge[q],polorder)<<"
2924  // sum="<<yinterp[q]<<endl;
2925  polorder--;
2926  }
2927  }
2928 
2929  // transf yprime into tx,ty:
2930  for (int n = 0; n < npoints; n++)
2931  {
2932  // ABS???!!
2933  txQedge[n] = 0;
2934  txQedge[n] = cos((atan((yprime[n]))));
2935  tyQedge[n] = sin((atan((yprime[n]))));
2936  cout << xcQedge[n] << " " << yinterp[n] << " " << yprime[n]
2937  << " " << txQedge[n] << " " << tyQedge[n] << endl;
2938  }
2939 }
2940 
2941 void PolyInterp(Array<OneD, NekDouble> xpol, Array<OneD, NekDouble> ypol,
2942  Array<OneD, NekDouble> &coeffsinterp,
2943  Array<OneD, NekDouble> &xcout, Array<OneD, NekDouble> &ycout,
2944  int edge, int npedge)
2945 {
2946  int np_pol = xpol.size();
2947  int N = np_pol;
2948  Array<OneD, NekDouble> A(N * N, 1.0);
2949  Array<OneD, NekDouble> b(N);
2950  int row = 0;
2951  // fill column by column
2952  for (int e = 0; e < N; e++)
2953  {
2954  row = 0;
2955  for (int w = 0; w < N; w++)
2956  {
2957  A[N * e + row] = std::pow(xpol[w], N - 1 - e);
2958  row++;
2959  }
2960  }
2961  row = 0;
2962  for (int r = 0; r < np_pol; r++)
2963  {
2964  b[row] = ypol[r];
2965  // cout<<"b="<<b[row]<<" y_c="<<y_c[r]<<endl;
2966  row++;
2967  }
2968 
2969  // cout<<"A elements="<<A.size()<<endl;
2970  Array<OneD, int> ipivot(N);
2971  int info = 0;
2972  // Lapack::Dgesv( N, 1, A.get(), N, ipivot.get(), b.get(), N, info);
2973  Lapack::Dgetrf(N, N, A.get(), N, ipivot.get(), info);
2974  if (info < 0)
2975  {
2976  std::string message =
2977  "ERROR: The " + boost::lexical_cast<std::string>(-info) +
2978  "th parameter had an illegal parameter for dgetrf";
2979  ASSERTL0(false, message.c_str());
2980  }
2981  else if (info > 0)
2982  {
2983  std::string message =
2984  "ERROR: Element u_" + boost::lexical_cast<std::string>(info) +
2985  boost::lexical_cast<std::string>(info) + " is 0 from dgetrf";
2986  ASSERTL0(false, message.c_str());
2987  }
2988 
2989  // N means no transponse (direct matrix)
2990  int ncolumns_b = 1;
2991  Lapack::Dgetrs('N', N, ncolumns_b, A.get(), N, ipivot.get(), b.get(), N,
2992  info);
2993  if (info < 0)
2994  {
2995  std::string message =
2996  "ERROR: The " + boost::lexical_cast<std::string>(-info) +
2997  "th parameter had an illegal parameter for dgetrf";
2998  ASSERTL0(false, message.c_str());
2999  }
3000  else if (info > 0)
3001  {
3002  std::string message =
3003  "ERROR: Element u_" + boost::lexical_cast<std::string>(info) +
3004  boost::lexical_cast<std::string>(info) + " is 0 from dgetrf";
3005  ASSERTL0(false, message.c_str());
3006  }
3007  /*
3008  for(int h=0; h<np_pol; h++)
3009  {
3010  cout<<"coeff:"<<b[h]<<endl;
3011  }
3012  */
3013 
3014  // ovewrite the ycPhysMOD:
3015  int polorder;
3016  for (int c = 0; c < npedge; c++)
3017  {
3018  polorder = np_pol - 1;
3019  // put ycPhysMOD to zero
3020  ycout[edge * (npedge) + c + 1] = 0;
3021  for (int d = 0; d < np_pol; d++)
3022  {
3023  ycout[edge * (npedge) + c + 1] +=
3024  b[d] * std::pow(xcout[edge * (npedge) + c + 1], polorder);
3025  polorder--;
3026  }
3027  }
3028 
3029  // write coeffs
3030  Vmath::Vcopy(np_pol, b, 1, coeffsinterp, 1);
3031 }
3032 
3033 void PolyFit(int polyorder, int npoints, Array<OneD, NekDouble> xin,
3034  Array<OneD, NekDouble> fin, Array<OneD, NekDouble> &coeffsinterp,
3035  Array<OneD, NekDouble> &xout, Array<OneD, NekDouble> &fout,
3036  int npout)
3037 {
3038 
3039  int N = polyorder + 1;
3040  Array<OneD, NekDouble> A(N * N, 0.0);
3041  Array<OneD, NekDouble> b(N, 0.0);
3042  cout << npoints << endl;
3043  for (int u = 0; u < npoints; u++)
3044  {
3045  cout << "c=" << xin[u] << " " << fin[u] << endl;
3046  }
3047  // fill column by column
3048  // e counts cols
3049  for (int e = 0; e < N; e++)
3050  {
3051 
3052  for (int row = 0; row < N; row++)
3053  {
3054  for (int w = 0; w < npoints; w++)
3055  {
3056  A[N * e + row] += std::pow(xin[w], e + row);
3057  }
3058  }
3059  }
3060 
3061  for (int row = 0; row < N; row++)
3062  {
3063  for (int h = 0; h < npoints; h++)
3064  {
3065  b[row] += fin[h] * std::pow(xin[h], row);
3066  // cout<<b="<<b[row]<<" y_c="<<y_c[r]<<endl;
3067  }
3068  }
3069 
3070  // cout<<"A elements="<<A.size()<<endl;
3071  Array<OneD, int> ipivot(N);
3072  int info = 0;
3073  // Lapack::Dgesv( N, 1, A.get(), N, ipivot.get(), b.get(), N, info);
3074  Lapack::Dgetrf(N, N, A.get(), N, ipivot.get(), info);
3075 
3076  if (info < 0)
3077  {
3078  std::string message =
3079  "ERROR: The " + boost::lexical_cast<std::string>(-info) +
3080  "th parameter had an illegal parameter for dgetrf";
3081  ASSERTL0(false, message.c_str());
3082  }
3083  else if (info > 0)
3084  {
3085  std::string message =
3086  "ERROR: Element u_" + boost::lexical_cast<std::string>(info) +
3087  boost::lexical_cast<std::string>(info) + " is 0 from dgetrf";
3088  ASSERTL0(false, message.c_str());
3089  }
3090  // N means no transponse (direct matrix)
3091  int ncolumns_b = 1;
3092  Lapack::Dgetrs('N', N, ncolumns_b, A.get(), N, ipivot.get(), b.get(), N,
3093  info);
3094  if (info < 0)
3095  {
3096  std::string message =
3097  "ERROR: The " + boost::lexical_cast<std::string>(-info) +
3098  "th parameter had an illegal parameter for dgetrf";
3099  ASSERTL0(false, message.c_str());
3100  }
3101  else if (info > 0)
3102  {
3103  std::string message =
3104  "ERROR: Element u_" + boost::lexical_cast<std::string>(info) +
3105  boost::lexical_cast<std::string>(info) + " is 0 from dgetrf";
3106  ASSERTL0(false, message.c_str());
3107  }
3108 
3109  // the lower coeff the lower is the grad
3110  // reverse:
3111  Array<OneD, NekDouble> tmp(N);
3112  Vmath::Vcopy(N, b, 1, tmp, 1);
3113  int cnt = N;
3114  for (int j = 0; j < N; j++)
3115  {
3116  b[j] = tmp[cnt - 1];
3117  cnt--;
3118  }
3119 
3120  for (int h = 0; h < N; h++)
3121  {
3122  cout << "coeff:" << b[h] << endl;
3123  }
3124 
3125  // ovewrite the ycPhysMOD:
3126  int polorder;
3127  for (int c = 0; c < npout; c++)
3128  {
3129  polorder = polyorder;
3130  // put ycPhysMOD to zero
3131  fout[c] = 0;
3132  for (int d = 0; d < N; d++)
3133  {
3134 
3135  fout[c] += b[d] * std::pow(xout[c], polorder);
3136  polorder--;
3137  }
3138  }
3139 
3140  // write coeffs
3141  Vmath::Vcopy(N, b, 1, coeffsinterp, 1);
3142 }
3143 
3144 void Orderfunctionx(Array<OneD, NekDouble> inarray_x,
3145  Array<OneD, NekDouble> inarray_y,
3146  Array<OneD, NekDouble> &outarray_x,
3147  Array<OneD, NekDouble> &outarray_y)
3148 {
3149  Array<OneD, NekDouble> tmpx(inarray_x.size());
3150  Array<OneD, NekDouble> tmpy(inarray_x.size());
3151  // local copy to prevent overwriting
3152  Vmath::Vcopy(inarray_x.size(), inarray_x, 1, tmpx, 1);
3153  Vmath::Vcopy(inarray_x.size(), inarray_y, 1, tmpy, 1);
3154 
3155  // order function with respect to x
3156  int index;
3157  NekDouble max = Vmath::Vmax(tmpx.size(), tmpx, 1);
3158  for (int w = 0; w < tmpx.size(); w++)
3159  {
3160  index = Vmath::Imin(tmpx.size(), tmpx, 1);
3161  outarray_x[w] = tmpx[index];
3162  outarray_y[w] = tmpy[index];
3163  if (w < tmpx.size() - 1) // case of repetitions
3164  {
3165  if (tmpx[index] == tmpx[index + 1])
3166  {
3167  outarray_x[w + 1] = tmpx[index + 1];
3168  outarray_y[w + 1] = tmpy[index + 1];
3169  tmpx[index + 1] = max + 1000;
3170  w++;
3171  }
3172  }
3173  /*
3174  if(w>0)//case of repetitions
3175  {
3176  if(inarray_x[index] == tmpx[index-1])
3177  {
3178  outarray_x[w+1]= tmpx[index-1];
3179  outarray_y[w+1]= tmpy[index-1];
3180  tmpx[index-1] = max+1000;
3181  w++;
3182  }
3183  }
3184  */
3185  tmpx[index] = max + 1000;
3186  }
3187 }
3188 
3189 void MoveLayersvertically(int nlays, int nvertl, int cntlow, int cntup,
3190  Array<OneD, Array<OneD, int>> lay_Vids,
3191  Array<OneD, NekDouble> xc, Array<OneD, NekDouble> yc,
3192  Array<OneD, int> Down, Array<OneD, int> Up,
3193  Array<OneD, NekDouble> &xnew,
3194  Array<OneD, NekDouble> &ynew,
3195  Array<OneD, Array<OneD, NekDouble>> &layers_x,
3196  Array<OneD, Array<OneD, NekDouble>> &layers_y)
3197 {
3198  int np_lay = layers_y[0].size();
3199  // 0<h<nlays-1 to fill only the 'internal' layers(no up,low);
3200  for (int h = 1; h < nlays - 1; h++)
3201  {
3202  layers_x[h] = Array<OneD, NekDouble>(np_lay);
3203  for (int s = 0; s < nvertl; s++)
3204  {
3205  // check if ynew is still empty
3206  ASSERTL0(ynew[lay_Vids[h][s]] == -20, "ynew layers not empty");
3207  if (h < cntlow + 1) // layers under the crit lay
3208  {
3209  // y= ylow+delta
3210  ynew[lay_Vids[h][s]] =
3211  ynew[Down[s]] +
3212  h * abs(ynew[Down[s]] - yc[s]) / (cntlow + 1);
3213  // put the layer vertical
3214  xnew[lay_Vids[h][s]] = xc[s];
3215  // cout<<"ynew="<<ynew[ lay_Vids[h][s] ]<<"
3216  // ydown="<<ynew[Down[s]]<< " delta="<<abs(ynew[Down[s]] -
3217  // y_c[s])/(cntlow+1)<<endl; until now layers_y=yold
3218  layers_y[h][s] = ynew[lay_Vids[h][s]];
3219  layers_x[h][s] = xnew[lay_Vids[h][s]];
3220  }
3221  else
3222  {
3223  // y = yc+delta
3224  ynew[lay_Vids[h][s]] = yc[s] + (h - cntlow) *
3225  abs(ynew[Up[s]] - yc[s]) /
3226  (cntup + 1);
3227  // put the layer vertical
3228  xnew[lay_Vids[h][s]] = xc[s];
3229  // until now layers_y=yold
3230  layers_y[h][s] = ynew[lay_Vids[h][s]];
3231  layers_x[h][s] = xnew[lay_Vids[h][s]];
3232  }
3233  }
3234  }
3235 }
3236 
3237 void MoveLayerNfixedxpos(int nvertl, int npedge, Array<OneD, NekDouble> xcPhys,
3238  Array<OneD, NekDouble> tmpx_lay,
3239  Array<OneD, NekDouble> tmpy_lay, Array<OneD, int> Vids,
3240  Array<OneD, NekDouble> &xlay,
3241  Array<OneD, NekDouble> &ylay,
3242  Array<OneD, NekDouble> &xnew,
3243  Array<OneD, NekDouble> &ynew)
3244 {
3245  int np_lay = xcPhys.size();
3246  int nedges = nvertl - 1;
3247  Array<OneD, NekDouble> tmpx(np_lay - (nedges - 1));
3248  Array<OneD, NekDouble> tmpy(np_lay - (nedges - 1));
3249  Cutrepetitions(nedges, tmpx_lay, tmpx);
3250  Cutrepetitions(nedges, tmpy_lay, tmpy);
3251  // order points in x:
3252  int index;
3253  int closepoints = 4;
3254  Array<OneD, NekDouble> Pxinterp(closepoints);
3255  Array<OneD, NekDouble> Pyinterp(closepoints);
3256  Orderfunctionx(tmpx, tmpy, tmpx, tmpy);
3257  // determine the neighbour points (-3;+3)
3258  for (int g = 0; g < nedges; g++)
3259  {
3260  // write vert coords
3261  // v1
3262  index = DetermineclosePointxindex(xcPhys[g * npedge + 0], tmpx);
3263  // generate neighbour arrays:
3264  GenerateNeighbourArrays(index, closepoints, tmpx, tmpy, Pxinterp,
3265  Pyinterp);
3266  ynew[Vids[g]] = LagrangeInterpolant(xcPhys[g * npedge + 0], closepoints,
3267  Pxinterp, Pyinterp);
3268  xnew[Vids[g]] = xcPhys[g * npedge + 0];
3269  ylay[g * npedge + 0] = ynew[Vids[g]];
3270  xlay[g * npedge + 0] = xnew[Vids[g]];
3271 
3272  // v2
3273  // determine closest index:
3274  index =
3275  DetermineclosePointxindex(xcPhys[g * npedge + npedge - 1], tmpx);
3276  // generate neighbour arrays:
3277  GenerateNeighbourArrays(index, closepoints, tmpx, tmpy, Pxinterp,
3278  Pyinterp);
3279  ynew[Vids[g + 1]] = LagrangeInterpolant(
3280  xcPhys[g * npedge + npedge - 1], closepoints, Pxinterp, Pyinterp);
3281  xnew[Vids[g + 1]] = xcPhys[g * npedge + npedge - 1];
3282  ylay[g * npedge + npedge - 1] = ynew[Vids[g + 1]];
3283  xlay[g * npedge + npedge - 1] = xnew[Vids[g + 1]];
3284 
3285  // middle points
3286  for (int r = 0; r < npedge - 2; r++)
3287  {
3288 
3289  // determine closest point index:
3290  index = DetermineclosePointxindex(xcPhys[g * npedge + r + 1], tmpx);
3291  // cout<<" vert+"<<index<<endl;
3292 
3293  ASSERTL0(index <= tmpy.size() - 1, " index wrong");
3294  // generate neighbour arrays Pyinterp,Pxinterp
3295  GenerateNeighbourArrays(index, closepoints, tmpx, tmpy, Pxinterp,
3296  Pyinterp);
3297 
3298  ylay[g * npedge + r + 1] = LagrangeInterpolant(
3299  xcPhys[g * npedge + r + 1], closepoints, Pxinterp, Pyinterp);
3300  // cout<<"x value="<<xcPhysMOD[g*npedge +r+1]<<endl;
3301  xlay[g * npedge + r + 1] = xcPhys[g * npedge + r + 1];
3302 
3303  /*
3304  for(int t=0; t<6; t++)
3305  {
3306  cout<<"Px="<<Pxinterp[t]<<" "<<Pyinterp[t]<<endl;
3307  }
3308  */
3309  }
3310 
3311  } // edge closed g
3312 }
3313 
3314 void MoveLayerNnormpos(int nvertl, int npedge, Array<OneD, NekDouble> xcPhys,
3315  Array<OneD, NekDouble> tmpx_lay,
3316  Array<OneD, NekDouble> tmpy_lay, Array<OneD, int> Vids,
3317  Array<OneD, NekDouble> &xlay,
3318  Array<OneD, NekDouble> &ylay,
3319  Array<OneD, NekDouble> &xnew,
3320  Array<OneD, NekDouble> &ynew)
3321 {
3322  int np_lay = xcPhys.size();
3323  int nedges = nvertl - 1;
3324  NekDouble x0, x1, xtmp;
3325  Array<OneD, NekDouble> tmpx(np_lay - (nedges - 1));
3326  Array<OneD, NekDouble> tmpy(np_lay - (nedges - 1));
3327  Cutrepetitions(nedges, tmpx_lay, tmpx);
3328  Cutrepetitions(nedges, tmpy_lay, tmpy);
3329  // order points in x:
3330  int index;
3331  int closepoints = 4;
3332  Array<OneD, NekDouble> Pxinterp(closepoints);
3333  Array<OneD, NekDouble> Pyinterp(closepoints);
3334  Orderfunctionx(tmpx, tmpy, tmpx, tmpy);
3335  // determine the neighbour points (-3;+3)
3336  for (int g = 0; g < nedges; g++)
3337  {
3338  // write vert coords
3339  // v1
3340  ynew[Vids[g]] = tmpy_lay[g * npedge + 0];
3341  xnew[Vids[g]] = tmpx_lay[g * npedge + 0];
3342 
3343  ylay[g * npedge + 0] = ynew[Vids[g]];
3344  xlay[g * npedge + 0] = xnew[Vids[g]];
3345 
3346  // v2
3347  ynew[Vids[g + 1]] = tmpy_lay[g * npedge + npedge - 1];
3348  xnew[Vids[g + 1]] = tmpx_lay[g * npedge + npedge - 1];
3349  ylay[g * npedge + npedge - 1] = ynew[Vids[g + 1]];
3350  xlay[g * npedge + npedge - 1] = xnew[Vids[g + 1]];
3351 
3352  // middle points
3353  for (int r = 0; r < npedge - 2; r++)
3354  {
3355  x0 = xlay[g * npedge + 0];
3356  x1 = xlay[g * npedge + npedge - 1];
3357  xtmp = x0 + r * (x1 - x0) / (npedge - 1);
3358  // determine closest point index:
3359  index = DetermineclosePointxindex(xtmp, tmpx);
3360  // cout<<" vert+"<<index<<endl;
3361 
3362  ASSERTL0(index <= tmpy.size() - 1, " index wrong");
3363  // generate neighbour arrays Pyinterp,Pxinterp
3364  GenerateNeighbourArrays(index, closepoints, tmpx, tmpy, Pxinterp,
3365  Pyinterp);
3366 
3367  ylay[g * npedge + r + 1] =
3368  LagrangeInterpolant(xtmp, closepoints, Pxinterp, Pyinterp);
3369  // cout<<"x value="<<xcPhysMOD[g*npedge +r+1]<<endl;
3370  xlay[g * npedge + r + 1] = xtmp;
3371  }
3372 
3373  } // edge closed g
3374 }
3375 
3376 void MoveOutsidePointsfixedxpos(
3377  int npedge, SpatialDomains::MeshGraphSharedPtr mesh,
3378  Array<OneD, NekDouble> xcold, Array<OneD, NekDouble> ycold,
3379  Array<OneD, NekDouble> xolddown, Array<OneD, NekDouble> yolddown,
3380  Array<OneD, NekDouble> xoldup, Array<OneD, NekDouble> yoldup,
3381  Array<OneD, NekDouble> ylaydown, Array<OneD, NekDouble> ylayup,
3382  Array<OneD, NekDouble> &xnew, Array<OneD, NekDouble> &ynew)
3383 {
3384  boost::ignore_unused(xolddown, xoldup);
3385 
3386  // update vertices coords outside layers region
3387  int nvertl = ycold.size();
3388  int nVertTot = mesh->GetNvertices();
3389  for (int n = 0; n < nVertTot; n++)
3390  {
3391  NekDouble ratio;
3392  SpatialDomains::PointGeomSharedPtr vertex = mesh->GetVertex(n);
3393  NekDouble x, y, z;
3394  vertex->GetCoords(x, y, z);
3395  int qp_closer = 0;
3396  // determine the closer xold_up
3397  NekDouble tmp = 1000;
3398 
3399  for (int k = 0; k < nvertl; k++)
3400  {
3401  if (abs(x - xcold[k]) < tmp)
3402  {
3403  tmp = abs(x - xcold[k]);
3404  qp_closer = k;
3405  }
3406  }
3407  // find nplay_closer
3408  int nplay_closer;
3409  if (qp_closer == 0)
3410  {
3411  nplay_closer = 0; // first vert
3412  }
3413  else
3414  {
3415  nplay_closer = (qp_closer - 1) * npedge + npedge - 1;
3416  }
3417 
3418  if (y > yoldup[qp_closer] && y < 1) // nlays-1 is layerup
3419  {
3420 
3421  // ratio =
3422  // (1-layers_y[nlays-1][qp_closer])*(1-y_c[qp_closer])/
3423  // ( (1-yold_up[n])*(1-yold_c[qp_closer]) );
3424  ratio = (1 - ylayup[nplay_closer]) / ((1 - yoldup[qp_closer]));
3425  // distance prop to layerup
3426  ynew[n] = ylayup[nplay_closer] + (y - yoldup[qp_closer]) * ratio;
3427  xnew[n] = x;
3428  }
3429  else if (y < yolddown[qp_closer] && y > -1) // 0 is layerdown
3430  {
3431 
3432  ratio = (1 + ylaydown[nplay_closer]) / ((1 + yolddown[qp_closer]));
3433  // distance prop to layerlow
3434  ynew[n] =
3435  ylaydown[nplay_closer] + (y - yolddown[qp_closer]) * ratio;
3436  xnew[n] = x;
3437  }
3438  }
3439 }
3440 
3441 void MoveOutsidePointsNnormpos(
3442  int npedge, SpatialDomains::MeshGraphSharedPtr mesh,
3443  Array<OneD, NekDouble> xcold, Array<OneD, NekDouble> ycold,
3444  Array<OneD, NekDouble> xolddown, Array<OneD, NekDouble> yolddown,
3445  Array<OneD, NekDouble> xoldup, Array<OneD, NekDouble> yoldup,
3446  Array<OneD, NekDouble> xlaydown, Array<OneD, NekDouble> ylaydown,
3447  Array<OneD, NekDouble> xlayup, Array<OneD, NekDouble> ylayup,
3448  Array<OneD, NekDouble> nxPhys, Array<OneD, NekDouble> nyPhys,
3449  Array<OneD, NekDouble> &xnew, Array<OneD, NekDouble> &ynew)
3450 {
3451  boost::ignore_unused(xcold, ycold);
3452  /*
3453  int nq1D =bndfieldup->GetTotPoints();
3454  Array<OneD, NekDouble> xlayoldup(nq1D);
3455  Array<OneD, NekDouble> xlayolddown(nq1D);
3456  Array<OneD, NekDouble> ylayoldup(nq1D);
3457  Array<OneD, NekDouble> ylayolddown(nq1D);
3458  Array<OneD, NekDouble> zlayoldup(nq1D);
3459  Array<OneD, NekDouble> zlayolddown(nq1D);
3460  bndfielddown->GetCoords( xlayolddown, ylayolddown,zlayolddown);
3461  bndfieldup->GetCoords( xlayoldup, ylayoldup,zlayoldup);
3462 
3463  NekDouble xmax = Vmath::Vmax(nq1D, xlayoldup,1);
3464  NekDouble xmin = Vmath::Vmin(nq1D, xlayoldup,1);
3465  */
3466  // determine the new verts up/down pos:
3467  int nvertl = xoldup.size();
3468  int nedges = nvertl - 1;
3469  Array<OneD, NekDouble> xnew_down(nvertl);
3470  Array<OneD, NekDouble> ynew_down(nvertl);
3471  Array<OneD, NekDouble> xnew_up(nvertl);
3472  Array<OneD, NekDouble> ynew_up(nvertl);
3473  Array<OneD, NekDouble> nxvert(nvertl);
3474  Array<OneD, NekDouble> nyvert(nvertl);
3475  Array<OneD, NekDouble> norm(nvertl);
3476  Array<OneD, NekDouble> tmp(nvertl);
3477  for (int a = 0; a < nedges; a++)
3478  {
3479  if (a == 0)
3480  {
3481  // v1
3482  xnew_down[a] = xlaydown[a * npedge + 0];
3483  ynew_down[a] = ylaydown[a * npedge + 0];
3484  xnew_up[a] = xlayup[a * npedge + 0];
3485  ynew_up[a] = ylayup[a * npedge + 0];
3486  nxvert[a] = nxPhys[a * npedge + 0];
3487  nyvert[a] = nyPhys[a * npedge + 0];
3488  // v2
3489  xnew_down[a + 1] = xlaydown[a * npedge + npedge - 1];
3490  ynew_down[a + 1] = ylaydown[a * npedge + npedge - 1];
3491  xnew_up[a + 1] = xlayup[a * npedge + npedge - 1];
3492  ynew_up[a + 1] = ylayup[a * npedge + npedge - 1];
3493  nxvert[a + 1] = nxPhys[a * npedge + npedge - 1];
3494  nyvert[a + 1] = nyPhys[a * npedge + npedge - 1];
3495  }
3496  else
3497  {
3498  // v2
3499  xnew_down[a + 1] = xlaydown[a * npedge + npedge - 1];
3500  ynew_down[a + 1] = ylaydown[a * npedge + npedge - 1];
3501  xnew_up[a + 1] = xlayup[a * npedge + npedge - 1];
3502  ynew_up[a + 1] = ylayup[a * npedge + npedge - 1];
3503  nxvert[a + 1] = nxPhys[a * npedge + npedge - 1];
3504  nyvert[a + 1] = nyPhys[a * npedge + npedge - 1];
3505  }
3506  }
3507 
3508  NekDouble xmax = Vmath::Vmax(nvertl, xoldup, 1);
3509  NekDouble xmin = Vmath::Vmin(nvertl, xoldup, 1);
3510 
3511  // update vertices coords outside layers region
3512  NekDouble ratiox;
3513  // NekDouble ratioy;
3514 
3515  int nVertTot = mesh->GetNvertices();
3516  for (int n = 0; n < nVertTot; n++)
3517  {
3518  NekDouble ratio;
3519  SpatialDomains::PointGeomSharedPtr vertex = mesh->GetVertex(n);
3520  NekDouble x, y, z;
3521  vertex->GetCoords(x, y, z);
3522  int qp_closeroldup = 0, qp_closerolddown = 0;
3523  NekDouble diffup, diffdown;
3524  // determine the closer xold_up,down
3525  diffdown = 1000;
3526  diffup = 1000;
3527 
3528  for (int k = 0; k < nvertl; k++)
3529  {
3530  if (abs(x - xolddown[k]) < diffdown)
3531  {
3532  diffdown = abs(x - xolddown[k]);
3533  qp_closerolddown = k;
3534  }
3535  if (abs(x - xoldup[k]) < diffup)
3536  {
3537  diffup = abs(x - xoldup[k]);
3538  qp_closeroldup = k;
3539  }
3540  }
3541 
3542  // find nplay_closer
3543  diffdown = 1000;
3544  diffup = 1000;
3545 
3546  int qp_closerup = 0, qp_closerdown = 0;
3547 
3548  for (int f = 0; f < nvertl; f++)
3549  {
3550  if (abs(x - xnew_down[f]) < diffdown)
3551  {
3552  diffdown = abs(x - xnew_down[f]);
3553  qp_closerdown = f;
3554  }
3555  if (abs(x - xnew_up[f]) < diffup)
3556  {
3557  diffup = abs(x - xnew_up[f]);
3558  qp_closerup = f;
3559  }
3560  }
3561 
3562  // int qp_closernormoldup;
3563  Vmath::Sadd(nvertl, -x, xoldup, 1, tmp, 1);
3564  Vmath::Vmul(nvertl, tmp, 1, tmp, 1, tmp, 1);
3565  Vmath::Sadd(nvertl, -y, yoldup, 1, norm, 1);
3566  Vmath::Vmul(nvertl, norm, 1, norm, 1, norm, 1);
3567  Vmath::Vadd(nvertl, norm, 1, tmp, 1, norm, 1);
3568  // qp_closernormoldup = Vmath::Imin(nvertl, norm,1);
3569 
3570  Vmath::Zero(nvertl, norm, 1);
3571  Vmath::Zero(nvertl, tmp, 1);
3572 
3573  // int qp_closernormolddown;
3574  Vmath::Sadd(nvertl, -x, xolddown, 1, tmp, 1);
3575  Vmath::Vmul(nvertl, tmp, 1, tmp, 1, tmp, 1);
3576  Vmath::Sadd(nvertl, -y, yolddown, 1, norm, 1);
3577  Vmath::Vmul(nvertl, norm, 1, norm, 1, norm, 1);
3578  Vmath::Vadd(nvertl, norm, 1, tmp, 1, norm, 1);
3579  // qp_closernormolddown = Vmath::Imin(nvertl, norm,1);
3580 
3581  Vmath::Zero(nvertl, norm, 1);
3582  Vmath::Zero(nvertl, tmp, 1);
3583 
3584  int qp_closernormup;
3585  Vmath::Sadd(nvertl, -x, xnew_up, 1, tmp, 1);
3586  Vmath::Vmul(nvertl, tmp, 1, tmp, 1, tmp, 1);
3587  Vmath::Sadd(nvertl, -y, ynew_up, 1, norm, 1);
3588  Vmath::Vmul(nvertl, norm, 1, norm, 1, norm, 1);
3589  Vmath::Vadd(nvertl, norm, 1, tmp, 1, norm, 1);
3590  qp_closernormup = Vmath::Imin(nvertl, norm, 1);
3591 
3592  Vmath::Zero(nvertl, norm, 1);
3593  Vmath::Zero(nvertl, tmp, 1);
3594 
3595  int qp_closernormdown;
3596  Vmath::Sadd(nvertl, -x, xnew_down, 1, tmp, 1);
3597  Vmath::Vmul(nvertl, tmp, 1, tmp, 1, tmp, 1);
3598  Vmath::Sadd(nvertl, -y, ynew_down, 1, norm, 1);
3599  Vmath::Vmul(nvertl, norm, 1, norm, 1, norm, 1);
3600  Vmath::Vadd(nvertl, norm, 1, tmp, 1, norm, 1);
3601  qp_closernormdown = Vmath::Imin(nvertl, norm, 1);
3602 
3603  if (y > yoldup[qp_closeroldup] && y < 1)
3604  {
3605 
3606  // ratio =
3607  // (1-layers_y[nlays-1][qp_closer])*(1-y_c[qp_closer])/
3608  // ( (1-yold_up[n])*(1-yold_c[qp_closer]) );
3609  ratio = (1 - ynew_up[qp_closerup]) / ((1 - yoldup[qp_closeroldup]));
3610 
3611  // ratioy = (1-ynew_up[qp_closernormup])/
3612  // ( (1-yoldup[qp_closernormoldup]) );
3613  // distance prop to layerup
3614  ynew[n] =
3615  ynew_up[qp_closerup] + (y - yoldup[qp_closeroldup]) * ratio;
3616  // ynew[n] = y
3617  // +abs(nyvert[qp_closernormup])*(ynew_up[qp_closeroldup]-yoldup[qp_closeroldup])*ratioy;
3618 
3619  // ynew[n] = y + 0.3*(ynew_up[qp_closerup]-yoldup[qp_closerup]);
3620  // xnew[n] = x +
3621  // abs(nxvert[qp_closeroldup])*(xnew_up[qp_closeroldup]-xoldup[qp_closeroldup]);
3622 
3623  if (x > (xmax - xmin) / 2. && x < xmax)
3624  {
3625  ratiox = (xmax - xnew_up[qp_closernormup]) /
3626  (xmax - xoldup[qp_closernormup]);
3627  if ((xmax - xoldup[qp_closernormup]) == 0)
3628  {
3629  ratiox = 1.0;
3630  }
3631 
3632  // xnew[n] = xnew_up[qp_closerup]
3633  // + (x-xoldup[qp_closerup])*ratiox;
3634  xnew[n] =
3635  x + abs(nxvert[qp_closernormup]) *
3636  (xnew_up[qp_closeroldup] - xoldup[qp_closeroldup]) *
3637  ratiox;
3638  ASSERTL0(x > xmin, " x value <xmin up second half");
3639  ASSERTL0(x < xmax, " x value >xmax up second half");
3640  }
3641  else if (x > xmin && x <= (xmax - xmin) / 2.)
3642  {
3643  // cout<<"up close normold="<<qp_closernormoldup<<"
3644  // closenorm="<<qp_closernormup<<endl;
3645  ratiox = (xnew_up[qp_closernormup] - xmin) /
3646  ((xoldup[qp_closernormup] - xmin));
3647  if ((xoldup[qp_closernormup] - xmin) == 0)
3648  {
3649  ratiox = 1.0;
3650  }
3651  // xnew[n] = xnew_up[qp_closerup]
3652  // + (x-xoldup[qp_closeroldup])*ratiox;
3653  xnew[n] =
3654  x + abs(nxvert[qp_closernormup]) *
3655  (xnew_up[qp_closeroldup] - xoldup[qp_closeroldup]) *
3656  ratiox;
3657  // cout<<"up xold="<<x<<" xnew="<<xnew[n]<<endl;
3658  ASSERTL0(x > xmin, " x value <xmin up first half");
3659  ASSERTL0(x < xmax, " x value >xmax up first half");
3660  }
3661  }
3662  else if (y < yolddown[qp_closerolddown] && y > -1)
3663  {
3664 
3665  ratio = (1 + ynew_down[qp_closerdown]) /
3666  ((1 + yolddown[qp_closerolddown]));
3667 
3668  // ratioy = (1-ynew_down[qp_closernormdown])/
3669  // ( (1-yolddown[qp_closernormolddown]) );
3670 
3671  // distance prop to layerlow
3672  ynew[n] = ynew_down[qp_closerdown] +
3673  (y - yolddown[qp_closerolddown]) * ratio;
3674  // ynew[n] = y +abs(nyvert[qp_closernormdown])*
3675  // (ynew_down[qp_closerolddown]-yolddown[qp_closerolddown])*ratioy;
3676  // ynew[n] = y +
3677  // 0.3*(ynew_down[qp_closerdown]-yolddown[qp_closerdown]); xnew[n] =
3678  // x +
3679  // abs(nxvert[qp_closerolddown])*(xnew_down[qp_closerolddown]-xolddown[qp_closerolddown]);
3680  /*
3681  if(n==74)
3682  {
3683  cout<<qp_closerolddown<<" nplaydown="<<qp_closerdown<<endl;
3684  cout<<"xolddown="<<xolddown[qp_closerolddown]<<"
3685  xnewdown="<<xnew_down[qp_closerdown]<<endl; cout<<"xold+"<<x<<"
3686  xnew+"<<xnew[n]<<endl;
3687  }
3688  */
3689 
3690  if (x > (xmax - xmin) / 2. && x < xmax)
3691  {
3692  ratiox = (xmax - xnew_down[qp_closernormdown]) /
3693  ((xmax - xolddown[qp_closernormdown]));
3694  if ((xmax - xolddown[qp_closernormdown]) == 0)
3695  {
3696  ratiox = 1.0;
3697  }
3698  // xnew[n] = xnew_down[qp_closerdown]
3699  // + (x-xolddown[qp_closerolddown])*ratiox;
3700  xnew[n] = x + abs(nxvert[qp_closernormdown]) *
3701  (xnew_down[qp_closerolddown] -
3702  xolddown[qp_closerolddown]) *
3703  ratiox;
3704  ASSERTL0(x > xmin, " x value <xmin down second half");
3705  ASSERTL0(x < xmax, " x value >xmax down second half");
3706  }
3707  else if (x > xmin && x <= (xmax - xmin) / 2.)
3708  {
3709  ratiox = (xnew_down[qp_closernormdown] - xmin) /
3710  ((xolddown[qp_closernormdown] - xmin));
3711  if ((xolddown[qp_closernormdown] - xmin) == 0)
3712  {
3713  ratiox = 1.0;
3714  }
3715  // xnew[n] = xnew_down[qp_closerdown]
3716  // + (x-xolddown[qp_closerolddown])*ratiox;
3717  xnew[n] = x + abs(nxvert[qp_closernormdown]) *
3718  (xnew_down[qp_closerolddown] -
3719  xolddown[qp_closerolddown]) *
3720  ratiox;
3721  ASSERTL0(x > xmin, " x value <xmin down first half");
3722  ASSERTL0(x < xmax, " x value >xmax down first half");
3723  }
3724  }
3725 
3726  cout << "xold" << x << " xnew=" << xnew[n] << endl;
3727  ASSERTL0(xnew[n] >= xmin, "newx < xmin");
3728  ASSERTL0(xnew[n] <= xmax, "newx > xmax");
3729 
3730  } // verts closed
3731 }
3732 
3733 void CheckSingularQuads(MultiRegions::ExpListSharedPtr Exp, Array<OneD, int> V1,
3734  Array<OneD, int> V2, Array<OneD, NekDouble> &xnew,
3735  Array<OneD, NekDouble> &ynew)
3736 {
3737  const std::shared_ptr<LocalRegions::ExpansionVector> exp2D = Exp->GetExp();
3738  int nel = exp2D->size();
3739  LocalRegions::QuadExpSharedPtr locQuadExp;
3740  LocalRegions::TriExpSharedPtr locTriExp;
3741  SpatialDomains::Geometry1DSharedPtr SegGeom;
3742  int idbef, idnext;
3743  NekDouble xV1, yV1, xV2, yV2;
3744  NekDouble slopebef = 0.0, slopenext = 0.0, slopenew = 0.0;
3745  Array<OneD, int> locEids(4);
3746  for (int i = 0; i < nel; i++)
3747  {
3748  if ((locQuadExp = (*exp2D)[i]->as<LocalRegions::QuadExp>()))
3749  {
3750  SegGeom = (locQuadExp->GetGeom2D())->GetEdge(0);
3751  idbef = SegGeom->GetGlobalID();
3752  if (xnew[V1[idbef]] <= xnew[V2[idbef]])
3753  {
3754  xV1 = xnew[V1[idbef]];
3755  yV1 = ynew[V1[idbef]];
3756  xV2 = xnew[V2[idbef]];
3757  yV2 = ynew[V2[idbef]];
3758  slopebef = (yV2 - yV1) / (xV2 - xV1);
3759  }
3760  else
3761  {
3762  xV1 = xnew[V2[idbef]];
3763  yV1 = ynew[V2[idbef]];
3764  xV2 = xnew[V1[idbef]];
3765  yV2 = ynew[V1[idbef]];
3766  slopebef = (yV2 - yV1) / (xV2 - xV1);
3767  }
3768  // cout<<"00 V1 x="<<xnew[ V1[idbef] ]<<" y="<<ynew[ V1[idbef]
3769  // ]<<endl; cout<<"00 V2 x="<<xnew[ V2[idbef] ]<<" y="<<ynew[
3770  // V2[idbef] ]<<endl;
3771  for (int j = 1; j < locQuadExp->GetNtraces(); ++j)
3772  {
3773  SegGeom = (locQuadExp->GetGeom2D())->GetEdge(j);
3774  idnext = SegGeom->GetGlobalID();
3775  // cout<<"id="<<idnext<<" locid="<<j<<endl;
3776  // cout<<" V1 x="<<xnew[ V1[idnext] ]<<" y="<<ynew[
3777  // V1[idnext] ]<<endl; cout<<" V2 x="<<xnew[ V2[idnext] ]<<"
3778  // y="<<ynew[ V2[idnext] ]<<endl;
3779  if (xV1 == xnew[V1[idnext]] && yV1 == ynew[V1[idnext]])
3780  {
3781  xV1 = xnew[V1[idnext]];
3782  yV1 = ynew[V1[idnext]];
3783  xV2 = xnew[V2[idnext]];
3784  yV2 = ynew[V2[idnext]];
3785  slopenext = (yV2 - yV1) / (xV2 - xV1);
3786  if (i == 23)
3787  {
3788  cout << "case1 x0=" << xV1 << " x1=" << xV2 << endl;
3789  cout << idnext << " 11slope bef =" << slopebef
3790  << " slopenext=" << slopenext << endl;
3791  }
3792  // compare with slope before
3793  if (slopebef / slopenext > 0.84 &&
3794  slopebef / slopenext < 1.18)
3795  {
3796  xnew[V1[idnext]] = xnew[V1[idnext]] - 0.01;
3797  slopenew = (yV2 - yV1) / (xV2 - xnew[V1[idnext]]);
3798 
3799  if (abs(slopebef - slopenew) <
3800  abs(slopebef - slopenext))
3801  {
3802  xnew[V1[idnext]] = xnew[V1[idnext]] + 0.02;
3803  slopenew = (yV2 - yV1) / (xV2 - xnew[V1[idnext]]);
3804  }
3805  slopenext = slopenew;
3806  cout << "slopenew=" << slopenew << endl;
3807  cout << "moved x=" << xnew[V1[idnext]] << endl;
3808  }
3809  }
3810  else if (xV2 == xnew[V2[idnext]] && yV2 == ynew[V2[idnext]])
3811  {
3812  xV1 = xnew[V2[idnext]];
3813  yV1 = ynew[V2[idnext]];
3814  xV2 = xnew[V1[idnext]];
3815  yV2 = ynew[V1[idnext]];
3816  slopenext = (yV2 - yV1) / (xV2 - xV1);
3817  if (i == 23)
3818  {
3819  cout << "case2 x0=" << xV1 << " x1=" << xV2 << endl;
3820  cout << idnext << " 22slope bef =" << slopebef
3821  << " slopenext=" << slopenext << endl;
3822  }
3823  // compare with slope before
3824  if (slopebef / slopenext > 0.84 &&
3825  slopebef / slopenext < 1.18)
3826  {
3827  xnew[V2[idnext]] = xnew[V2[idnext]] - 0.01;
3828  slopenew = (yV2 - yV1) / (xV2 - xnew[V2[idnext]]);
3829 
3830  if (abs(slopebef - slopenew) <
3831  abs(slopebef - slopenext))
3832  {
3833  xnew[V2[idnext]] = xnew[V2[idnext]] + 0.02;
3834  slopenew = (yV2 - yV1) / (xV2 - xnew[V2[idnext]]);
3835  }
3836  slopenext = slopenew;
3837  cout << "slopenew=" << slopenew << endl;
3838  cout << "moved x=" << xnew[V2[idnext]] << endl;
3839  }
3840  }
3841  else if (xV1 == xnew[V2[idnext]] && yV1 == ynew[V2[idnext]])
3842  {
3843  xV1 = xnew[V2[idnext]];
3844  yV1 = ynew[V2[idnext]];
3845  xV2 = xnew[V1[idnext]];
3846  yV2 = ynew[V1[idnext]];
3847  slopenext = (yV2 - yV1) / (xV2 - xV1);
3848  if (i == 23)
3849  {
3850  cout << "case3 x0=" << xV1 << " x1=" << xV2 << endl;
3851  cout << idnext << " 22slope bef =" << slopebef
3852  << " slopenext=" << slopenext << endl;
3853  }
3854  // compare with slope before
3855  if (slopebef / slopenext > 0.84 &&
3856  slopebef / slopenext < 1.18)
3857  {
3858  xnew[V2[idnext]] = xnew[V2[idnext]] - 0.01;
3859  slopenew = (yV2 - yV1) / (xV2 - xnew[V2[idnext]]);
3860 
3861  if (abs(slopebef - slopenew) <
3862  abs(slopebef - slopenext))
3863  {
3864  xnew[V2[idnext]] = xnew[V2[idnext]] + 0.02;
3865  slopenew = (yV2 - yV1) / (xV2 - xnew[V2[idnext]]);
3866  }
3867  slopenext = slopenew;
3868  cout << "slopenew=" << slopenew << endl;
3869  cout << "moved x=" << xnew[V2[idnext]] << endl;
3870  }
3871  }
3872  else if (xV2 == xnew[V1[idnext]] && yV2 == ynew[V1[idnext]])
3873  {
3874  xV1 = xnew[V1[idnext]];
3875  yV1 = ynew[V1[idnext]];
3876  xV2 = xnew[V2[idnext]];
3877  yV2 = ynew[V2[idnext]];
3878  slopenext = (yV2 - yV1) / (xV2 - xV1);
3879  if (i == 23)
3880  {
3881  cout << "case4 x0=" << xV1 << " x1=" << xV2 << endl;
3882  cout << idnext << " 22slope bef =" << slopebef
3883  << " slopenext=" << slopenext << endl;
3884  }
3885  // compare with slope before
3886  if (slopebef / slopenext > 0.84 &&
3887  slopebef / slopenext < 1.18)
3888  {
3889  xnew[V1[idnext]] = xnew[V1[idnext]] - 0.01;
3890  slopenew = (yV2 - yV1) / (xV2 - xnew[V1[idnext]]);
3891 
3892  if (abs(slopebef - slopenew) <
3893  abs(slopebef - slopenext))
3894  {
3895  xnew[V1[idnext]] = xnew[V1[idnext]] + 0.02;
3896  slopenew = (yV2 - yV1) / (xV2 - xnew[V1[idnext]]);
3897  }
3898  slopenext = slopenew;
3899  cout << "slopenew=" << slopenew << endl;
3900  cout << "moved x=" << xnew[V1[idnext]] << endl;
3901  }
3902  }
3903  else
3904  {
3905  ASSERTL0(false, "edge not connected");
3906  }
3907  slopebef = slopenext;
3908  }
3909  }
3910  }
3911 }
3912 
3913 void Replacevertices(string filename, Array<OneD, NekDouble> newx,
3914  Array<OneD, NekDouble> newy, Array<OneD, NekDouble> xcPhys,
3915  Array<OneD, NekDouble> ycPhys, Array<OneD, int> Eids,
3916  int Npoints, string s_alp,
3917  Array<OneD, Array<OneD, NekDouble>> x_lay,
3918  Array<OneD, Array<OneD, NekDouble>> y_lay,
3919  Array<OneD, Array<OneD, int>> lay_eids, bool curv_lay)
3920 {
3921  // load existing file
3922  string newfile;
3923  TiXmlDocument doc(filename);
3924  // load xscale parameter (if exists)
3925  TiXmlElement *master = doc.FirstChildElement("NEKTAR");
3926  TiXmlElement *mesh = master->FirstChildElement("GEOMETRY");
3927  TiXmlElement *element = mesh->FirstChildElement("VERTEX");
3928  NekDouble xscale = 1.0;
3929  LibUtilities::Interpreter expEvaluator;
3930  const char *xscal = element->Attribute("XSCALE");
3931  if (xscal)
3932  {
3933  std::string xscalstr = xscal;
3934  int expr_id = expEvaluator.DefineFunction("", xscalstr);
3935  xscale = expEvaluator.Evaluate(expr_id);
3936  }
3937 
3938  // Save a new XML file.
3939  newfile = filename.substr(0, filename.find_last_of(".")) + "_moved.xml";
3940  doc.SaveFile(newfile);
3941 
3942  // write the new vertices
3943  TiXmlDocument docnew(newfile);
3944  bool loadOkaynew = docnew.LoadFile();
3945 
3946  std::string errstr = "Unable to load file: ";
3947  errstr += newfile;
3948  ASSERTL0(loadOkaynew, errstr.c_str());
3949 
3950  TiXmlHandle docHandlenew(&docnew);
3951  TiXmlElement *meshnew = NULL;
3952  TiXmlElement *masternew = NULL;
3953  TiXmlElement *condnew = NULL;
3954  TiXmlElement *Parsnew = NULL;
3955  TiXmlElement *parnew = NULL;
3956 
3957  // Master tag within which all data is contained.
3958 
3959  masternew = docnew.FirstChildElement("NEKTAR");
3960  ASSERTL0(masternew, "Unable to find NEKTAR tag in file.");
3961 
3962  // set the alpha value
3963  string alphastring;
3964  condnew = masternew->FirstChildElement("CONDITIONS");
3965  Parsnew = condnew->FirstChildElement("PARAMETERS");
3966  cout << "alpha=" << s_alp << endl;
3967  parnew = Parsnew->FirstChildElement("P");
3968  while (parnew)
3969  {
3970  TiXmlNode *node = parnew->FirstChild();
3971  if (node)
3972  {
3973  // Format is "paramName = value"
3974  std::string line = node->ToText()->Value();
3975  std::string lhs;
3976  std::string rhs;
3977  /// Pull out lhs and rhs and eliminate any spaces.
3978  int beg = line.find_first_not_of(" ");
3979  int end = line.find_first_of("=");
3980  // Check for no parameter name
3981  if (beg == end)
3982  throw 1;
3983  // Check for no parameter value
3984  if (end != line.find_last_of("="))
3985  throw 1;
3986  // Check for no equals sign
3987  if (end == std::string::npos)
3988  throw 1;
3989  lhs = line.substr(line.find_first_not_of(" "), end - beg);
3990  lhs = lhs.substr(0, lhs.find_last_not_of(" ") + 1);
3991 
3992  // rhs = line.substr(line.find_last_of("=")+1);
3993  // rhs = rhs.substr(rhs.find_first_not_of(" "));
3994  // rhs = rhs.substr(0, rhs.find_last_not_of(" ")+1);
3995 
3996  boost::to_upper(lhs);
3997  if (lhs == "ALPHA")
3998  {
3999  alphastring = "Alpha = " + s_alp;
4000  parnew->RemoveChild(node);
4001  parnew->LinkEndChild(new TiXmlText(alphastring));
4002  }
4003  }
4004 
4005  parnew = parnew->NextSiblingElement("P");
4006  }
4007 
4008  // Find the Mesh tag and same the dim and space attributes
4009  meshnew = masternew->FirstChildElement("GEOMETRY");
4010 
4011  ASSERTL0(meshnew, "Unable to find GEOMETRY tag in file.");
4012  // Now read the vertices
4013  TiXmlElement *elementnew = meshnew->FirstChildElement("VERTEX");
4014  ASSERTL0(elementnew, "Unable to find mesh VERTEX tag in file.");
4015  // set xscale 1!!
4016  if (xscale != 1.0)
4017  {
4018  elementnew->SetAttribute("XSCALE", 1.0);
4019  }
4020  TiXmlElement *vertexnew = elementnew->FirstChildElement("V");
4021 
4022  int indx;
4023  int err, numPts;
4024  int nextVertexNumber = -1;
4025 
4026  while (vertexnew)
4027  {
4028  nextVertexNumber++;
4029  // delete the old one
4030  TiXmlAttribute *vertexAttr = vertexnew->FirstAttribute();
4031  std::string attrName(vertexAttr->Name());
4032  ASSERTL0(attrName == "ID",
4033  (std::string("Unknown attribute name: ") + attrName).c_str());
4034 
4035  err = vertexAttr->QueryIntValue(&indx);
4036  ASSERTL0(err == TIXML_SUCCESS, "Unable to read attribute ID.");
4037  ASSERTL0(indx == nextVertexNumber,
4038  "Vertex IDs must begin with zero and be sequential.");
4039 
4040  std::string vertexBodyStr;
4041  // Now read body of vertex
4042  TiXmlNode *vertexBody = vertexnew->FirstChild();
4043  // Accumulate all non-comment body data.
4044  if (vertexBody->Type() == TiXmlNode::TINYXML_TEXT)
4045  {
4046  vertexBodyStr += vertexBody->ToText()->Value();
4047  vertexBodyStr += " ";
4048  }
4049  ASSERTL0(!vertexBodyStr.empty(),
4050  "Vertex definitions must contain vertex data.");
4051  // remove the old coordinates
4052  vertexnew->RemoveChild(vertexBody);
4053  // write the new one
4054  // cout<<"writing.. v:"<<nextVertexNumber<<endl;
4055  stringstream s;
4056  // we need at least 5 digits (setprecision 5) to get the streak position
4057  // with
4058  // precision 10^-10
4059  s << std::scientific << std::setprecision(8) << newx[nextVertexNumber]
4060  << " " << newy[nextVertexNumber] << " " << 0.0;
4061  vertexnew->LinkEndChild(new TiXmlText(s.str()));
4062  // TiXmlNode *newvertexBody = vertexnew->FirstChild();
4063  // string newvertexbodystr= newvertexBody->SetValue(s.str());
4064  // vertexnew->ReplaceChild(vertexBody,new TiXmlText(newvertexbodystr));
4065 
4066  vertexnew = vertexnew->NextSiblingElement("V");
4067  }
4068 
4069  // read the curved tag
4070  TiXmlElement *curvednew = meshnew->FirstChildElement("CURVED");
4071  ASSERTL0(curvednew, "Unable to find mesh CURVED tag in file.");
4072  TiXmlElement *edgenew = curvednew->FirstChildElement("E");
4073  int cnt = -1;
4074  // ID is different from index...
4075  std::string charindex;
4076  int eid;
4077  int index;
4078  int indexeid;
4079  int neids_lay = lay_eids[0].size();
4080  // if edgenew belongs to the crit lay replace it, else delete it.
4081  while (edgenew)
4082  {
4083  indexeid = -1;
4084  cnt++;
4085  // get the index...
4086  TiXmlAttribute *edgeAttr = edgenew->FirstAttribute();
4087  std::string attrName(edgeAttr->Name());
4088  charindex = edgeAttr->Value();
4089  std::istringstream iss(charindex);
4090  iss >> std::dec >> index;
4091  // get the eid
4092  edgenew->QueryIntAttribute("EDGEID", &eid);
4093  // cout<<"eid="<<eid<<" neid="<<Eids.size()<<endl;
4094  // find the corresponding index curve point
4095  for (int u = 0; u < Eids.size(); u++)
4096  {
4097  // cout<<"Eids="<<Eids[u]<<" eid="<<eid<<endl;
4098  if (Eids[u] == eid)
4099  {
4100  indexeid = u;
4101  }
4102  }
4103  if (indexeid == -1)
4104  {
4105  curvednew->RemoveChild(edgenew);
4106  // ASSERTL0(false, "edge to update not found");
4107  }
4108  else
4109  {
4110 
4111  std::string edgeBodyStr;
4112  // read the body of the edge
4113  TiXmlNode *edgeBody = edgenew->FirstChild();
4114  if (edgeBody->Type() == TiXmlNode::TINYXML_TEXT)
4115  {
4116  edgeBodyStr += edgeBody->ToText()->Value();
4117  edgeBodyStr += " ";
4118  }
4119  ASSERTL0(!edgeBodyStr.empty(),
4120  "Edge definitions must contain edge data");
4121  // remove the old coordinates
4122  edgenew->RemoveChild(edgeBody);
4123  // write the new points coordinates
4124  // we need at least 5 digits (setprecision 5) to get the streak
4125  // position with
4126  // precision 10^-10
4127 
4128  // Determine the number of points
4129  err = edgenew->QueryIntAttribute("NUMPOINTS", &numPts);
4130  ASSERTL0(err == TIXML_SUCCESS,
4131  "Unable to read curve attribute NUMPOINTS.");
4132 
4133  stringstream st;
4134  edgenew->SetAttribute("NUMPOINTS", Npoints);
4135  for (int u = 0; u < Npoints; u++)
4136  {
4137  st << std::scientific << std::setprecision(8)
4138  << xcPhys[cnt * Npoints + u] << " "
4139  << ycPhys[cnt * Npoints + u] << " " << 0.000 << " ";
4140  }
4141 
4142  edgenew->LinkEndChild(new TiXmlText(st.str()));
4143 
4144  /*
4145  st << std::scientific <<
4146  std::setprecision(8) << x_crit[v1] << " "
4147  << y_crit[v1] << " " << 0.000<<" ";
4148  for(int a=0; a< Npoints-2; a++)
4149  {
4150  st << std::scientific <<
4151  std::setprecision(8) << " "<<Pcurvx[indexeid*(Npoints-2)
4152  +a]<<" "<<Pcurvy[indexeid*(Npoints-2) +a]
4153  <<" "<<0.000<<" ";
4154 
4155  }
4156  st << std::scientific <<
4157  std::setprecision(8) << " "<<x_crit[v2]<<" "<< y_crit[v2]
4158  <<" "<< 0.000; edgenew->LinkEndChild(new TiXmlText(st.str()));
4159 
4160  */
4161  }
4162 
4163  edgenew = edgenew->NextSiblingElement("E");
4164  }
4165 
4166  // write also the others layers curve points
4167  if (curv_lay == true)
4168  {
4169  cout << "write other curved edges" << endl;
4170  TiXmlElement *curved = meshnew->FirstChildElement("CURVED");
4171  int idcnt = 300;
4172  int nlays = lay_eids.size();
4173 
4174  // TiXmlComment * comment = new TiXmlComment();
4175  // comment->SetValue(" new edges ");
4176  // curved->LinkEndChild(comment);
4177  for (int g = 0; g < nlays; ++g)
4178  {
4179  for (int p = 0; p < neids_lay; p++)
4180  {
4181  stringstream st;
4182  TiXmlElement *e = new TiXmlElement("E");
4183  e->SetAttribute("ID", idcnt++);
4184  e->SetAttribute("EDGEID", lay_eids[g][p]);
4185  e->SetAttribute("NUMPOINTS", Npoints);
4186  e->SetAttribute("TYPE", "PolyEvenlySpaced");
4187  for (int c = 0; c < Npoints; c++)
4188  {
4189  st << std::scientific << std::setprecision(8)
4190  << x_lay[g][p * Npoints + c] << " "
4191  << y_lay[g][p * Npoints + c] << " " << 0.000 << " ";
4192  }
4193 
4194  TiXmlText *t0 = new TiXmlText(st.str());
4195  e->LinkEndChild(t0);
4196  curved->LinkEndChild(e);
4197  }
4198  }
4199  }
4200 
4201  docnew.SaveFile(newfile);
4202 
4203  cout << "new file: " << newfile << endl;
4204 }
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:2320
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:3314
EvaluateTangent
void EvaluateTangent(int npoints, Array< OneD, NekDouble > xcQedge, Array< OneD, NekDouble > coeffsinterp, Array< OneD, NekDouble > &txQedge, Array< OneD, NekDouble > &tyQedge)
Definition: MeshMove.cpp:2895
LagrangeInterpolant
NekDouble LagrangeInterpolant(NekDouble x, int npts, Array< OneD, NekDouble > xpts, Array< OneD, NekDouble > funcvals)
Definition: MeshMove.cpp:2867
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:2104
Cutrepetitions
void Cutrepetitions(int nedges, Array< OneD, NekDouble > inarray, Array< OneD, NekDouble > &outarray)
Definition: MeshMove.cpp:2770
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:3441
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:3733
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:3144
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:3237
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:2009
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:3033
GenerateAddPointsNewtonIt
void GenerateAddPointsNewtonIt(NekDouble xi, NekDouble yi, NekDouble &xout, NekDouble &yout, MultiRegions::ExpListSharedPtr function, Array< OneD, NekDouble > derfunction, NekDouble cr)
Definition: MeshMove.cpp:2230
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:2815
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:3913
checkcommonvert
bool checkcommonvert(Array< OneD, int > Vids_laybefore, Array< OneD, int > Vids_c, int Vid)
Definition: MeshMove.cpp:2755
DetermineclosePointxindex
int DetermineclosePointxindex(NekDouble x, Array< OneD, NekDouble > xArray)
Definition: MeshMove.cpp:2797
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:3189
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:3376
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:2941
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:2383
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/Interpreter.cpp:2121
Nektar::LibUtilities::Interpreter::Evaluate
NekDouble Evaluate(const int id)
Evaluate a function which depends only on constants and/or parameters.
Definition: Interpreter/Interpreter.cpp:2127
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:585
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:290
Nektar::LibUtilities::SessionReaderSharedPtr
std::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: SessionReader.h:115
Nektar::LocalRegions::QuadExpSharedPtr
std::shared_ptr< QuadExp > QuadExpSharedPtr
Definition: QuadExp.h:257
Nektar::LocalRegions::SegExpSharedPtr
std::shared_ptr< SegExp > SegExpSharedPtr
Definition: SegExp.h:251
Nektar::LocalRegions::TriExpSharedPtr
std::shared_ptr< TriExp > TriExpSharedPtr
Definition: TriExp.h:253
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:289
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:2
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 scalar 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:298
tinysimd::sqrt
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:294
Generated on Sat Mar 25 2023 12:35:51 for Nektar++ by doxygen 1.9.1