Nektar++
Geometry3D.cpp
Go to the documentation of this file.
1////////////////////////////////////////////////////////////////////////////////
2//
3// File: Geometry3D.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: 3D geometry information.
32//
33////////////////////////////////////////////////////////////////////////////////
34
35#include <iomanip>
36
37#include <SpatialDomains/GeomFactors.h>
38#include <SpatialDomains/Geometry2D.h>
39#include <SpatialDomains/Geometry3D.h>
40#include <SpatialDomains/PointGeom.h>
41#include <SpatialDomains/SegGeom.h>
42#include <iomanip>
43
44namespace Nektar::SpatialDomains
45{
46
47Geometry3D::Geometry3D()
48{
49}
50
51Geometry3D::Geometry3D(const int coordim) : Geometry(coordim)
52{
53 ASSERTL0(m_coordim > 2,
54 "Coordinate dimension should be at least 3 for a 3D geometry.");
55}
56
57//---------------------------------------
58// Helper functions
59//---------------------------------------
60
61//---------------------------------------
62// 3D Geometry Methods
63//---------------------------------------
64void Geometry3D::NewtonIterationForLocCoord(
65 const Array<OneD, const NekDouble> &coords, Array<OneD, NekDouble> &Lcoords)
66{
67 // maximum iterations for convergence
68 const int MaxIterations = 51;
69 // |x-xp|^2 < EPSILON error tolerance
70 const NekDouble Tol = 1.e-8;
71 // |r,s| > LcoordDIV stop the search
72 const NekDouble LcoordDiv = 15.0;
73
74 Array<OneD, const NekDouble> Jac =
75 m_geomFactors->GetJac(m_xmap->GetPointsKeys());
76
77 NekDouble ScaledTol =
78 Vmath::Vsum(Jac.size(), Jac, 1) / ((NekDouble)Jac.size());
79 ScaledTol *= Tol * Tol;
80
81 NekDouble xmap, ymap, zmap, res;
82 int cnt = 0;
83 Array<OneD, NekDouble> var(8, 1.);
84 Array<OneD, NekDouble> deriv(9);
85 Array<OneD, NekDouble> tmp(3);
86 while (cnt++ < MaxIterations)
87 {
88 var[1] = Lcoords[0];
89 var[2] = Lcoords[1];
90 var[3] = Lcoords[2];
91 var[4] = Lcoords[0] * Lcoords[1];
92 var[5] = Lcoords[1] * Lcoords[2];
93 var[6] = Lcoords[0] * Lcoords[2];
94 var[7] = var[4] * Lcoords[2];
95 // calculate the global point corresponding to Lcoords
96 xmap =
97 Vmath::Dot(m_isoParameter[0].size(), var, 1, m_isoParameter[0], 1);
98 ymap =
99 Vmath::Dot(m_isoParameter[0].size(), var, 1, m_isoParameter[1], 1);
100 zmap =
101 Vmath::Dot(m_isoParameter[0].size(), var, 1, m_isoParameter[2], 1);
102
103 tmp[0] = coords[0] - xmap;
104 tmp[1] = coords[1] - ymap;
105 tmp[2] = coords[2] - zmap;
106
107 res = tmp[0] * tmp[0] + tmp[1] * tmp[1] + tmp[2] * tmp[2];
108 if (res < ScaledTol)
109 {
110 break;
111 }
112
113 // Interpolate derivative metric at Lcoords (ddx1, ddx2, ddx3)
114 deriv[0] = m_isoParameter[0][1] + m_isoParameter[0][4] * Lcoords[1];
115 deriv[1] = m_isoParameter[0][2] + m_isoParameter[0][4] * Lcoords[0];
116 deriv[2] = m_isoParameter[0][3];
117 deriv[3] = m_isoParameter[1][1] + m_isoParameter[1][4] * Lcoords[1];
118 deriv[4] = m_isoParameter[1][2] + m_isoParameter[1][4] * Lcoords[0];
119 deriv[5] = m_isoParameter[1][3];
120 deriv[6] = m_isoParameter[2][1] + m_isoParameter[2][4] * Lcoords[1];
121 deriv[7] = m_isoParameter[2][2] + m_isoParameter[2][4] * Lcoords[0];
122 deriv[8] = m_isoParameter[2][3];
123 if (m_isoParameter[0].size() >= 6)
124 {
125 deriv[1] += m_isoParameter[0][5] * Lcoords[2];
126 deriv[2] += m_isoParameter[0][5] * Lcoords[1];
127 deriv[4] += m_isoParameter[1][5] * Lcoords[2];
128 deriv[5] += m_isoParameter[1][5] * Lcoords[1];
129 deriv[7] += m_isoParameter[2][5] * Lcoords[2];
130 deriv[8] += m_isoParameter[2][5] * Lcoords[1];
131 }
132 if (m_isoParameter[0].size() >= 8)
133 {
134 deriv[0] += m_isoParameter[0][6] * Lcoords[2] +
135 m_isoParameter[0][7] * var[5];
136 deriv[1] += m_isoParameter[0][7] * var[6];
137 deriv[2] += m_isoParameter[0][6] * Lcoords[0] +
138 m_isoParameter[0][7] * var[4];
139 deriv[3] += m_isoParameter[1][6] * Lcoords[2] +
140 m_isoParameter[1][7] * var[5];
141 deriv[4] += m_isoParameter[1][7] * var[6];
142 deriv[5] += m_isoParameter[1][6] * Lcoords[0] +
143 m_isoParameter[1][7] * var[4];
144 deriv[6] += m_isoParameter[2][6] * Lcoords[2] +
145 m_isoParameter[2][7] * var[5];
146 deriv[7] += m_isoParameter[2][7] * var[6];
147 deriv[8] += m_isoParameter[2][6] * Lcoords[0] +
148 m_isoParameter[2][7] * var[4];
149 }
150 DNekMatSharedPtr mat =
151 MemoryManager<DNekMat>::AllocateSharedPtr(3, 3, 0, eFULL);
152 Vmath::Vcopy(9, deriv, 1, mat->GetPtr(), 1);
153 mat->Invert();
154 Lcoords[0] += Vmath::Dot(3, mat->GetPtr(), 1, tmp, 1);
155 Lcoords[1] += Vmath::Dot(3, mat->GetPtr() + 3, 1, tmp, 1);
156 Lcoords[2] += Vmath::Dot(3, mat->GetPtr() + 6, 1, tmp, 1);
157
158 if (!(std::isfinite(Lcoords[0]) && std::isfinite(Lcoords[1]) &&
159 std::isfinite(Lcoords[2])) ||
160 fabs(Lcoords[0]) > LcoordDiv || fabs(Lcoords[1]) > LcoordDiv ||
161 fabs(Lcoords[2]) > LcoordDiv)
162 {
163 std::ostringstream ss;
164 ss << "Iteration has diverged in NewtonIterationForLocCoord in "
165 "element "
166 << GetGlobalID();
167 WARNINGL1(false, ss.str());
168 return;
169 }
170 }
171
172 if (cnt >= MaxIterations)
173 {
174 std::ostringstream ss;
175
176 ss << "Reached MaxIterations (" << MaxIterations
177 << ") in Newton iteration ";
178
179 WARNINGL1(cnt < MaxIterations, ss.str());
180 }
181}
182
183int Geometry3D::v_AllLeftCheck(
184 [[maybe_unused]] const Array<OneD, const NekDouble> &gloCoord)
185{
186 // todo: check only plane surfaces
187 return 0;
188}
189
190void Geometry3D::NewtonIterationForLocCoord(
191 const Array<OneD, const NekDouble> &coords,
192 const Array<OneD, const NekDouble> &ptsx,
193 const Array<OneD, const NekDouble> &ptsy,
194 const Array<OneD, const NekDouble> &ptsz, Array<OneD, NekDouble> &Lcoords,
195 NekDouble &dist)
196{
197 // maximum iterations for convergence
198 const int MaxIterations = NekConstants::kNewtonIterations;
199 // |x-xp|^2 < EPSILON error tolerance
200 const NekDouble Tol = 1.e-8;
201 // |r,s| > LcoordDIV stop the search
202 const NekDouble LcoordDiv = 15.0;
203
204 Array<OneD, const NekDouble> Jac =
205 m_geomFactors->GetJac(m_xmap->GetPointsKeys());
206
207 NekDouble ScaledTol =
208 Vmath::Vsum(Jac.size(), Jac, 1) / ((NekDouble)Jac.size());
209 ScaledTol *= Tol;
210
211 NekDouble xmap, ymap, zmap, F1, F2, F3;
212
213 NekDouble derx_1, derx_2, derx_3, dery_1, dery_2, dery_3, derz_1, derz_2,
214 derz_3, jac;
215
216 // save intiial guess for later reference if required.
217 NekDouble init0 = Lcoords[0], init1 = Lcoords[1], init2 = Lcoords[2];
218
219 Array<OneD, NekDouble> DxD1(ptsx.size());
220 Array<OneD, NekDouble> DxD2(ptsx.size());
221 Array<OneD, NekDouble> DxD3(ptsx.size());
222 Array<OneD, NekDouble> DyD1(ptsx.size());
223 Array<OneD, NekDouble> DyD2(ptsx.size());
224 Array<OneD, NekDouble> DyD3(ptsx.size());
225 Array<OneD, NekDouble> DzD1(ptsx.size());
226 Array<OneD, NekDouble> DzD2(ptsx.size());
227 Array<OneD, NekDouble> DzD3(ptsx.size());
228
229 // Ideally this will be stored in m_geomfactors
230 m_xmap->PhysDeriv(ptsx, DxD1, DxD2, DxD3);
231 m_xmap->PhysDeriv(ptsy, DyD1, DyD2, DyD3);
232 m_xmap->PhysDeriv(ptsz, DzD1, DzD2, DzD3);
233
234 int cnt = 0;
235 Array<OneD, DNekMatSharedPtr> I(3);
236 Array<OneD, NekDouble> eta(3);
237
238 F1 = F2 = F3 = 2000; // Starting value of Function
239 NekDouble resid = sqrt(F1 * F1 + F2 * F2 + F3 * F3);
240 while (cnt++ < MaxIterations)
241 {
242 // evaluate lagrange interpolant at Lcoords
243 m_xmap->LocCoordToLocCollapsed(Lcoords, eta);
244 I[0] = m_xmap->GetBasis(0)->GetI(eta);
245 I[1] = m_xmap->GetBasis(1)->GetI(eta + 1);
246 I[2] = m_xmap->GetBasis(2)->GetI(eta + 2);
247
248 // calculate the global point `corresponding to Lcoords
249 xmap = m_xmap->PhysEvaluate(I, ptsx);
250 ymap = m_xmap->PhysEvaluate(I, ptsy);
251 zmap = m_xmap->PhysEvaluate(I, ptsz);
252
253 F1 = coords[0] - xmap;
254 F2 = coords[1] - ymap;
255 F3 = coords[2] - zmap;
256
257 if (F1 * F1 + F2 * F2 + F3 * F3 < ScaledTol)
258 {
259 resid = sqrt(F1 * F1 + F2 * F2 + F3 * F3);
260 break;
261 }
262
263 // Interpolate derivative metric at Lcoords
264 derx_1 = m_xmap->PhysEvaluate(I, DxD1);
265 derx_2 = m_xmap->PhysEvaluate(I, DxD2);
266 derx_3 = m_xmap->PhysEvaluate(I, DxD3);
267 dery_1 = m_xmap->PhysEvaluate(I, DyD1);
268 dery_2 = m_xmap->PhysEvaluate(I, DyD2);
269 dery_3 = m_xmap->PhysEvaluate(I, DyD3);
270 derz_1 = m_xmap->PhysEvaluate(I, DzD1);
271 derz_2 = m_xmap->PhysEvaluate(I, DzD2);
272 derz_3 = m_xmap->PhysEvaluate(I, DzD3);
273
274 jac = derx_1 * (dery_2 * derz_3 - dery_3 * derz_2) -
275 derx_2 * (dery_1 * derz_3 - dery_3 * derz_1) +
276 derx_3 * (dery_1 * derz_2 - dery_2 * derz_1);
277
278 // use analytical inverse of derivitives which are also similar to
279 // those of metric factors.
280 Lcoords[0] =
281 Lcoords[0] +
282 ((dery_2 * derz_3 - dery_3 * derz_2) * (coords[0] - xmap) -
283 (derx_2 * derz_3 - derx_3 * derz_2) * (coords[1] - ymap) +
284 (derx_2 * dery_3 - derx_3 * dery_2) * (coords[2] - zmap)) /
285 jac;
286
287 Lcoords[1] =
288 Lcoords[1] -
289 ((dery_1 * derz_3 - dery_3 * derz_1) * (coords[0] - xmap) -
290 (derx_1 * derz_3 - derx_3 * derz_1) * (coords[1] - ymap) +
291 (derx_1 * dery_3 - derx_3 * dery_1) * (coords[2] - zmap)) /
292 jac;
293
294 Lcoords[2] =
295 Lcoords[2] +
296 ((dery_1 * derz_2 - dery_2 * derz_1) * (coords[0] - xmap) -
297 (derx_1 * derz_2 - derx_2 * derz_1) * (coords[1] - ymap) +
298 (derx_1 * dery_2 - derx_2 * dery_1) * (coords[2] - zmap)) /
299 jac;
300
301 if (!(std::isfinite(Lcoords[0]) && std::isfinite(Lcoords[1]) &&
302 std::isfinite(Lcoords[2])))
303 {
304 dist = 1e16;
305 std::ostringstream ss;
306 ss << "nan or inf found in NewtonIterationForLocCoord in element "
307 << GetGlobalID();
308 WARNINGL1(false, ss.str());
309 return;
310 }
311 if (fabs(Lcoords[0]) > LcoordDiv || fabs(Lcoords[1]) > LcoordDiv ||
312 fabs(Lcoords[2]) > LcoordDiv)
313 {
314 break; // lcoords have diverged so stop iteration
315 }
316 }
317
318 m_xmap->LocCoordToLocCollapsed(Lcoords, eta);
319 if (ClampLocCoords(eta, 0.))
320 {
321 I[0] = m_xmap->GetBasis(0)->GetI(eta);
322 I[1] = m_xmap->GetBasis(1)->GetI(eta + 1);
323 I[2] = m_xmap->GetBasis(2)->GetI(eta + 2);
324 // calculate the global point corresponding to Lcoords
325 xmap = m_xmap->PhysEvaluate(I, ptsx);
326 ymap = m_xmap->PhysEvaluate(I, ptsy);
327 zmap = m_xmap->PhysEvaluate(I, ptsz);
328 F1 = coords[0] - xmap;
329 F2 = coords[1] - ymap;
330 F3 = coords[2] - zmap;
331 dist = sqrt(F1 * F1 + F2 * F2 + F3 * F3);
332 }
333 else
334 {
335 dist = 0.;
336 }
337
338 if (cnt >= MaxIterations)
339 {
340 Array<OneD, NekDouble> collCoords(3);
341 m_xmap->LocCoordToLocCollapsed(Lcoords, collCoords);
342
343 // if coordinate is inside element dump error!
344 if ((collCoords[0] >= -1.0 && collCoords[0] <= 1.0) &&
345 (collCoords[1] >= -1.0 && collCoords[1] <= 1.0) &&
346 (collCoords[2] >= -1.0 && collCoords[2] <= 1.0))
347 {
348 std::ostringstream ss;
349
350 ss << "Reached MaxIterations (" << MaxIterations
351 << ") in Newton iteration ";
352 ss << "Init value (" << std::setprecision(4) << init0 << ","
353 << init1 << "," << init2 << ") ";
354 ss << "Fin value (" << Lcoords[0] << "," << Lcoords[1] << ","
355 << Lcoords[2] << ") ";
356 ss << "Resid = " << resid
357 << " Tolerance = " << std::sqrt(ScaledTol);
358
359 WARNINGL1(cnt < MaxIterations, ss.str());
360 }
361 }
362}
363
364NekDouble Geometry3D::v_GetLocCoords(const Array<OneD, const NekDouble> &coords,
365 Array<OneD, NekDouble> &Lcoords)
366{
367 NekDouble dist = std::numeric_limits<double>::max();
368 Array<OneD, NekDouble> tmpcoords(3);
369 tmpcoords[0] = coords[0];
370 tmpcoords[1] = coords[1];
371 tmpcoords[2] = coords[2];
372 if (GetMetricInfo()->GetGtype() == eRegular)
373 {
374 tmpcoords[0] -= m_isoParameter[0][0];
375 tmpcoords[1] -= m_isoParameter[1][0];
376 tmpcoords[2] -= m_isoParameter[2][0];
377 Lcoords[0] = m_invIsoParam[0][0] * tmpcoords[0] +
378 m_invIsoParam[0][1] * tmpcoords[1] +
379 m_invIsoParam[0][2] * tmpcoords[2];
380 Lcoords[1] = m_invIsoParam[1][0] * tmpcoords[0] +
381 m_invIsoParam[1][1] * tmpcoords[1] +
382 m_invIsoParam[1][2] * tmpcoords[2];
383 Lcoords[2] = m_invIsoParam[2][0] * tmpcoords[0] +
384 m_invIsoParam[2][1] * tmpcoords[1] +
385 m_invIsoParam[2][2] * tmpcoords[2];
386 }
387 else if (m_straightEdge)
388 {
389 ClampLocCoords(Lcoords, 0.);
390 NewtonIterationForLocCoord(tmpcoords, Lcoords);
391 }
392 else
393 {
394 v_FillGeom();
395 // Determine nearest point of coords to values in m_xmap
396 int npts = m_xmap->GetTotPoints();
397 Array<OneD, NekDouble> ptsx(npts), ptsy(npts), ptsz(npts);
398 Array<OneD, NekDouble> tmp1(npts), tmp2(npts);
399
400 m_xmap->BwdTrans(m_coeffs[0], ptsx);
401 m_xmap->BwdTrans(m_coeffs[1], ptsy);
402 m_xmap->BwdTrans(m_coeffs[2], ptsz);
403
404 const Array<OneD, const NekDouble> za = m_xmap->GetPoints(0);
405 const Array<OneD, const NekDouble> zb = m_xmap->GetPoints(1);
406 const Array<OneD, const NekDouble> zc = m_xmap->GetPoints(2);
407
408 // guess the first local coords based on nearest point
409 Vmath::Sadd(npts, -coords[0], ptsx, 1, tmp1, 1);
410 Vmath::Vmul(npts, tmp1, 1, tmp1, 1, tmp1, 1);
411 Vmath::Sadd(npts, -coords[1], ptsy, 1, tmp2, 1);
412 Vmath::Vvtvp(npts, tmp2, 1, tmp2, 1, tmp1, 1, tmp1, 1);
413 Vmath::Sadd(npts, -coords[2], ptsz, 1, tmp2, 1);
414 Vmath::Vvtvp(npts, tmp2, 1, tmp2, 1, tmp1, 1, tmp1, 1);
415
416 int min_i = Vmath::Imin(npts, tmp1, 1);
417
418 // Get Local coordinates
419 int qa = za.size(), qb = zb.size();
420
421 Array<OneD, NekDouble> eta(3, 0.);
422 eta[2] = zc[min_i / (qa * qb)];
423 min_i = min_i % (qa * qb);
424 eta[1] = zb[min_i / qa];
425 eta[0] = za[min_i % qa];
426 m_xmap->LocCollapsedToLocCoord(eta, Lcoords);
427
428 // Perform newton iteration to find local coordinates
429 NewtonIterationForLocCoord(coords, ptsx, ptsy, ptsz, Lcoords, dist);
430 }
431 return dist;
432}
433
434/**
435 * @brief Put all quadrature information into face/edge structure and
436 * backward transform.
437 *
438 * Note verts, edges, and faces are listed according to anticlockwise
439 * convention but points in _coeffs have to be in array format from left
440 * to right.
441 */
442void Geometry3D::v_FillGeom()
443{
444 if (m_state == ePtsFilled)
445 {
446 return;
447 }
448
449 int i, j, k;
450
451 for (i = 0; i < m_forient.size(); i++)
452 {
453 m_faces[i]->FillGeom();
454
455 int nFaceCoeffs = m_faces[i]->GetXmap()->GetNcoeffs();
456
457 Array<OneD, unsigned int> mapArray(nFaceCoeffs);
458 Array<OneD, int> signArray(nFaceCoeffs);
459
460 if (m_forient[i] < 9)
461 {
462 m_xmap->GetTraceToElementMap(
463 i, mapArray, signArray, m_forient[i],
464 m_faces[i]->GetXmap()->GetTraceNcoeffs(0),
465 m_faces[i]->GetXmap()->GetTraceNcoeffs(1));
466 }
467 else
468 {
469 m_xmap->GetTraceToElementMap(
470 i, mapArray, signArray, m_forient[i],
471 m_faces[i]->GetXmap()->GetTraceNcoeffs(1),
472 m_faces[i]->GetXmap()->GetTraceNcoeffs(0));
473 }
474
475 for (j = 0; j < m_coordim; j++)
476 {
477 const Array<OneD, const NekDouble> &coeffs =
478 m_faces[i]->GetCoeffs(j);
479
480 for (k = 0; k < nFaceCoeffs; k++)
481 {
482 NekDouble v = signArray[k] * coeffs[k];
483 m_coeffs[j][mapArray[k]] = v;
484 }
485 }
486 }
487
488 m_state = ePtsFilled;
489}
490
491/**
492 * @brief Given local collapsed coordinate Lcoord return the value of
493 * physical coordinate in direction i.
494 */
495NekDouble Geometry3D::v_GetCoord(const int i,
496 const Array<OneD, const NekDouble> &Lcoord)
497{
498 ASSERTL1(m_state == ePtsFilled, "Geometry is not in physical space");
499
500 Array<OneD, NekDouble> tmp(m_xmap->GetTotPoints());
501 m_xmap->BwdTrans(m_coeffs[i], tmp);
502
503 return m_xmap->PhysEvaluate(Lcoord, tmp);
504}
505
506//---------------------------------------
507// Helper functions
508//---------------------------------------
509
510int Geometry3D::v_GetShapeDim() const
511{
512 return 3;
513}
514
515int Geometry3D::v_GetNumVerts() const
516{
517 return m_verts.size();
518}
519
520int Geometry3D::v_GetNumEdges() const
521{
522 return m_edges.size();
523}
524
525int Geometry3D::v_GetNumFaces() const
526{
527 return m_faces.size();
528}
529
530PointGeomSharedPtr Geometry3D::v_GetVertex(int i) const
531{
532 return m_verts[i];
533}
534
535Geometry1DSharedPtr Geometry3D::v_GetEdge(int i) const
536{
537 ASSERTL2(i >= 0 && i <= m_edges.size() - 1,
538 "Edge ID must be between 0 and " +
539 std::to_string(m_edges.size() - 1));
540 return m_edges[i];
541}
542
543Geometry2DSharedPtr Geometry3D::v_GetFace(int i) const
544{
545 ASSERTL2((i >= 0) && (i <= 5), "Edge id must be between 0 and 4");
546 return m_faces[i];
547}
548
549inline StdRegions::Orientation Geometry3D::v_GetEorient(const int i) const
550{
551 ASSERTL2(i >= 0 && i <= m_edges.size() - 1,
552 "Edge ID must be between 0 and " +
553 std::to_string(m_edges.size() - 1));
554 return m_eorient[i];
555}
556
557StdRegions::Orientation Geometry3D::v_GetForient(const int i) const
558{
559 ASSERTL2(i >= 0 && i <= m_faces.size() - 1,
560 "Face ID must be between 0 and " +
561 std::to_string(m_faces.size() - 1));
562 return m_forient[i];
563}
564
565void Geometry3D::v_CalculateInverseIsoParam()
566{
567 DNekMatSharedPtr mat =
568 MemoryManager<DNekMat>::AllocateSharedPtr(3, 3, 0, eFULL);
569 for (int i = 0; i < 3; ++i)
570 {
571 for (int j = 0; j < 3; ++j)
572 {
573
574 mat->SetValue(i, j, m_isoParameter[i][j + 1]);
575 }
576 }
577 mat->Invert();
578 m_invIsoParam = Array<OneD, Array<OneD, NekDouble>>(3);
579 for (int i = 0; i < 3; ++i)
580 {
581 m_invIsoParam[i] = Array<OneD, NekDouble>(3);
582 for (int j = 0; j < 3; ++j)
583 {
584 m_invIsoParam[i][j] = mat->GetValue(i, j);
585 }
586 }
587}
588
589} // namespace Nektar::SpatialDomains
WARNINGL1
#define WARNINGL1(condition, msg)
Definition: ErrorUtil.hpp:243
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
ASSERTL1
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:242
ASSERTL2
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed to...
Definition: ErrorUtil.hpp:265
Nektar::MemoryManager::AllocateSharedPtr
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
Definition: NekMemoryManager.hpp:166
Nektar::SpatialDomains::Geometry3D::v_GetFace
Geometry2DSharedPtr v_GetFace(int i) const override
Returns face i of this object.
Definition: Geometry3D.cpp:543
Nektar::SpatialDomains::Geometry3D::v_GetNumFaces
int v_GetNumFaces() const override
Get the number of faces of this object.
Definition: Geometry3D.cpp:525
Nektar::SpatialDomains::Geometry3D::v_GetLocCoords
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: Geometry3D.cpp:364
Nektar::SpatialDomains::Geometry3D::v_AllLeftCheck
int v_AllLeftCheck(const Array< OneD, const NekDouble > &gloCoord) override
Definition: Geometry3D.cpp:183
Nektar::SpatialDomains::Geometry3D::v_GetNumEdges
int v_GetNumEdges() const override
Get the number of edges of this object.
Definition: Geometry3D.cpp:520
Nektar::SpatialDomains::Geometry3D::m_forient
std::vector< StdRegions::Orientation > m_forient
Definition: Geometry3D.h:81
Nektar::SpatialDomains::Geometry3D::v_CalculateInverseIsoParam
void v_CalculateInverseIsoParam() override
Definition: Geometry3D.cpp:565
Nektar::SpatialDomains::Geometry3D::v_FillGeom
void v_FillGeom() override
Put all quadrature information into face/edge structure and backward transform.
Definition: Geometry3D.cpp:442
Nektar::SpatialDomains::Geometry3D::v_GetCoord
NekDouble v_GetCoord(const int i, const Array< OneD, const NekDouble > &Lcoord) override
Given local collapsed coordinate Lcoord return the value of physical coordinate in direction i.
Definition: Geometry3D.cpp:495
Nektar::SpatialDomains::Geometry3D::v_GetVertex
PointGeomSharedPtr v_GetVertex(int i) const override
Definition: Geometry3D.cpp:530
Nektar::SpatialDomains::Geometry3D::NewtonIterationForLocCoord
void NewtonIterationForLocCoord(const Array< OneD, const NekDouble > &coords, const Array< OneD, const NekDouble > &ptsx, const Array< OneD, const NekDouble > &ptsy, const Array< OneD, const NekDouble > &ptsz, Array< OneD, NekDouble > &Lcoords, NekDouble &dist)
Definition: Geometry3D.cpp:190
Nektar::SpatialDomains::Geometry3D::v_GetEdge
Geometry1DSharedPtr v_GetEdge(int i) const override
Returns edge i of this object.
Definition: Geometry3D.cpp:535
Nektar::SpatialDomains::Geometry3D::v_GetForient
StdRegions::Orientation v_GetForient(const int i) const override
Returns the orientation of face i with respect to the ordering of faces in the standard element.
Definition: Geometry3D.cpp:557
Nektar::SpatialDomains::Geometry3D::v_GetShapeDim
int v_GetShapeDim() const override
Get the object's shape dimension.
Definition: Geometry3D.cpp:510
Nektar::SpatialDomains::Geometry3D::v_GetEorient
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: Geometry3D.cpp:549
Nektar::SpatialDomains::Geometry3D::m_eorient
std::vector< StdRegions::Orientation > m_eorient
Definition: Geometry3D.h:80
Nektar::SpatialDomains::Geometry3D::v_GetNumVerts
int v_GetNumVerts() const override
Get the number of vertices of this object.
Definition: Geometry3D.cpp:515
Nektar::SpatialDomains::Geometry
Base class for shape geometry information.
Definition: Geometry.h:79
Nektar::SpatialDomains::Geometry::m_state
GeomState m_state
Enumeration to dictate whether coefficients are filled.
Definition: Geometry.h:196
Nektar::SpatialDomains::Geometry::m_isoParameter
Array< OneD, Array< OneD, NekDouble > > m_isoParameter
Definition: Geometry.h:209
Nektar::SpatialDomains::Geometry::m_invIsoParam
Array< OneD, Array< OneD, NekDouble > > m_invIsoParam
Definition: Geometry.h:210
Nektar::SpatialDomains::Geometry::ClampLocCoords
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
Nektar::SpatialDomains::Geometry::GetGlobalID
int GetGlobalID(void) const
Get the ID of this object.
Definition: Geometry.h:326
Nektar::SpatialDomains::Geometry::m_coeffs
Array< OneD, Array< OneD, NekDouble > > m_coeffs
Array containing expansion coefficients of m_xmap.
Definition: Geometry.h:206
Nektar::SpatialDomains::Geometry::GetMetricInfo
GeomFactorsSharedPtr GetMetricInfo()
Get the geometric factors for this object.
Definition: Geometry.h:310
Nektar::SpatialDomains::Geometry::m_xmap
StdRegions::StdExpansionSharedPtr m_xmap
mapping containing isoparametric transformation.
Definition: Geometry.h:194
Nektar::SpatialDomains::Geometry::GetXmap
StdRegions::StdExpansionSharedPtr GetXmap() const
Return the mapping object Geometry::m_xmap that represents the coordinate transformation from standar...
Definition: Geometry.h:435
Nektar::SpatialDomains::Geometry::m_geomFactors
GeomFactorsSharedPtr m_geomFactors
Geometric factors.
Definition: Geometry.h:190
Nektar::SpatialDomains::Geometry::m_coordim
int m_coordim
Coordinate dimension of this geometry object.
Definition: Geometry.h:188
Nektar::NekConstants::kNewtonIterations
static const unsigned int kNewtonIterations
Definition: NektarUnivConsts.hpp:64
Nektar::SpatialDomains::eRegular
@ eRegular
Geometry is straight-sided with constant geometric factors.
Definition: SpatialDomains.hpp:52
Nektar::SpatialDomains::PointGeomSharedPtr
std::shared_ptr< PointGeom > PointGeomSharedPtr
Definition: Geometry.h:57
Nektar::SpatialDomains::Geometry2DSharedPtr
std::shared_ptr< Geometry2D > Geometry2DSharedPtr
Definition: Geometry.h:62
Nektar::SpatialDomains::ePtsFilled
@ ePtsFilled
Geometric information has been generated.
Definition: SpatialDomains.hpp:65
Nektar::SpatialDomains::Geometry1DSharedPtr
std::shared_ptr< Geometry1D > Geometry1DSharedPtr
Definition: Geometry.h:61
Nektar::DNekMatSharedPtr
std::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:75
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Vmath::Vmul
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.hpp:72
Vmath::Vvtvp
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.hpp:366
Vmath::Vsum
T Vsum(int n, const T *x, const int incx)
Subtract return sum(x)
Definition: Vmath.hpp:608
Vmath::Dot
T Dot(int n, const T *w, const T *x)
dot product
Definition: Vmath.hpp:761
Vmath::Imin
int Imin(int n, const T *x, const int incx)
Return the index of the minimum element in x.
Definition: Vmath.hpp:704
Vmath::Sadd
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add vector y = alpha + x.
Definition: Vmath.hpp:194
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:825
tinysimd::sqrt
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:285
Generated on Sun Jun 15 2025 07:52:33 for Nektar++ by doxygen 1.9.4