Nektar++
Geometry2D.cpp
Go to the documentation of this file.
1////////////////////////////////////////////////////////////////////////////////
2//
3// File: Geometry2D.cpp
4//
5// For more information, please see: http://www.nektar.info/
6//
7// The MIT License
8//
9// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10// Department of Aeronautics, Imperial College London (UK), and Scientific
11// Computing and Imaging Institute, University of Utah (USA).
12//
13// Permission is hereby granted, free of charge, to any person obtaining a
14// copy of this software and associated documentation files (the "Software"),
15// to deal in the Software without restriction, including without limitation
16// the rights to use, copy, modify, merge, publish, distribute, sublicense,
17// and/or sell copies of the Software, and to permit persons to whom the
18// Software is furnished to do so, subject to the following conditions:
19//
20// The above copyright notice and this permission notice shall be included
21// in all copies or substantial portions of the Software.
22//
23// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29// DEALINGS IN THE SOFTWARE.
30//
31// Description: 2D geometry information
32//
33//
34////////////////////////////////////////////////////////////////////////////////
35
36#include <iomanip>
37
40
41#include <iomanip>
42
44{
45
47{
48}
49
50Geometry2D::Geometry2D(const int coordim, CurveSharedPtr curve)
51 : Geometry(coordim), m_curve(curve)
52{
54 "Coordinate dimension should be at least 2 for a 2D geometry");
55}
56
58{
59 int nc = 1, d0 = m_manifold[0], d1 = m_manifold[1];
60 if (0 == m_edgeNormal.size())
61 {
64 x[0] = Array<OneD, NekDouble>(3);
65 x[1] = Array<OneD, NekDouble>(3);
66 m_verts[0]->GetCoords(x[0]);
67 int i0 = 1, i1 = 0, direction = 1;
68 for (size_t i = 0; i < m_verts.size(); ++i)
69 {
70 i0 ^= 1;
71 i1 ^= 1;
72 m_verts[(i + 1) % m_verts.size()]->GetCoords(x[i1]);
73 if (m_edges[i]->GetXmap()->GetBasis(0)->GetNumModes() > 2)
74 {
75 continue;
76 }
78 m_edgeNormal[i][0] = x[i0][d1] - x[i1][d1];
79 m_edgeNormal[i][1] = x[i1][d0] - x[i0][d0];
80 }
81 if (m_coordim == 3)
82 {
83 for (size_t i = 0; i < m_verts.size(); ++i)
84 {
85 if (m_edgeNormal[i].size() == 2)
86 {
87 m_verts[i]->GetCoords(x[0]);
88 m_verts[(i + 2) % m_verts.size()]->GetCoords(x[1]);
89 if (m_edgeNormal[i][0] * (x[1][d0] - x[0][d0]) <
90 m_edgeNormal[i][1] * (x[0][d1] - x[1][d1]))
91 {
92 direction = -1;
93 }
94 break;
95 }
96 }
97 }
98 if (direction == -1)
99 {
100 for (size_t i = 0; i < m_verts.size(); ++i)
101 {
102 if (m_edgeNormal[i].size() == 2)
103 {
104 m_edgeNormal[i][0] = -m_edgeNormal[i][0];
105 m_edgeNormal[i][1] = -m_edgeNormal[i][1];
106 }
107 }
108 }
109 }
110
111 Array<OneD, NekDouble> vertex(3);
112 for (size_t i = 0; i < m_verts.size(); ++i)
113 {
114 int i1 = (i + 1) % m_verts.size();
115 if (m_verts[i]->GetVid() < m_verts[i1]->GetVid())
116 {
117 m_verts[i]->GetCoords(vertex);
118 }
119 else
120 {
121 m_verts[i1]->GetCoords(vertex);
122 }
123 if (m_edgeNormal[i].size() == 0)
124 {
125 nc = 0; // not sure
126 continue;
127 }
128 if (m_edgeNormal[i][0] * (gloCoord[d0] - vertex[d0]) <
129 m_edgeNormal[i][1] * (vertex[d1] - gloCoord[d1]))
130 {
131 return -1; // outside
132 }
133 }
134 // 3D manifold needs to check the distance
135 if (m_coordim == 3)
136 {
137 nc = 0;
138 }
139 // nc: 1 (side element), 0 (maybe inside), -1 (outside)
140 return nc;
141}
142
145{
146 int i0 = 0, i1 = 1, j1 = 0, j2 = 1;
147 if (m_straightEdge & 2)
148 {
149 i0 = 1;
150 i1 = 0;
151 }
152 if (m_straightEdge & 4)
153 {
154 j1 = 1;
155 j2 = 0;
156 }
158 NekDouble gamma = m_isoParameter[1][3];
159 NekDouble tty = (coords[i1] - gamma * coords[i0] - m_isoParameter[1][0]) *
160 m_isoParameter[1][1];
161 NekDouble denom = 1. / (m_isoParameter[0][2] + m_isoParameter[0][3] * tty);
162 NekDouble epsilon = -m_isoParameter[0][3] * beta * denom;
163 NekDouble h =
164 (m_isoParameter[0][0] + m_isoParameter[0][1] * tty - coords[i0]) *
165 denom;
166 Lcoords[j2] = -h / (0.5 + sqrt(0.25 - epsilon * h));
167 Lcoords[j1] = -beta * Lcoords[j2] + tty;
168}
169
171 const Array<OneD, const NekDouble> &coords,
174 NekDouble &dist)
175{
176 // Maximum iterations for convergence
177 const int MaxIterations = NekConstants::kNewtonIterations;
178 // |x-xp|^2 < EPSILON error tolerance
179 const NekDouble Tol = 1.e-8;
180 // |r,s| > LcoordDIV stop the search
181 const NekDouble LcoordDiv = 15.0;
182
184 m_geomFactors->GetJac(m_xmap->GetPointsKeys());
185
186 NekDouble ScaledTol =
187 Vmath::Vsum(Jac.size(), Jac, 1) / ((NekDouble)Jac.size());
188 ScaledTol *= Tol;
189
190 NekDouble xmap, ymap, F1, F2;
191 NekDouble derx_1, derx_2, dery_1, dery_2, jac;
192
193 // save intiial guess for later reference if required.
194 NekDouble init0 = Lcoords[0], init1 = Lcoords[1];
195
196 Array<OneD, NekDouble> DxD1(ptsx.size());
197 Array<OneD, NekDouble> DxD2(ptsx.size());
198 Array<OneD, NekDouble> DyD1(ptsx.size());
199 Array<OneD, NekDouble> DyD2(ptsx.size());
200
201 // Ideally this will be stored in m_geomfactors
202 m_xmap->PhysDeriv(ptsx, DxD1, DxD2);
203 m_xmap->PhysDeriv(ptsy, DyD1, DyD2);
204
205 int cnt = 0;
208
209 F1 = F2 = 2000; // Starting value of Function
210 NekDouble resid = sqrt(F1 * F1 + F2 * F2);
211 while (cnt++ < MaxIterations)
212 {
213 // evaluate lagrange interpolant at Lcoords
214 m_xmap->LocCoordToLocCollapsed(Lcoords, eta);
215 I[0] = m_xmap->GetBasis(0)->GetI(eta);
216 I[1] = m_xmap->GetBasis(1)->GetI(eta + 1);
217
218 // calculate the global point `corresponding to Lcoords
219 xmap = m_xmap->PhysEvaluate(I, ptsx);
220 ymap = m_xmap->PhysEvaluate(I, ptsy);
221
222 F1 = coords[0] - xmap;
223 F2 = coords[1] - ymap;
224
225 if (F1 * F1 + F2 * F2 < ScaledTol)
226 {
227 resid = sqrt(F1 * F1 + F2 * F2);
228 break;
229 }
230
231 // Interpolate derivative metric at Lcoords
232 derx_1 = m_xmap->PhysEvaluate(I, DxD1);
233 derx_2 = m_xmap->PhysEvaluate(I, DxD2);
234 dery_1 = m_xmap->PhysEvaluate(I, DyD1);
235 dery_2 = m_xmap->PhysEvaluate(I, DyD2);
236
237 jac = dery_2 * derx_1 - dery_1 * derx_2;
238
239 // use analytical inverse of derivitives which are
240 // also similar to those of metric factors.
241 Lcoords[0] =
242 Lcoords[0] +
243 (dery_2 * (coords[0] - xmap) - derx_2 * (coords[1] - ymap)) / jac;
244
245 Lcoords[1] =
246 Lcoords[1] +
247 (-dery_1 * (coords[0] - xmap) + derx_1 * (coords[1] - ymap)) / jac;
248
249 if (!(std::isfinite(Lcoords[0]) && std::isfinite(Lcoords[1])))
250 {
251 dist = 1e16;
252 std::ostringstream ss;
253 ss << "nan or inf found in NewtonIterationForLocCoord in element "
254 << GetGlobalID();
255 WARNINGL1(false, ss.str());
256 return;
257 }
258 if (fabs(Lcoords[0]) > LcoordDiv || fabs(Lcoords[1]) > LcoordDiv)
259 {
260 break; // lcoords have diverged so stop iteration
261 }
262 }
263
264 m_xmap->LocCoordToLocCollapsed(Lcoords, eta);
265 if (ClampLocCoords(eta, 0.))
266 {
267 I[0] = m_xmap->GetBasis(0)->GetI(eta);
268 I[1] = m_xmap->GetBasis(1)->GetI(eta + 1);
269 // calculate the global point corresponding to Lcoords
270 xmap = m_xmap->PhysEvaluate(I, ptsx);
271 ymap = m_xmap->PhysEvaluate(I, ptsy);
272 F1 = coords[0] - xmap;
273 F2 = coords[1] - ymap;
274 dist = sqrt(F1 * F1 + F2 * F2);
275 }
276 else
277 {
278 dist = 0.;
279 }
280
281 if (cnt >= MaxIterations)
282 {
283 Array<OneD, NekDouble> collCoords(2);
284 m_xmap->LocCoordToLocCollapsed(Lcoords, collCoords);
285
286 // if coordinate is inside element dump error!
287 if ((collCoords[0] >= -1.0 && collCoords[0] <= 1.0) &&
288 (collCoords[1] >= -1.0 && collCoords[1] <= 1.0))
289 {
290 std::ostringstream ss;
291
292 ss << "Reached MaxIterations (" << MaxIterations
293 << ") in Newton iteration ";
294 ss << "Init value (" << std::setprecision(4) << init0 << ","
295 << init1 << ","
296 << ") ";
297 ss << "Fin value (" << Lcoords[0] << "," << Lcoords[1] << ","
298 << ") ";
299 ss << "Resid = " << resid
300 << " Tolerance = " << std::sqrt(ScaledTol);
301
302 WARNINGL1(cnt < MaxIterations, ss.str());
303 }
304 }
305}
306
308 Array<OneD, NekDouble> &Lcoords)
309{
310 NekDouble dist = std::numeric_limits<double>::max();
311 Array<OneD, NekDouble> tmpcoords(2);
312 tmpcoords[0] = coords[m_manifold[0]];
313 tmpcoords[1] = coords[m_manifold[1]];
314 if (GetMetricInfo()->GetGtype() == eRegular)
315 {
316 tmpcoords[0] -= m_isoParameter[0][0];
317 tmpcoords[1] -= m_isoParameter[1][0];
318 Lcoords[0] = m_invIsoParam[0][0] * tmpcoords[0] +
319 m_invIsoParam[0][1] * tmpcoords[1];
320 Lcoords[1] = m_invIsoParam[1][0] * tmpcoords[0] +
321 m_invIsoParam[1][1] * tmpcoords[1];
322 }
323 else if (m_straightEdge)
324 {
325 SolveStraightEdgeQuad(tmpcoords, Lcoords);
326 }
327 else if (GetMetricInfo()->GetGtype() == eDeformed)
328 {
329 v_FillGeom();
330 // Determine nearest point of coords to values in m_xmap
331 int npts = m_xmap->GetTotPoints();
332 Array<OneD, NekDouble> ptsx(npts), ptsy(npts);
333 Array<OneD, NekDouble> tmpx(npts), tmpy(npts);
334
335 // Determine 3D manifold orientation
336 m_xmap->BwdTrans(m_coeffs[m_manifold[0]], ptsx);
337 m_xmap->BwdTrans(m_coeffs[m_manifold[1]], ptsy);
338
339 Array<OneD, NekDouble> eta(2, 0.);
340 m_xmap->LocCoordToLocCollapsed(Lcoords, eta);
341 ClampLocCoords(eta, 0.);
342
343 m_xmap->LocCollapsedToLocCoord(eta, Lcoords);
344
345 // Perform newton iteration to find local coordinates
346 NewtonIterationForLocCoord(tmpcoords, ptsx, ptsy, Lcoords, dist);
347 }
348 if (m_coordim == 3)
349 {
350 Array<OneD, NekDouble> eta(2, 0.), xi(2, 0.);
351 m_xmap->LocCoordToLocCollapsed(Lcoords, eta);
352 ClampLocCoords(eta, 0.);
353 m_xmap->LocCollapsedToLocCoord(eta, xi);
354 int npts = m_xmap->GetTotPoints();
355 Array<OneD, NekDouble> ptsz(npts);
356 m_xmap->BwdTrans(m_coeffs[m_manifold[2]], ptsz);
357 NekDouble z = m_xmap->PhysEvaluate(xi, ptsz) - coords[m_manifold[2]];
358 if (GetMetricInfo()->GetGtype() == eDeformed)
359 {
360 dist = sqrt(z * z + dist * dist);
361 }
362 else
363 {
364 dist = fabs(z);
365 }
366 }
367 return dist;
368}
369
371{
372 return m_verts.size();
373}
374
376{
377 return m_edges.size();
378}
379
381{
382 ASSERTL2(i >= 0 && i < m_verts.size(), "Index out of range");
383 return m_verts[i];
384}
385
387{
388 ASSERTL2(i >= 0 && i < m_edges.size(), "Index out of range");
389 return m_edges[i];
390}
391
393{
394 ASSERTL2(i >= 0 && i < m_eorient.size(), "Index out of range");
395 return m_eorient[i];
396}
397
399{
400 return 2;
401}
402
405{
406 if (m_geomFactors->GetGtype() == eRegular)
407 {
408 xiOut = Array<OneD, NekDouble>(2, 0.0);
409
410 GetLocCoords(xs, xiOut);
411 ClampLocCoords(xiOut);
412
413 Array<OneD, NekDouble> gloCoord(3);
414 gloCoord[0] = GetCoord(0, xiOut);
415 gloCoord[1] = GetCoord(1, xiOut);
416 gloCoord[2] = GetCoord(2, xiOut);
417
418 return sqrt((xs[0] - gloCoord[0]) * (xs[0] - gloCoord[0]) +
419 (xs[1] - gloCoord[1]) * (xs[1] - gloCoord[1]) +
420 (xs[2] - gloCoord[2]) * (xs[2] - gloCoord[2]));
421 }
422 // If deformed edge then the inverse mapping is non-linear so need to
423 // numerically solve for the local coordinate
424 else if (m_geomFactors->GetGtype() == eDeformed)
425 {
426 // Choose starting based on closest quad
427 Array<OneD, NekDouble> xi(2, 0.0), eta(2, 0.0);
428 m_xmap->LocCollapsedToLocCoord(eta, xi);
429
430 // Armijo constants:
431 // https://en.wikipedia.org/wiki/Backtracking_line_search
432 const NekDouble c1 = 1e-4, c2 = 0.9;
433
434 int nq = m_xmap->GetTotPoints();
435
436 Array<OneD, NekDouble> x(nq), y(nq), z(nq);
437 m_xmap->BwdTrans(m_coeffs[0], x);
438 m_xmap->BwdTrans(m_coeffs[1], y);
439 m_xmap->BwdTrans(m_coeffs[2], z);
440
441 Array<OneD, NekDouble> xderxi1(nq, 0.0), yderxi1(nq, 0.0),
442 zderxi1(nq, 0.0), xderxi2(nq, 0.0), yderxi2(nq, 0.0),
443 zderxi2(nq, 0.0), xderxi1xi1(nq, 0.0), yderxi1xi1(nq, 0.0),
444 zderxi1xi1(nq, 0.0), xderxi1xi2(nq, 0.0), yderxi1xi2(nq, 0.0),
445 zderxi1xi2(nq, 0.0), xderxi2xi1(nq, 0.0), yderxi2xi1(nq, 0.0),
446 zderxi2xi1(nq, 0.0), xderxi2xi2(nq, 0.0), yderxi2xi2(nq, 0.0),
447 zderxi2xi2(nq, 0.0);
448
449 // Get first & second derivatives & partial derivatives of x,y,z values
450 std::array<NekDouble, 3> xc_derxi, yc_derxi, zc_derxi;
451
452 m_xmap->PhysDeriv(x, xderxi1, xderxi2);
453 m_xmap->PhysDeriv(y, yderxi1, yderxi2);
454 m_xmap->PhysDeriv(z, zderxi1, zderxi2);
455
456 m_xmap->PhysDeriv(xderxi1, xderxi1xi1, xderxi1xi2);
457 m_xmap->PhysDeriv(yderxi1, yderxi1xi1, yderxi1xi2);
458 m_xmap->PhysDeriv(zderxi1, zderxi1xi1, zderxi1xi2);
459
460 m_xmap->PhysDeriv(yderxi2, yderxi2xi1, yderxi2xi2);
461 m_xmap->PhysDeriv(xderxi2, xderxi2xi1, xderxi2xi2);
462 m_xmap->PhysDeriv(zderxi2, zderxi2xi1, zderxi2xi2);
463
464 // Minimisation loop (Quasi-newton method)
465 NekDouble fx_prev = std::numeric_limits<NekDouble>::max();
466 for (int i = 0; i < NekConstants::kNewtonIterations; ++i)
467 {
468 // Compute the objective function, f(x_k) and its derivatives
469 NekDouble xc = m_xmap->PhysEvaluate(xi, x, xc_derxi);
470 NekDouble yc = m_xmap->PhysEvaluate(xi, y, yc_derxi);
471 NekDouble zc = m_xmap->PhysEvaluate(xi, z, zc_derxi);
472
473 NekDouble xc_derxi1xi1 = m_xmap->PhysEvaluate(xi, xderxi1xi1);
474 NekDouble yc_derxi1xi1 = m_xmap->PhysEvaluate(xi, yderxi1xi1);
475 NekDouble zc_derxi1xi1 = m_xmap->PhysEvaluate(xi, zderxi1xi1);
476
477 NekDouble xc_derxi1xi2 = m_xmap->PhysEvaluate(xi, xderxi1xi2);
478 NekDouble yc_derxi1xi2 = m_xmap->PhysEvaluate(xi, yderxi1xi2);
479 NekDouble zc_derxi1xi2 = m_xmap->PhysEvaluate(xi, zderxi1xi2);
480
481 NekDouble xc_derxi2xi2 = m_xmap->PhysEvaluate(xi, xderxi2xi2);
482 NekDouble yc_derxi2xi2 = m_xmap->PhysEvaluate(xi, yderxi2xi2);
483 NekDouble zc_derxi2xi2 = m_xmap->PhysEvaluate(xi, zderxi2xi2);
484
485 // Objective function is the distance to the search point
486 NekDouble xdiff = xc - xs[0];
487 NekDouble ydiff = yc - xs[1];
488 NekDouble zdiff = zc - xs[2];
489
490 NekDouble fx = xdiff * xdiff + ydiff * ydiff + zdiff * zdiff;
491
492 NekDouble fx_derxi1 = 2.0 * xdiff * xc_derxi[0] +
493 2.0 * ydiff * yc_derxi[0] +
494 2.0 * zdiff * zc_derxi[0];
495
496 NekDouble fx_derxi2 = 2.0 * xdiff * xc_derxi[1] +
497 2.0 * ydiff * yc_derxi[1] +
498 2.0 * zdiff * zc_derxi[1];
499
500 NekDouble fx_derxi1xi1 =
501 2.0 * xdiff * xc_derxi1xi1 + 2.0 * xc_derxi[0] * xc_derxi[0] +
502 2.0 * ydiff * yc_derxi1xi1 + 2.0 * yc_derxi[0] * yc_derxi[0] +
503 2.0 * zdiff * zc_derxi1xi1 + 2.0 * zc_derxi[0] * zc_derxi[0];
504
505 NekDouble fx_derxi1xi2 =
506 2.0 * xdiff * xc_derxi1xi2 + 2.0 * xc_derxi[1] * xc_derxi[0] +
507 2.0 * ydiff * yc_derxi1xi2 + 2.0 * yc_derxi[1] * yc_derxi[0] +
508 2.0 * zdiff * zc_derxi1xi2 + 2.0 * zc_derxi[1] * zc_derxi[0];
509
510 NekDouble fx_derxi2xi2 =
511 2.0 * xdiff * xc_derxi2xi2 + 2.0 * xc_derxi[1] * xc_derxi[1] +
512 2.0 * ydiff * yc_derxi2xi2 + 2.0 * yc_derxi[1] * yc_derxi[1] +
513 2.0 * zdiff * zc_derxi2xi2 + 2.0 * zc_derxi[1] * zc_derxi[1];
514
515 // Jacobian
516 NekDouble jac[2];
517 jac[0] = fx_derxi1;
518 jac[1] = fx_derxi2;
519
520 // Inverse of 2x2 hessian
521 NekDouble hessInv[2][2];
522
523 NekDouble det =
524 1 / (fx_derxi1xi1 * fx_derxi2xi2 - fx_derxi1xi2 * fx_derxi1xi2);
525 hessInv[0][0] = det * fx_derxi2xi2;
526 hessInv[0][1] = det * -fx_derxi1xi2;
527 hessInv[1][0] = det * -fx_derxi1xi2;
528 hessInv[1][1] = det * fx_derxi1xi1;
529
530 // Check for convergence
531 if (abs(fx - fx_prev) < 1e-12)
532 {
533 fx_prev = fx;
534 break;
535 }
536 else
537 {
538 fx_prev = fx;
539 }
540
541 NekDouble gamma = 1.0;
542 bool conv = false;
543
544 // Search direction: Newton's method
545 NekDouble pk[2];
546 pk[0] = -(hessInv[0][0] * jac[0] + hessInv[1][0] * jac[1]);
547 pk[1] = -(hessInv[0][1] * jac[0] + hessInv[1][1] * jac[1]);
548
549 // Backtracking line search
550 while (gamma > 1e-10)
551 {
552 Array<OneD, NekDouble> xi_pk(2);
553 xi_pk[0] = xi[0] + pk[0] * gamma;
554 xi_pk[1] = xi[1] + pk[1] * gamma;
555
556 Array<OneD, NekDouble> eta_pk(2, 0.0);
557 m_xmap->LocCoordToLocCollapsed(xi_pk, eta_pk);
558
559 if (eta_pk[0] <
560 (-1 - std::numeric_limits<NekDouble>::epsilon()) ||
561 eta_pk[0] >
562 (1 + std::numeric_limits<NekDouble>::epsilon()) ||
563 eta_pk[1] <
564 (-1 - std::numeric_limits<NekDouble>::epsilon()) ||
565 eta_pk[1] > (1 + std::numeric_limits<NekDouble>::epsilon()))
566 {
567 gamma /= 2.0;
568 continue;
569 }
570
571 std::array<NekDouble, 3> xc_pk_derxi, yc_pk_derxi, zc_pk_derxi;
572
573 NekDouble xc_pk = m_xmap->PhysEvaluate(xi_pk, x, xc_pk_derxi);
574 NekDouble yc_pk = m_xmap->PhysEvaluate(xi_pk, y, yc_pk_derxi);
575 NekDouble zc_pk = m_xmap->PhysEvaluate(xi_pk, z, zc_pk_derxi);
576
577 NekDouble xc_pk_diff = xc_pk - xs[0];
578 NekDouble yc_pk_diff = yc_pk - xs[1];
579 NekDouble zc_pk_diff = zc_pk - xs[2];
580
581 NekDouble fx_pk = xc_pk_diff * xc_pk_diff +
582 yc_pk_diff * yc_pk_diff +
583 zc_pk_diff * zc_pk_diff;
584
585 NekDouble fx_pk_derxi1 = 2.0 * xc_pk_diff * xc_pk_derxi[0] +
586 2.0 * yc_pk_diff * yc_pk_derxi[0] +
587 2.0 * zc_pk_diff * zc_pk_derxi[0];
588
589 NekDouble fx_pk_derxi2 = 2.0 * xc_pk_diff * xc_pk_derxi[1] +
590 2.0 * yc_pk_diff * yc_pk_derxi[1] +
591 2.0 * zc_pk_diff * zc_pk_derxi[1];
592
593 // Check Wolfe conditions using Armijo constants
594 // https://en.wikipedia.org/wiki/Wolfe_conditions
595 NekDouble tmp = pk[0] * fx_derxi1 + pk[1] * fx_derxi2;
596 NekDouble tmp2 = pk[0] * fx_pk_derxi1 + pk[1] * fx_pk_derxi2;
597 if ((fx_pk - (fx + c1 * gamma * tmp)) <
598 std::numeric_limits<NekDouble>::epsilon() &&
599 (-tmp2 - (-c2 * tmp)) <
600 std::numeric_limits<NekDouble>::epsilon())
601 {
602 conv = true;
603 break;
604 }
605
606 gamma /= 2.0;
607 }
608
609 if (!conv)
610 {
611 break;
612 }
613
614 xi[0] += gamma * pk[0];
615 xi[1] += gamma * pk[1];
616 }
617
618 xiOut = xi;
619 return sqrt(fx_prev);
620 }
621 else
622 {
623 ASSERTL0(false, "Geometry type unknown")
624 }
625
626 return -1.0;
627}
628
630{
631 NekDouble Jac = m_isoParameter[0][1] * m_isoParameter[1][2] -
632 m_isoParameter[1][1] * m_isoParameter[0][2];
633 Jac = 1. / Jac;
634 // a12, -a02, -a11, a01
638 m_invIsoParam[0][0] = m_isoParameter[1][2] * Jac;
639 m_invIsoParam[0][1] = -m_isoParameter[0][2] * Jac;
640 m_invIsoParam[1][0] = -m_isoParameter[1][1] * Jac;
641 m_invIsoParam[1][1] = m_isoParameter[0][1] * Jac;
642}
643
644} // namespace Nektar::SpatialDomains
#define WARNINGL1(condition, msg)
Definition: ErrorUtil.hpp:243
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed to...
Definition: ErrorUtil.hpp:265
int v_GetShapeDim() const override
Get the object's shape dimension.
Definition: Geometry2D.cpp:398
void NewtonIterationForLocCoord(const Array< OneD, const NekDouble > &coords, const Array< OneD, const NekDouble > &ptsx, const Array< OneD, const NekDouble > &ptsy, Array< OneD, NekDouble > &Lcoords, NekDouble &dist)
Definition: Geometry2D.cpp:170
Array< OneD, int > m_manifold
Definition: Geometry2D.h:83
Array< OneD, Array< OneD, NekDouble > > m_edgeNormal
Definition: Geometry2D.h:84
void SolveStraightEdgeQuad(const Array< OneD, const NekDouble > &coords, Array< OneD, NekDouble > &Lcoords)
Definition: Geometry2D.cpp:143
PointGeomSharedPtr v_GetVertex(int i) const override
Definition: Geometry2D.cpp:380
Geometry1DSharedPtr v_GetEdge(int i) const override
Returns edge i of this object.
Definition: Geometry2D.cpp:386
std::vector< StdRegions::Orientation > m_eorient
Definition: Geometry2D.h:81
NekDouble v_GetLocCoords(const Array< OneD, const NekDouble > &coords, Array< OneD, NekDouble > &Lcoords) override
Determine the local collapsed coordinates that correspond to a given Cartesian coordinate for this ge...
Definition: Geometry2D.cpp:307
int v_AllLeftCheck(const Array< OneD, const NekDouble > &gloCoord) override
Definition: Geometry2D.cpp:57
void v_CalculateInverseIsoParam() override
Definition: Geometry2D.cpp:629
StdRegions::Orientation v_GetEorient(const int i) const override
Returns the orientation of edge i with respect to the ordering of edges in the standard element.
Definition: Geometry2D.cpp:392
NekDouble v_FindDistance(const Array< OneD, const NekDouble > &xs, Array< OneD, NekDouble > &xi) override
Definition: Geometry2D.cpp:403
int v_GetNumEdges() const override
Get the number of edges of this object.
Definition: Geometry2D.cpp:375
int v_GetNumVerts() const override
Get the number of vertices of this object.
Definition: Geometry2D.cpp:370
Base class for shape geometry information.
Definition: Geometry.h:79
NekDouble GetCoord(const int i, const Array< OneD, const NekDouble > &Lcoord)
Given local collapsed coordinate Lcoord, return the value of physical coordinate in direction i.
Definition: Geometry.h:554
NekDouble GetLocCoords(const Array< OneD, const NekDouble > &coords, Array< OneD, NekDouble > &Lcoords)
Determine the local collapsed coordinates that correspond to a given Cartesian coordinate for this ge...
Definition: Geometry.h:544
virtual void v_FillGeom()
Populate the coordinate mapping Geometry::m_coeffs information from any children geometry elements.
Definition: Geometry.cpp:355
int GetVid(int i) const
Get the ID of vertex i of this object.
Definition: Geometry.cpp:104
Array< OneD, Array< OneD, NekDouble > > m_isoParameter
Definition: Geometry.h:209
Array< OneD, Array< OneD, NekDouble > > m_invIsoParam
Definition: Geometry.h:210
bool ClampLocCoords(Array< OneD, NekDouble > &locCoord, NekDouble tol=std::numeric_limits< NekDouble >::epsilon())
Clamp local coords to be within standard regions [-1, 1]^dim.
Definition: Geometry.cpp:530
int GetGlobalID(void) const
Get the ID of this object.
Definition: Geometry.h:326
Array< OneD, Array< OneD, NekDouble > > m_coeffs
Array containing expansion coefficients of m_xmap.
Definition: Geometry.h:206
GeomFactorsSharedPtr GetMetricInfo()
Get the geometric factors for this object.
Definition: Geometry.h:310
StdRegions::StdExpansionSharedPtr m_xmap
mapping containing isoparametric transformation.
Definition: Geometry.h:194
StdRegions::StdExpansionSharedPtr GetXmap() const
Return the mapping object Geometry::m_xmap that represents the coordinate transformation from standar...
Definition: Geometry.h:435
GeomFactorsSharedPtr m_geomFactors
Geometric factors.
Definition: Geometry.h:190
int m_coordim
Coordinate dimension of this geometry object.
Definition: Geometry.h:188
@ beta
Gauss Radau pinned at x=-1,.
Definition: PointsType.h:59
static const unsigned int kNewtonIterations
std::shared_ptr< Curve > CurveSharedPtr
Definition: Curve.hpp:58
@ eRegular
Geometry is straight-sided with constant geometric factors.
@ eDeformed
Geometry is curved or has non-constant factors.
std::shared_ptr< PointGeom > PointGeomSharedPtr
Definition: Geometry.h:57
std::shared_ptr< Geometry1D > Geometry1DSharedPtr
Definition: Geometry.h:61
std::vector< double > z(NPUPPER)
double NekDouble
T Vsum(int n, const T *x, const int incx)
Subtract return sum(x)
Definition: Vmath.hpp:608
scalarT< T > abs(scalarT< T > in)
Definition: scalar.hpp:289
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:285