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