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