MFEM  v4.0 Finite element discretization library
intrules.hpp
Go to the documentation of this file.
1 // Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
2 // the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
3 // reserved. See file COPYRIGHT for details.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability see http://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
10 // Software Foundation) version 2.1 dated February 1999.
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
30  void Init() { x = y = z = weight = 0.0; }
31
32  void Set(const double *p, const int dim)
33  {
34  MFEM_ASSERT(1 <= dim && dim <= 3, "invalid dim: " << dim);
35  x = p[0];
36  if (dim > 1)
37  {
38  y = p[1];
39  if (dim > 2)
40  {
41  z = p[2];
42  }
43  }
44  }
45
46  void Get(double *p, const int dim) const
47  {
48  MFEM_ASSERT(1 <= dim && dim <= 3, "invalid dim: " << dim);
49  p[0] = x;
50  if (dim > 1)
51  {
52  p[1] = y;
53  if (dim > 2)
54  {
55  p[2] = z;
56  }
57  }
58  }
59
60  void Set(const double x1, const double x2, const double x3, const double w)
61  { x = x1; y = x2; z = x3; weight = w; }
62
63  void Set3w(const double *p) { x = p[0]; y = p[1]; z = p[2]; weight = p[3]; }
64
65  void Set3(const double x1, const double x2, const double x3)
66  { x = x1; y = x2; z = x3; }
67
68  void Set3(const double *p) { x = p[0]; y = p[1]; z = p[2]; }
69
70  void Set2w(const double x1, const double x2, const double w)
71  { x = x1; y = x2; weight = w; }
72
73  void Set2w(const double *p) { x = p[0]; y = p[1]; weight = p[2]; }
74
75  void Set2(const double x1, const double x2) { x = x1; y = x2; }
76
77  void Set2(const double *p) { x = p[0]; y = p[1]; }
78
79  void Set1w(const double x1, const double w) { x = x1; weight = w; }
80
81  void Set1w(const double *p) { x = p[0]; weight = p[1]; }
82 };
83
84 /// Class for an integration rule - an Array of IntegrationPoint.
85 class IntegrationRule : public Array<IntegrationPoint>
86 {
87 private:
88  friend class IntegrationRules;
89  int Order;
90  /** @brief The quadrature weights gathered as a contiguous array. Created
91  by request with the method GetWeights(). */
92  mutable Array<double> weights;
93
94  /// Define n-simplex rule (triangle/tetrahedron for n=2/3) of order (2s+1)
95  void GrundmannMollerSimplexRule(int s, int n = 3);
96
97  void AddTriMidPoint(const int off, const double weight)
98  { IntPoint(off).Set2w(1./3., 1./3., weight); }
99
100  void AddTriPoints3(const int off, const double a, const double b,
101  const double weight)
102  {
103  IntPoint(off + 0).Set2w(a, a, weight);
104  IntPoint(off + 1).Set2w(a, b, weight);
105  IntPoint(off + 2).Set2w(b, a, weight);
106  }
107
108  void AddTriPoints3(const int off, const double a, const double weight)
109  { AddTriPoints3(off, a, 1. - 2.*a, weight); }
110
111  void AddTriPoints3b(const int off, const double b, const double weight)
112  { AddTriPoints3(off, (1. - b)/2., b, weight); }
113
114  void AddTriPoints3R(const int off, const double a, const double b,
115  const double c, const double weight)
116  {
117  IntPoint(off + 0).Set2w(a, b, weight);
118  IntPoint(off + 1).Set2w(c, a, weight);
119  IntPoint(off + 2).Set2w(b, c, weight);
120  }
121
122  void AddTriPoints3R(const int off, const double a, const double b,
123  const double weight)
124  { AddTriPoints3R(off, a, b, 1. - a - b, weight); }
125
126  void AddTriPoints6(const int off, const double a, const double b,
127  const double c, const double weight)
128  {
129  IntPoint(off + 0).Set2w(a, b, weight);
130  IntPoint(off + 1).Set2w(b, a, weight);
131  IntPoint(off + 2).Set2w(a, c, weight);
132  IntPoint(off + 3).Set2w(c, a, weight);
133  IntPoint(off + 4).Set2w(b, c, weight);
134  IntPoint(off + 5).Set2w(c, b, weight);
135  }
136
137  void AddTriPoints6(const int off, const double a, const double b,
138  const double weight)
139  { AddTriPoints6(off, a, b, 1. - a - b, weight); }
140
141  // add the permutations of (a,a,b)
142  void AddTetPoints3(const int off, const double a, const double b,
143  const double weight)
144  {
145  IntPoint(off + 0).Set(a, a, b, weight);
146  IntPoint(off + 1).Set(a, b, a, weight);
147  IntPoint(off + 2).Set(b, a, a, weight);
148  }
149
150  // add the permutations of (a,b,c)
151  void AddTetPoints6(const int off, const double a, const double b,
152  const double c, const double weight)
153  {
154  IntPoint(off + 0).Set(a, b, c, weight);
155  IntPoint(off + 1).Set(a, c, b, weight);
156  IntPoint(off + 2).Set(b, c, a, weight);
157  IntPoint(off + 3).Set(b, a, c, weight);
158  IntPoint(off + 4).Set(c, a, b, weight);
159  IntPoint(off + 5).Set(c, b, a, weight);
160  }
161
162  void AddTetMidPoint(const int off, const double weight)
163  { IntPoint(off).Set(0.25, 0.25, 0.25, weight); }
164
165  // given a, add the permutations of (a,a,a,b), where 3*a + b = 1
166  void AddTetPoints4(const int off, const double a, const double weight)
167  {
168  IntPoint(off).Set(a, a, a, weight);
169  AddTetPoints3(off + 1, a, 1. - 3.*a, weight);
170  }
171
172  // given b, add the permutations of (a,a,a,b), where 3*a + b = 1
173  void AddTetPoints4b(const int off, const double b, const double weight)
174  {
175  const double a = (1. - b)/3.;
176  IntPoint(off).Set(a, a, a, weight);
177  AddTetPoints3(off + 1, a, b, weight);
178  }
179
180  // add the permutations of (a,a,b,b), 2*(a + b) = 1
181  void AddTetPoints6(const int off, const double a, const double weight)
182  {
183  const double b = 0.5 - a;
185  AddTetPoints3(off + 3, b, a, weight);
186  }
187
188  // given (a,b) or (a,c), add the permutations of (a,a,b,c), 2*a + b + c = 1
189  void AddTetPoints12(const int off, const double a, const double bc,
190  const double weight)
191  {
192  const double cb = 1. - 2*a - bc;
194  AddTetPoints3(off + 3, a, cb, weight);
195  AddTetPoints6(off + 6, a, bc, cb, weight);
196  }
197
198  // given (b,c), add the permutations of (a,a,b,c), 2*a + b + c = 1
199  void AddTetPoints12bc(const int off, const double b, const double c,
200  const double weight)
201  {
202  const double a = (1. - b - c)/2.;
204  AddTetPoints3(off + 3, a, c, weight);
205  AddTetPoints6(off + 6, a, b, c, weight);
206  }
207
208 public:
210  Array<IntegrationPoint>(), Order(0) { }
211
212  /// Construct an integration rule with given number of points
213  explicit IntegrationRule(int NP) :
214  Array<IntegrationPoint>(NP), Order(0)
215  {
216  for (int i = 0; i < this->Size(); i++)
217  {
218  (*this)[i].Init();
219  }
220  }
221
222  /// Tensor product of two 1D integration rules
224
225  /// Tensor product of three 1D integration rules
227  IntegrationRule &irz);
228
229  /// Returns the order of the integration rule
230  int GetOrder() const { return Order; }
231
232  /** @brief Sets the order of the integration rule. This is only for keeping
233  order information, it does not alter any data in the IntegrationRule. */
234  void SetOrder(const int order) { Order = order; }
235
236  /// Returns the number of the points in the integration rule
237  int GetNPoints() const { return Size(); }
238
239  /// Returns a reference to the i-th integration point
240  IntegrationPoint &IntPoint(int i) { return (*this)[i]; }
241
242  /// Returns a const reference to the i-th integration point
243  const IntegrationPoint &IntPoint(int i) const { return (*this)[i]; }
244
245  /// Return the quadrature weights in a contiguous array.
246  /** If a contiguous array is not required, the weights can be accessed with
247  a call like this: IntPoint(i).weight. */
248  const Array<double> &GetWeights() const;
249
250  /// Destroys an IntegrationRule object
252 };
253
254 /// A Class that defines 1-D numerical quadrature rules on [0,1].
256 {
257 public:
258  /** @name Methods for calculating quadrature rules.
259  These methods calculate the actual points and weights for the different
260  types of quadrature rules. */
261  ///@{
262  void GaussLegendre(const int np, IntegrationRule* ir);
263  void GaussLobatto(const int np, IntegrationRule *ir);
264  void OpenUniform(const int np, IntegrationRule *ir);
265  void ClosedUniform(const int np, IntegrationRule *ir);
266  void OpenHalfUniform(const int np, IntegrationRule *ir);
267  ///@}
268
269  /// A helper function that will play nice with Poly_1D::OpenPoints and
270  /// Poly_1D::ClosedPoints
271  void GivePolyPoints(const int np, double *pts, const int type);
272
273 private:
274  void CalculateUniformWeights(IntegrationRule *ir, const int type);
275 };
276
277 /// A class container for 1D quadrature type constants.
279 {
280 public:
281  enum
282  {
283  Invalid = -1,
286  OpenUniform = 2, ///< aka open Newton-Cotes
287  ClosedUniform = 3, ///< aka closed Newton-Cotes
288  OpenHalfUniform = 4 ///< aka "open half" Newton-Cotes
289  };
290  /** @brief If the Quadrature1D type is not closed return Invalid; otherwise
291  return type. */
292  static int CheckClosed(int type);
293  /** @brief If the Quadrature1D type is not open return Invalid; otherwise
294  return type. */
295  static int CheckOpen(int type);
296 };
297
298 /// Container class for integration rules
300 {
301 private:
302  /// Taken from the Quadrature1D class anonymous enum
303  /// Determines the type of numerical quadrature used for
304  /// segment, square, and cube geometries
306
307  int own_rules, refined;
308
309  /// Function that generates quadrature points and weights on [0,1]
311
312  Array<IntegrationRule *> PointIntRules;
313  Array<IntegrationRule *> SegmentIntRules;
314  Array<IntegrationRule *> TriangleIntRules;
315  Array<IntegrationRule *> SquareIntRules;
316  Array<IntegrationRule *> TetrahedronIntRules;
317  Array<IntegrationRule *> PrismIntRules;
318  Array<IntegrationRule *> CubeIntRules;
319
320  void AllocIntRule(Array<IntegrationRule *> &ir_array, int Order)
321  {
322  if (ir_array.Size() <= Order)
323  {
324  ir_array.SetSize(Order + 1, NULL);
325  }
326  }
327  bool HaveIntRule(Array<IntegrationRule *> &ir_array, int Order)
328  {
329  return (ir_array.Size() > Order && ir_array[Order] != NULL);
330  }
331  int GetSegmentRealOrder(int Order) const
332  {
333  return Order | 1; // valid for all quad_type's
334  }
335
336  /// The following methods allocate new IntegrationRule objects without
337  /// checking if they already exist. To avoid memory leaks use
338  /// IntegrationRules::Get(int GeomType, int Order) instead.
339  IntegrationRule *GenerateIntegrationRule(int GeomType, int Order);
340  IntegrationRule *PointIntegrationRule(int Order);
341  IntegrationRule *SegmentIntegrationRule(int Order);
342  IntegrationRule *TriangleIntegrationRule(int Order);
343  IntegrationRule *SquareIntegrationRule(int Order);
344  IntegrationRule *TetrahedronIntegrationRule(int Order);
345  IntegrationRule *PrismIntegrationRule(int Order);
346  IntegrationRule *CubeIntegrationRule(int Order);
347
348  void DeleteIntRuleArray(Array<IntegrationRule *> &ir_array);
349
350 public:
351  /// Sets initial sizes for the integration rule arrays, but rules
352  /// are defined the first time they are requested with the Get method.
353  explicit IntegrationRules(int Ref = 0,
355
356  /// Returns an integration rule for given GeomType and Order.
357  const IntegrationRule &Get(int GeomType, int Order);
358
359  void Set(int GeomType, int Order, IntegrationRule &IntRule);
360
361  void SetOwnRules(int o) { own_rules = o; }
362
363  /// Destroys an IntegrationRules object
365 };
366
367 /// A global object with all integration rules (defined in intrules.cpp)
369
370 /// A global object with all refined integration rules
372
373 }
374
375 #endif
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:237
void Set(const double *p, const int dim)
Definition: intrules.hpp:32
int Size() const
Logical size of the array.
Definition: array.hpp:118
void Get(double *p, const int dim) const
Definition: intrules.hpp:46
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:85
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Definition: intrules.cpp:877
void ClosedUniform(const int np, IntegrationRule *ir)
Definition: intrules.cpp:578
void Set1w(const double *p)
Definition: intrules.hpp:81
void GivePolyPoints(const int np, double *pts, const int type)
Definition: intrules.cpp:609
~IntegrationRules()
Destroys an IntegrationRules object.
Definition: intrules.cpp:968
aka closed Newton-Cotes
Definition: intrules.hpp:287
Container class for integration rules.
Definition: intrules.hpp:299
void Set(const double x1, const double x2, const double x3, const double w)
Definition: intrules.hpp:60
const Array< double > & GetWeights() const
Return the quadrature weights in a contiguous array.
Definition: intrules.cpp:81
void OpenUniform(const int np, IntegrationRule *ir)
Definition: intrules.cpp:564
const IntegrationPoint & IntPoint(int i) const
Returns a const reference to the i-th integration point.
Definition: intrules.hpp:243
A class container for 1D quadrature type constants.
Definition: intrules.hpp:278
aka open Newton-Cotes
Definition: intrules.hpp:286
void GaussLobatto(const int np, IntegrationRule *ir)
Definition: intrules.cpp:442
void Set3(const double *p)
Definition: intrules.hpp:68
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:240
void SetOwnRules(int o)
Definition: intrules.hpp:361
int dim
Definition: ex3.cpp:48
~IntegrationRule()
Destroys an IntegrationRule object.
Definition: intrules.hpp:251
void Set3w(const double *p)
Definition: intrules.hpp:63
void Set2(const double x1, const double x2)
Definition: intrules.hpp:75
void Set3(const double x1, const double x2, const double x3)
Definition: intrules.hpp:65
void Set2w(const double *p)
Definition: intrules.hpp:73
static int CheckClosed(int type)
If the Quadrature1D type is not closed return Invalid; otherwise return type.
Definition: intrules.cpp:812
void OpenHalfUniform(const int np, IntegrationRule *ir)
Definition: intrules.cpp:596
int GetOrder() const
Returns the order of the integration rule.
Definition: intrules.hpp:230
void SetSize(int nsize)
Change logical size of the array, keep existing entries.
Definition: array.hpp:618
static int CheckOpen(int type)
If the Quadrature1D type is not open return Invalid; otherwise return type.
Definition: intrules.cpp:824
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:234
Definition: intrules.cpp:844
Class for integration point with weight.
Definition: intrules.hpp:25
aka &quot;open half&quot; Newton-Cotes
Definition: intrules.hpp:288
void Set2(const double *p)
Definition: intrules.hpp:77
A global object with all refined integration rules.
Definition: intrules.hpp:371
void Set2w(const double x1, const double x2, const double w)
Definition: intrules.hpp:70
void pts(int iphi, int t, double x[])
IntegrationRule(int NP)
Construct an integration rule with given number of points.
Definition: intrules.hpp:213
void Set(int GeomType, int Order, IntegrationRule &IntRule)
Definition: intrules.cpp:923