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