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