MFEM  v4.1.0 Finite element discretization library
geom.hpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2020, 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"
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 */
31 class Geometry
32 {
33 public:
34  enum Type
35  {
36  INVALID = -1,
39  };
40
41  static const int NumGeom = NUM_GEOMETRIES;
42  static const int MaxDim = 3;
43  static const int NumBdrArray[NumGeom];
44  static const char *Name[NumGeom];
45  static const double Volume[NumGeom];
46  static const int Dimension[NumGeom];
47  static const int DimStart[MaxDim+2]; // including MaxDim+1
48  static const int NumVerts[NumGeom];
49  static const int NumEdges[NumGeom];
50  static const int NumFaces[NumGeom];
51
52  // Structure that holds constants describing the Geometries.
53  template <Type Geom> struct Constants;
54
55 private:
56  IntegrationRule *GeomVert[NumGeom];
57  IntegrationPoint GeomCenter[NumGeom];
58  DenseMatrix *GeomToPerfGeomJac[NumGeom];
59  DenseMatrix *PerfGeomToGeomJac[NumGeom];
60
61 public:
62  Geometry();
63  ~Geometry();
64
65  /** @brief Return an IntegrationRule consisting of all vertices of the given
66  Geometry::Type, @a GeomType. */
67  const IntegrationRule *GetVertices(int GeomType);
68
69  /// Return the center of the given Geometry::Type, @a GeomType.
70  const IntegrationPoint &GetCenter(int GeomType)
71  { return GeomCenter[GeomType]; }
72
73  /// Get a random point in the reference element specified by @a GeomType.
74  /** This method uses the function rand() for random number generation. */
75  static void GetRandomPoint(int GeomType, IntegrationPoint &ip);
76
77  /// Check if the given point is inside the given reference element.
78  static bool CheckPoint(int GeomType, const IntegrationPoint &ip);
79  /** @brief Check if the given point is inside the given reference element.
80  Overload for fuzzy tolerance. */
81  static bool CheckPoint(int GeomType, const IntegrationPoint &ip, double eps);
82
83  /// Project a point @a end, onto the given Geometry::Type, @a GeomType.
84  /** Check if the @a end point is inside the reference element, if not
85  overwrite it with the point on the boundary that lies on the line segment
86  between @a beg and @a end (@a beg must be inside the element). Return
87  true if @a end is inside the element, and false otherwise. */
88  static bool ProjectPoint(int GeomType, const IntegrationPoint &beg,
89  IntegrationPoint &end);
90
91  /// Project a point @a ip, onto the given Geometry::Type, @a GeomType.
92  /** If @a ip is outside the element, replace it with the point on the
93  boundary that is closest to the original @a ip and return false;
94  otherwise, return true without changing @a ip. */
95  static bool ProjectPoint(int GeomType, IntegrationPoint &ip);
96
97  const DenseMatrix &GetGeomToPerfGeomJac(int GeomType) const
98  { return *GeomToPerfGeomJac[GeomType]; }
100  { return PerfGeomToGeomJac[GeomType]; }
101  void GetPerfPointMat(int GeomType, DenseMatrix &pm);
102  void JacToPerfJac(int GeomType, const DenseMatrix &J,
103  DenseMatrix &PJ) const;
104
105  /// Return the number of boundary "faces" of a given Geometry::Type.
106  int NumBdr(int GeomType) { return NumBdrArray[GeomType]; }
107 };
108
109 template <> struct Geometry::Constants<Geometry::POINT>
110 {
111  static const int Dimension = 0;
112  static const int NumVert = 1;
113
114  static const int NumOrient = 1;
115  static const int Orient[NumOrient][NumVert];
116  static const int InvOrient[NumOrient];
117 };
118
119 template <> struct Geometry::Constants<Geometry::SEGMENT>
120 {
121  static const int Dimension = 1;
122  static const int NumVert = 2;
123  static const int NumEdges = 1;
124  static const int Edges[NumEdges][2];
125
126  static const int NumOrient = 2;
127  static const int Orient[NumOrient][NumVert];
128  static const int InvOrient[NumOrient];
129 };
130
131 template <> struct Geometry::Constants<Geometry::TRIANGLE>
132 {
133  static const int Dimension = 2;
134  static const int NumVert = 3;
135  static const int NumEdges = 3;
136  static const int Edges[NumEdges][2];
137  // Upper-triangular part of the local vertex-to-vertex graph.
138  struct VertToVert
139  {
140  static const int I[NumVert];
141  static const int J[NumEdges][2]; // {end,edge_idx}
142  };
143  static const int NumFaces = 1;
144  static const int FaceVert[NumFaces][NumVert];
145
146  // For a given base tuple v={v0,v1,v2}, the orientation of a permutation
147  // u={u0,u1,u2} of v, is an index 'j' such that u[i]=v[Orient[j][i]].
148  // The static method Mesh::GetTriOrientation, computes the index 'j' of the
149  // permutation that maps the second argument 'test' to the first argument
150  // 'base': test[Orient[j][i]]=base[i].
151  static const int NumOrient = 6;
152  static const int Orient[NumOrient][NumVert];
153  // The inverse of orientation 'j' is InvOrient[j].
154  static const int InvOrient[NumOrient];
155 };
156
157 template <> struct Geometry::Constants<Geometry::SQUARE>
158 {
159  static const int Dimension = 2;
160  static const int NumVert = 4;
161  static const int NumEdges = 4;
162  static const int Edges[NumEdges][2];
163  // Upper-triangular part of the local vertex-to-vertex graph.
164  struct VertToVert
165  {
166  static const int I[NumVert];
167  static const int J[NumEdges][2]; // {end,edge_idx}
168  };
169  static const int NumFaces = 1;
170  static const int FaceVert[NumFaces][NumVert];
171
172  static const int NumOrient = 8;
173  static const int Orient[NumOrient][NumVert];
174  static const int InvOrient[NumOrient];
175 };
176
177 template <> struct Geometry::Constants<Geometry::TETRAHEDRON>
178 {
179  static const int Dimension = 3;
180  static const int NumVert = 4;
181  static const int NumEdges = 6;
182  static const int Edges[NumEdges][2];
183  static const int NumFaces = 4;
184  static const int FaceTypes[NumFaces];
185  static const int MaxFaceVert = 3;
186  static const int FaceVert[NumFaces][MaxFaceVert];
187  // Upper-triangular part of the local vertex-to-vertex graph.
188  struct VertToVert
189  {
190  static const int I[NumVert];
191  static const int J[NumEdges][2]; // {end,edge_idx}
192  };
193 };
194
195 template <> struct Geometry::Constants<Geometry::CUBE>
196 {
197  static const int Dimension = 3;
198  static const int NumVert = 8;
199  static const int NumEdges = 12;
200  static const int Edges[NumEdges][2];
201  static const int NumFaces = 6;
202  static const int FaceTypes[NumFaces];
203  static const int MaxFaceVert = 4;
204  static const int FaceVert[NumFaces][MaxFaceVert];
205  // Upper-triangular part of the local vertex-to-vertex graph.
206  struct VertToVert
207  {
208  static const int I[NumVert];
209  static const int J[NumEdges][2]; // {end,edge_idx}
210  };
211 };
212
213 template <> struct Geometry::Constants<Geometry::PRISM>
214 {
215  static const int Dimension = 3;
216  static const int NumVert = 6;
217  static const int NumEdges = 9;
218  static const int Edges[NumEdges][2];
219  static const int NumFaces = 5;
220  static const int FaceTypes[NumFaces];
221  static const int MaxFaceVert = 4;
222  static const int FaceVert[NumFaces][MaxFaceVert];
223  // Upper-triangular part of the local vertex-to-vertex graph.
224  struct VertToVert
225  {
226  static const int I[NumVert];
227  static const int J[NumEdges][2]; // {end,edge_idx}
228  };
229 };
230
231 // Defined in fe.cpp to ensure construction after 'mfem::WedgeFE'.
232 extern Geometry Geometries;
233
234
236 {
237 public:
238  int Times, ETimes;
241  int NumBdrEdges; // at the beginning of RefEdges
242  int Type;
243
244  RefinedGeometry(int NPts, int NRefG, int NRefE, int NBdrE = 0) :
245  RefPts(NPts), RefGeoms(NRefG), RefEdges(NRefE), NumBdrEdges(NBdrE) { }
246 };
247
249 {
250 private:
251  int type; // Quadrature1D type (ClosedUniform is default)
254
255  RefinedGeometry *FindInRGeom(Geometry::Type Geom, int Times, int ETimes,
256  int Type);
257  IntegrationRule *FindInIntPts(Geometry::Type Geom, int NPts);
258
259 public:
260  GeometryRefiner();
261
262  /// Set the Quadrature1D type of points to use for subdivision.
263  void SetType(const int t) { type = t; }
264  /// Get the Quadrature1D type of points used for subdivision.
265  int GetType() const { return type; }
266
267  RefinedGeometry *Refine(Geometry::Type Geom, int Times, int ETimes = 1);
268
269  /// @note This method always uses Quadrature1D::OpenUniform points.
270  const IntegrationRule *RefineInterior(Geometry::Type Geom, int Times);
271
273 };
274
275 extern GeometryRefiner GlobGeometryRefiner;
276
277 }
278
279 #endif
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
RefinedGeometry(int NPts, int NRefG, int NRefE, int NBdrE=0)
Definition: geom.hpp:244
static const int NumGeom
Definition: geom.hpp:41
void JacToPerfJac(int GeomType, const DenseMatrix &J, DenseMatrix &PJ) const
Definition: geom.cpp:737
static void GetRandomPoint(int GeomType, IntegrationPoint &ip)
Get a random point in the reference element specified by GeomType.
Definition: geom.cpp:247
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
static const double Volume[NumGeom]
Definition: geom.hpp:45
static const int NumEdges[NumGeom]
Definition: geom.hpp:49
Array< int > RefEdges
Definition: geom.hpp:240
const IntegrationPoint & GetCenter(int GeomType)
Return the center of the given Geometry::Type, GeomType.
Definition: geom.hpp:70
const IntegrationRule * GetVertices(int GeomType)
Return an IntegrationRule consisting of all vertices of the given Geometry::Type, GeomType...
Definition: geom.cpp:228
const IntegrationRule * RefineInterior(Geometry::Type Geom, int Times)
Definition: geom.cpp:1316
const DenseMatrix & GetGeomToPerfGeomJac(int GeomType) const
Definition: geom.hpp:97
static const int NumFaces[NumGeom]
Definition: geom.hpp:50
Geometry Geometries
Definition: fe.cpp:12638
static const int Dimension[NumGeom]
Definition: geom.hpp:46
void SetType(const int t)
Set the Quadrature1D type of points to use for subdivision.
Definition: geom.hpp:263
static const int NumVerts[NumGeom]
Definition: geom.hpp:48
GeometryRefiner GlobGeometryRefiner
Definition: geom.cpp:1395
IntegrationRule RefPts
Definition: geom.hpp:239
static const int NumBdrArray[NumGeom]
Definition: geom.hpp:43
RefinedGeometry * Refine(Geometry::Type Geom, int Times, int ETimes=1)
Definition: geom.cpp:925
static const char * Name[NumGeom]
Definition: geom.hpp:44
static const int MaxDim
Definition: geom.hpp:42
static bool ProjectPoint(int GeomType, const IntegrationPoint &beg, IntegrationPoint &end)
Project a point end, onto the given Geometry::Type, GeomType.
Definition: geom.cpp:511
void GetPerfPointMat(int GeomType, DenseMatrix &pm)
Definition: geom.cpp:664
int NumBdr(int GeomType)
Return the number of boundary &quot;faces&quot; of a given Geometry::Type.
Definition: geom.hpp:106
DenseMatrix * GetPerfGeomToGeomJac(int GeomType)
Definition: geom.hpp:99
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:265
static const int DimStart[MaxDim+2]
Definition: geom.hpp:47
static bool CheckPoint(int GeomType, const IntegrationPoint &ip)
Check if the given point is inside the given reference element.
Definition: geom.cpp:348
Array< int > RefGeoms
Definition: geom.hpp:240