Nektar++
TriGeom.cpp
Go to the documentation of this file.
1////////////////////////////////////////////////////////////////////////////////
2//
3// File: TriGeom.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:
32//
33//
34////////////////////////////////////////////////////////////////////////////////
35
40
43
45{
46
48{
50}
51
52TriGeom::TriGeom(const int id, const SegGeomSharedPtr edges[],
53 const CurveSharedPtr curve)
54 : Geometry2D(edges[0]->GetVertex(0)->GetCoordim(), curve)
55{
56 int j;
57
59 m_globalID = id;
60
61 // Copy the edge shared pointers.
62 m_edges.insert(m_edges.begin(), edges, edges + TriGeom::kNedges);
63 m_eorient.resize(kNedges);
64
65 for (j = 0; j < kNedges; ++j)
66 {
67 m_eorient[j] =
68 SegGeom::GetEdgeOrientation(*edges[j], *edges[(j + 1) % kNedges]);
69 m_verts.push_back(
70 edges[j]->GetVertex(m_eorient[j] == StdRegions::eForwards ? 0 : 1));
71 }
72
76
77 m_coordim = edges[0]->GetVertex(0)->GetCoordim();
78 ASSERTL0(m_coordim > 1, "Cannot call function with dim == 1");
79}
80
82{
83 // From Geometry
85
86 // From TriFaceComponent
88
89 // From TriGeom
90 m_verts = in.m_verts;
91 m_edges = in.m_edges;
92 for (int i = 0; i < kNedges; i++)
93 {
94 m_eorient[i] = in.m_eorient[i];
95 }
96}
97
99 const Array<OneD, const NekDouble> &Lcoord)
100{
101 ASSERTL1(m_state == ePtsFilled, "Geometry is not in physical space");
102
103 Array<OneD, NekDouble> tmp(m_xmap->GetTotPoints());
104 m_xmap->BwdTrans(m_coeffs[i], tmp);
105
106 return m_xmap->PhysEvaluate(Lcoord, tmp);
107}
108
110 const TriGeom &face2,
111 bool doRot, int dir,
112 NekDouble angle,
113 NekDouble tol)
114{
115 return GetFaceOrientation(face1.m_verts, face2.m_verts, doRot, dir, angle,
116 tol);
117}
118
120 const PointGeomVector &face1, const PointGeomVector &face2, bool doRot,
121 int dir, NekDouble angle, NekDouble tol)
122{
123 int i, j, vmap[3] = {-1, -1, -1};
124
125 if (doRot)
126 {
127 PointGeom rotPt;
128
129 for (i = 0; i < 3; ++i)
130 {
131 rotPt.Rotate((*face1[i]), dir, angle);
132 for (j = 0; j < 3; ++j)
133 {
134 if (rotPt.dist(*face2[j]) < tol)
135 {
136 vmap[j] = i;
137 break;
138 }
139 }
140 }
141 }
142 else
143 {
144
145 NekDouble x, y, z, x1, y1, z1, cx = 0.0, cy = 0.0, cz = 0.0;
146
147 // For periodic faces, we calculate the vector between the centre
148 // points of the two faces. (For connected faces this will be
149 // zero). We can then use this to determine alignment later in the
150 // algorithm.
151 for (i = 0; i < 3; ++i)
152 {
153 cx += (*face2[i])(0) - (*face1[i])(0);
154 cy += (*face2[i])(1) - (*face1[i])(1);
155 cz += (*face2[i])(2) - (*face1[i])(2);
156 }
157 cx /= 3;
158 cy /= 3;
159 cz /= 3;
160
161 // Now construct a mapping which takes us from the vertices of one
162 // face to the other. That is, vertex j of face2 corresponds to
163 // vertex vmap[j] of face1.
164 for (i = 0; i < 3; ++i)
165 {
166 x = (*face1[i])(0);
167 y = (*face1[i])(1);
168 z = (*face1[i])(2);
169 for (j = 0; j < 3; ++j)
170 {
171 x1 = (*face2[j])(0) - cx;
172 y1 = (*face2[j])(1) - cy;
173 z1 = (*face2[j])(2) - cz;
174 if (sqrt((x1 - x) * (x1 - x) + (y1 - y) * (y1 - y) +
175 (z1 - z) * (z1 - z)) < 1e-8)
176 {
177 vmap[j] = i;
178 break;
179 }
180 }
181 }
182 }
183
184 if (vmap[1] == (vmap[0] + 1) % 3)
185 {
186 switch (vmap[0])
187 {
188 case 0:
190 break;
191 case 1:
193 break;
194 case 2:
196 break;
197 }
198 }
199 else
200 {
201 switch (vmap[0])
202 {
203 case 0:
205 break;
206 case 1:
208 break;
209 case 2:
211 break;
212 }
213 }
214
215 ASSERTL0(false, "Unable to determine triangle orientation");
217}
218
220{
221 if (!m_setupState)
222 {
224 }
225
227 {
228 GeomType Gtype = eRegular;
229
231
232 // check to see if expansions are linear
233 m_straightEdge = 1;
234 if (m_xmap->GetBasisNumModes(0) != 2 ||
235 m_xmap->GetBasisNumModes(1) != 2)
236 {
237 Gtype = eDeformed;
238 m_straightEdge = 0;
239 }
240
242 m_manifold[0] = 0;
243 m_manifold[1] = 1;
244 if (m_coordim == 3)
245 {
246 PointGeom e01, e21, norm;
247 e01.Sub(*m_verts[0], *m_verts[1]);
248 e21.Sub(*m_verts[2], *m_verts[1]);
249 norm.Mult(e01, e21);
250 int tmpi = 0;
251 double tmp = std::fabs(norm[0]);
252 if (tmp < fabs(norm[1]))
253 {
254 tmp = fabs(norm[1]);
255 tmpi = 1;
256 }
257 if (tmp < fabs(norm[2]))
258 {
259 tmpi = 2;
260 }
261 m_manifold[0] = (tmpi + 1) % 3;
262 m_manifold[1] = (tmpi + 2) % 3;
263 m_manifold[2] = (tmpi + 3) % 3;
264 }
265 if (Gtype == eRegular)
266 {
268 for (int i = 0; i < m_verts.size(); ++i)
269 {
270 verts[i] = Array<OneD, NekDouble>(3);
271 m_verts[i]->GetCoords(verts[i]);
272 }
273 // a00 + a01 xi1 + a02 xi2
274 // a10 + a11 xi1 + a12 xi2
276 for (int i = 0; i < 2; ++i)
277 {
278 unsigned int d = m_manifold[i];
280 NekDouble A = verts[0][d];
281 NekDouble B = verts[1][d];
282 NekDouble C = verts[2][d];
283 m_isoParameter[i][0] = 0.5 * (B + C); // 1
284 m_isoParameter[i][1] = 0.5 * (-A + B); // xi1
285 m_isoParameter[i][2] = 0.5 * (-A + C); // xi2
286 }
288 }
289
291 Gtype, m_coordim, m_xmap, m_coeffs);
292
294 }
295}
296
297/**
298 * Note verts and edges are listed according to anticlockwise
299 * convention but points in _coeffs have to be in array format from
300 * left to right.
301 */
303{
304 // check to see if geometry structure is already filled
305 if (m_state == ePtsFilled)
306 {
307 return;
308 }
309
310 int i, j, k;
311 int nEdgeCoeffs = m_xmap->GetTraceNcoeffs(0);
312
313 if (m_curve)
314 {
316 2, m_curve->m_ptype)]
317 ->GetPointsDim();
318
319 // Deal with 2D points type separately
320 // (e.g. electrostatic or Fekete points) to 1D tensor
321 // product.
322 if (pdim == 2)
323 {
324 int N = m_curve->m_points.size();
325 int nEdgePts =
326 (-1 + (int)sqrt(static_cast<NekDouble>(8 * N + 1))) / 2;
327
328 ASSERTL0(nEdgePts * (nEdgePts + 1) / 2 == N,
329 "NUMPOINTS should be a triangle number for"
330 " triangle curved face " +
331 std::to_string(m_globalID));
332
333 // Sanity check 1: are curved vertices consistent with
334 // triangle vertices?
335 for (i = 0; i < 3; ++i)
336 {
337 NekDouble dist = m_verts[i]->dist(*(m_curve->m_points[i]));
339 {
340 std::stringstream ss;
341 ss << "Curved vertex " << i << " of triangle " << m_globalID
342 << " is separated from expansion vertex by"
343 << " more than " << NekConstants::kVertexTheSameDouble
344 << " (dist = " << dist << ")";
345 NEKERROR(ErrorUtil::ewarning, ss.str().c_str());
346 }
347 }
348
349 // Sanity check 2: are curved edges from the face curvature
350 // consistent with curved edges?
351 for (i = 0; i < kNedges; ++i)
352 {
353 CurveSharedPtr edgeCurve = m_edges[i]->GetCurve();
354
355 ASSERTL0(edgeCurve->m_points.size() == nEdgePts,
356 "Number of edge points does not correspond "
357 "to number of face points in triangle " +
358 std::to_string(m_globalID));
359
360 const int offset = 3 + i * (nEdgePts - 2);
361 NekDouble maxDist = 0.0;
362
363 // Account for different ordering of nodal coordinates
364 // vs. Cartesian ordering of element.
366
367 if (i == 2)
368 {
369 orient = orient == StdRegions::eForwards
372 }
373
374 if (orient == StdRegions::eForwards)
375 {
376 for (j = 0; j < nEdgePts - 2; ++j)
377 {
378 NekDouble dist = m_curve->m_points[offset + j]->dist(
379 *(edgeCurve->m_points[j + 1]));
380 maxDist = dist > maxDist ? dist : maxDist;
381 }
382 }
383 else
384 {
385 for (j = 0; j < nEdgePts - 2; ++j)
386 {
387 NekDouble dist = m_curve->m_points[offset + j]->dist(
388 *(edgeCurve->m_points[nEdgePts - 2 - j]));
389 maxDist = dist > maxDist ? dist : maxDist;
390 }
391 }
392
394 {
395 std::stringstream ss;
396 ss << "Curved edge " << i << " of triangle " << m_globalID
397 << " has a point separated from edge interior"
398 << " points by more than "
400 << " (maxdist = " << maxDist << ")";
401 NEKERROR(ErrorUtil::ewarning, ss.str().c_str());
402 }
403 }
404
408 nEdgePts, LibUtilities::eGaussRadauMAlpha1Beta0);
410 P0);
412 P1);
414 std::max(nEdgePts * nEdgePts, m_xmap->GetTotPoints()));
415 Array<OneD, NekDouble> tmp(nEdgePts * nEdgePts);
416
417 for (i = 0; i < m_coordim; ++i)
418 {
419 // Create a StdNodalTriExp.
422 AllocateSharedPtr(T0, T1, m_curve->m_ptype);
423
424 for (j = 0; j < N; ++j)
425 {
426 phys[j] = (m_curve->m_points[j]->GetPtr())[i];
427 }
428
429 t->BwdTrans(phys, tmp);
430
431 // Interpolate points to standard region.
433 P0, P1, tmp, m_xmap->GetBasis(0)->GetPointsKey(),
434 m_xmap->GetBasis(1)->GetPointsKey(), phys);
435
436 // Forwards transform to get coefficient space.
437 m_xmap->FwdTrans(phys, m_coeffs[i]);
438 }
439 }
440 else if (pdim == 1)
441 {
442 int npts = m_curve->m_points.size();
443 int nEdgePts = (int)sqrt(static_cast<NekDouble>(npts));
444 Array<OneD, NekDouble> tmp(npts);
445 Array<OneD, NekDouble> phys(m_xmap->GetTotPoints());
446 LibUtilities::PointsKey curveKey(nEdgePts, m_curve->m_ptype);
447
448 // Sanity checks:
449 // - Curved faces should have square number of points;
450 // - Each edge should have sqrt(npts) points.
451 ASSERTL0(nEdgePts * nEdgePts == npts,
452 "NUMPOINTS should be a square number for"
453 " triangle " +
454 std::to_string(m_globalID));
455
456 for (i = 0; i < kNedges; ++i)
457 {
458 ASSERTL0(m_edges[i]->GetXmap()->GetNcoeffs() == nEdgePts,
459 "Number of edge points does not correspond to "
460 "number of face points in triangle " +
461 std::to_string(m_globalID));
462 }
463
464 for (i = 0; i < m_coordim; ++i)
465 {
466 for (j = 0; j < npts; ++j)
467 {
468 tmp[j] = (m_curve->m_points[j]->GetPtr())[i];
469 }
470
471 // Interpolate curve points to standard triangle
472 // points.
473 LibUtilities::Interp2D(curveKey, curveKey, tmp,
474 m_xmap->GetBasis(0)->GetPointsKey(),
475 m_xmap->GetBasis(1)->GetPointsKey(),
476 phys);
477
478 // Forwards transform to get coefficient space.
479 m_xmap->FwdTrans(phys, m_coeffs[i]);
480 }
481 }
482 else
483 {
484 ASSERTL0(false, "Only 1D/2D points distributions "
485 "supported.");
486 }
487 }
488
489 Array<OneD, unsigned int> mapArray(nEdgeCoeffs);
490 Array<OneD, int> signArray(nEdgeCoeffs);
491
492 for (i = 0; i < kNedges; i++)
493 {
494 m_edges[i]->FillGeom();
495 m_xmap->GetTraceToElementMap(i, mapArray, signArray, m_eorient[i]);
496
497 nEdgeCoeffs = m_edges[i]->GetXmap()->GetNcoeffs();
498
499 for (j = 0; j < m_coordim; j++)
500 {
501 for (k = 0; k < nEdgeCoeffs; k++)
502 {
503 m_coeffs[j][mapArray[k]] =
504 signArray[k] * m_edges[i]->GetCoeffs(j)[k];
505 }
506 }
507 }
508
510}
511
512int TriGeom::v_GetDir(const int i, [[maybe_unused]] const int j) const
513{
514 return i == 0 ? 0 : 1;
515}
516
517void TriGeom::v_Reset(CurveMap &curvedEdges, CurveMap &curvedFaces)
518{
519 Geometry::v_Reset(curvedEdges, curvedFaces);
520 CurveMap::iterator it = curvedFaces.find(m_globalID);
521
522 if (it != curvedFaces.end())
523 {
524 m_curve = it->second;
525 }
526
527 for (int i = 0; i < 3; ++i)
528 {
529 m_edges[i]->Reset(curvedEdges, curvedFaces);
530 }
531
532 SetUpXmap();
533 SetUpCoeffs(m_xmap->GetNcoeffs());
534}
535
537{
538 if (!m_setupState)
539 {
540 for (int i = 0; i < 3; ++i)
541 {
542 m_edges[i]->Setup();
543 }
544 SetUpXmap();
545 SetUpCoeffs(m_xmap->GetNcoeffs());
546 m_setupState = true;
547 }
548}
549
551{
552 int order0 = m_edges[0]->GetXmap()->GetBasis(0)->GetNumModes();
553 int order1 = std::max(
554 order0, std::max(m_edges[1]->GetXmap()->GetBasis(0)->GetNumModes(),
555 m_edges[2]->GetXmap()->GetBasis(0)->GetNumModes()));
556
557 const LibUtilities::BasisKey B0(
559 LibUtilities::PointsKey(order0 + 1,
561 const LibUtilities::BasisKey B1(
563 LibUtilities::PointsKey(order1, LibUtilities::eGaussRadauMAlpha1Beta0));
564
566}
567
568} // namespace Nektar::SpatialDomains
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mode...
Definition: ErrorUtil.hpp:202
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:242
Describes the specification for a Basis.
Definition: Basis.h:45
Defines a specification for a set of points.
Definition: Points.h:50
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
2D geometry information
Definition: Geometry2D.h:65
Array< OneD, int > m_manifold
Definition: Geometry2D.h:83
std::vector< StdRegions::Orientation > m_eorient
Definition: Geometry2D.h:81
void v_CalculateInverseIsoParam() override
Definition: Geometry2D.cpp:629
bool m_setupState
Wether or not the setup routines have been run.
Definition: Geometry.h:198
PointGeomSharedPtr GetVertex(int i) const
Returns vertex i of this object.
Definition: Geometry.h:357
GeomState m_state
Enumeration to dictate whether coefficients are filled.
Definition: Geometry.h:196
void SetUpCoeffs(const int nCoeffs)
Initialise the Geometry::m_coeffs array.
Definition: Geometry.h:700
Array< OneD, Array< OneD, NekDouble > > m_isoParameter
Definition: Geometry.h:209
virtual void v_Reset(CurveMap &curvedEdges, CurveMap &curvedFaces)
Reset this geometry object: unset the current state, zero Geometry::m_coeffs and remove allocated Geo...
Definition: Geometry.cpp:364
LibUtilities::ShapeType m_shapeType
Type of shape.
Definition: Geometry.h:202
Array< OneD, Array< OneD, NekDouble > > m_coeffs
Array containing expansion coefficients of m_xmap.
Definition: Geometry.h:206
GeomState m_geomFactorsState
State of the geometric factors.
Definition: Geometry.h:192
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
void Sub(PointGeom &a, PointGeom &b)
Definition: PointGeom.cpp:119
void Mult(PointGeom &a, PointGeom &b)
_this = a x b
Definition: PointGeom.cpp:128
void Rotate(PointGeom &a, int dir, NekDouble angle)
_this = rotation of a by angle 'angle' around axis dir
Definition: PointGeom.cpp:137
NekDouble dist(PointGeom &a)
return distance between this and input a
Definition: PointGeom.cpp:175
static StdRegions::Orientation GetEdgeOrientation(const SegGeom &edge1, const SegGeom &edge2)
Get the orientation of edge1.
Definition: SegGeom.cpp:190
static const int kNedges
Get the orientation of face1.
Definition: TriGeom.h:70
void v_GenGeomFactors() override
Definition: TriGeom.cpp:219
int v_GetDir(const int faceidx, const int facedir) const override
Returns the element coordinate direction corresponding to a given face coordinate direction.
Definition: TriGeom.cpp:512
static StdRegions::Orientation GetFaceOrientation(const TriGeom &face1, const TriGeom &face2, bool doRot, int dir, NekDouble angle, NekDouble tol)
Definition: TriGeom.cpp:109
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: TriGeom.cpp:98
void v_Reset(CurveMap &curvedEdges, CurveMap &curvedFaces) override
Reset this geometry object: unset the current state, zero Geometry::m_coeffs and remove allocated Geo...
Definition: TriGeom.cpp:517
PointsManagerT & PointsManager(void)
void Interp2D(const BasisKey &fbasis0, const BasisKey &fbasis1, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, const BasisKey &tbasis1, Array< OneD, NekDouble > &to)
this function interpolates a 2D function evaluated at the quadrature points of the 2D basis,...
Definition: Interp.cpp:101
@ eGaussLobattoLegendre
1D Gauss-Lobatto-Legendre quadrature points
Definition: PointsType.h:51
@ eModified_B
Principle Modified Functions .
Definition: BasisType.h:49
@ eOrtho_A
Principle Orthogonal Functions .
Definition: BasisType.h:42
@ eOrtho_B
Principle Orthogonal Functions .
Definition: BasisType.h:44
@ eModified_A
Principle Modified Functions .
Definition: BasisType.h:48
static const NekDouble kVertexTheSameDouble
std::vector< PointGeomSharedPtr > PointGeomVector
Definition: Geometry2D.h:60
std::shared_ptr< Curve > CurveSharedPtr
Definition: Curve.hpp:58
std::unordered_map< int, CurveSharedPtr > CurveMap
Definition: Curve.hpp:59
GeomType
Indicates the type of element geometry.
@ eRegular
Geometry is straight-sided with constant geometric factors.
@ eDeformed
Geometry is curved or has non-constant factors.
std::shared_ptr< SegGeom > SegGeomSharedPtr
Definition: Geometry2D.h:58
@ ePtsFilled
Geometric information has been generated.
std::shared_ptr< StdNodalTriExp > StdNodalTriExpSharedPtr
std::vector< double > z(NPUPPER)
std::vector< double > d(NPUPPER *NPUPPER)
double NekDouble
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:285