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