1 // Copyright (c) 2010-2021, 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 */
31 class Geometry
32 {
33 public:
34    enum Type
35    {
36       INVALID = -1,
37       POINT = 0, SEGMENT, TRIANGLE, SQUARE, TETRAHEDRON, CUBE, PRISM,
38       NUM_GEOMETRIES
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.
GetCenter(int 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 
GetGeomToPerfGeomJac(int GeomType) const97    const DenseMatrix &GetGeomToPerfGeomJac(int GeomType) const
98    { return *GeomToPerfGeomJac[GeomType]; }
GetPerfGeomToGeomJac(int GeomType)99    DenseMatrix *GetPerfGeomToGeomJac(int 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    /// Returns true if the given @a geom is of tensor-product type (i.e. if geom
106    /// is a segment, quadrilateral, or hexahedron), returns false otherwise.
IsTensorProduct(Type geom)107    static bool IsTensorProduct(Type geom)
108    { return geom == SEGMENT || geom == SQUARE || geom == CUBE; }
109 
110    /// Returns the Geometry::Type corresponding to a tensor-product of the
111    /// given dimension.
TensorProductGeometry(int dim)112    static Type TensorProductGeometry(int dim)
113    {
114       switch (dim)
115       {
116          case 0: return POINT;
117          case 1: return SEGMENT;
118          case 2: return SQUARE;
119          case 3: return CUBE;
120          default: MFEM_ABORT("Invalid dimension."); return INVALID;
121       }
122    }
123 
124    /// Return the number of boundary "faces" of a given Geometry::Type.
NumBdr(int GeomType)125    int NumBdr(int GeomType) { return NumBdrArray[GeomType]; }
126 };
127 
128 template <> struct Geometry::Constants<Geometry::POINT>
129 {
130    static const int Dimension = 0;
131    static const int NumVert = 1;
132 
133    static const int NumOrient = 1;
134    static const int Orient[NumOrient][NumVert];
135    static const int InvOrient[NumOrient];
136 };
137 
138 template <> struct Geometry::Constants<Geometry::SEGMENT>
139 {
140    static const int Dimension = 1;
141    static const int NumVert = 2;
142    static const int NumEdges = 1;
143    static const int Edges[NumEdges][2];
144 
145    static const int NumOrient = 2;
146    static const int Orient[NumOrient][NumVert];
147    static const int InvOrient[NumOrient];
148 };
149 
150 template <> struct Geometry::Constants<Geometry::TRIANGLE>
151 {
152    static const int Dimension = 2;
153    static const int NumVert = 3;
154    static const int NumEdges = 3;
155    static const int Edges[NumEdges][2];
156    // Upper-triangular part of the local vertex-to-vertex graph.
157    struct VertToVert
158    {
159       static const int I[NumVert];
160       static const int J[NumEdges][2]; // {end,edge_idx}
161    };
162    static const int NumFaces = 1;
163    static const int FaceVert[NumFaces][NumVert];
164 
165    // For a given base tuple v={v0,v1,v2}, the orientation of a permutation
166    // u={u0,u1,u2} of v, is an index 'j' such that u[i]=v[Orient[j][i]].
167    // The static method Mesh::GetTriOrientation, computes the index 'j' of the
168    // permutation that maps the second argument 'test' to the first argument
169    // 'base': test[Orient[j][i]]=base[i].
170    static const int NumOrient = 6;
171    static const int Orient[NumOrient][NumVert];
172    // The inverse of orientation 'j' is InvOrient[j].
173    static const int InvOrient[NumOrient];
174 };
175 
176 template <> struct Geometry::Constants<Geometry::SQUARE>
177 {
178    static const int Dimension = 2;
179    static const int NumVert = 4;
180    static const int NumEdges = 4;
181    static const int Edges[NumEdges][2];
182    // Upper-triangular part of the local vertex-to-vertex graph.
183    struct VertToVert
184    {
185       static const int I[NumVert];
186       static const int J[NumEdges][2]; // {end,edge_idx}
187    };
188    static const int NumFaces = 1;
189    static const int FaceVert[NumFaces][NumVert];
190 
191    static const int NumOrient = 8;
192    static const int Orient[NumOrient][NumVert];
193    static const int InvOrient[NumOrient];
194 };
195 
196 template <> struct Geometry::Constants<Geometry::TETRAHEDRON>
197 {
198    static const int Dimension = 3;
199    static const int NumVert = 4;
200    static const int NumEdges = 6;
201    static const int Edges[NumEdges][2];
202    static const int NumFaces = 4;
203    static const int FaceTypes[NumFaces];
204    static const int MaxFaceVert = 3;
205    static const int FaceVert[NumFaces][MaxFaceVert];
206    // Upper-triangular part of the local vertex-to-vertex graph.
207    struct VertToVert
208    {
209       static const int I[NumVert];
210       static const int J[NumEdges][2]; // {end,edge_idx}
211    };
212 
213    static const int NumOrient = 24;
214    static const int Orient[NumOrient][NumVert];
215    static const int InvOrient[NumOrient];
216 };
217 
218 template <> struct Geometry::Constants<Geometry::CUBE>
219 {
220    static const int Dimension = 3;
221    static const int NumVert = 8;
222    static const int NumEdges = 12;
223    static const int Edges[NumEdges][2];
224    static const int NumFaces = 6;
225    static const int FaceTypes[NumFaces];
226    static const int MaxFaceVert = 4;
227    static const int FaceVert[NumFaces][MaxFaceVert];
228    // Upper-triangular part of the local vertex-to-vertex graph.
229    struct VertToVert
230    {
231       static const int I[NumVert];
232       static const int J[NumEdges][2]; // {end,edge_idx}
233    };
234 };
235 
236 template <> struct Geometry::Constants<Geometry::PRISM>
237 {
238    static const int Dimension = 3;
239    static const int NumVert = 6;
240    static const int NumEdges = 9;
241    static const int Edges[NumEdges][2];
242    static const int NumFaces = 5;
243    static const int FaceTypes[NumFaces];
244    static const int MaxFaceVert = 4;
245    static const int FaceVert[NumFaces][MaxFaceVert];
246    // Upper-triangular part of the local vertex-to-vertex graph.
247    struct VertToVert
248    {
249       static const int I[NumVert];
250       static const int J[NumEdges][2]; // {end,edge_idx}
251    };
252 };
253 
254 // Defined in fe.cpp to ensure construction after 'mfem::WedgeFE'.
255 extern Geometry Geometries;
256 
257 
258 class RefinedGeometry
259 {
260 public:
261    int Times, ETimes;
262    IntegrationRule RefPts;
263    Array<int> RefGeoms, RefEdges;
264    int NumBdrEdges; // at the beginning of RefEdges
265    int Type;
266 
RefinedGeometry(int NPts,int NRefG,int NRefE,int NBdrE=0)267    RefinedGeometry(int NPts, int NRefG, int NRefE, int NBdrE = 0) :
268       RefPts(NPts), RefGeoms(NRefG), RefEdges(NRefE), NumBdrEdges(NBdrE) { }
269 };
270 
271 class GeometryRefiner
272 {
273 private:
274    int type; // Quadrature1D type (ClosedUniform is default)
275    Array<RefinedGeometry *> RGeom[Geometry::NumGeom];
276    Array<IntegrationRule *> IntPts[Geometry::NumGeom];
277 
278    RefinedGeometry *FindInRGeom(Geometry::Type Geom, int Times, int ETimes,
279                                 int Type);
280    IntegrationRule *FindInIntPts(Geometry::Type Geom, int NPts);
281 
282 public:
283    GeometryRefiner();
284 
285    /// Set the Quadrature1D type of points to use for subdivision.
SetType(const int t)286    void SetType(const int t) { type = t; }
287    /// Get the Quadrature1D type of points used for subdivision.
GetType() const288    int GetType() const { return type; }
289 
290    RefinedGeometry *Refine(Geometry::Type Geom, int Times, int ETimes = 1);
291 
292    /// @note This method always uses Quadrature1D::OpenUniform points.
293    const IntegrationRule *RefineInterior(Geometry::Type Geom, int Times);
294 
295    /// Get the Refinement level based on number of points
296    virtual int GetRefinementLevelFromPoints(Geometry::Type Geom, int Npts);
297 
298    /// Get the Refinement level based on number of elements
299    virtual int GetRefinementLevelFromElems(Geometry::Type geom, int Npts);
300 
301    ~GeometryRefiner();
302 };
303 
304 extern GeometryRefiner GlobGeometryRefiner;
305 
306 }
307 
308 #endif
309