MFEM v4.8.0
Finite element discretization library
Loading...
Searching...
No Matches
geom.hpp
Go to the documentation of this file.
1// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
2// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3// LICENSE and NOTICE for details. LLNL-CODE-806117.
4//
5// This file is part of the MFEM library. For more information and source code
6// availability visit https://mfem.org.
7//
8// MFEM is free software; you can redistribute it and/or modify it under the
9// terms of the BSD-3 license. We welcome feedback and contributions, see file
10// CONTRIBUTING.md for details.
11
12#ifndef MFEM_GEOM
13#define MFEM_GEOM
14
15#include "../config/config.hpp"
17#include "intrules.hpp"
18#include "../general/hash.hpp"
19
20#include <memory>
21#include <unordered_map>
22
23namespace mfem
24{
25
26/** Types of domains for integration rules and reference finite elements:
27 Geometry::POINT - a point
28 Geometry::SEGMENT - the interval [0,1]
29 Geometry::TRIANGLE - triangle with vertices (0,0), (1,0), (0,1)
30 Geometry::SQUARE - the unit square (0,1)x(0,1)
31 Geometry::TETRAHEDRON - w/ vert. (0,0,0),(1,0,0),(0,1,0),(0,0,1)
32 Geometry::CUBE - the unit cube
33 Geometry::PRISM - w/ vert. (0,0,0),(1,0,0),(0,1,0),(0,0,1),(1,0,1),(0,1,1)
34 Geometry::PYRAMID - w/ vert. (0,0,0),(1,0,0),(1,1,0),(0,1,0),(0,0,1)
35*/
36class MFEM_EXPORT Geometry
37{
38public:
39 enum Type
40 {
41 INVALID = -1,
42 POINT = 0, SEGMENT, TRIANGLE, SQUARE, TETRAHEDRON, CUBE, PRISM, PYRAMID,
43 NUM_GEOMETRIES
44 };
45
46 static const int NumGeom = NUM_GEOMETRIES;
47 static const int MaxDim = 3;
48 static const int NumBdrArray[NumGeom];
49 static const char *Name[NumGeom];
50 static const real_t Volume[NumGeom];
51 static const int Dimension[NumGeom];
52 static const int DimStart[MaxDim+2]; // including MaxDim+1
53 static const int NumVerts[NumGeom];
54 static const int NumEdges[NumGeom];
55 static const int NumFaces[NumGeom];
56
57 // Structure that holds constants describing the Geometries.
58 template <Type Geom> struct Constants;
59
60private:
61 IntegrationRule *GeomVert[NumGeom];
62 IntegrationPoint GeomCenter[NumGeom];
63 DenseMatrix *GeomToPerfGeomJac[NumGeom];
64 DenseMatrix *PerfGeomToGeomJac[NumGeom];
65
66public:
67 Geometry();
68 ~Geometry();
69
70 /** @brief Return an IntegrationRule consisting of all vertices of the given
71 Geometry::Type, @a GeomType. */
72 const IntegrationRule *GetVertices(int GeomType) const;
73
74 /// Return the center of the given Geometry::Type, @a GeomType.
75 const IntegrationPoint &GetCenter(int GeomType) const
76 { return GeomCenter[GeomType]; }
77
78 /// Get a random point in the reference element specified by @a GeomType.
79 /** This method uses the function rand() for random number generation. */
80 static void GetRandomPoint(int GeomType, IntegrationPoint &ip);
81
82 /// Check if the given point is inside the given reference element.
83 static bool CheckPoint(int GeomType, const IntegrationPoint &ip);
84 /** @brief Check if the given point is inside the given reference element.
85 Overload for fuzzy tolerance. */
86 static bool CheckPoint(int GeomType, const IntegrationPoint &ip, real_t eps);
87
88 /// Project a point @a end, onto the given Geometry::Type, @a GeomType.
89 /** Check if the @a end point is inside the reference element, if not
90 overwrite it with the point on the boundary that lies on the line segment
91 between @a beg and @a end (@a beg must be inside the element). Return
92 true if @a end is inside the element, and false otherwise. */
93 static bool ProjectPoint(int GeomType, const IntegrationPoint &beg,
94 IntegrationPoint &end);
95
96 /// Project a point @a ip, onto the given Geometry::Type, @a GeomType.
97 /** If @a ip is outside the element, replace it with the point on the
98 boundary that is closest to the original @a ip and return false;
99 otherwise, return true without changing @a ip. */
100 static bool ProjectPoint(int GeomType, IntegrationPoint &ip);
101
102 const DenseMatrix &GetGeomToPerfGeomJac(int GeomType) const
103 { return *GeomToPerfGeomJac[GeomType]; }
104 const DenseMatrix *GetPerfGeomToGeomJac(int GeomType) const
105 { return PerfGeomToGeomJac[GeomType]; }
106 void GetPerfPointMat(int GeomType, DenseMatrix &pm) const;
107 void JacToPerfJac(int GeomType, const DenseMatrix &J,
108 DenseMatrix &PJ) const;
109
110 /// Returns true if the given @a geom is of tensor-product type (i.e. if geom
111 /// is a segment, quadrilateral, or hexahedron), returns false otherwise.
112 static bool IsTensorProduct(Type geom)
113 { return geom == SEGMENT || geom == SQUARE || geom == CUBE; }
114
115 /// Returns the Geometry::Type corresponding to a tensor-product of the
116 /// given dimension.
118 {
119 switch (dim)
120 {
121 case 0: return POINT;
122 case 1: return SEGMENT;
123 case 2: return SQUARE;
124 case 3: return CUBE;
125 default: MFEM_ABORT("Invalid dimension."); return INVALID;
126 }
127 }
128
129 /// Return the inverse of the given orientation for the specified geometry type.
130 static int GetInverseOrientation(Type geom_type, int orientation);
131
132 /// Return the number of boundary "faces" of a given Geometry::Type.
133 int NumBdr(int GeomType) const { return NumBdrArray[GeomType]; }
134};
135
136template <> struct
137/// @cond Suppress_Doxygen_warnings
138 MFEM_EXPORT
139/// @endcond
141{
142 static const int Dimension = 0;
143 static const int NumVert = 1;
144
145 static const int NumOrient = 1;
146 static const int Orient[NumOrient][NumVert];
147 static const int InvOrient[NumOrient];
148};
149
150template <> struct
151/// @cond Suppress_Doxygen_warnings
152 MFEM_EXPORT
153/// @endcond
155{
156 static const int Dimension = 1;
157 static const int NumVert = 2;
158 static const int NumEdges = 1;
159 static const int Edges[NumEdges][2];
160
161 static const int NumOrient = 2;
162 static const int Orient[NumOrient][NumVert];
163 static const int InvOrient[NumOrient];
164};
165
166template <> struct
167/// @cond Suppress_Doxygen_warnings
168 MFEM_EXPORT
169/// @endcond
171{
172 static const int Dimension = 2;
173 static const int NumVert = 3;
174 static const int NumEdges = 3;
175 static const int Edges[NumEdges][2];
176 // Upper-triangular part of the local vertex-to-vertex graph.
178 {
179 static const int I[NumVert];
180 static const int J[NumEdges][2]; // {end,edge_idx}
181 };
182 static const int NumFaces = 1;
183 static const int FaceVert[NumFaces][NumVert];
184
185 // For a given base tuple v={v0,v1,v2}, the orientation of a permutation
186 // u={u0,u1,u2} of v, is an index 'j' such that u[i]=v[Orient[j][i]].
187 // The static method Mesh::GetTriOrientation, computes the index 'j' of the
188 // permutation that maps the second argument 'test' to the first argument
189 // 'base': test[Orient[j][i]]=base[i].
190 static const int NumOrient = 6;
191 static const int Orient[NumOrient][NumVert];
192 // The inverse of orientation 'j' is InvOrient[j].
193 static const int InvOrient[NumOrient];
194};
195
196template <> struct
197/// @cond Suppress_Doxygen_warnings
198 MFEM_EXPORT
199/// @endcond
201{
202 static const int Dimension = 2;
203 static const int NumVert = 4;
204 static const int NumEdges = 4;
205 static const int Edges[NumEdges][2];
206 // Upper-triangular part of the local vertex-to-vertex graph.
208 {
209 static const int I[NumVert];
210 static const int J[NumEdges][2]; // {end,edge_idx}
211 };
212 static const int NumFaces = 1;
213 static const int FaceVert[NumFaces][NumVert];
214
215 static const int NumOrient = 8;
216 static const int Orient[NumOrient][NumVert];
217 static const int InvOrient[NumOrient];
218};
219
220template <> struct
221/// @cond Suppress_Doxygen_warnings
222 MFEM_EXPORT
223/// @endcond
225{
226 static const int Dimension = 3;
227 static const int NumVert = 4;
228 static const int NumEdges = 6;
229 static const int Edges[NumEdges][2];
230 static const int NumFaces = 4;
231 static const int FaceTypes[NumFaces];
232 static const int MaxFaceVert = 3;
233 static const int FaceVert[NumFaces][MaxFaceVert];
234 // Upper-triangular part of the local vertex-to-vertex graph.
236 {
237 static const int I[NumVert];
238 static const int J[NumEdges][2]; // {end,edge_idx}
239 };
240
241 static const int NumOrient = 24;
242 static const int Orient[NumOrient][NumVert];
243 static const int InvOrient[NumOrient];
244};
245
246template <> struct
247/// @cond Suppress_Doxygen_warnings
248 MFEM_EXPORT
249/// @endcond
251{
252 static const int Dimension = 3;
253 static const int NumVert = 8;
254 static const int NumEdges = 12;
255 static const int Edges[NumEdges][2];
256 static const int NumFaces = 6;
257 static const int FaceTypes[NumFaces];
258 static const int MaxFaceVert = 4;
259 static const int FaceVert[NumFaces][MaxFaceVert];
260 // Upper-triangular part of the local vertex-to-vertex graph.
262 {
263 static const int I[NumVert];
264 static const int J[NumEdges][2]; // {end,edge_idx}
265 };
266};
267
268template <> struct
269/// @cond Suppress_Doxygen_warnings
270 MFEM_EXPORT
271/// @endcond
273{
274 static const int Dimension = 3;
275 static const int NumVert = 6;
276 static const int NumEdges = 9;
277 static const int Edges[NumEdges][2];
278 static const int NumFaces = 5;
279 static const int FaceTypes[NumFaces];
280 static const int MaxFaceVert = 4;
281 static const int FaceVert[NumFaces][MaxFaceVert];
282 // Upper-triangular part of the local vertex-to-vertex graph.
284 {
285 static const int I[NumVert];
286 static const int J[NumEdges][2]; // {end,edge_idx}
287 };
288};
289
290template <> struct
291/// @cond Suppress_Doxygen_warnings
292 MFEM_EXPORT
293/// @endcond
295{
296 static const int Dimension = 3;
297 static const int NumVert = 5;
298 static const int NumEdges = 8;
299 static const int Edges[NumEdges][2];
300 static const int NumFaces = 5;
301 static const int FaceTypes[NumFaces];
302 static const int MaxFaceVert = 4;
303 static const int FaceVert[NumFaces][MaxFaceVert];
304 // Upper-triangular part of the local vertex-to-vertex graph.
306 {
307 static const int I[NumVert];
308 static const int J[NumEdges][2]; // {end,edge_idx}
309 };
310};
311
312// Defined in fe.cpp to ensure construction after 'mfem::TriangleFE' and
313// `mfem::TetrahedronFE`.
314extern MFEM_EXPORT Geometry Geometries;
315
316
318{
319public:
323 int NumBdrEdges; // at the beginning of RefEdges
324 int Type;
325
326 RefinedGeometry(int NPts, int NRefG, int NRefE, int NBdrE = 0) :
327 RefPts(NPts), RefGeoms(NRefG), RefEdges(NRefE), NumBdrEdges(NBdrE) {}
328};
329
331{
332private:
333 int Type; // Quadrature1D type (ClosedUniform is default)
334 /// Cache of RefinedGeometry for Refine
336 /// Cache of integration rules for EdgeScan
337 /// key: (type, geom, times)
338 std::unordered_map<std::array<int, 3>, std::unique_ptr<IntegrationRule>,
340 SGeom;
342
343 RefinedGeometry *FindInRGeom(Geometry::Type Geom, int Times,
344 int ETimes) const;
345 IntegrationRule *FindInIntPts(Geometry::Type Geom, int NPts) const;
346
347public:
349
350 /// Set the Quadrature1D type of points to use for subdivision.
351 void SetType(int t) { Type = t; }
352 /// Get the Quadrature1D type of points used for subdivision.
353 int GetType() const { return Type; }
354
355 RefinedGeometry *Refine(Geometry::Type Geom, int Times, int ETimes = 1);
356
357 /// Get an integration rule which scans along the r/s/t=0 edges of the element.
358 const IntegrationRule *EdgeScan(Geometry::Type Geom, int NPts1d);
359
360 /// @note This method always uses Quadrature1D::OpenUniform points.
361 const IntegrationRule *RefineInterior(Geometry::Type Geom, int Times);
362
363 /// Get the Refinement level based on number of points
364 static int GetRefinementLevelFromPoints(Geometry::Type Geom, int Npts);
365
366 /// Get the Refinement level based on number of elements
367 static int GetRefinementLevelFromElems(Geometry::Type geom, int Npts);
368
370};
371
372extern MFEM_EXPORT GeometryRefiner GlobGeometryRefiner;
373
374}
375
376#endif
Data type dense matrix using column-major storage.
Definition densemat.hpp:24
RefinedGeometry * Refine(Geometry::Type Geom, int Times, int ETimes=1)
Definition geom.cpp:1136
int GetType() const
Get the Quadrature1D type of points used for subdivision.
Definition geom.hpp:353
void SetType(int t)
Set the Quadrature1D type of points to use for subdivision.
Definition geom.hpp:351
const IntegrationRule * RefineInterior(Geometry::Type Geom, int Times)
Definition geom.cpp:1793
const IntegrationRule * EdgeScan(Geometry::Type Geom, int NPts1d)
Get an integration rule which scans along the r/s/t=0 edges of the element.
Definition geom.cpp:1670
static int GetRefinementLevelFromElems(Geometry::Type geom, int Npts)
Get the Refinement level based on number of elements.
Definition geom.cpp:1972
static int GetRefinementLevelFromPoints(Geometry::Type Geom, int Npts)
Get the Refinement level based on number of points.
Definition geom.cpp:1904
GeometryRefiner(int t=Quadrature1D::ClosedUniform)
Definition geom.hpp:348
static const int NumGeom
Definition geom.hpp:46
const IntegrationPoint & GetCenter(int GeomType) const
Return the center of the given Geometry::Type, GeomType.
Definition geom.hpp:75
static bool IsTensorProduct(Type geom)
Definition geom.hpp:112
int NumBdr(int GeomType) const
Return the number of boundary "faces" of a given Geometry::Type.
Definition geom.hpp:133
const DenseMatrix * GetPerfGeomToGeomJac(int GeomType) const
Definition geom.hpp:104
static Type TensorProductGeometry(int dim)
Definition geom.hpp:117
const DenseMatrix & GetGeomToPerfGeomJac(int GeomType) const
Definition geom.hpp:102
Class for integration point with weight.
Definition intrules.hpp:35
Class for an integration rule - an Array of IntegrationPoint.
Definition intrules.hpp:100
@ ClosedUniform
aka closed Newton-Cotes
Definition intrules.hpp:408
IntegrationRule RefPts
Definition geom.hpp:321
RefinedGeometry(int NPts, int NRefG, int NRefE, int NBdrE=0)
Definition geom.hpp:326
Array< int > RefGeoms
Definition geom.hpp:322
Array< int > RefEdges
Definition geom.hpp:322
int dim
Definition ex24.cpp:53
real_t Volume()
Analytic volume integral over subdomain with positive level-set.
Definition ex38.cpp:109
GeometryRefiner GlobGeometryRefiner
Definition geom.cpp:2014
Geometry Geometries
Definition fe.cpp:49
float real_t
Definition config.hpp:43
Helper class for hashing std::array. Usable in place of std::hash<std::array<T,N>>
Definition hash.hpp:511