MFEM  v4.5.2
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 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
26 {
27 public:
28  double x, y, z, weight;
29  int index;
30 
31  void Init(int const i)
32  {
33  x = y = z = weight = 0.0;
34  index = i;
35  }
36 
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 
51  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 
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 
68  void Set3w(const double *p) { x = p[0]; y = p[1]; z = p[2]; weight = p[3]; }
69 
70  void Set3(const double x1, const double x2, const double x3)
71  { x = x1; y = x2; z = x3; }
72 
73  void Set3(const double *p) { x = p[0]; y = p[1]; z = p[2]; }
74 
75  void Set2w(const double x1, const double x2, const double w)
76  { x = x1; y = x2; weight = w; }
77 
78  void Set2w(const double *p) { x = p[0]; y = p[1]; weight = p[2]; }
79 
80  void Set2(const double x1, const double x2) { x = x1; y = x2; }
81 
82  void Set2(const double *p) { x = p[0]; y = p[1]; }
83 
84  void Set1w(const double x1, const double w) { x = x1; weight = w; }
85 
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 
102  void AddTriMidPoint(const int off, const double weight)
103  { IntPoint(off).Set2w(1./3., 1./3., weight); }
104 
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 
113  void AddTriPoints3(const int off, const double a, const double weight)
114  { AddTriPoints3(off, a, 1. - 2.*a, weight); }
115 
116  void AddTriPoints3b(const int off, const double b, const double weight)
117  { AddTriPoints3(off, (1. - b)/2., b, weight); }
118 
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 
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 
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 
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)
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)
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 
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
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
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
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
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
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:
215  Array<IntegrationPoint>(), Order(0) { }
216 
217  /// Construct an integration rule with given number of points
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
234 
235  /// Tensor product of three 1D integration rules
237  IntegrationRule &irz);
238 
239  /// Returns the order of the integration rule
240  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. */
244  void SetOrder(const int order) { Order = order; }
245 
246  /// Returns the number of the points in the integration rule
247  int GetNPoints() const { return Size(); }
248 
249  /// Returns a reference to the i-th integration point
250  IntegrationPoint &IntPoint(int i) { return (*this)[i]; }
251 
252  /// Returns a const reference to the i-th integration point
253  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
262 };
263 
264 /// A Class that defines 1-D numerical quadrature rules on [0,1].
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.
290 {
291 public:
292  enum
293  {
294  Invalid = -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
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 *> PyramidIntRules;
327  Array<IntegrationRule *> PrismIntRules;
328  Array<IntegrationRule *> CubeIntRules;
329 
330  void AllocIntRule(Array<IntegrationRule *> &ir_array, int Order)
331  {
332  if (ir_array.Size() <= Order)
333  {
334  ir_array.SetSize(Order + 1, NULL);
335  }
336  }
337  bool HaveIntRule(Array<IntegrationRule *> &ir_array, int Order)
338  {
339  return (ir_array.Size() > Order && ir_array[Order] != NULL);
340  }
341  int GetSegmentRealOrder(int Order) const
342  {
343  return Order | 1; // valid for all quad_type's
344  }
345 
346  /// The following methods allocate new IntegrationRule objects without
347  /// checking if they already exist. To avoid memory leaks use
348  /// IntegrationRules::Get(int GeomType, int Order) instead.
349  IntegrationRule *GenerateIntegrationRule(int GeomType, int Order);
350  IntegrationRule *PointIntegrationRule(int Order);
351  IntegrationRule *SegmentIntegrationRule(int Order);
352  IntegrationRule *TriangleIntegrationRule(int Order);
353  IntegrationRule *SquareIntegrationRule(int Order);
354  IntegrationRule *TetrahedronIntegrationRule(int Order);
355  IntegrationRule *PyramidIntegrationRule(int Order);
356  IntegrationRule *PrismIntegrationRule(int Order);
357  IntegrationRule *CubeIntegrationRule(int Order);
358 
359  void DeleteIntRuleArray(Array<IntegrationRule *> &ir_array);
360 
361 public:
362  /// Sets initial sizes for the integration rule arrays, but rules
363  /// are defined the first time they are requested with the Get method.
364  explicit IntegrationRules(int Ref = 0,
365  int type = Quadrature1D::GaussLegendre);
366 
367  /// Returns an integration rule for given GeomType and Order.
368  const IntegrationRule &Get(int GeomType, int Order);
369 
370  void Set(int GeomType, int Order, IntegrationRule &IntRule);
371 
372  void SetOwnRules(int o) { own_rules = o; }
373 
374  /// Destroys an IntegrationRules object
376 };
377 
378 /// A global object with all integration rules (defined in intrules.cpp)
379 extern MFEM_EXPORT IntegrationRules IntRules;
380 
381 /// A global object with all refined integration rules
383 
384 }
385 
386 #endif
void Set(const double *p, const int dim)
Definition: intrules.hpp:37
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:247
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Definition: intrules.cpp:923
static void ClosedUniform(const int np, IntegrationRule *ir)
Definition: intrules.cpp:592
void Set1w(const double *p)
Definition: intrules.hpp:86
static void GivePolyPoints(const int np, double *pts, const int type)
Definition: intrules.cpp:646
~IntegrationRules()
Destroys an IntegrationRules object.
Definition: intrules.cpp:1016
aka "open half" Newton-Cotes
Definition: intrules.hpp:299
Container class for integration rules.
Definition: intrules.hpp:311
void Set(const double x1, const double x2, const double x3, const double w)
Definition: intrules.hpp:65
IntegrationRules RefinedIntRules(1, Quadrature1D::GaussLegendre)
A global object with all refined integration rules.
Definition: intrules.hpp:382
void Init(int const i)
Definition: intrules.hpp:31
static void OpenUniform(const int np, IntegrationRule *ir)
Definition: intrules.cpp:577
A class container for 1D quadrature type constants.
Definition: intrules.hpp:289
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Definition: intrules.hpp:379
static void GaussLobatto(const int np, IntegrationRule *ir)
Definition: intrules.cpp:454
void Set3(const double *p)
Definition: intrules.hpp:73
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:250
void SetOwnRules(int o)
Definition: intrules.hpp:372
const Array< double > & GetWeights() const
Return the quadrature weights in a contiguous array.
Definition: intrules.cpp:83
double b
Definition: lissajous.cpp:42
int GetOrder() const
Returns the order of the integration rule.
Definition: intrules.hpp:240
void Get(double *p, const int dim) const
Definition: intrules.hpp:51
~IntegrationRule()
Destroys an IntegrationRule object.
Definition: intrules.hpp:261
void Set3w(const double *p)
Definition: intrules.hpp:68
void Set2(const double x1, const double x2)
Definition: intrules.hpp:80
const IntegrationPoint & IntPoint(int i) const
Returns a const reference to the i-th integration point.
Definition: intrules.hpp:253
void SetPointIndices()
Sets the indices of each quadrature point on initialization.
Definition: intrules.cpp:96
void Set3(const double x1, const double x2, const double x3)
Definition: intrules.hpp:70
void Set2w(const double *p)
Definition: intrules.hpp:78
static int CheckClosed(int type)
If the Quadrature1D type is not closed return Invalid; otherwise return type.
Definition: intrules.cpp:854
static void OpenHalfUniform(const int np, IntegrationRule *ir)
Definition: intrules.cpp:611
static void ClosedGL(const int np, IntegrationRule *ir)
Definition: intrules.cpp:625
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:684
static int CheckOpen(int type)
If the Quadrature1D type is not open return Invalid; otherwise return type.
Definition: intrules.cpp:866
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:244
aka closed Gauss Legendre
Definition: intrules.hpp:300
double a
Definition: lissajous.cpp:41
aka closed Newton-Cotes
Definition: intrules.hpp:298
IntegrationRules(int Ref=0, int type=Quadrature1D::GaussLegendre)
Definition: intrules.cpp:886
Class for integration point with weight.
Definition: intrules.hpp:25
int dim
Definition: ex24.cpp:53
aka open Newton-Cotes
Definition: intrules.hpp:297
int Size() const
Return the logical size of the array.
Definition: array.hpp:141
void Set2(const double *p)
Definition: intrules.hpp:82
void Set2w(const double x1, const double x2, const double w)
Definition: intrules.hpp:75
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:218
void Set(int GeomType, int Order, IntegrationRule &IntRule)
Definition: intrules.cpp:970
A Class that defines 1-D numerical quadrature rules on [0,1].
Definition: intrules.hpp:265
void Set1w(const double x1, const double w)
Definition: intrules.hpp:84
static void GaussLegendre(const int np, IntegrationRule *ir)
Definition: intrules.cpp:374