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