MFEM  v4.5.1 Finite element discretization library
geom.hpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2022, 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 Geometry
33 {
34 public:
35  enum Type
36  {
37  INVALID = -1,
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 Geometry::Constants<Geometry::POINT>
130 {
131  static const int Dimension = 0;
132  static const int NumVert = 1;
133
134  static const int NumOrient = 1;
135  static const int Orient[NumOrient][NumVert];
136  static const int InvOrient[NumOrient];
137 };
138
139 template <> struct Geometry::Constants<Geometry::SEGMENT>
140 {
141  static const int Dimension = 1;
142  static const int NumVert = 2;
143  static const int NumEdges = 1;
144  static const int Edges[NumEdges][2];
145
146  static const int NumOrient = 2;
147  static const int Orient[NumOrient][NumVert];
148  static const int InvOrient[NumOrient];
149 };
150
151 template <> struct Geometry::Constants<Geometry::TRIANGLE>
152 {
153  static const int Dimension = 2;
154  static const int NumVert = 3;
155  static const int NumEdges = 3;
156  static const int Edges[NumEdges][2];
157  // Upper-triangular part of the local vertex-to-vertex graph.
158  struct VertToVert
159  {
160  static const int I[NumVert];
161  static const int J[NumEdges][2]; // {end,edge_idx}
162  };
163  static const int NumFaces = 1;
164  static const int FaceVert[NumFaces][NumVert];
165
166  // For a given base tuple v={v0,v1,v2}, the orientation of a permutation
167  // u={u0,u1,u2} of v, is an index 'j' such that u[i]=v[Orient[j][i]].
168  // The static method Mesh::GetTriOrientation, computes the index 'j' of the
169  // permutation that maps the second argument 'test' to the first argument
170  // 'base': test[Orient[j][i]]=base[i].
171  static const int NumOrient = 6;
172  static const int Orient[NumOrient][NumVert];
173  // The inverse of orientation 'j' is InvOrient[j].
174  static const int InvOrient[NumOrient];
175 };
176
177 template <> struct Geometry::Constants<Geometry::SQUARE>
178 {
179  static const int Dimension = 2;
180  static const int NumVert = 4;
181  static const int NumEdges = 4;
182  static const int Edges[NumEdges][2];
183  // Upper-triangular part of the local vertex-to-vertex graph.
184  struct VertToVert
185  {
186  static const int I[NumVert];
187  static const int J[NumEdges][2]; // {end,edge_idx}
188  };
189  static const int NumFaces = 1;
190  static const int FaceVert[NumFaces][NumVert];
191
192  static const int NumOrient = 8;
193  static const int Orient[NumOrient][NumVert];
194  static const int InvOrient[NumOrient];
195 };
196
197 template <> struct Geometry::Constants<Geometry::TETRAHEDRON>
198 {
199  static const int Dimension = 3;
200  static const int NumVert = 4;
201  static const int NumEdges = 6;
202  static const int Edges[NumEdges][2];
203  static const int NumFaces = 4;
204  static const int FaceTypes[NumFaces];
205  static const int MaxFaceVert = 3;
206  static const int FaceVert[NumFaces][MaxFaceVert];
207  // Upper-triangular part of the local vertex-to-vertex graph.
208  struct VertToVert
209  {
210  static const int I[NumVert];
211  static const int J[NumEdges][2]; // {end,edge_idx}
212  };
213
214  static const int NumOrient = 24;
215  static const int Orient[NumOrient][NumVert];
216  static const int InvOrient[NumOrient];
217 };
218
219 template <> struct Geometry::Constants<Geometry::CUBE>
220 {
221  static const int Dimension = 3;
222  static const int NumVert = 8;
223  static const int NumEdges = 12;
224  static const int Edges[NumEdges][2];
225  static const int NumFaces = 6;
226  static const int FaceTypes[NumFaces];
227  static const int MaxFaceVert = 4;
228  static const int FaceVert[NumFaces][MaxFaceVert];
229  // Upper-triangular part of the local vertex-to-vertex graph.
230  struct VertToVert
231  {
232  static const int I[NumVert];
233  static const int J[NumEdges][2]; // {end,edge_idx}
234  };
235 };
236
237 template <> struct Geometry::Constants<Geometry::PRISM>
238 {
239  static const int Dimension = 3;
240  static const int NumVert = 6;
241  static const int NumEdges = 9;
242  static const int Edges[NumEdges][2];
243  static const int NumFaces = 5;
244  static const int FaceTypes[NumFaces];
245  static const int MaxFaceVert = 4;
246  static const int FaceVert[NumFaces][MaxFaceVert];
247  // Upper-triangular part of the local vertex-to-vertex graph.
248  struct VertToVert
249  {
250  static const int I[NumVert];
251  static const int J[NumEdges][2]; // {end,edge_idx}
252  };
253 };
254
255 template <> struct Geometry::Constants<Geometry::PYRAMID>
256 {
257  static const int Dimension = 3;
258  static const int NumVert = 5;
259  static const int NumEdges = 8;
260  static const int Edges[NumEdges][2];
261  static const int NumFaces = 5;
262  static const int FaceTypes[NumFaces];
263  static const int MaxFaceVert = 4;
264  static const int FaceVert[NumFaces][MaxFaceVert];
265  // Upper-triangular part of the local vertex-to-vertex graph.
266  struct VertToVert
267  {
268  static const int I[NumVert];
269  static const int J[NumEdges][2]; // {end,edge_idx}
270  };
271 };
272
273 // Defined in fe.cpp to ensure construction after 'mfem::TriangleFE' and
274 // mfem::TetrahedronFE.
275 extern Geometry Geometries;
276
277
279 {
280 public:
281  int Times, ETimes;
284  int NumBdrEdges; // at the beginning of RefEdges
285  int Type;
286
287  RefinedGeometry(int NPts, int NRefG, int NRefE, int NBdrE = 0) :
288  RefPts(NPts), RefGeoms(NRefG), RefEdges(NRefE), NumBdrEdges(NBdrE) { }
289 };
290
292 {
293 private:
294  int type; // Quadrature1D type (ClosedUniform is default)
297
298  RefinedGeometry *FindInRGeom(Geometry::Type Geom, int Times, int ETimes,
299  int Type);
300  IntegrationRule *FindInIntPts(Geometry::Type Geom, int NPts);
301
302 public:
303  GeometryRefiner();
304
305  /// Set the Quadrature1D type of points to use for subdivision.
306  void SetType(const int t) { type = t; }
307  /// Get the Quadrature1D type of points used for subdivision.
308  int GetType() const { return type; }
309
310  RefinedGeometry *Refine(Geometry::Type Geom, int Times, int ETimes = 1);
311
312  /// @note This method always uses Quadrature1D::OpenUniform points.
313  const IntegrationRule *RefineInterior(Geometry::Type Geom, int Times);
314
315  /// Get the Refinement level based on number of points
316  virtual int GetRefinementLevelFromPoints(Geometry::Type Geom, int Npts);
317
318  /// Get the Refinement level based on number of elements
319  virtual int GetRefinementLevelFromElems(Geometry::Type geom, int Npts);
320
322 };
323
324 extern GeometryRefiner GlobGeometryRefiner;
325
326 }
327
328 #endif
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
static Type TensorProductGeometry(int dim)
Definition: geom.hpp:113
RefinedGeometry(int NPts, int NRefG, int NRefE, int NBdrE=0)
Definition: geom.hpp:287
static const int NumGeom
Definition: geom.hpp:42
void JacToPerfJac(int GeomType, const DenseMatrix &J, DenseMatrix &PJ) const
Definition: geom.cpp:867
static void GetRandomPoint(int GeomType, IntegrationPoint &ip)
Get a random point in the reference element specified by GeomType.
Definition: geom.cpp:285
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
static const double Volume[NumGeom]
Definition: geom.hpp:46
static const int NumEdges[NumGeom]
Definition: geom.hpp:50
Array< int > RefEdges
Definition: geom.hpp:283
const IntegrationPoint & GetCenter(int GeomType)
Return the center of the given Geometry::Type, GeomType.
Definition: geom.hpp:71
const IntegrationRule * GetVertices(int GeomType)
Return an IntegrationRule consisting of all vertices of the given Geometry::Type, GeomType...
Definition: geom.cpp:265
const IntegrationRule * RefineInterior(Geometry::Type Geom, int Times)
Definition: geom.cpp:1586
const DenseMatrix & GetGeomToPerfGeomJac(int GeomType) const
Definition: geom.hpp:98
static const int NumFaces[NumGeom]
Definition: geom.hpp:51
Geometry Geometries
Definition: fe.cpp:49
static const int Dimension[NumGeom]
Definition: geom.hpp:47
void SetType(const int t)
Set the Quadrature1D type of points to use for subdivision.
Definition: geom.hpp:306
static const int NumVerts[NumGeom]
Definition: geom.hpp:49
GeometryRefiner GlobGeometryRefiner
Definition: geom.cpp:1773
IntegrationRule RefPts
Definition: geom.hpp:282
static const int NumBdrArray[NumGeom]
Definition: geom.hpp:44
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
static const char * Name[NumGeom]
Definition: geom.hpp:45
static const int MaxDim
Definition: geom.hpp:43
static bool ProjectPoint(int GeomType, const IntegrationPoint &beg, IntegrationPoint &end)
Project a point end, onto the given Geometry::Type, GeomType.
Definition: geom.cpp:583
void GetPerfPointMat(int GeomType, DenseMatrix &pm)
Definition: geom.cpp:783
int NumBdr(int GeomType)
Return the number of boundary &quot;faces&quot; 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:25
int GetType() const
Get the Quadrature1D type of points used for subdivision.
Definition: geom.hpp:308
virtual int GetRefinementLevelFromPoints(Geometry::Type Geom, int Npts)
Get the Refinement level based on number of points.
Definition: geom.cpp:1665
static const int DimStart[MaxDim+2]
Definition: geom.hpp:48
int dim
Definition: ex24.cpp:53
RefCoord t[3]
static bool CheckPoint(int GeomType, const IntegrationPoint &ip)
Check if the given point is inside the given reference element.
Definition: geom.cpp:405
Array< int > RefGeoms
Definition: geom.hpp:283