MFEM  v3.4 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  x = p[0];
35  if (dim > 1)
36  {
37  y = p[1];
38  if (dim > 2)
39  {
40  z = p[2];
41  }
42  }
43  }
44
45  void Get(double *p, const int dim) const
46  {
47  p[0] = x;
48  if (dim > 1)
49  {
50  p[1] = y;
51  if (dim > 2)
52  {
53  p[2] = z;
54  }
55  }
56  }
57
58  void Set(const double x1, const double x2, const double x3, const double w)
59  { x = x1; y = x2; z = x3; weight = w; }
60
61  void Set3w(const double *p) { x = p[0]; y = p[1]; z = p[2]; weight = p[3]; }
62
63  void Set3(const double x1, const double x2, const double x3)
64  { x = x1; y = x2; z = x3; }
65
66  void Set3(const double *p) { x = p[0]; y = p[1]; z = p[2]; }
67
68  void Set2w(const double x1, const double x2, const double w)
69  { x = x1; y = x2; weight = w; }
70
71  void Set2w(const double *p) { x = p[0]; y = p[1]; weight = p[2]; }
72
73  void Set2(const double x1, const double x2) { x = x1; y = x2; }
74
75  void Set2(const double *p) { x = p[0]; y = p[1]; }
76
77  void Set1w(const double x1, const double w) { x = x1; weight = w; }
78
79  void Set1w(const double *p) { x = p[0]; weight = p[1]; }
80 };
81
82 /// Class for an integration rule - an Array of IntegrationPoint.
83 class IntegrationRule : public Array<IntegrationPoint>
84 {
85 private:
86  friend class IntegrationRules;
87  int Order;
88
89  /// Define n-simplex rule (triangle/tetrahedron for n=2/3) of order (2s+1)
90  void GrundmannMollerSimplexRule(int s, int n = 3);
91
92  void AddTriMidPoint(const int off, const double weight)
93  { IntPoint(off).Set2w(1./3., 1./3., weight); }
94
95  void AddTriPoints3(const int off, const double a, const double b,
96  const double weight)
97  {
98  IntPoint(off + 0).Set2w(a, a, weight);
99  IntPoint(off + 1).Set2w(a, b, weight);
100  IntPoint(off + 2).Set2w(b, a, weight);
101  }
102
103  void AddTriPoints3(const int off, const double a, const double weight)
104  { AddTriPoints3(off, a, 1. - 2.*a, weight); }
105
106  void AddTriPoints3b(const int off, const double b, const double weight)
107  { AddTriPoints3(off, (1. - b)/2., b, weight); }
108
109  void AddTriPoints3R(const int off, const double a, const double b,
110  const double c, const double weight)
111  {
112  IntPoint(off + 0).Set2w(a, b, weight);
113  IntPoint(off + 1).Set2w(c, a, weight);
114  IntPoint(off + 2).Set2w(b, c, weight);
115  }
116
117  void AddTriPoints3R(const int off, const double a, const double b,
118  const double weight)
119  { AddTriPoints3R(off, a, b, 1. - a - b, weight); }
120
121  void AddTriPoints6(const int off, const double a, const double b,
122  const double c, const double weight)
123  {
124  IntPoint(off + 0).Set2w(a, b, weight);
125  IntPoint(off + 1).Set2w(b, a, weight);
126  IntPoint(off + 2).Set2w(a, c, weight);
127  IntPoint(off + 3).Set2w(c, a, weight);
128  IntPoint(off + 4).Set2w(b, c, weight);
129  IntPoint(off + 5).Set2w(c, b, weight);
130  }
131
132  void AddTriPoints6(const int off, const double a, const double b,
133  const double weight)
134  { AddTriPoints6(off, a, b, 1. - a - b, weight); }
135
136  // add the permutations of (a,a,b)
137  void AddTetPoints3(const int off, const double a, const double b,
138  const double weight)
139  {
140  IntPoint(off + 0).Set(a, a, b, weight);
141  IntPoint(off + 1).Set(a, b, a, weight);
142  IntPoint(off + 2).Set(b, a, a, weight);
143  }
144
145  // add the permutations of (a,b,c)
146  void AddTetPoints6(const int off, const double a, const double b,
147  const double c, const double weight)
148  {
149  IntPoint(off + 0).Set(a, b, c, weight);
150  IntPoint(off + 1).Set(a, c, b, weight);
151  IntPoint(off + 2).Set(b, c, a, weight);
152  IntPoint(off + 3).Set(b, a, c, weight);
153  IntPoint(off + 4).Set(c, a, b, weight);
154  IntPoint(off + 5).Set(c, b, a, weight);
155  }
156
157  void AddTetMidPoint(const int off, const double weight)
158  { IntPoint(off).Set(0.25, 0.25, 0.25, weight); }
159
160  // given a, add the permutations of (a,a,a,b), where 3*a + b = 1
161  void AddTetPoints4(const int off, const double a, const double weight)
162  {
163  IntPoint(off).Set(a, a, a, weight);
164  AddTetPoints3(off + 1, a, 1. - 3.*a, weight);
165  }
166
167  // given b, add the permutations of (a,a,a,b), where 3*a + b = 1
168  void AddTetPoints4b(const int off, const double b, const double weight)
169  {
170  const double a = (1. - b)/3.;
171  IntPoint(off).Set(a, a, a, weight);
172  AddTetPoints3(off + 1, a, b, weight);
173  }
174
175  // add the permutations of (a,a,b,b), 2*(a + b) = 1
176  void AddTetPoints6(const int off, const double a, const double weight)
177  {
178  const double b = 0.5 - a;
180  AddTetPoints3(off + 3, b, a, weight);
181  }
182
183  // given (a,b) or (a,c), add the permutations of (a,a,b,c), 2*a + b + c = 1
184  void AddTetPoints12(const int off, const double a, const double bc,
185  const double weight)
186  {
187  const double cb = 1. - 2*a - bc;
189  AddTetPoints3(off + 3, a, cb, weight);
190  AddTetPoints6(off + 6, a, bc, cb, weight);
191  }
192
193  // given (b,c), add the permutations of (a,a,b,c), 2*a + b + c = 1
194  void AddTetPoints12bc(const int off, const double b, const double c,
195  const double weight)
196  {
197  const double a = (1. - b - c)/2.;
199  AddTetPoints3(off + 3, a, c, weight);
200  AddTetPoints6(off + 6, a, b, c, weight);
201  }
202
203 public:
205  Array<IntegrationPoint>(), Order(0) { }
206
207  /// Construct an integration rule with given number of points
208  explicit IntegrationRule(int NP) :
209  Array<IntegrationPoint>(NP), Order(0)
210  {
211  for (int i = 0; i < this->Size(); i++)
212  {
213  (*this)[i].Init();
214  }
215  }
216
217  /// Tensor product of two 1D integration rules
219
220  /// Tensor product of three 1D integration rules
222  IntegrationRule &irz);
223
224  /// Returns the order of the integration rule
225  int GetOrder() const { return Order; }
226
227  /** @brief Sets the order of the integration rule. This is only for keeping
228  order information, it does not alter any data in the IntegrationRule. */
229  void SetOrder(const int order) { Order = order; }
230
231  /// Returns the number of the points in the integration rule
232  int GetNPoints() const { return Size(); }
233
234  /// Returns a reference to the i-th integration point
235  IntegrationPoint &IntPoint(int i) { return (*this)[i]; }
236
237  /// Returns a const reference to the i-th integration point
238  const IntegrationPoint &IntPoint(int i) const { return (*this)[i]; }
239
240  /// Destroys an IntegrationRule object
242 };
243
244 /// A Class that defines 1-D numerical quadrature rules on [0,1].
246 {
247 public:
248  /** @name Methods for calculating quadrature rules.
249  These methods calculate the actual points and weights for the different
250  types of quadrature rules. */
251  ///@{
252  void GaussLegendre(const int np, IntegrationRule* ir);
253  void GaussLobatto(const int np, IntegrationRule *ir);
254  void OpenUniform(const int np, IntegrationRule *ir);
255  void ClosedUniform(const int np, IntegrationRule *ir);
256  void OpenHalfUniform(const int np, IntegrationRule *ir);
257  ///@}
258
259  /// A helper function that will play nice with Poly_1D::OpenPoints and
260  /// Poly_1D::ClosedPoints
261  void GivePolyPoints(const int np, double *pts, const int type);
262
263 private:
264  void CalculateUniformWeights(IntegrationRule *ir, const int type);
265 };
266
267 /// A class container for 1D quadrature type constants.
269 {
270 public:
271  enum
272  {
273  Invalid = -1,
276  OpenUniform = 2, ///< aka open Newton-Cotes
277  ClosedUniform = 3, ///< aka closed Newton-Cotes
278  OpenHalfUniform = 4 ///< aka "open half" Newton-Cotes
279  };
280  /** @brief If the Quadrature1D type is not closed return Invalid; otherwise
281  return type. */
282  static int CheckClosed(int type);
283  /** @brief If the Quadrature1D type is not open return Invalid; otherwise
284  return type. */
285  static int CheckOpen(int type);
286 };
287
288 /// Container class for integration rules
290 {
291 private:
292  /// Taken from the Quadrature1D class anonymous enum
293  /// Determines the type of numerical quadrature used for
294  /// segment, square, and cube geometries
296
297  int own_rules, refined;
298
299  /// Function that generates quadrature points and weights on [0,1]
301
302  Array<IntegrationRule *> PointIntRules;
303  Array<IntegrationRule *> SegmentIntRules;
304  Array<IntegrationRule *> TriangleIntRules;
305  Array<IntegrationRule *> SquareIntRules;
306  Array<IntegrationRule *> TetrahedronIntRules;
307  Array<IntegrationRule *> CubeIntRules;
308
309  void AllocIntRule(Array<IntegrationRule *> &ir_array, int Order)
310  {
311  if (ir_array.Size() <= Order)
312  {
313  ir_array.SetSize(Order + 1, NULL);
314  }
315  }
316  bool HaveIntRule(Array<IntegrationRule *> &ir_array, int Order)
317  {
318  return (ir_array.Size() > Order && ir_array[Order] != NULL);
319  }
320  int GetSegmentRealOrder(int Order) const
321  {
322  return Order | 1; // valid for all quad_type's
323  }
324
325  IntegrationRule *GenerateIntegrationRule(int GeomType, int Order);
326  IntegrationRule *PointIntegrationRule(int Order);
327  IntegrationRule *SegmentIntegrationRule(int Order);
328  IntegrationRule *TriangleIntegrationRule(int Order);
329  IntegrationRule *SquareIntegrationRule(int Order);
330  IntegrationRule *TetrahedronIntegrationRule(int Order);
331  IntegrationRule *CubeIntegrationRule(int Order);
332
333  void DeleteIntRuleArray(Array<IntegrationRule *> &ir_array);
334
335 public:
336  /// Sets initial sizes for the integration rule arrays, but rules
337  /// are defined the first time they are requested with the Get method.
338  explicit IntegrationRules(int Ref = 0,
340
341  /// Returns an integration rule for given GeomType and Order.
342  const IntegrationRule &Get(int GeomType, int Order);
343
344  void Set(int GeomType, int Order, IntegrationRule &IntRule);
345
346  void SetOwnRules(int o) { own_rules = o; }
347
348  /// Destroys an IntegrationRules object
350 };
351
352 /// A global object with all integration rules (defined in intrules.cpp)
354
355 /// A global object with all refined integration rules
357
358 }
359
360 #endif
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:232
void Set(const double *p, const int dim)
Definition: intrules.hpp:32
int Size() const
Logical size of the array.
Definition: array.hpp:133
void Get(double *p, const int dim) const
Definition: intrules.hpp:45
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:83
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Definition: intrules.cpp:861
void ClosedUniform(const int np, IntegrationRule *ir)
Definition: intrules.cpp:565
void Set1w(const double *p)
Definition: intrules.hpp:79
void GivePolyPoints(const int np, double *pts, const int type)
Definition: intrules.cpp:596
~IntegrationRules()
Destroys an IntegrationRules object.
Definition: intrules.cpp:950
Container class for integration rules.
Definition: intrules.hpp:289
void Set(const double x1, const double x2, const double x3, const double w)
Definition: intrules.hpp:58
void OpenUniform(const int np, IntegrationRule *ir)
Definition: intrules.cpp:551
const IntegrationPoint & IntPoint(int i) const
Returns a const reference to the i-th integration point.
Definition: intrules.hpp:238
A class container for 1D quadrature type constants.
Definition: intrules.hpp:268
void GaussLobatto(const int np, IntegrationRule *ir)
Definition: intrules.cpp:429
void Set3(const double *p)
Definition: intrules.hpp:66
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:235
void SetOwnRules(int o)
Definition: intrules.hpp:346
int dim
Definition: ex3.cpp:47
~IntegrationRule()
Destroys an IntegrationRule object.
Definition: intrules.hpp:241
void Set3w(const double *p)
Definition: intrules.hpp:61
void Set2(const double x1, const double x2)
Definition: intrules.hpp:73
aka &quot;open half&quot; Newton-Cotes
Definition: intrules.hpp:278
void Set3(const double x1, const double x2, const double x3)
Definition: intrules.hpp:63
void Set2w(const double *p)
Definition: intrules.hpp:71
static int CheckClosed(int type)
If the Quadrature1D type is not closed return Invalid; otherwise return type.
Definition: intrules.cpp:799
void OpenHalfUniform(const int np, IntegrationRule *ir)
Definition: intrules.cpp:583
int GetOrder() const
Returns the order of the integration rule.
Definition: intrules.hpp:225
void SetSize(int nsize)
Change logical size of the array, keep existing entries.
Definition: array.hpp:569
static int CheckOpen(int type)
If the Quadrature1D type is not open return Invalid; otherwise return type.
Definition: intrules.cpp:811
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:229
Definition: intrules.cpp:831
Class for integration point with weight.
Definition: intrules.hpp:25
aka open Newton-Cotes
Definition: intrules.hpp:276
void Set2(const double *p)
Definition: intrules.hpp:75
A global object with all refined integration rules.
Definition: intrules.hpp:356
void Set2w(const double x1, const double x2, const double w)
Definition: intrules.hpp:68
IntegrationRule(int NP)
Construct an integration rule with given number of points.
Definition: intrules.hpp:208
void Set(int GeomType, int Order, IntegrationRule &IntRule)
Definition: intrules.cpp:906