Nektar++
Geometry.cpp
Go to the documentation of this file.
1////////////////////////////////////////////////////////////////////////////////
2//
3// File: Geometry.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: This file contains the base class implementation for the
32// Geometry class.
33//
34//
35////////////////////////////////////////////////////////////////////////////////
36
40
42{
43
44// static class property
46
47/**
48 * @brief Default constructor.
49 */
51 : m_coordim(0), m_geomFactorsState(eNotFilled), m_state(eNotFilled),
52 m_setupState(false), m_shapeType(LibUtilities::eNoShapeType),
53 m_globalID(-1), m_straightEdge(0)
54{
55}
56
57/**
58 * @brief Constructor when supplied a coordinate dimension.
59 */
60Geometry::Geometry(const int coordim)
61 : m_coordim(coordim), m_geomFactorsState(eNotFilled), m_state(eNotFilled),
62 m_setupState(false), m_shapeType(LibUtilities::eNoShapeType),
63 m_globalID(-1), m_straightEdge(0)
64{
65}
66
67/**
68 * @brief Check to see if a geometric factor has already been created that
69 * contains the same regular information.
70 *
71 * The principle behind this is that many regular (i.e. constant Jacobian)
72 * elements have identicial geometric factors. Memory may therefore be reduced
73 * by storing only the unique factors.
74 *
75 * @param geomFactor The GeomFactor to check.
76 *
77 * @return Either the cached GeomFactor or @p geomFactor.
78 *
79 * @todo Currently this method is disabled since the lookup is very expensive.
80 */
82 GeomFactorsSharedPtr geomFactor)
83{
84 GeomFactorsSharedPtr returnval = geomFactor;
85
86 return returnval;
87}
88
89bool SortByGlobalId(const std::shared_ptr<Geometry> &lhs,
90 const std::shared_ptr<Geometry> &rhs)
91{
92 return lhs->GetGlobalID() < rhs->GetGlobalID();
93}
94
95bool GlobalIdEquality(const std::shared_ptr<Geometry> &lhs,
96 const std::shared_ptr<Geometry> &rhs)
97{
98 return lhs->GetGlobalID() == rhs->GetGlobalID();
99}
100
101/**
102 * @brief Get the ID of vertex @p i of this object.
103 */
104int Geometry::GetVid(int i) const
105{
106 return GetVertex(i)->GetGlobalID();
107}
108
109/**
110 * @brief Get the ID of edge @p i of this object.
111 */
112int Geometry::GetEid(int i) const
113{
114 return GetEdge(i)->GetGlobalID();
115}
116
117/**
118 * @brief Get the ID of face @p i of this object.
119 */
120int Geometry::GetFid(int i) const
121{
122 return GetFace(i)->GetGlobalID();
123}
124
125/**
126 * @copydoc Geometry::GetEdge()
127 */
128Geometry1DSharedPtr Geometry::v_GetEdge([[maybe_unused]] int i) const
129{
131 "This function is only valid for shape type geometries");
132 return Geometry1DSharedPtr();
133}
134
135/**
136 * @copydoc Geometry::GetFace()
137 */
138Geometry2DSharedPtr Geometry::v_GetFace([[maybe_unused]] int i) const
139{
141 "This function is only valid for shape type geometries");
142 return Geometry2DSharedPtr();
143}
144
145/**
146 * @copydoc Geometry::GetNumVerts()
147 */
149{
151 "This function is only valid for shape type geometries");
152 return 0;
153}
154
155/**
156 * @copydoc Geometry::GetEorient()
157 */
159 [[maybe_unused]] const int i) const
160{
162 "This function is not valid for this geometry.");
164}
165
166/**
167 * @copydoc Geometry::GetForient()
168 */
170 [[maybe_unused]] const int i) const
171{
173 "This function is not valid for this geometry.");
174 return StdRegions::eFwd;
175}
176
177/**
178 * @copydoc Geometry::GetNumEdges()
179 */
181{
182 return 0;
183}
184
185/**
186 * @copydoc Geometry::GetNumFaces()
187 */
189{
190 return 0;
191}
192
193/**
194 * @copydoc Geometry::GetShapeDim()
195 */
197{
199 "This function is only valid for shape type geometries");
200 return 0;
201}
202
204 [[maybe_unused]] const Array<OneD, const NekDouble> &gloCoord)
205{
206 return 0;
207}
208
210{
212 "This function is only valid for shape type geometries");
213}
214
215/**
216 * @copydoc Geometry::GetXmap()
217 */
219{
220 return m_xmap;
221}
222
223/**
224 * @copydoc Geometry::ContainsPoint(
225 * const Array<OneD, const NekDouble> &, Array<OneD, NekDouble> &,
226 * NekDouble, NekDouble&)
227 * dist is assigned value for curved elements
228 */
230 Array<OneD, NekDouble> &locCoord, NekDouble tol,
231 NekDouble &dist)
232{
233 int inside = PreliminaryCheck(gloCoord);
234 if (inside == -1)
235 {
236 dist = std::numeric_limits<double>::max();
237 return false;
238 }
239 dist = GetLocCoords(gloCoord, locCoord);
240 if (inside == 1)
241 {
242 dist = 0.;
243 return true;
244 }
245 else
246 {
248 m_xmap->LocCoordToLocCollapsed(locCoord, eta);
249 if (ClampLocCoords(eta, tol))
250 {
251 if (GetMetricInfo()->GetGtype() == eRegular)
252 {
253 dist = std::numeric_limits<double>::max();
254 }
255 return false;
256 }
257 return 3 != m_coordim ||
260 dist <= tol;
261 }
262}
263
265 [[maybe_unused]] const Array<OneD, const NekDouble> &xs,
266 [[maybe_unused]] Array<OneD, NekDouble> &xi)
267{
269 "This function has not been defined for this geometry");
270 return false;
271}
272
273/**
274 * @copydoc Geometry::GetVertexEdgeMap()
275 */
276int Geometry::v_GetVertexEdgeMap([[maybe_unused]] const int i,
277 [[maybe_unused]] const int j) const
278{
280 "This function has not been defined for this geometry");
281 return 0;
282}
283
284/**
285 * @copydoc Geometry::GetVertexFaceMap()
286 */
287int Geometry::v_GetVertexFaceMap([[maybe_unused]] const int i,
288 [[maybe_unused]] const int j) const
289{
291 "This function has not been defined for this geometry");
292 return 0;
293}
294
295/**
296 * @copydoc Geometry::GetEdgeFaceMap()
297 */
298int Geometry::v_GetEdgeFaceMap([[maybe_unused]] const int i,
299 [[maybe_unused]] const int j) const
300{
302 "This function has not been defined for this geometry");
303 return 0;
304}
305
306/**
307 * @copydoc Geometry::GetEdgeNormalToFaceVert()
308 */
309int Geometry::v_GetEdgeNormalToFaceVert([[maybe_unused]] const int i,
310 [[maybe_unused]] const int j) const
311{
313 "This function has not been defined for this geometry");
314 return 0;
315}
316
317/**
318 * @copydoc Geometry::GetDir()
319 */
320int Geometry::v_GetDir([[maybe_unused]] const int i,
321 [[maybe_unused]] const int j) const
322{
324 "This function has not been defined for this geometry");
325 return 0;
326}
327
328/**
329 * @copydoc Geometry::GetCoord()
330 */
332 [[maybe_unused]] const int i,
333 [[maybe_unused]] const Array<OneD, const NekDouble> &Lcoord)
334{
336 "This function is only valid for expansion type geometries");
337 return 0.0;
338}
339
340/**
341 * @copydoc Geometry::GetLocCoords()
342 */
344 [[maybe_unused]] const Array<OneD, const NekDouble> &coords,
345 [[maybe_unused]] Array<OneD, NekDouble> &Lcoords)
346{
348 "This function is only valid for expansion type geometries");
349 return 0.0;
350}
351
352/**
353 * @copydoc Geometry::FillGeom()
354 */
356{
358 "This function is only valid for expansion type geometries");
359}
360
361/**
362 * @copydoc Geometry::Reset()
363 */
364void Geometry::v_Reset([[maybe_unused]] CurveMap &curvedEdges,
365 [[maybe_unused]] CurveMap &curvedFaces)
366{
367
368 // Reset state
371
372 // Junk geometric factors
374}
375
377{
379 "This function is only valid for expansion type geometries");
380}
381
382/**
383 * @brief Generates the bounding box for the element.
384 *
385 * For regular elements, the vertices are sufficient to define the extent of
386 * the bounding box. For non-regular elements, the extremes of the quadrature
387 * point coordinates are used. A 10% margin is added around this computed
388 * region to account for convex hull elements where the true extent of the
389 * element may extend slightly beyond the quadrature points.
390 */
391std::array<NekDouble, 6> Geometry::GetBoundingBox()
392{
393 if (m_boundingBox.size() == 6)
394 {
395 return {{m_boundingBox[0], m_boundingBox[1], m_boundingBox[2],
397 }
398 // NekDouble minx, miny, minz, maxx, maxy, maxz;
399 Array<OneD, NekDouble> min(3), max(3);
400
401 // Always get vertexes min/max
403 Array<OneD, NekDouble> x(3, 0.0);
404 p->GetCoords(x[0], x[1], x[2]);
405 for (int j = 0; j < 3; ++j)
406 {
407 min[j] = x[j];
408 max[j] = x[j];
409 }
410 for (int i = 1; i < GetNumVerts(); ++i)
411 {
412 p = GetVertex(i);
413 p->GetCoords(x[0], x[1], x[2]);
414 for (int j = 0; j < 3; ++j)
415 {
416 min[j] = (x[j] < min[j] ? x[j] : min[j]);
417 max[j] = (x[j] > max[j] ? x[j] : max[j]);
418 }
419 }
420 // If element is deformed loop over quadrature points
422 if (GetGeomFactors()->GetGtype() != eRegular)
423 {
424 marginFactor = 0.1;
425 const int nq = GetXmap()->GetTotPoints();
427 for (int j = 0; j < 3; ++j)
428 {
429 xvec[j] = Array<OneD, NekDouble>(nq, 0.0);
430 }
431 for (int j = 0; j < GetCoordim(); ++j)
432 {
433 GetXmap()->BwdTrans(m_coeffs[j], xvec[j]);
434 }
435 for (int j = 0; j < 3; ++j)
436 {
437 for (int i = 0; i < nq; ++i)
438 {
439 min[j] = (xvec[j][i] < min[j] ? xvec[j][i] : min[j]);
440 max[j] = (xvec[j][i] > max[j] ? xvec[j][i] : max[j]);
441 }
442 }
443 }
444 // Add margin to bounding box, in order to
445 // return the nearest element
446 for (int j = 0; j < 3; ++j)
447 {
448 NekDouble margin =
449 marginFactor * (max[j] - min[j]) + NekConstants::kFindDistanceMin;
450 min[j] -= margin;
451 max[j] += margin;
452 }
453
454 // save bounding box
456 for (int j = 0; j < 3; ++j)
457 {
458 m_boundingBox[j] = min[j];
459 m_boundingBox[j + 3] = max[j];
460 }
461 // Return bounding box
462 return {{min[0], min[1], min[2], max[0], max[1], max[2]}};
463}
464
466{
467 m_boundingBox = {};
468}
469
470/**
471 * @brief A fast and robust check if a given global coord is
472 * outside of a deformed element. For regular elements, this
473 * check is unnecessary
474 *
475 * @param coords Input Cartesian global coordinates
476 *
477 * @return 1 is inside of the element.
478 * 0 maybe inside
479 * -1 outside of the element
480 */
482{
483 // bounding box check
484 if (!MinMaxCheck(gloCoord))
485 {
486 return -1;
487 }
488
489 // regular element check
490 if (GetMetricInfo()->GetGtype() == eRegular)
491 {
492 return 0;
493 }
494
495 // All left check for straight edges/plane surfaces
496 return v_AllLeftCheck(gloCoord);
497}
498
499/**
500 * @brief Check if given global coord is within the BoundingBox
501 * of the element.
502 *
503 * @param coords Input Cartesian global coordinates
504 *
505 * @return True if within distance or False otherwise.
506 */
508{
509 // Validation checks
510 ASSERTL1(gloCoord.size() >= m_coordim,
511 "Expects number of global coordinates supplied to be greater than "
512 "or equal to the mesh dimension.");
513
514 std::array<NekDouble, 6> minMax = GetBoundingBox();
515 for (int i = 0; i < m_coordim; ++i)
516 {
517 if ((gloCoord[i] < minMax[i]) || (gloCoord[i] > minMax[i + 3]))
518 {
519 return false;
520 }
521 }
522 return true;
523}
524
525/**
526 * @brief Clamp local coords to be within standard regions [-1, 1]^dim.
527 *
528 * @param Lcoords Corresponding local coordinates
529 */
531{
532 // Validation checks
533 ASSERTL1(locCoord.size() >= GetShapeDim(),
534 "Expects local coordinates to be same or "
535 "larger than shape dimension.");
536
537 // If out of range clamp locCoord to be within [-1,1]^dim
538 // since any larger value will be very oscillatory if
539 // called by 'returnNearestElmt' option in
540 // ExpList::GetExpIndex
541 bool clamp = false;
542 for (int i = 0; i < GetShapeDim(); ++i)
543 {
544 if (!std::isfinite(locCoord[i]))
545 {
546 locCoord[i] = 0.;
547 clamp = true;
548 }
549 else if (locCoord[i] < -(1. + tol))
550 {
551 locCoord[i] = -(1. + tol);
552 clamp = true;
553 }
554 else if (locCoord[i] > (1. + tol))
555 {
556 locCoord[i] = 1. + tol;
557 clamp = true;
558 }
559 }
560 return clamp;
561}
562
563} // namespace Nektar::SpatialDomains
#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
virtual int v_GetNumEdges() const
Get the number of edges of this object.
Definition: Geometry.cpp:180
virtual StdRegions::Orientation v_GetForient(const int i) const
Returns the orientation of face i with respect to the ordering of faces in the standard element.
Definition: Geometry.cpp:169
virtual NekDouble v_GetCoord(const int i, const Array< OneD, const NekDouble > &Lcoord)
Given local collapsed coordinate Lcoord, return the value of physical coordinate in direction i.
Definition: Geometry.cpp:331
Array< OneD, NekDouble > m_boundingBox
Array containing bounding box.
Definition: Geometry.h:208
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
virtual void v_CalculateInverseIsoParam()
Definition: Geometry.cpp:209
virtual int v_AllLeftCheck(const Array< OneD, const NekDouble > &gloCoord)
Definition: Geometry.cpp:203
NekDouble GetLocCoords(const Array< OneD, const NekDouble > &coords, Array< OneD, NekDouble > &Lcoords)
Determine the local collapsed coordinates that correspond to a given Cartesian coordinate for this ge...
Definition: Geometry.h:544
int PreliminaryCheck(const Array< OneD, const NekDouble > &gloCoord)
A fast and robust check if a given global coord is outside of a deformed element. For regular element...
Definition: Geometry.cpp:481
virtual int v_GetEdgeNormalToFaceVert(const int i, const int j) const
Returns the standard lement edge IDs that are normal to a given face vertex.
Definition: Geometry.cpp:309
Geometry1DSharedPtr GetEdge(int i) const
Returns edge i of this object.
Definition: Geometry.h:365
int GetShapeDim() const
Get the object's shape dimension.
Definition: Geometry.h:426
virtual void v_FillGeom()
Populate the coordinate mapping Geometry::m_coeffs information from any children geometry elements.
Definition: Geometry.cpp:355
Geometry()
Default constructor.
Definition: Geometry.cpp:50
virtual NekDouble v_GetLocCoords(const Array< OneD, const NekDouble > &coords, Array< OneD, NekDouble > &Lcoords)
Determine the local collapsed coordinates that correspond to a given Cartesian coordinate for this ge...
Definition: Geometry.cpp:343
virtual Geometry2DSharedPtr v_GetFace(int i) const
Returns face i of this object.
Definition: Geometry.cpp:138
int GetVid(int i) const
Get the ID of vertex i of this object.
Definition: Geometry.cpp:104
virtual int v_GetNumFaces() const
Get the number of faces of this object.
Definition: Geometry.cpp:188
virtual StdRegions::StdExpansionSharedPtr v_GetXmap() const
Return the mapping object Geometry::m_xmap that represents the coordinate transformation from standar...
Definition: Geometry.cpp:218
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
bool ClampLocCoords(Array< OneD, NekDouble > &locCoord, NekDouble tol=std::numeric_limits< NekDouble >::epsilon())
Clamp local coords to be within standard regions [-1, 1]^dim.
Definition: Geometry.cpp:530
int GetCoordim() const
Return the coordinate dimension of this object (i.e. the dimension of the space in which this object ...
Definition: Geometry.h:283
int GetFid(int i) const
Get the ID of face i of this object.
Definition: Geometry.cpp:120
virtual int v_GetEdgeFaceMap(int i, int j) const
Returns the standard element edge IDs that are connected to a given face.
Definition: Geometry.cpp:298
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
GeomFactorsSharedPtr GetMetricInfo()
Get the geometric factors for this object.
Definition: Geometry.h:310
StdRegions::StdExpansionSharedPtr m_xmap
mapping containing isoparametric transformation.
Definition: Geometry.h:194
virtual Geometry1DSharedPtr v_GetEdge(int i) const
Returns edge i of this object.
Definition: Geometry.cpp:128
virtual StdRegions::Orientation v_GetEorient(const int i) const
Returns the orientation of edge i with respect to the ordering of edges in the standard element.
Definition: Geometry.cpp:158
StdRegions::StdExpansionSharedPtr GetXmap() const
Return the mapping object Geometry::m_xmap that represents the coordinate transformation from standar...
Definition: Geometry.h:435
static GeomFactorsVector m_regGeomFactorsManager
Definition: Geometry.h:185
std::array< NekDouble, 6 > GetBoundingBox()
Generates the bounding box for the element.
Definition: Geometry.cpp:391
int GetNumVerts() const
Get the number of vertices of this object.
Definition: Geometry.h:399
Geometry2DSharedPtr GetFace(int i) const
Returns face i of this object.
Definition: Geometry.h:373
bool MinMaxCheck(const Array< OneD, const NekDouble > &gloCoord)
Check if given global coord is within the BoundingBox of the element.
Definition: Geometry.cpp:507
virtual int v_GetVertexEdgeMap(int i, int j) const
Returns the standard element edge IDs that are connected to a given vertex.
Definition: Geometry.cpp:276
virtual NekDouble v_FindDistance(const Array< OneD, const NekDouble > &xs, Array< OneD, NekDouble > &xi)
Definition: Geometry.cpp:264
GeomFactorsSharedPtr GetGeomFactors()
Get the geometric factors for this object, generating them if required.
Definition: Geometry.h:301
static GeomFactorsSharedPtr ValidateRegGeomFactor(GeomFactorsSharedPtr geomFactor)
Check to see if a geometric factor has already been created that contains the same regular informatio...
Definition: Geometry.cpp:81
virtual int v_GetNumVerts() const
Get the number of vertices of this object.
Definition: Geometry.cpp:148
GeomFactorsSharedPtr m_geomFactors
Geometric factors.
Definition: Geometry.h:190
int m_coordim
Coordinate dimension of this geometry object.
Definition: Geometry.h:188
virtual int v_GetDir(const int faceidx, const int facedir) const
Returns the element coordinate direction corresponding to a given face coordinate direction.
Definition: Geometry.cpp:320
virtual int v_GetShapeDim() const
Get the object's shape dimension.
Definition: Geometry.cpp:196
virtual int v_GetVertexFaceMap(int i, int j) const
Returns the standard element face IDs that are connected to a given vertex.
Definition: Geometry.cpp:287
virtual bool v_ContainsPoint(const Array< OneD, const NekDouble > &gloCoord, Array< OneD, NekDouble > &locCoord, NekDouble tol, NekDouble &dist)
Determine whether an element contains a particular Cartesian coordinate .
Definition: Geometry.cpp:229
int GetEid(int i) const
Get the ID of edge i of this object.
Definition: Geometry.cpp:112
static const NekDouble kGeomFactorsTol
static const NekDouble kFindDistanceMin
std::vector< GeomFactorsSharedPtr > GeomFactorsVector
A vector of GeomFactor pointers.
Definition: GeomFactors.h:62
bool GlobalIdEquality(const std::shared_ptr< Geometry > &lhs, const std::shared_ptr< Geometry > &rhs)
Definition: Geometry.cpp:95
std::shared_ptr< GeomFactors > GeomFactorsSharedPtr
Pointer to a GeomFactors object.
Definition: GeomFactors.h:60
bool SortByGlobalId(const std::shared_ptr< Geometry > &lhs, const std::shared_ptr< Geometry > &rhs)
Less than operator to sort Geometry objects by global id when sorting STL containers.
Definition: Geometry.cpp:89
std::unordered_map< int, CurveSharedPtr > CurveMap
Definition: Curve.hpp:59
@ eRegular
Geometry is straight-sided with constant geometric factors.
std::shared_ptr< PointGeom > PointGeomSharedPtr
Definition: Geometry.h:57
std::shared_ptr< Geometry2D > Geometry2DSharedPtr
Definition: Geometry.h:62
@ eNotFilled
Geometric information has not been generated.
std::shared_ptr< Geometry1D > Geometry1DSharedPtr
Definition: Geometry.h:61
std::shared_ptr< StdExpansion > StdExpansionSharedPtr
double NekDouble