quadrature_gauss.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2018 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 // libMesh includes
20 #include "libmesh/quadrature_gauss.h"
21 #include "libmesh/enum_quadrature_type.h"
22 
23 namespace libMesh
24 {
25 
26 // See the files:
27 // quadrature_gauss_1D.C
28 // quadrature_gauss_2D.C
29 // quadrature_gauss_3D.C
30 // for implementation of specific element types.
31 
32 
33 QuadratureType QGauss::type() const
34 {
35  return QGAUSS;
36 }
37 
38 void QGauss::keast_rule(const Real rule_data[][4],
39  const unsigned int n_pts)
40 {
41  // Like the Dunavant rule, the input data should have 4 columns. These columns
42  // have the following format and implied permutations (w=weight).
43  // {a, 0, 0, w} = 1-permutation (a,a,a)
44  // {a, b, 0, w} = 4-permutation (a,b,b), (b,a,b), (b,b,a), (b,b,b)
45  // {a, 0, b, w} = 6-permutation (a,a,b), (a,b,b), (b,b,a), (b,a,b), (b,a,a), (a,b,a)
46  // {a, b, c, w} = 12-permutation (a,a,b), (a,a,c), (b,a,a), (c,a,a), (a,b,a), (a,c,a)
47  // (a,b,c), (a,c,b), (b,a,c), (b,c,a), (c,a,b), (c,b,a)
48 
49  // Always insert into the points & weights vector relative to the offset
50  unsigned int offset=0;
51 
52 
53  for (unsigned int p=0; p<n_pts; ++p)
54  {
55 
56  // There must always be a non-zero entry to start the row
57  libmesh_assert_not_equal_to (rule_data[p][0], static_cast<Real>(0.0));
58 
59  // A zero weight may imply you did not set up the raw data correctly
60  libmesh_assert_not_equal_to (rule_data[p][3], static_cast<Real>(0.0));
61 
62  // What kind of point is this?
63  // One non-zero entry in first 3 cols ? 1-perm (centroid) point = 1
64  // Two non-zero entries in first 3 cols ? 3-perm point = 3
65  // Three non-zero entries ? 6-perm point = 6
66  unsigned int pointtype=1;
67 
68  if (rule_data[p][1] != static_cast<Real>(0.0))
69  {
70  if (rule_data[p][2] != static_cast<Real>(0.0))
71  pointtype = 12;
72  else
73  pointtype = 4;
74  }
75  else
76  {
77  // The second entry is zero. What about the third?
78  if (rule_data[p][2] != static_cast<Real>(0.0))
79  pointtype = 6;
80  }
81 
82 
83  switch (pointtype)
84  {
85  case 1:
86  {
87  // Be sure we have enough space to insert this point
88  libmesh_assert_less (offset + 0, _points.size());
89 
90  const Real a = rule_data[p][0];
91 
92  // The point has only a single permutation (the centroid!)
93  _points[offset + 0] = Point(a,a,a);
94 
95  // The weight is always the last entry in the row.
96  _weights[offset + 0] = rule_data[p][3];
97 
98  offset += pointtype;
99  break;
100  }
101 
102  case 4:
103  {
104  // Be sure we have enough space to insert these points
105  libmesh_assert_less (offset + 3, _points.size());
106 
107  const Real a = rule_data[p][0];
108  const Real b = rule_data[p][1];
109  const Real wt = rule_data[p][3];
110 
111  // Here it's understood the second entry is to be used twice, and
112  // thus there are three possible permutations.
113  _points[offset + 0] = Point(a,b,b);
114  _points[offset + 1] = Point(b,a,b);
115  _points[offset + 2] = Point(b,b,a);
116  _points[offset + 3] = Point(b,b,b);
117 
118  for (unsigned int j=0; j<pointtype; ++j)
119  _weights[offset + j] = wt;
120 
121  offset += pointtype;
122  break;
123  }
124 
125  case 6:
126  {
127  // Be sure we have enough space to insert these points
128  libmesh_assert_less (offset + 5, _points.size());
129 
130  const Real a = rule_data[p][0];
131  const Real b = rule_data[p][2];
132  const Real wt = rule_data[p][3];
133 
134  // Three individual entries with six permutations.
135  _points[offset + 0] = Point(a,a,b);
136  _points[offset + 1] = Point(a,b,b);
137  _points[offset + 2] = Point(b,b,a);
138  _points[offset + 3] = Point(b,a,b);
139  _points[offset + 4] = Point(b,a,a);
140  _points[offset + 5] = Point(a,b,a);
141 
142  for (unsigned int j=0; j<pointtype; ++j)
143  _weights[offset + j] = wt;
144 
145  offset += pointtype;
146  break;
147  }
148 
149 
150  case 12:
151  {
152  // Be sure we have enough space to insert these points
153  libmesh_assert_less (offset + 11, _points.size());
154 
155  const Real a = rule_data[p][0];
156  const Real b = rule_data[p][1];
157  const Real c = rule_data[p][2];
158  const Real wt = rule_data[p][3];
159 
160  // Three individual entries with six permutations.
161  _points[offset + 0] = Point(a,a,b); _points[offset + 6] = Point(a,b,c);
162  _points[offset + 1] = Point(a,a,c); _points[offset + 7] = Point(a,c,b);
163  _points[offset + 2] = Point(b,a,a); _points[offset + 8] = Point(b,a,c);
164  _points[offset + 3] = Point(c,a,a); _points[offset + 9] = Point(b,c,a);
165  _points[offset + 4] = Point(a,b,a); _points[offset + 10] = Point(c,a,b);
166  _points[offset + 5] = Point(a,c,a); _points[offset + 11] = Point(c,b,a);
167 
168  for (unsigned int j=0; j<pointtype; ++j)
169  _weights[offset + j] = wt;
170 
171  offset += pointtype;
172  break;
173  }
174 
175  default:
176  libmesh_error_msg("Don't know what to do with this many permutation points!");
177  }
178 
179  }
180 
181 }
182 
183 
184 // A number of different rules for triangles can be described by
185 // permutations of the following types of points:
186 // I: "1"-permutation, (1/3,1/3) (single point only)
187 // II: 3-permutation, (a,a,1-2a)
188 // III: 6-permutation, (a,b,1-a-b)
189 // The weights for a given set of permutations are all the same.
190 void QGauss::dunavant_rule2(const Real * wts,
191  const Real * a,
192  const Real * b,
193  const unsigned int * permutation_ids,
194  unsigned int n_wts)
195 {
196  // Figure out how many total points by summing up the entries
197  // in the permutation_ids array, and resize the _points and _weights
198  // vectors appropriately.
199  unsigned int total_pts = 0;
200  for (unsigned int p=0; p<n_wts; ++p)
201  total_pts += permutation_ids[p];
202 
203  // Resize point and weight vectors appropriately.
204  _points.resize(total_pts);
205  _weights.resize(total_pts);
206 
207  // Always insert into the points & weights vector relative to the offset
208  unsigned int offset=0;
209 
210  for (unsigned int p=0; p<n_wts; ++p)
211  {
212  switch (permutation_ids[p])
213  {
214  case 1:
215  {
216  // The point has only a single permutation (the centroid!)
217  // So we don't even need to look in the a or b arrays.
218  _points [offset + 0] = Point(Real(1)/3, Real(1)/3);
219  _weights[offset + 0] = wts[p];
220 
221  offset += 1;
222  break;
223  }
224 
225 
226  case 3:
227  {
228  // For this type of rule, don't need to look in the b array.
229  _points[offset + 0] = Point(a[p], a[p]); // (a,a)
230  _points[offset + 1] = Point(a[p], 1-2*a[p]); // (a,1-2a)
231  _points[offset + 2] = Point(1-2*a[p], a[p]); // (1-2a,a)
232 
233  for (unsigned int j=0; j<3; ++j)
234  _weights[offset + j] = wts[p];
235 
236  offset += 3;
237  break;
238  }
239 
240  case 6:
241  {
242  // This type of point uses all 3 arrays...
243  _points[offset + 0] = Point(a[p], b[p]);
244  _points[offset + 1] = Point(b[p], a[p]);
245  _points[offset + 2] = Point(a[p], 1-a[p]-b[p]);
246  _points[offset + 3] = Point(1-a[p]-b[p], a[p]);
247  _points[offset + 4] = Point(b[p], 1-a[p]-b[p]);
248  _points[offset + 5] = Point(1-a[p]-b[p], b[p]);
249 
250  for (unsigned int j=0; j<6; ++j)
251  _weights[offset + j] = wts[p];
252 
253  offset += 6;
254  break;
255  }
256 
257  default:
258  libmesh_error_msg("Unknown permutation id: " << permutation_ids[p] << "!");
259  }
260  }
261 
262 }
263 
264 
265 void QGauss::dunavant_rule(const Real rule_data[][4],
266  const unsigned int n_pts)
267 {
268  // The input data array has 4 columns. The first 3 are the permutation points.
269  // The last column is the weights for a given set of permutation points. A zero
270  // in two of the first 3 columns implies the point is a 1-permutation (centroid).
271  // A zero in one of the first 3 columns implies the point is a 3-permutation.
272  // Otherwise each point is assumed to be a 6-permutation.
273 
274  // Always insert into the points & weights vector relative to the offset
275  unsigned int offset=0;
276 
277 
278  for (unsigned int p=0; p<n_pts; ++p)
279  {
280 
281  // There must always be a non-zero entry to start the row
282  libmesh_assert_not_equal_to ( rule_data[p][0], static_cast<Real>(0.0) );
283 
284  // A zero weight may imply you did not set up the raw data correctly
285  libmesh_assert_not_equal_to ( rule_data[p][3], static_cast<Real>(0.0) );
286 
287  // What kind of point is this?
288  // One non-zero entry in first 3 cols ? 1-perm (centroid) point = 1
289  // Two non-zero entries in first 3 cols ? 3-perm point = 3
290  // Three non-zero entries ? 6-perm point = 6
291  unsigned int pointtype=1;
292 
293  if (rule_data[p][1] != static_cast<Real>(0.0))
294  {
295  if (rule_data[p][2] != static_cast<Real>(0.0))
296  pointtype = 6;
297  else
298  pointtype = 3;
299  }
300 
301  switch (pointtype)
302  {
303  case 1:
304  {
305  // Be sure we have enough space to insert this point
306  libmesh_assert_less (offset + 0, _points.size());
307 
308  // The point has only a single permutation (the centroid!)
309  _points[offset + 0] = Point(rule_data[p][0], rule_data[p][0]);
310 
311  // The weight is always the last entry in the row.
312  _weights[offset + 0] = rule_data[p][3];
313 
314  offset += 1;
315  break;
316  }
317 
318  case 3:
319  {
320  // Be sure we have enough space to insert these points
321  libmesh_assert_less (offset + 2, _points.size());
322 
323  // Here it's understood the second entry is to be used twice, and
324  // thus there are three possible permutations.
325  _points[offset + 0] = Point(rule_data[p][0], rule_data[p][1]);
326  _points[offset + 1] = Point(rule_data[p][1], rule_data[p][0]);
327  _points[offset + 2] = Point(rule_data[p][1], rule_data[p][1]);
328 
329  for (unsigned int j=0; j<3; ++j)
330  _weights[offset + j] = rule_data[p][3];
331 
332  offset += 3;
333  break;
334  }
335 
336  case 6:
337  {
338  // Be sure we have enough space to insert these points
339  libmesh_assert_less (offset + 5, _points.size());
340 
341  // Three individual entries with six permutations.
342  _points[offset + 0] = Point(rule_data[p][0], rule_data[p][1]);
343  _points[offset + 1] = Point(rule_data[p][0], rule_data[p][2]);
344  _points[offset + 2] = Point(rule_data[p][1], rule_data[p][0]);
345  _points[offset + 3] = Point(rule_data[p][1], rule_data[p][2]);
346  _points[offset + 4] = Point(rule_data[p][2], rule_data[p][0]);
347  _points[offset + 5] = Point(rule_data[p][2], rule_data[p][1]);
348 
349  for (unsigned int j=0; j<6; ++j)
350  _weights[offset + j] = rule_data[p][3];
351 
352  offset += 6;
353  break;
354  }
355 
356  default:
357  libmesh_error_msg("Don't know what to do with this many permutation points!");
358  }
359  }
360 }
361 
362 } // namespace libMesh
enum_quadrature_type.h
libMesh::QGauss::type
virtual QuadratureType type() const override
Definition: quadrature_gauss.C:33
libMesh::QBase::_points
std::vector< Point > _points
Definition: quadrature.h:349
libMesh::QBase::_weights
std::vector< Real > _weights
Definition: quadrature.h:355
libMesh::QGauss::dunavant_rule2
void dunavant_rule2(const Real *wts, const Real *a, const Real *b, const unsigned int *permutation_ids, const unsigned int n_wts)
Definition: quadrature_gauss.C:190
libMesh::QGauss::dunavant_rule
void dunavant_rule(const Real rule_data[][4], const unsigned int n_pts)
Definition: quadrature_gauss.C:265
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:131
libMesh::QGauss::keast_rule
void keast_rule(const Real rule_data[][4], const unsigned int n_pts)
Definition: quadrature_gauss.C:38
quadrature_gauss.h
libMesh::Point
A geometric point in (x,y,z) space.
Definition: point.h:38
generated by doxygen