MFEM  v4.6.0
Finite element discretization library
geom.hpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, 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"
16 #include "../linalg/densemat.hpp"
17 #include "intrules.hpp"
18 
19 namespace 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 */
32 class MFEM_EXPORT Geometry
33 {
34 public:
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 double 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 
56 private:
57  IntegrationRule *GeomVert[NumGeom];
58  IntegrationPoint GeomCenter[NumGeom];
59  DenseMatrix *GeomToPerfGeomJac[NumGeom];
60  DenseMatrix *PerfGeomToGeomJac[NumGeom];
61 
62 public:
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);
69 
70  /// Return the center of the given Geometry::Type, @a GeomType.
71  const IntegrationPoint &GetCenter(int GeomType)
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, double 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]; }
101  { return PerfGeomToGeomJac[GeomType]; }
102  void GetPerfPointMat(int GeomType, DenseMatrix &pm);
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 number of boundary "faces" of a given Geometry::Type.
126  int NumBdr(int GeomType) { return NumBdrArray[GeomType]; }
127 };
128 
129 template <> struct
130 /// @cond Suppress_Doxygen_warnings
131  MFEM_EXPORT
132 /// @endcond
134 {
135  static const int Dimension = 0;
136  static const int NumVert = 1;
137 
138  static const int NumOrient = 1;
139  static const int Orient[NumOrient][NumVert];
140  static const int InvOrient[NumOrient];
141 };
142 
143 template <> struct
144 /// @cond Suppress_Doxygen_warnings
145  MFEM_EXPORT
146 /// @endcond
148 {
149  static const int Dimension = 1;
150  static const int NumVert = 2;
151  static const int NumEdges = 1;
152  static const int Edges[NumEdges][2];
153 
154  static const int NumOrient = 2;
155  static const int Orient[NumOrient][NumVert];
156  static const int InvOrient[NumOrient];
157 };
158 
159 template <> struct
160 /// @cond Suppress_Doxygen_warnings
161  MFEM_EXPORT
162 /// @endcond
164 {
165  static const int Dimension = 2;
166  static const int NumVert = 3;
167  static const int NumEdges = 3;
168  static const int Edges[NumEdges][2];
169  // Upper-triangular part of the local vertex-to-vertex graph.
170  struct VertToVert
171  {
172  static const int I[NumVert];
173  static const int J[NumEdges][2]; // {end,edge_idx}
174  };
175  static const int NumFaces = 1;
176  static const int FaceVert[NumFaces][NumVert];
177 
178  // For a given base tuple v={v0,v1,v2}, the orientation of a permutation
179  // u={u0,u1,u2} of v, is an index 'j' such that u[i]=v[Orient[j][i]].
180  // The static method Mesh::GetTriOrientation, computes the index 'j' of the
181  // permutation that maps the second argument 'test' to the first argument
182  // 'base': test[Orient[j][i]]=base[i].
183  static const int NumOrient = 6;
184  static const int Orient[NumOrient][NumVert];
185  // The inverse of orientation 'j' is InvOrient[j].
186  static const int InvOrient[NumOrient];
187 };
188 
189 template <> struct
190 /// @cond Suppress_Doxygen_warnings
191  MFEM_EXPORT
192 /// @endcond
194 {
195  static const int Dimension = 2;
196  static const int NumVert = 4;
197  static const int NumEdges = 4;
198  static const int Edges[NumEdges][2];
199  // Upper-triangular part of the local vertex-to-vertex graph.
200  struct VertToVert
201  {
202  static const int I[NumVert];
203  static const int J[NumEdges][2]; // {end,edge_idx}
204  };
205  static const int NumFaces = 1;
206  static const int FaceVert[NumFaces][NumVert];
207 
208  static const int NumOrient = 8;
209  static const int Orient[NumOrient][NumVert];
210  static const int InvOrient[NumOrient];
211 };
212 
213 template <> struct
214 /// @cond Suppress_Doxygen_warnings
215  MFEM_EXPORT
216 /// @endcond
218 {
219  static const int Dimension = 3;
220  static const int NumVert = 4;
221  static const int NumEdges = 6;
222  static const int Edges[NumEdges][2];
223  static const int NumFaces = 4;
224  static const int FaceTypes[NumFaces];
225  static const int MaxFaceVert = 3;
226  static const int FaceVert[NumFaces][MaxFaceVert];
227  // Upper-triangular part of the local vertex-to-vertex graph.
228  struct VertToVert
229  {
230  static const int I[NumVert];
231  static const int J[NumEdges][2]; // {end,edge_idx}
232  };
233 
234  static const int NumOrient = 24;
235  static const int Orient[NumOrient][NumVert];
236  static const int InvOrient[NumOrient];
237 };
238 
239 template <> struct
240 /// @cond Suppress_Doxygen_warnings
241  MFEM_EXPORT
242 /// @endcond
244 {
245  static const int Dimension = 3;
246  static const int NumVert = 8;
247  static const int NumEdges = 12;
248  static const int Edges[NumEdges][2];
249  static const int NumFaces = 6;
250  static const int FaceTypes[NumFaces];
251  static const int MaxFaceVert = 4;
252  static const int FaceVert[NumFaces][MaxFaceVert];
253  // Upper-triangular part of the local vertex-to-vertex graph.
254  struct VertToVert
255  {
256  static const int I[NumVert];
257  static const int J[NumEdges][2]; // {end,edge_idx}
258  };
259 };
260 
261 template <> struct
262 /// @cond Suppress_Doxygen_warnings
263  MFEM_EXPORT
264 /// @endcond
266 {
267  static const int Dimension = 3;
268  static const int NumVert = 6;
269  static const int NumEdges = 9;
270  static const int Edges[NumEdges][2];
271  static const int NumFaces = 5;
272  static const int FaceTypes[NumFaces];
273  static const int MaxFaceVert = 4;
274  static const int FaceVert[NumFaces][MaxFaceVert];
275  // Upper-triangular part of the local vertex-to-vertex graph.
276  struct VertToVert
277  {
278  static const int I[NumVert];
279  static const int J[NumEdges][2]; // {end,edge_idx}
280  };
281 };
282 
283 template <> struct
284 /// @cond Suppress_Doxygen_warnings
285  MFEM_EXPORT
286 /// @endcond
288 {
289  static const int Dimension = 3;
290  static const int NumVert = 5;
291  static const int NumEdges = 8;
292  static const int Edges[NumEdges][2];
293  static const int NumFaces = 5;
294  static const int FaceTypes[NumFaces];
295  static const int MaxFaceVert = 4;
296  static const int FaceVert[NumFaces][MaxFaceVert];
297  // Upper-triangular part of the local vertex-to-vertex graph.
298  struct VertToVert
299  {
300  static const int I[NumVert];
301  static const int J[NumEdges][2]; // {end,edge_idx}
302  };
303 };
304 
305 // Defined in fe.cpp to ensure construction after 'mfem::TriangleFE' and
306 // `mfem::TetrahedronFE`.
307 extern MFEM_EXPORT Geometry Geometries;
308 
309 
311 {
312 public:
313  int Times, ETimes;
316  int NumBdrEdges; // at the beginning of RefEdges
317  int Type;
318 
319  RefinedGeometry(int NPts, int NRefG, int NRefE, int NBdrE = 0) :
320  RefPts(NPts), RefGeoms(NRefG), RefEdges(NRefE), NumBdrEdges(NBdrE) { }
321 };
322 
324 {
325 private:
326  int type; // Quadrature1D type (ClosedUniform is default)
329 
330  RefinedGeometry *FindInRGeom(Geometry::Type Geom, int Times, int ETimes,
331  int Type);
332  IntegrationRule *FindInIntPts(Geometry::Type Geom, int NPts);
333 
334 public:
335  GeometryRefiner();
336 
337  /// Set the Quadrature1D type of points to use for subdivision.
338  void SetType(const int t) { type = t; }
339  /// Get the Quadrature1D type of points used for subdivision.
340  int GetType() const { return type; }
341 
342  RefinedGeometry *Refine(Geometry::Type Geom, int Times, int ETimes = 1);
343 
344  /// @note This method always uses Quadrature1D::OpenUniform points.
345  const IntegrationRule *RefineInterior(Geometry::Type Geom, int Times);
346 
347  /// Get the Refinement level based on number of points
348  virtual int GetRefinementLevelFromPoints(Geometry::Type Geom, int Npts);
349 
350  /// Get the Refinement level based on number of elements
351  virtual int GetRefinementLevelFromElems(Geometry::Type geom, int Npts);
352 
354 };
355 
356 extern MFEM_EXPORT GeometryRefiner GlobGeometryRefiner;
357 
358 }
359 
360 #endif
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:96
static Type TensorProductGeometry(int dim)
Definition: geom.hpp:113
RefinedGeometry(int NPts, int NRefG, int NRefE, int NBdrE=0)
Definition: geom.hpp:319
static const int NumGeom
Definition: geom.hpp:42
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
Array< int > RefEdges
Definition: geom.hpp:315
const IntegrationPoint & GetCenter(int GeomType)
Return the center of the given Geometry::Type, GeomType.
Definition: geom.hpp:71
const IntegrationRule * RefineInterior(Geometry::Type Geom, int Times)
Definition: geom.cpp:1586
Geometry Geometries
Definition: fe.cpp:49
GeometryRefiner GlobGeometryRefiner
Definition: geom.cpp:1773
void SetType(const int t)
Set the Quadrature1D type of points to use for subdivision.
Definition: geom.hpp:338
IntegrationRule RefPts
Definition: geom.hpp:314
RefinedGeometry * Refine(Geometry::Type Geom, int Times, int ETimes=1)
Definition: geom.cpp:1099
virtual int GetRefinementLevelFromElems(Geometry::Type geom, int Npts)
Get the Refinement level based on number of elements.
Definition: geom.cpp:1732
int NumBdr(int GeomType)
Return the number of boundary "faces" of a given Geometry::Type.
Definition: geom.hpp:126
DenseMatrix * GetPerfGeomToGeomJac(int GeomType)
Definition: geom.hpp:100
static bool IsTensorProduct(Type geom)
Definition: geom.hpp:108
Class for integration point with weight.
Definition: intrules.hpp:31
virtual int GetRefinementLevelFromPoints(Geometry::Type Geom, int Npts)
Get the Refinement level based on number of points.
Definition: geom.cpp:1665
int dim
Definition: ex24.cpp:53
const DenseMatrix & GetGeomToPerfGeomJac(int GeomType) const
Definition: geom.hpp:98
RefCoord t[3]
int GetType() const
Get the Quadrature1D type of points used for subdivision.
Definition: geom.hpp:340
Array< int > RefGeoms
Definition: geom.hpp:315