MFEM  v4.6.0
Finite element discretization library
intrules.hpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, 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 #include <vector>
19 #include <map>
20 
21 namespace mfem
22 {
23 
24 class KnotVector;
25 class Mesh;
26 
27 /* Classes for IntegrationPoint, IntegrationRule, and container class
28  IntegrationRules. Declares the global variable IntRules */
29 
30 /// Class for integration point with weight
32 {
33 public:
34  double x, y, z, weight;
35  int index;
36 
37  void Init(int const i)
38  {
39  x = y = z = weight = 0.0;
40  index = i;
41  }
42 
43  void Set(const double *p, const int dim)
44  {
45  MFEM_ASSERT(1 <= dim && dim <= 3, "invalid dim: " << dim);
46  x = p[0];
47  if (dim > 1)
48  {
49  y = p[1];
50  if (dim > 2)
51  {
52  z = p[2];
53  }
54  }
55  }
56 
57  void Get(double *p, const int dim) const
58  {
59  MFEM_ASSERT(1 <= dim && dim <= 3, "invalid dim: " << dim);
60  p[0] = x;
61  if (dim > 1)
62  {
63  p[1] = y;
64  if (dim > 2)
65  {
66  p[2] = z;
67  }
68  }
69  }
70 
71  void Set(const double x1, const double x2, const double x3, const double w)
72  { x = x1; y = x2; z = x3; weight = w; }
73 
74  void Set3w(const double *p) { x = p[0]; y = p[1]; z = p[2]; weight = p[3]; }
75 
76  void Set3(const double x1, const double x2, const double x3)
77  { x = x1; y = x2; z = x3; }
78 
79  void Set3(const double *p) { x = p[0]; y = p[1]; z = p[2]; }
80 
81  void Set2w(const double x1, const double x2, const double w)
82  { x = x1; y = x2; weight = w; }
83 
84  void Set2w(const double *p) { x = p[0]; y = p[1]; weight = p[2]; }
85 
86  void Set2(const double x1, const double x2) { x = x1; y = x2; }
87 
88  void Set2(const double *p) { x = p[0]; y = p[1]; }
89 
90  void Set1w(const double x1, const double w) { x = x1; weight = w; }
91 
92  void Set1w(const double *p) { x = p[0]; weight = p[1]; }
93 };
94 
95 /// Class for an integration rule - an Array of IntegrationPoint.
96 class IntegrationRule : public Array<IntegrationPoint>
97 {
98 private:
99  friend class IntegrationRules;
100  int Order = 0;
101  /** @brief The quadrature weights gathered as a contiguous array. Created
102  by request with the method GetWeights(). */
103  mutable Array<double> weights;
104 
105  /// Define n-simplex rule (triangle/tetrahedron for n=2/3) of order (2s+1)
106  void GrundmannMollerSimplexRule(int s, int n = 3);
107 
108  void AddTriMidPoint(const int off, const double weight)
109  { IntPoint(off).Set2w(1./3., 1./3., weight); }
110 
111  void AddTriPoints3(const int off, const double a, const double b,
112  const double weight)
113  {
114  IntPoint(off + 0).Set2w(a, a, weight);
115  IntPoint(off + 1).Set2w(a, b, weight);
116  IntPoint(off + 2).Set2w(b, a, weight);
117  }
118 
119  void AddTriPoints3(const int off, const double a, const double weight)
120  { AddTriPoints3(off, a, 1. - 2.*a, weight); }
121 
122  void AddTriPoints3b(const int off, const double b, const double weight)
123  { AddTriPoints3(off, (1. - b)/2., b, weight); }
124 
125  void AddTriPoints3R(const int off, const double a, const double b,
126  const double c, const double weight)
127  {
128  IntPoint(off + 0).Set2w(a, b, weight);
129  IntPoint(off + 1).Set2w(c, a, weight);
130  IntPoint(off + 2).Set2w(b, c, weight);
131  }
132 
133  void AddTriPoints3R(const int off, const double a, const double b,
134  const double weight)
135  { AddTriPoints3R(off, a, b, 1. - a - b, weight); }
136 
137  void AddTriPoints6(const int off, const double a, const double b,
138  const double c, const double weight)
139  {
140  IntPoint(off + 0).Set2w(a, b, weight);
141  IntPoint(off + 1).Set2w(b, a, weight);
142  IntPoint(off + 2).Set2w(a, c, weight);
143  IntPoint(off + 3).Set2w(c, a, weight);
144  IntPoint(off + 4).Set2w(b, c, weight);
145  IntPoint(off + 5).Set2w(c, b, weight);
146  }
147 
148  void AddTriPoints6(const int off, const double a, const double b,
149  const double weight)
150  { AddTriPoints6(off, a, b, 1. - a - b, weight); }
151 
152  // add the permutations of (a,a,b)
153  void AddTetPoints3(const int off, const double a, const double b,
154  const double weight)
155  {
156  IntPoint(off + 0).Set(a, a, b, weight);
157  IntPoint(off + 1).Set(a, b, a, weight);
158  IntPoint(off + 2).Set(b, a, a, weight);
159  }
160 
161  // add the permutations of (a,b,c)
162  void AddTetPoints6(const int off, const double a, const double b,
163  const double c, const double weight)
164  {
165  IntPoint(off + 0).Set(a, b, c, weight);
166  IntPoint(off + 1).Set(a, c, b, weight);
167  IntPoint(off + 2).Set(b, c, a, weight);
168  IntPoint(off + 3).Set(b, a, c, weight);
169  IntPoint(off + 4).Set(c, a, b, weight);
170  IntPoint(off + 5).Set(c, b, a, weight);
171  }
172 
173  void AddTetMidPoint(const int off, const double weight)
174  { IntPoint(off).Set(0.25, 0.25, 0.25, weight); }
175 
176  // given a, add the permutations of (a,a,a,b), where 3*a + b = 1
177  void AddTetPoints4(const int off, const double a, const double weight)
178  {
179  IntPoint(off).Set(a, a, a, weight);
180  AddTetPoints3(off + 1, a, 1. - 3.*a, weight);
181  }
182 
183  // given b, add the permutations of (a,a,a,b), where 3*a + b = 1
184  void AddTetPoints4b(const int off, const double b, const double weight)
185  {
186  const double a = (1. - b)/3.;
187  IntPoint(off).Set(a, a, a, weight);
188  AddTetPoints3(off + 1, a, b, weight);
189  }
190 
191  // add the permutations of (a,a,b,b), 2*(a + b) = 1
192  void AddTetPoints6(const int off, const double a, const double weight)
193  {
194  const double b = 0.5 - a;
195  AddTetPoints3(off, a, b, weight);
196  AddTetPoints3(off + 3, b, a, weight);
197  }
198 
199  // given (a,b) or (a,c), add the permutations of (a,a,b,c), 2*a + b + c = 1
200  void AddTetPoints12(const int off, const double a, const double bc,
201  const double weight)
202  {
203  const double cb = 1. - 2*a - bc;
204  AddTetPoints3(off, a, bc, weight);
205  AddTetPoints3(off + 3, a, cb, weight);
206  AddTetPoints6(off + 6, a, bc, cb, weight);
207  }
208 
209  // given (b,c), add the permutations of (a,a,b,c), 2*a + b + c = 1
210  void AddTetPoints12bc(const int off, const double b, const double c,
211  const double weight)
212  {
213  const double a = (1. - b - c)/2.;
214  AddTetPoints3(off, a, b, weight);
215  AddTetPoints3(off + 3, a, c, weight);
216  AddTetPoints6(off + 6, a, b, c, weight);
217  }
218 
219 public:
221  Array<IntegrationPoint>() { }
222 
223  /// Construct an integration rule with given number of points
224  explicit IntegrationRule(int NP) :
226  {
227  for (int i = 0; i < this->Size(); i++)
228  {
229  (*this)[i].Init(i);
230  }
231  }
232 
233  /// Sets the indices of each quadrature point on initialization.
234  /** Note that most calls to IntegrationRule::SetSize should be paired with a
235  call to SetPointIndices in order for the indices to be set correctly. */
236  void SetPointIndices();
237 
238  /// Tensor product of two 1D integration rules
240 
241  /// Tensor product of three 1D integration rules
243  IntegrationRule &irz);
244 
245  /// Returns the order of the integration rule
246  int GetOrder() const { return Order; }
247 
248  /** @brief Sets the order of the integration rule. This is only for keeping
249  order information, it does not alter any data in the IntegrationRule. */
250  void SetOrder(const int order) { Order = order; }
251 
252  /// Returns the number of the points in the integration rule
253  int GetNPoints() const { return Size(); }
254 
255  /// Returns a reference to the i-th integration point
256  IntegrationPoint &IntPoint(int i) { return (*this)[i]; }
257 
258  /// Returns a const reference to the i-th integration point
259  const IntegrationPoint &IntPoint(int i) const { return (*this)[i]; }
260 
261  /// Return the quadrature weights in a contiguous array.
262  /** If a contiguous array is not required, the weights can be accessed with
263  a call like this: `IntPoint(i).weight`. */
264  const Array<double> &GetWeights() const;
265 
266  /// @brief Return an integration rule for KnotVector @a kv, defined by
267  /// applying this rule on each knot interval.
269 
270  /// Destroys an IntegrationRule object
272 };
273 
274 /// Class for defining different integration rules on each NURBS patch.
276 {
277 public:
278  /// Construct a rule for each patch, using SetPatchRules1D.
279  NURBSMeshRules(const int numPatches, const int dim_) :
280  patchRules1D(numPatches, dim_),
281  npatches(numPatches), dim(dim_) { }
282 
283  /// Returns a rule for the element.
284  IntegrationRule &GetElementRule(const int elem, const int patch,
285  const int *ijk,
286  Array<const KnotVector*> const& kv,
287  bool & deleteRule) const;
288 
289  /// Add a rule to be used for individual elements. Returns the rule index.
290  std::size_t AddElementRule(IntegrationRule *ir_element)
291  {
292  elementRule.push_back(ir_element);
293  return elementRule.size() - 1;
294  }
295 
296  /// @brief Set the integration rule for the element of the given index. This
297  /// rule is used instead of the rule for the patch containing the element.
298  void SetElementRule(const std::size_t element,
299  const std::size_t elementRuleIndex)
300  {
301  elementToRule[element] = elementRuleIndex;
302  }
303 
304  /// @brief Set 1D integration rules to be used as a tensor product rule on
305  /// the patch with index @a patch. This class takes ownership of these rules.
306  void SetPatchRules1D(const int patch,
307  std::vector<const IntegrationRule*> & ir1D);
308 
309  /// @brief For tensor product rules defined on each patch by
310  /// SetPatchRules1D(), return a pointer to the 1D rule in the specified
311  /// @a dimension.
312  const IntegrationRule* GetPatchRule1D(const int patch,
313  const int dimension) const
314  {
315  return patchRules1D(patch, dimension);
316  }
317 
318  /// @brief For tensor product rules defined on each patch by
319  /// SetPatchRules1D(), return the integration point with index (i,j,k).
320  void GetIntegrationPointFrom1D(const int patch, int i, int j, int k,
321  IntegrationPoint & ip);
322 
323  /// @brief Finalize() must be called before this class can be used for
324  /// assembly. In particular, it defines data used by GetPointElement().
325  void Finalize(Mesh const& mesh);
326 
327  /// @brief For tensor product rules defined on each patch by
328  /// SetPatchRules1D(), returns the index of the element containing
329  /// integration point (i,j,k) for patch index @a patch. Finalize() must be
330  /// called first.
331  int GetPointElement(int patch, int i, int j, int k) const
332  {
333  return pointToElem[patch](i,j,k);
334  }
335 
336  int GetDim() const { return dim; }
337 
338  /// @brief For tensor product rules defined on each patch by
339  /// SetPatchRules1D(), returns an array of knot span indices for each
340  /// integration point in the specified @a dimension.
341  const Array<int>& GetPatchRule1D_KnotSpan(const int patch,
342  const int dimension) const
343  {
344  return patchRules1D_KnotSpan[patch][dimension];
345  }
346 
347  ~NURBSMeshRules();
348 
349 private:
350  /// Tensor-product rules defined on all patches independently.
351  Array2D<const IntegrationRule*> patchRules1D;
352 
353  /// Integration rules defined on elements.
354  std::vector<IntegrationRule*> elementRule;
355 
356  std::map<std::size_t, std::size_t> elementToRule;
357 
358  std::vector<Array3D<int>> pointToElem;
359  std::vector<std::vector<Array<int>>> patchRules1D_KnotSpan;
360 
361  const int npatches;
362  const int dim;
363 };
364 
365 /// A Class that defines 1-D numerical quadrature rules on [0,1].
367 {
368 public:
369  /** @name Methods for calculating quadrature rules.
370  These methods calculate the actual points and weights for the different
371  types of quadrature rules. */
372  ///@{
373  static void GaussLegendre(const int np, IntegrationRule* ir);
374  static void GaussLobatto(const int np, IntegrationRule *ir);
375  static void OpenUniform(const int np, IntegrationRule *ir);
376  static void ClosedUniform(const int np, IntegrationRule *ir);
377  static void OpenHalfUniform(const int np, IntegrationRule *ir);
378  static void ClosedGL(const int np, IntegrationRule *ir);
379  ///@}
380 
381  /// A helper function that will play nice with Poly_1D::OpenPoints and
382  /// Poly_1D::ClosedPoints
383  static void GivePolyPoints(const int np, double *pts, const int type);
384 
385 private:
386  static void CalculateUniformWeights(IntegrationRule *ir, const int type);
387 };
388 
389 /// A class container for 1D quadrature type constants.
391 {
392 public:
393  enum
394  {
395  Invalid = -1,
398  OpenUniform = 2, ///< aka open Newton-Cotes
399  ClosedUniform = 3, ///< aka closed Newton-Cotes
400  OpenHalfUniform = 4, ///< aka "open half" Newton-Cotes
401  ClosedGL = 5 ///< aka closed Gauss Legendre
402  };
403  /** @brief If the Quadrature1D type is not closed return Invalid; otherwise
404  return type. */
405  static int CheckClosed(int type);
406  /** @brief If the Quadrature1D type is not open return Invalid; otherwise
407  return type. */
408  static int CheckOpen(int type);
409 };
410 
411 /// Container class for integration rules
413 {
414 private:
415  /// Taken from the Quadrature1D class anonymous enum
416  /// Determines the type of numerical quadrature used for
417  /// segment, square, and cube geometries
418  const int quad_type;
419 
420  int own_rules, refined;
421 
422  Array<IntegrationRule *> PointIntRules;
423  Array<IntegrationRule *> SegmentIntRules;
424  Array<IntegrationRule *> TriangleIntRules;
425  Array<IntegrationRule *> SquareIntRules;
426  Array<IntegrationRule *> TetrahedronIntRules;
427  Array<IntegrationRule *> PyramidIntRules;
428  Array<IntegrationRule *> PrismIntRules;
429  Array<IntegrationRule *> CubeIntRules;
430 
431  void AllocIntRule(Array<IntegrationRule *> &ir_array, int Order)
432  {
433  if (ir_array.Size() <= Order)
434  {
435  ir_array.SetSize(Order + 1, NULL);
436  }
437  }
438  bool HaveIntRule(Array<IntegrationRule *> &ir_array, int Order)
439  {
440  return (ir_array.Size() > Order && ir_array[Order] != NULL);
441  }
442  int GetSegmentRealOrder(int Order) const
443  {
444  return Order | 1; // valid for all quad_type's
445  }
446 
447  /// The following methods allocate new IntegrationRule objects without
448  /// checking if they already exist. To avoid memory leaks use
449  /// IntegrationRules::Get(int GeomType, int Order) instead.
450  IntegrationRule *GenerateIntegrationRule(int GeomType, int Order);
451  IntegrationRule *PointIntegrationRule(int Order);
452  IntegrationRule *SegmentIntegrationRule(int Order);
453  IntegrationRule *TriangleIntegrationRule(int Order);
454  IntegrationRule *SquareIntegrationRule(int Order);
455  IntegrationRule *TetrahedronIntegrationRule(int Order);
456  IntegrationRule *PyramidIntegrationRule(int Order);
457  IntegrationRule *PrismIntegrationRule(int Order);
458  IntegrationRule *CubeIntegrationRule(int Order);
459 
460  void DeleteIntRuleArray(Array<IntegrationRule *> &ir_array);
461 
462 public:
463  /// Sets initial sizes for the integration rule arrays, but rules
464  /// are defined the first time they are requested with the Get method.
465  explicit IntegrationRules(int Ref = 0,
466  int type = Quadrature1D::GaussLegendre);
467 
468  /// Returns an integration rule for given GeomType and Order.
469  const IntegrationRule &Get(int GeomType, int Order);
470 
471  void Set(int GeomType, int Order, IntegrationRule &IntRule);
472 
473  void SetOwnRules(int o) { own_rules = o; }
474 
475  /// Destroys an IntegrationRules object
477 };
478 
479 /// A global object with all integration rules (defined in intrules.cpp)
480 extern MFEM_EXPORT IntegrationRules IntRules;
481 
482 /// A global object with all refined integration rules
484 
485 }
486 
487 #endif
void Set(const double *p, const int dim)
Definition: intrules.hpp:43
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:253
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:96
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Definition: intrules.cpp:980
static void ClosedUniform(const int np, IntegrationRule *ir)
Definition: intrules.cpp:646
void Set1w(const double *p)
Definition: intrules.hpp:92
static void GivePolyPoints(const int np, double *pts, const int type)
Definition: intrules.cpp:703
~IntegrationRules()
Destroys an IntegrationRules object.
Definition: intrules.cpp:1077
void GetIntegrationPointFrom1D(const int patch, int i, int j, int k, IntegrationPoint &ip)
For tensor product rules defined on each patch by SetPatchRules1D(), return the integration point wit...
Definition: intrules.cpp:1903
aka "open half" Newton-Cotes
Definition: intrules.hpp:400
Container class for integration rules.
Definition: intrules.hpp:412
void Set(const double x1, const double x2, const double x3, const double w)
Definition: intrules.hpp:71
IntegrationRules RefinedIntRules(1, Quadrature1D::GaussLegendre)
A global object with all refined integration rules.
Definition: intrules.hpp:483
void Init(int const i)
Definition: intrules.hpp:37
static void OpenUniform(const int np, IntegrationRule *ir)
Definition: intrules.cpp:630
void Finalize(Mesh const &mesh)
Finalize() must be called before this class can be used for assembly. In particular, it defines data used by GetPointElement().
Definition: intrules.cpp:1925
A class container for 1D quadrature type constants.
Definition: intrules.hpp:390
void SetElementRule(const std::size_t element, const std::size_t elementRuleIndex)
Set the integration rule for the element of the given index. This rule is used instead of the rule fo...
Definition: intrules.hpp:298
IntegrationRule * ApplyToKnotIntervals(KnotVector const &kv) const
Return an integration rule for KnotVector kv, defined by applying this rule on each knot interval...
Definition: intrules.cpp:181
int GetDim() const
Definition: intrules.hpp:336
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Definition: intrules.hpp:480
static void GaussLobatto(const int np, IntegrationRule *ir)
Definition: intrules.cpp:505
const Array< int > & GetPatchRule1D_KnotSpan(const int patch, const int dimension) const
For tensor product rules defined on each patch by SetPatchRules1D(), returns an array of knot span in...
Definition: intrules.hpp:341
void Set3(const double *p)
Definition: intrules.hpp:79
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:256
void SetPatchRules1D(const int patch, std::vector< const IntegrationRule *> &ir1D)
Set 1D integration rules to be used as a tensor product rule on the patch with index patch...
Definition: intrules.cpp:2031
void SetOwnRules(int o)
Definition: intrules.hpp:473
const Array< double > & GetWeights() const
Return the quadrature weights in a contiguous array.
Definition: intrules.cpp:86
double b
Definition: lissajous.cpp:42
int GetOrder() const
Returns the order of the integration rule.
Definition: intrules.hpp:246
void Get(double *p, const int dim) const
Definition: intrules.hpp:57
const IntegrationRule * GetPatchRule1D(const int patch, const int dimension) const
For tensor product rules defined on each patch by SetPatchRules1D(), return a pointer to the 1D rule ...
Definition: intrules.hpp:312
~IntegrationRule()
Destroys an IntegrationRule object.
Definition: intrules.hpp:271
void Set3w(const double *p)
Definition: intrules.hpp:74
void Set2(const double x1, const double x2)
Definition: intrules.hpp:86
const IntegrationPoint & IntPoint(int i) const
Returns a const reference to the i-th integration point.
Definition: intrules.hpp:259
Dynamic 2D array using row-major layout.
Definition: array.hpp:354
void SetPointIndices()
Sets the indices of each quadrature point on initialization.
Definition: intrules.cpp:99
void Set3(const double x1, const double x2, const double x3)
Definition: intrules.hpp:76
NURBSMeshRules(const int numPatches, const int dim_)
Construct a rule for each patch, using SetPatchRules1D.
Definition: intrules.hpp:279
void Set2w(const double *p)
Definition: intrules.hpp:84
static int CheckClosed(int type)
If the Quadrature1D type is not closed return Invalid; otherwise return type.
Definition: intrules.cpp:911
static void OpenHalfUniform(const int np, IntegrationRule *ir)
Definition: intrules.cpp:666
static void ClosedGL(const int np, IntegrationRule *ir)
Definition: intrules.cpp:681
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:687
static int CheckOpen(int type)
If the Quadrature1D type is not open return Invalid; otherwise return type.
Definition: intrules.cpp:923
void SetOrder(const int order)
Sets the order of the integration rule. This is only for keeping order information, it does not alter any data in the IntegrationRule.
Definition: intrules.hpp:250
aka closed Gauss Legendre
Definition: intrules.hpp:401
double a
Definition: lissajous.cpp:41
aka closed Newton-Cotes
Definition: intrules.hpp:399
IntegrationRules(int Ref=0, int type=Quadrature1D::GaussLegendre)
Definition: intrules.cpp:943
Class for integration point with weight.
Definition: intrules.hpp:31
std::size_t AddElementRule(IntegrationRule *ir_element)
Add a rule to be used for individual elements. Returns the rule index.
Definition: intrules.hpp:290
int GetPointElement(int patch, int i, int j, int k) const
For tensor product rules defined on each patch by SetPatchRules1D(), returns the index of the element...
Definition: intrules.hpp:331
int dim
Definition: ex24.cpp:53
aka open Newton-Cotes
Definition: intrules.hpp:398
constexpr int dimension
This example only works in 3D. Kernels for 2D are not implemented.
Definition: hooke.cpp:45
int Size() const
Return the logical size of the array.
Definition: array.hpp:141
void Set2(const double *p)
Definition: intrules.hpp:88
void Set2w(const double x1, const double x2, const double w)
Definition: intrules.hpp:81
void pts(int iphi, int t, double x[])
RefCoord s[3]
IntegrationRule(int NP)
Construct an integration rule with given number of points.
Definition: intrules.hpp:224
IntegrationRule & GetElementRule(const int elem, const int patch, const int *ijk, Array< const KnotVector *> const &kv, bool &deleteRule) const
Returns a rule for the element.
Definition: intrules.cpp:1820
void Set(int GeomType, int Order, IntegrationRule &IntRule)
Definition: intrules.cpp:1031
A Class that defines 1-D numerical quadrature rules on [0,1].
Definition: intrules.hpp:366
void Set1w(const double x1, const double w)
Definition: intrules.hpp:90
static void GaussLegendre(const int np, IntegrationRule *ir)
Definition: intrules.cpp:424
Class for defining different integration rules on each NURBS patch.
Definition: intrules.hpp:275