MFEM v4.7.0 Finite element discretization library
Searching...
No Matches
geom.hpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
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
19namespace mfem
20{
21
22/** Types of domains for integration rules and reference finite elements:
23 Geometry::POINT - a point
24 Geometry::SEGMENT - the interval [0,1]
25 Geometry::TRIANGLE - triangle with vertices (0,0), (1,0), (0,1)
26 Geometry::SQUARE - the unit square (0,1)x(0,1)
27 Geometry::TETRAHEDRON - w/ vert. (0,0,0),(1,0,0),(0,1,0),(0,0,1)
28 Geometry::CUBE - the unit cube
29 Geometry::PRISM - w/ vert. (0,0,0),(1,0,0),(0,1,0),(0,0,1),(1,0,1),(0,1,1)
30 Geometry::PYRAMID - w/ vert. (0,0,0),(1,0,0),(1,1,0),(0,1,0),(0,0,1)
31*/
32class MFEM_EXPORT Geometry
33{
34public:
35 enum Type
36 {
37 INVALID = -1,
38 POINT = 0, SEGMENT, TRIANGLE, SQUARE, TETRAHEDRON, CUBE, PRISM, PYRAMID,
39 NUM_GEOMETRIES
40 };
41
42 static const int NumGeom = NUM_GEOMETRIES;
43 static const int MaxDim = 3;
44 static const int NumBdrArray[NumGeom];
45 static const char *Name[NumGeom];
46 static const real_t Volume[NumGeom];
47 static const int Dimension[NumGeom];
48 static const int DimStart[MaxDim+2]; // including MaxDim+1
49 static const int NumVerts[NumGeom];
50 static const int NumEdges[NumGeom];
51 static const int NumFaces[NumGeom];
52
53 // Structure that holds constants describing the Geometries.
54 template <Type Geom> struct Constants;
55
56private:
57 IntegrationRule *GeomVert[NumGeom];
58 IntegrationPoint GeomCenter[NumGeom];
59 DenseMatrix *GeomToPerfGeomJac[NumGeom];
60 DenseMatrix *PerfGeomToGeomJac[NumGeom];
61
62public:
63 Geometry();
64 ~Geometry();
65
66 /** @brief Return an IntegrationRule consisting of all vertices of the given
67 Geometry::Type, @a GeomType. */
68 const IntegrationRule *GetVertices(int GeomType) const;
69
70 /// Return the center of the given Geometry::Type, @a GeomType.
71 const IntegrationPoint &GetCenter(int GeomType) const
72 { return GeomCenter[GeomType]; }
73
74 /// Get a random point in the reference element specified by @a GeomType.
75 /** This method uses the function rand() for random number generation. */
76 static void GetRandomPoint(int GeomType, IntegrationPoint &ip);
77
78 /// Check if the given point is inside the given reference element.
79 static bool CheckPoint(int GeomType, const IntegrationPoint &ip);
80 /** @brief Check if the given point is inside the given reference element.
81 Overload for fuzzy tolerance. */
82 static bool CheckPoint(int GeomType, const IntegrationPoint &ip, real_t eps);
83
84 /// Project a point @a end, onto the given Geometry::Type, @a GeomType.
85 /** Check if the @a end point is inside the reference element, if not
86 overwrite it with the point on the boundary that lies on the line segment
87 between @a beg and @a end (@a beg must be inside the element). Return
88 true if @a end is inside the element, and false otherwise. */
89 static bool ProjectPoint(int GeomType, const IntegrationPoint &beg,
90 IntegrationPoint &end);
91
92 /// Project a point @a ip, onto the given Geometry::Type, @a GeomType.
93 /** If @a ip is outside the element, replace it with the point on the
94 boundary that is closest to the original @a ip and return false;
95 otherwise, return true without changing @a ip. */
96 static bool ProjectPoint(int GeomType, IntegrationPoint &ip);
97
98 const DenseMatrix &GetGeomToPerfGeomJac(int GeomType) const
99 { return *GeomToPerfGeomJac[GeomType]; }
100 const DenseMatrix *GetPerfGeomToGeomJac(int GeomType) const
101 { return PerfGeomToGeomJac[GeomType]; }
102 void GetPerfPointMat(int GeomType, DenseMatrix &pm) const;
103 void JacToPerfJac(int GeomType, const DenseMatrix &J,
104 DenseMatrix &PJ) const;
105
106 /// Returns true if the given @a geom is of tensor-product type (i.e. if geom
107 /// is a segment, quadrilateral, or hexahedron), returns false otherwise.
108 static bool IsTensorProduct(Type geom)
109 { return geom == SEGMENT || geom == SQUARE || geom == CUBE; }
110
111 /// Returns the Geometry::Type corresponding to a tensor-product of the
112 /// given dimension.
114 {
115 switch (dim)
116 {
117 case 0: return POINT;
118 case 1: return SEGMENT;
119 case 2: return SQUARE;
120 case 3: return CUBE;
121 default: MFEM_ABORT("Invalid dimension."); return INVALID;
122 }
123 }
124
125 /// Return the inverse of the given orientation for the specified geometry type.
126 static int GetInverseOrientation(Type geom_type, int orientation);
127
128 /// Return the number of boundary "faces" of a given Geometry::Type.
129 int NumBdr(int GeomType) const { return NumBdrArray[GeomType]; }
130};
131
132template <> struct
133/// @cond Suppress_Doxygen_warnings
134 MFEM_EXPORT
135/// @endcond
137{
138 static const int Dimension = 0;
139 static const int NumVert = 1;
140
141 static const int NumOrient = 1;
142 static const int Orient[NumOrient][NumVert];
143 static const int InvOrient[NumOrient];
144};
145
146template <> struct
147/// @cond Suppress_Doxygen_warnings
148 MFEM_EXPORT
149/// @endcond
151{
152 static const int Dimension = 1;
153 static const int NumVert = 2;
154 static const int NumEdges = 1;
155 static const int Edges[NumEdges][2];
156
157 static const int NumOrient = 2;
158 static const int Orient[NumOrient][NumVert];
159 static const int InvOrient[NumOrient];
160};
161
162template <> struct
163/// @cond Suppress_Doxygen_warnings
164 MFEM_EXPORT
165/// @endcond
167{
168 static const int Dimension = 2;
169 static const int NumVert = 3;
170 static const int NumEdges = 3;
171 static const int Edges[NumEdges][2];
172 // Upper-triangular part of the local vertex-to-vertex graph.
174 {
175 static const int I[NumVert];
176 static const int J[NumEdges][2]; // {end,edge_idx}
177 };
178 static const int NumFaces = 1;
179 static const int FaceVert[NumFaces][NumVert];
180
181 // For a given base tuple v={v0,v1,v2}, the orientation of a permutation
182 // u={u0,u1,u2} of v, is an index 'j' such that u[i]=v[Orient[j][i]].
183 // The static method Mesh::GetTriOrientation, computes the index 'j' of the
184 // permutation that maps the second argument 'test' to the first argument
185 // 'base': test[Orient[j][i]]=base[i].
186 static const int NumOrient = 6;
187 static const int Orient[NumOrient][NumVert];
188 // The inverse of orientation 'j' is InvOrient[j].
189 static const int InvOrient[NumOrient];
190};
191
192template <> struct
193/// @cond Suppress_Doxygen_warnings
194 MFEM_EXPORT
195/// @endcond
197{
198 static const int Dimension = 2;
199 static const int NumVert = 4;
200 static const int NumEdges = 4;
201 static const int Edges[NumEdges][2];
202 // Upper-triangular part of the local vertex-to-vertex graph.
204 {
205 static const int I[NumVert];
206 static const int J[NumEdges][2]; // {end,edge_idx}
207 };
208 static const int NumFaces = 1;
209 static const int FaceVert[NumFaces][NumVert];
210
211 static const int NumOrient = 8;
212 static const int Orient[NumOrient][NumVert];
213 static const int InvOrient[NumOrient];
214};
215
216template <> struct
217/// @cond Suppress_Doxygen_warnings
218 MFEM_EXPORT
219/// @endcond
221{
222 static const int Dimension = 3;
223 static const int NumVert = 4;
224 static const int NumEdges = 6;
225 static const int Edges[NumEdges][2];
226 static const int NumFaces = 4;
227 static const int FaceTypes[NumFaces];
228 static const int MaxFaceVert = 3;
229 static const int FaceVert[NumFaces][MaxFaceVert];
230 // Upper-triangular part of the local vertex-to-vertex graph.
232 {
233 static const int I[NumVert];
234 static const int J[NumEdges][2]; // {end,edge_idx}
235 };
236
237 static const int NumOrient = 24;
238 static const int Orient[NumOrient][NumVert];
239 static const int InvOrient[NumOrient];
240};
241
242template <> struct
243/// @cond Suppress_Doxygen_warnings
244 MFEM_EXPORT
245/// @endcond
247{
248 static const int Dimension = 3;
249 static const int NumVert = 8;
250 static const int NumEdges = 12;
251 static const int Edges[NumEdges][2];
252 static const int NumFaces = 6;
253 static const int FaceTypes[NumFaces];
254 static const int MaxFaceVert = 4;
255 static const int FaceVert[NumFaces][MaxFaceVert];
256 // Upper-triangular part of the local vertex-to-vertex graph.
258 {
259 static const int I[NumVert];
260 static const int J[NumEdges][2]; // {end,edge_idx}
261 };
262};
263
264template <> struct
265/// @cond Suppress_Doxygen_warnings
266 MFEM_EXPORT
267/// @endcond
269{
270 static const int Dimension = 3;
271 static const int NumVert = 6;
272 static const int NumEdges = 9;
273 static const int Edges[NumEdges][2];
274 static const int NumFaces = 5;
275 static const int FaceTypes[NumFaces];
276 static const int MaxFaceVert = 4;
277 static const int FaceVert[NumFaces][MaxFaceVert];
278 // Upper-triangular part of the local vertex-to-vertex graph.
280 {
281 static const int I[NumVert];
282 static const int J[NumEdges][2]; // {end,edge_idx}
283 };
284};
285
286template <> struct
287/// @cond Suppress_Doxygen_warnings
288 MFEM_EXPORT
289/// @endcond
291{
292 static const int Dimension = 3;
293 static const int NumVert = 5;
294 static const int NumEdges = 8;
295 static const int Edges[NumEdges][2];
296 static const int NumFaces = 5;
297 static const int FaceTypes[NumFaces];
298 static const int MaxFaceVert = 4;
299 static const int FaceVert[NumFaces][MaxFaceVert];
300 // Upper-triangular part of the local vertex-to-vertex graph.
302 {
303 static const int I[NumVert];
304 static const int J[NumEdges][2]; // {end,edge_idx}
305 };
306};
307
308// Defined in fe.cpp to ensure construction after 'mfem::TriangleFE' and
309// mfem::TetrahedronFE.
310extern MFEM_EXPORT Geometry Geometries;
311
312
314{
315public:
319 int NumBdrEdges; // at the beginning of RefEdges
320 int Type;
321
322 RefinedGeometry(int NPts, int NRefG, int NRefE, int NBdrE = 0) :
323 RefPts(NPts), RefGeoms(NRefG), RefEdges(NRefE), NumBdrEdges(NBdrE) {}
324};
325
327{
328private:
329 int Type; // Quadrature1D type (ClosedUniform is default)
332
333 RefinedGeometry *FindInRGeom(Geometry::Type Geom, int Times,
334 int ETimes) const;
335 IntegrationRule *FindInIntPts(Geometry::Type Geom, int NPts) const;
336
337public:
339
340 /// Set the Quadrature1D type of points to use for subdivision.
341 void SetType(int t) { Type = t; }
342 /// Get the Quadrature1D type of points used for subdivision.
343 int GetType() const { return Type; }
344
345 RefinedGeometry *Refine(Geometry::Type Geom, int Times, int ETimes = 1);
346
347 /// @note This method always uses Quadrature1D::OpenUniform points.
348 const IntegrationRule *RefineInterior(Geometry::Type Geom, int Times);
349
350 /// Get the Refinement level based on number of points
351 static int GetRefinementLevelFromPoints(Geometry::Type Geom, int Npts);
352
353 /// Get the Refinement level based on number of elements
354 static int GetRefinementLevelFromElems(Geometry::Type geom, int Npts);
355
357};
358
359extern MFEM_EXPORT GeometryRefiner GlobGeometryRefiner;
360
361}
362
363#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:343
void SetType(int t)
Set the Quadrature1D type of points to use for subdivision.
Definition geom.hpp:341
const IntegrationRule * RefineInterior(Geometry::Type Geom, int Times)
Definition geom.cpp:1670
static int GetRefinementLevelFromElems(Geometry::Type geom, int Npts)
Get the Refinement level based on number of elements.
Definition geom.cpp:1849
static int GetRefinementLevelFromPoints(Geometry::Type Geom, int Npts)
Get the Refinement level based on number of points.
Definition geom.cpp:1781
Definition geom.hpp:338
static const int NumGeom
Definition geom.hpp:42
const IntegrationPoint & GetCenter(int GeomType) const
Return the center of the given Geometry::Type, GeomType.
Definition geom.hpp:71
static bool IsTensorProduct(Type geom)
Definition geom.hpp:108
int NumBdr(int GeomType) const
Return the number of boundary "faces" of a given Geometry::Type.
Definition geom.hpp:129
const DenseMatrix * GetPerfGeomToGeomJac(int GeomType) const
Definition geom.hpp:100
static Type TensorProductGeometry(int dim)
Definition geom.hpp:113
const DenseMatrix & GetGeomToPerfGeomJac(int GeomType) const
Definition geom.hpp:98
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:402
IntegrationRule RefPts
Definition geom.hpp:317
RefinedGeometry(int NPts, int NRefG, int NRefE, int NBdrE=0)
Definition geom.hpp:322
Array< int > RefGeoms
Definition geom.hpp:318
Array< int > RefEdges
Definition geom.hpp:318
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:1891
Geometry Geometries
Definition fe.cpp:49
float real_t
Definition config.hpp:43
RefCoord t[3]