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