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