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_INTRULES 13 #define MFEM_INTRULES 14 15 #include "../config/config.hpp" 16 #include "../general/array.hpp" 17 18 namespace mfem 19 { 20 21 /* Classes for IntegrationPoint, IntegrationRule, and container class 22 IntegrationRules. Declares the global variable IntRules */ 23 24 /// Class for integration point with weight 25 class IntegrationPoint 26 { 27 public: 28 double x, y, z, weight; 29 int index; 30 Init(int const i)31 void Init(int const i) 32 { 33 x = y = z = weight = 0.0; 34 index = i; 35 } 36 Set(const double * p,const int dim)37 void Set(const double *p, const int dim) 38 { 39 MFEM_ASSERT(1 <= dim && dim <= 3, "invalid dim: " << dim); 40 x = p[0]; 41 if (dim > 1) 42 { 43 y = p[1]; 44 if (dim > 2) 45 { 46 z = p[2]; 47 } 48 } 49 } 50 Get(double * p,const int dim) const51 void Get(double *p, const int dim) const 52 { 53 MFEM_ASSERT(1 <= dim && dim <= 3, "invalid dim: " << dim); 54 p[0] = x; 55 if (dim > 1) 56 { 57 p[1] = y; 58 if (dim > 2) 59 { 60 p[2] = z; 61 } 62 } 63 } 64 Set(const double x1,const double x2,const double x3,const double w)65 void Set(const double x1, const double x2, const double x3, const double w) 66 { x = x1; y = x2; z = x3; weight = w; } 67 Set3w(const double * p)68 void Set3w(const double *p) { x = p[0]; y = p[1]; z = p[2]; weight = p[3]; } 69 Set3(const double x1,const double x2,const double x3)70 void Set3(const double x1, const double x2, const double x3) 71 { x = x1; y = x2; z = x3; } 72 Set3(const double * p)73 void Set3(const double *p) { x = p[0]; y = p[1]; z = p[2]; } 74 Set2w(const double x1,const double x2,const double w)75 void Set2w(const double x1, const double x2, const double w) 76 { x = x1; y = x2; weight = w; } 77 Set2w(const double * p)78 void Set2w(const double *p) { x = p[0]; y = p[1]; weight = p[2]; } 79 Set2(const double x1,const double x2)80 void Set2(const double x1, const double x2) { x = x1; y = x2; } 81 Set2(const double * p)82 void Set2(const double *p) { x = p[0]; y = p[1]; } 83 Set1w(const double x1,const double w)84 void Set1w(const double x1, const double w) { x = x1; weight = w; } 85 Set1w(const double * p)86 void Set1w(const double *p) { x = p[0]; weight = p[1]; } 87 }; 88 89 /// Class for an integration rule - an Array of IntegrationPoint. 90 class IntegrationRule : public Array<IntegrationPoint> 91 { 92 private: 93 friend class IntegrationRules; 94 int Order; 95 /** @brief The quadrature weights gathered as a contiguous array. Created 96 by request with the method GetWeights(). */ 97 mutable Array<double> weights; 98 99 /// Define n-simplex rule (triangle/tetrahedron for n=2/3) of order (2s+1) 100 void GrundmannMollerSimplexRule(int s, int n = 3); 101 AddTriMidPoint(const int off,const double weight)102 void AddTriMidPoint(const int off, const double weight) 103 { IntPoint(off).Set2w(1./3., 1./3., weight); } 104 AddTriPoints3(const int off,const double a,const double b,const double weight)105 void AddTriPoints3(const int off, const double a, const double b, 106 const double weight) 107 { 108 IntPoint(off + 0).Set2w(a, a, weight); 109 IntPoint(off + 1).Set2w(a, b, weight); 110 IntPoint(off + 2).Set2w(b, a, weight); 111 } 112 AddTriPoints3(const int off,const double a,const double weight)113 void AddTriPoints3(const int off, const double a, const double weight) 114 { AddTriPoints3(off, a, 1. - 2.*a, weight); } 115 AddTriPoints3b(const int off,const double b,const double weight)116 void AddTriPoints3b(const int off, const double b, const double weight) 117 { AddTriPoints3(off, (1. - b)/2., b, weight); } 118 AddTriPoints3R(const int off,const double a,const double b,const double c,const double weight)119 void AddTriPoints3R(const int off, const double a, const double b, 120 const double c, const double weight) 121 { 122 IntPoint(off + 0).Set2w(a, b, weight); 123 IntPoint(off + 1).Set2w(c, a, weight); 124 IntPoint(off + 2).Set2w(b, c, weight); 125 } 126 AddTriPoints3R(const int off,const double a,const double b,const double weight)127 void AddTriPoints3R(const int off, const double a, const double b, 128 const double weight) 129 { AddTriPoints3R(off, a, b, 1. - a - b, weight); } 130 AddTriPoints6(const int off,const double a,const double b,const double c,const double weight)131 void AddTriPoints6(const int off, const double a, const double b, 132 const double c, const double weight) 133 { 134 IntPoint(off + 0).Set2w(a, b, weight); 135 IntPoint(off + 1).Set2w(b, a, weight); 136 IntPoint(off + 2).Set2w(a, c, weight); 137 IntPoint(off + 3).Set2w(c, a, weight); 138 IntPoint(off + 4).Set2w(b, c, weight); 139 IntPoint(off + 5).Set2w(c, b, weight); 140 } 141 AddTriPoints6(const int off,const double a,const double b,const double weight)142 void AddTriPoints6(const int off, const double a, const double b, 143 const double weight) 144 { AddTriPoints6(off, a, b, 1. - a - b, weight); } 145 146 // add the permutations of (a,a,b) AddTetPoints3(const int off,const double a,const double b,const double weight)147 void AddTetPoints3(const int off, const double a, const double b, 148 const double weight) 149 { 150 IntPoint(off + 0).Set(a, a, b, weight); 151 IntPoint(off + 1).Set(a, b, a, weight); 152 IntPoint(off + 2).Set(b, a, a, weight); 153 } 154 155 // add the permutations of (a,b,c) AddTetPoints6(const int off,const double a,const double b,const double c,const double weight)156 void AddTetPoints6(const int off, const double a, const double b, 157 const double c, const double weight) 158 { 159 IntPoint(off + 0).Set(a, b, c, weight); 160 IntPoint(off + 1).Set(a, c, b, weight); 161 IntPoint(off + 2).Set(b, c, a, weight); 162 IntPoint(off + 3).Set(b, a, c, weight); 163 IntPoint(off + 4).Set(c, a, b, weight); 164 IntPoint(off + 5).Set(c, b, a, weight); 165 } 166 AddTetMidPoint(const int off,const double weight)167 void AddTetMidPoint(const int off, const double weight) 168 { IntPoint(off).Set(0.25, 0.25, 0.25, weight); } 169 170 // given a, add the permutations of (a,a,a,b), where 3*a + b = 1 AddTetPoints4(const int off,const double a,const double weight)171 void AddTetPoints4(const int off, const double a, const double weight) 172 { 173 IntPoint(off).Set(a, a, a, weight); 174 AddTetPoints3(off + 1, a, 1. - 3.*a, weight); 175 } 176 177 // given b, add the permutations of (a,a,a,b), where 3*a + b = 1 AddTetPoints4b(const int off,const double b,const double weight)178 void AddTetPoints4b(const int off, const double b, const double weight) 179 { 180 const double a = (1. - b)/3.; 181 IntPoint(off).Set(a, a, a, weight); 182 AddTetPoints3(off + 1, a, b, weight); 183 } 184 185 // add the permutations of (a,a,b,b), 2*(a + b) = 1 AddTetPoints6(const int off,const double a,const double weight)186 void AddTetPoints6(const int off, const double a, const double weight) 187 { 188 const double b = 0.5 - a; 189 AddTetPoints3(off, a, b, weight); 190 AddTetPoints3(off + 3, b, a, weight); 191 } 192 193 // given (a,b) or (a,c), add the permutations of (a,a,b,c), 2*a + b + c = 1 AddTetPoints12(const int off,const double a,const double bc,const double weight)194 void AddTetPoints12(const int off, const double a, const double bc, 195 const double weight) 196 { 197 const double cb = 1. - 2*a - bc; 198 AddTetPoints3(off, a, bc, weight); 199 AddTetPoints3(off + 3, a, cb, weight); 200 AddTetPoints6(off + 6, a, bc, cb, weight); 201 } 202 203 // given (b,c), add the permutations of (a,a,b,c), 2*a + b + c = 1 AddTetPoints12bc(const int off,const double b,const double c,const double weight)204 void AddTetPoints12bc(const int off, const double b, const double c, 205 const double weight) 206 { 207 const double a = (1. - b - c)/2.; 208 AddTetPoints3(off, a, b, weight); 209 AddTetPoints3(off + 3, a, c, weight); 210 AddTetPoints6(off + 6, a, b, c, weight); 211 } 212 213 public: IntegrationRule()214 IntegrationRule() : 215 Array<IntegrationPoint>(), Order(0) { } 216 217 /// Construct an integration rule with given number of points IntegrationRule(int NP)218 explicit IntegrationRule(int NP) : 219 Array<IntegrationPoint>(NP), Order(0) 220 { 221 for (int i = 0; i < this->Size(); i++) 222 { 223 (*this)[i].Init(i); 224 } 225 } 226 227 /// Sets the indices of each quadrature point on initialization. 228 /** Note that most calls to IntegrationRule::SetSize should be paired with a 229 call to SetPointIndices in order for the indices to be set correctly. */ 230 void SetPointIndices(); 231 232 /// Tensor product of two 1D integration rules 233 IntegrationRule(IntegrationRule &irx, IntegrationRule &iry); 234 235 /// Tensor product of three 1D integration rules 236 IntegrationRule(IntegrationRule &irx, IntegrationRule &iry, 237 IntegrationRule &irz); 238 239 /// Returns the order of the integration rule GetOrder() const240 int GetOrder() const { return Order; } 241 242 /** @brief Sets the order of the integration rule. This is only for keeping 243 order information, it does not alter any data in the IntegrationRule. */ SetOrder(const int order)244 void SetOrder(const int order) { Order = order; } 245 246 /// Returns the number of the points in the integration rule GetNPoints() const247 int GetNPoints() const { return Size(); } 248 249 /// Returns a reference to the i-th integration point IntPoint(int i)250 IntegrationPoint &IntPoint(int i) { return (*this)[i]; } 251 252 /// Returns a const reference to the i-th integration point IntPoint(int i) const253 const IntegrationPoint &IntPoint(int i) const { return (*this)[i]; } 254 255 /// Return the quadrature weights in a contiguous array. 256 /** If a contiguous array is not required, the weights can be accessed with 257 a call like this: `IntPoint(i).weight`. */ 258 const Array<double> &GetWeights() const; 259 260 /// Destroys an IntegrationRule object ~IntegrationRule()261 ~IntegrationRule() { } 262 }; 263 264 /// A Class that defines 1-D numerical quadrature rules on [0,1]. 265 class QuadratureFunctions1D 266 { 267 public: 268 /** @name Methods for calculating quadrature rules. 269 These methods calculate the actual points and weights for the different 270 types of quadrature rules. */ 271 ///@{ 272 static void GaussLegendre(const int np, IntegrationRule* ir); 273 static void GaussLobatto(const int np, IntegrationRule *ir); 274 static void OpenUniform(const int np, IntegrationRule *ir); 275 static void ClosedUniform(const int np, IntegrationRule *ir); 276 static void OpenHalfUniform(const int np, IntegrationRule *ir); 277 static void ClosedGL(const int np, IntegrationRule *ir); 278 ///@} 279 280 /// A helper function that will play nice with Poly_1D::OpenPoints and 281 /// Poly_1D::ClosedPoints 282 static void GivePolyPoints(const int np, double *pts, const int type); 283 284 private: 285 static void CalculateUniformWeights(IntegrationRule *ir, const int type); 286 }; 287 288 /// A class container for 1D quadrature type constants. 289 class Quadrature1D 290 { 291 public: 292 enum 293 { 294 Invalid = -1, 295 GaussLegendre = 0, 296 GaussLobatto = 1, 297 OpenUniform = 2, ///< aka open Newton-Cotes 298 ClosedUniform = 3, ///< aka closed Newton-Cotes 299 OpenHalfUniform = 4, ///< aka "open half" Newton-Cotes 300 ClosedGL = 5 ///< aka closed Gauss Legendre 301 }; 302 /** @brief If the Quadrature1D type is not closed return Invalid; otherwise 303 return type. */ 304 static int CheckClosed(int type); 305 /** @brief If the Quadrature1D type is not open return Invalid; otherwise 306 return type. */ 307 static int CheckOpen(int type); 308 }; 309 310 /// Container class for integration rules 311 class IntegrationRules 312 { 313 private: 314 /// Taken from the Quadrature1D class anonymous enum 315 /// Determines the type of numerical quadrature used for 316 /// segment, square, and cube geometries 317 const int quad_type; 318 319 int own_rules, refined; 320 321 Array<IntegrationRule *> PointIntRules; 322 Array<IntegrationRule *> SegmentIntRules; 323 Array<IntegrationRule *> TriangleIntRules; 324 Array<IntegrationRule *> SquareIntRules; 325 Array<IntegrationRule *> TetrahedronIntRules; 326 Array<IntegrationRule *> PrismIntRules; 327 Array<IntegrationRule *> CubeIntRules; 328 AllocIntRule(Array<IntegrationRule * > & ir_array,int Order)329 void AllocIntRule(Array<IntegrationRule *> &ir_array, int Order) 330 { 331 if (ir_array.Size() <= Order) 332 { 333 ir_array.SetSize(Order + 1, NULL); 334 } 335 } HaveIntRule(Array<IntegrationRule * > & ir_array,int Order)336 bool HaveIntRule(Array<IntegrationRule *> &ir_array, int Order) 337 { 338 return (ir_array.Size() > Order && ir_array[Order] != NULL); 339 } GetSegmentRealOrder(int Order) const340 int GetSegmentRealOrder(int Order) const 341 { 342 return Order | 1; // valid for all quad_type's 343 } 344 345 /// The following methods allocate new IntegrationRule objects without 346 /// checking if they already exist. To avoid memory leaks use 347 /// IntegrationRules::Get(int GeomType, int Order) instead. 348 IntegrationRule *GenerateIntegrationRule(int GeomType, int Order); 349 IntegrationRule *PointIntegrationRule(int Order); 350 IntegrationRule *SegmentIntegrationRule(int Order); 351 IntegrationRule *TriangleIntegrationRule(int Order); 352 IntegrationRule *SquareIntegrationRule(int Order); 353 IntegrationRule *TetrahedronIntegrationRule(int Order); 354 IntegrationRule *PrismIntegrationRule(int Order); 355 IntegrationRule *CubeIntegrationRule(int Order); 356 357 void DeleteIntRuleArray(Array<IntegrationRule *> &ir_array); 358 359 public: 360 /// Sets initial sizes for the integration rule arrays, but rules 361 /// are defined the first time they are requested with the Get method. 362 explicit IntegrationRules(int Ref = 0, 363 int type = Quadrature1D::GaussLegendre); 364 365 /// Returns an integration rule for given GeomType and Order. 366 const IntegrationRule &Get(int GeomType, int Order); 367 368 void Set(int GeomType, int Order, IntegrationRule &IntRule); 369 SetOwnRules(int o)370 void SetOwnRules(int o) { own_rules = o; } 371 372 /// Destroys an IntegrationRules object 373 ~IntegrationRules(); 374 }; 375 376 /// A global object with all integration rules (defined in intrules.cpp) 377 extern IntegrationRules IntRules; 378 379 /// A global object with all refined integration rules 380 extern IntegrationRules RefinedIntRules; 381 382 } 383 384 #endif 385