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