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