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 :
20 : // Local includes
21 : #include "libmesh/quadrature_gauss.h"
22 : #include "libmesh/quadrature_conical.h"
23 : #include "libmesh/quadrature_gm.h"
24 : #include "libmesh/enum_to_string.h"
25 : #include "libmesh/cell_c0polyhedron.h"
26 :
27 : namespace libMesh
28 : {
29 :
30 171105 : void QGauss::init_3D()
31 : {
32 : #if LIBMESH_DIM == 3
33 :
34 : //-----------------------------------------------------------------------
35 : // 3D quadrature rules
36 171105 : switch (_type)
37 : {
38 : //---------------------------------------------
39 : // Hex quadrature rules
40 32730 : case HEX8:
41 : case HEX20:
42 : case HEX27:
43 : {
44 : // We compute the 3D quadrature rule as a tensor
45 : // product of the 1D quadrature rule.
46 35968 : QGauss q1D(1, get_order());
47 32730 : tensor_product_hex( q1D );
48 1619 : return;
49 : }
50 :
51 :
52 :
53 : //---------------------------------------------
54 : // Tetrahedral quadrature rules
55 92923 : case TET4:
56 : case TET10:
57 : case TET14:
58 : {
59 99291 : switch(get_order())
60 : {
61 : // Taken from pg. 222 of "The finite element method," vol. 1
62 : // ed. 5 by Zienkiewicz & Taylor
63 437 : case CONSTANT:
64 : case FIRST:
65 : {
66 : // Exact for linears
67 437 : _points.resize(1);
68 437 : _weights.resize(1);
69 :
70 :
71 437 : _points[0](0) = .25;
72 437 : _points[0](1) = .25;
73 437 : _points[0](2) = .25;
74 :
75 437 : _weights[0] = Real(1)/6;
76 :
77 437 : return;
78 : }
79 284 : case SECOND:
80 : {
81 : // Exact for quadratics
82 284 : _points.resize(4);
83 284 : _weights.resize(4);
84 :
85 :
86 : // Can't be constexpr with my version of Boost quad
87 : // precision
88 8 : const Real b = 0.25*(1-std::sqrt(Real(5))/5);
89 8 : const Real a = 1-3*b;
90 :
91 284 : _points[0](0) = a;
92 284 : _points[0](1) = b;
93 284 : _points[0](2) = b;
94 :
95 284 : _points[1](0) = b;
96 284 : _points[1](1) = a;
97 284 : _points[1](2) = b;
98 :
99 284 : _points[2](0) = b;
100 284 : _points[2](1) = b;
101 284 : _points[2](2) = a;
102 :
103 284 : _points[3](0) = b;
104 284 : _points[3](1) = b;
105 284 : _points[3](2) = b;
106 :
107 :
108 :
109 284 : _weights[0] = Real(1)/24;
110 284 : _weights[1] = _weights[0];
111 284 : _weights[2] = _weights[0];
112 284 : _weights[3] = _weights[0];
113 :
114 284 : return;
115 : }
116 :
117 :
118 :
119 : // Can be found in the class notes
120 : // http://www.cs.rpi.edu/~flaherje/FEM/fem6.ps
121 : // by Flaherty.
122 : //
123 : // Caution: this rule has a negative weight and may be
124 : // unsuitable for some problems.
125 : // Exact for cubics.
126 : //
127 : // Note: Keast (see ref. elsewhere in this file) also gives
128 : // a third-order rule with positive weights, but it contains points
129 : // on the ref. elt. boundary, making it less suitable for FEM calculations.
130 8550 : case THIRD:
131 : {
132 8550 : if (allow_rules_with_negative_weights)
133 : {
134 8408 : _points.resize(5);
135 8408 : _weights.resize(5);
136 :
137 :
138 8408 : _points[0](0) = .25;
139 8408 : _points[0](1) = .25;
140 8408 : _points[0](2) = .25;
141 :
142 8408 : _points[1](0) = .5;
143 8408 : _points[1](1) = Real(1)/6;
144 8408 : _points[1](2) = Real(1)/6;
145 :
146 8408 : _points[2](0) = Real(1)/6;
147 8408 : _points[2](1) = .5;
148 8408 : _points[2](2) = Real(1)/6;
149 :
150 8408 : _points[3](0) = Real(1)/6;
151 8408 : _points[3](1) = Real(1)/6;
152 8408 : _points[3](2) = .5;
153 :
154 8408 : _points[4](0) = Real(1)/6;
155 8408 : _points[4](1) = Real(1)/6;
156 8408 : _points[4](2) = Real(1)/6;
157 :
158 :
159 8408 : _weights[0] = Real(-2)/15;
160 8408 : _weights[1] = .075;
161 8408 : _weights[2] = _weights[1];
162 8408 : _weights[3] = _weights[1];
163 8408 : _weights[4] = _weights[1];
164 :
165 8408 : return;
166 : } // end if (allow_rules_with_negative_weights)
167 : else
168 : {
169 : // If a rule with positive weights is required, the 2x2x2 conical
170 : // product rule is third-order accurate and has less points than
171 : // the next-available positive-weight rule at FIFTH order.
172 4 : QConical conical_rule(3, _order);
173 142 : conical_rule.init(*this);
174 :
175 : // Swap points and weights with the about-to-be destroyed rule.
176 4 : _points.swap (conical_rule.get_points() );
177 4 : _weights.swap(conical_rule.get_weights());
178 :
179 4 : return;
180 : }
181 : // Note: if !allow_rules_with_negative_weights, fall through to next case.
182 : }
183 :
184 :
185 :
186 : // Originally a Keast rule,
187 : // Patrick Keast,
188 : // Moderate Degree Tetrahedral Quadrature Formulas,
189 : // Computer Methods in Applied Mechanics and Engineering,
190 : // Volume 55, Number 3, May 1986, pages 339-348.
191 : //
192 : // Can also be found the class notes
193 : // http://www.cs.rpi.edu/~flaherje/FEM/fem6.ps
194 : // by Flaherty.
195 : //
196 : // Caution: this rule has a negative weight and may be
197 : // unsuitable for some problems.
198 1207 : case FOURTH:
199 : {
200 1207 : if (allow_rules_with_negative_weights)
201 : {
202 213 : _points.resize(11);
203 213 : _weights.resize(11);
204 :
205 : // The raw data for the quadrature rule.
206 213 : const Real rule_data[3][4] = {
207 : {0.250000000000000000e+00_R, 0._R, 0._R, -0.131555555555555556e-01_R}, // 1
208 : {0.785714285714285714e+00_R, 0.714285714285714285e-01_R, 0._R, 0.762222222222222222e-02_R}, // 4
209 : {0.399403576166799219e+00_R, 0._R, 0.100596423833200785e+00_R, 0.248888888888888889e-01_R} // 6
210 : };
211 :
212 :
213 : // Now call the keast routine to generate _points and _weights
214 213 : keast_rule(rule_data, 3);
215 :
216 6 : return;
217 56 : } // end if (allow_rules_with_negative_weights)
218 : // Note: if !allow_rules_with_negative_weights, fall through to next case.
219 : }
220 :
221 : libmesh_fallthrough();
222 :
223 :
224 : // Walkington's fifth-order 14-point rule from
225 : // "Quadrature on Simplices of Arbitrary Dimension"
226 : //
227 : // We originally had a Keast rule here, but this rule had
228 : // more points than an equivalent rule by Walkington and
229 : // also contained points on the boundary of the ref. elt,
230 : // making it less suitable for FEM calculations.
231 : case FIFTH:
232 : {
233 40308 : _points.resize(14);
234 40308 : _weights.resize(14);
235 :
236 : // permutations of these points and suitably-modified versions of
237 : // these points are the quadrature point locations
238 40308 : const Real a[3] = {0.31088591926330060980_R, // a1 from the paper
239 : 0.092735250310891226402_R, // a2 from the paper
240 : 0.045503704125649649492_R}; // a3 from the paper
241 :
242 : // weights. a[] and wt[] are the only floating-point inputs required
243 : // for this rule.
244 40308 : const Real wt[3] = {0.018781320953002641800_R, // w1 from the paper
245 : 0.012248840519393658257_R, // w2 from the paper
246 : 0.0070910034628469110730_R}; // w3 from the paper
247 :
248 : // The first two sets of 4 points are formed in a similar manner
249 120924 : for (unsigned int i=0; i<2; ++i)
250 : {
251 : // Where we will insert values into _points and _weights
252 80616 : const unsigned int offset=4*i;
253 :
254 : // Stuff points and weights values into their arrays
255 80616 : const Real b = 1. - 3.*a[i];
256 :
257 : // Here are the permutations. Order of these is not important,
258 : // all have the same weight
259 80616 : _points[offset + 0] = Point(a[i], a[i], a[i]);
260 80616 : _points[offset + 1] = Point(a[i], b, a[i]);
261 80616 : _points[offset + 2] = Point( b, a[i], a[i]);
262 86426 : _points[offset + 3] = Point(a[i], a[i], b);
263 :
264 : // These 4 points all have the same weights
265 403080 : for (unsigned int j=0; j<4; ++j)
266 345704 : _weights[offset + j] = wt[i];
267 : } // end for
268 :
269 :
270 : {
271 : // The third set contains 6 points and is formed a little differently
272 2905 : const unsigned int offset = 8;
273 2905 : const Real b = 0.5*(1. - 2.*a[2]);
274 :
275 : // Here are the permutations. Order of these is not important,
276 : // all have the same weight
277 40308 : _points[offset + 0] = Point(b , b, a[2]);
278 40308 : _points[offset + 1] = Point(b , a[2], a[2]);
279 40308 : _points[offset + 2] = Point(a[2], a[2], b);
280 40308 : _points[offset + 3] = Point(a[2], b, a[2]);
281 40308 : _points[offset + 4] = Point( b, a[2], b);
282 40308 : _points[offset + 5] = Point(a[2], b, b);
283 :
284 : // These 6 points all have the same weights
285 282156 : for (unsigned int j=0; j<6; ++j)
286 259278 : _weights[offset + j] = wt[2];
287 : }
288 :
289 :
290 : // Original rule by Keast, unsuitable because it has points on the
291 : // reference element boundary!
292 : // _points.resize(15);
293 : // _weights.resize(15);
294 :
295 : // _points[0](0) = 0.25;
296 : // _points[0](1) = 0.25;
297 : // _points[0](2) = 0.25;
298 :
299 : // {
300 : // const Real a = 0.;
301 : // const Real b = Real(1)/3;
302 :
303 : // _points[1](0) = a;
304 : // _points[1](1) = b;
305 : // _points[1](2) = b;
306 :
307 : // _points[2](0) = b;
308 : // _points[2](1) = a;
309 : // _points[2](2) = b;
310 :
311 : // _points[3](0) = b;
312 : // _points[3](1) = b;
313 : // _points[3](2) = a;
314 :
315 : // _points[4](0) = b;
316 : // _points[4](1) = b;
317 : // _points[4](2) = b;
318 : // }
319 : // {
320 : // const Real a = Real(8)/11;
321 : // const Real b = Real(1)/11;
322 :
323 : // _points[5](0) = a;
324 : // _points[5](1) = b;
325 : // _points[5](2) = b;
326 :
327 : // _points[6](0) = b;
328 : // _points[6](1) = a;
329 : // _points[6](2) = b;
330 :
331 : // _points[7](0) = b;
332 : // _points[7](1) = b;
333 : // _points[7](2) = a;
334 :
335 : // _points[8](0) = b;
336 : // _points[8](1) = b;
337 : // _points[8](2) = b;
338 : // }
339 : // {
340 : // const Real a = 0.066550153573664;
341 : // const Real b = 0.433449846426336;
342 :
343 : // _points[9](0) = b;
344 : // _points[9](1) = a;
345 : // _points[9](2) = a;
346 :
347 : // _points[10](0) = a;
348 : // _points[10](1) = a;
349 : // _points[10](2) = b;
350 :
351 : // _points[11](0) = a;
352 : // _points[11](1) = b;
353 : // _points[11](2) = b;
354 :
355 : // _points[12](0) = b;
356 : // _points[12](1) = a;
357 : // _points[12](2) = b;
358 :
359 : // _points[13](0) = b;
360 : // _points[13](1) = b;
361 : // _points[13](2) = a;
362 :
363 : // _points[14](0) = a;
364 : // _points[14](1) = b;
365 : // _points[14](2) = a;
366 : // }
367 :
368 : // _weights[0] = 0.030283678097089;
369 : // _weights[1] = 0.006026785714286;
370 : // _weights[2] = _weights[1];
371 : // _weights[3] = _weights[1];
372 : // _weights[4] = _weights[1];
373 : // _weights[5] = 0.011645249086029;
374 : // _weights[6] = _weights[5];
375 : // _weights[7] = _weights[5];
376 : // _weights[8] = _weights[5];
377 : // _weights[9] = 0.010949141561386;
378 : // _weights[10] = _weights[9];
379 : // _weights[11] = _weights[9];
380 : // _weights[12] = _weights[9];
381 : // _weights[13] = _weights[9];
382 : // _weights[14] = _weights[9];
383 :
384 2905 : return;
385 : }
386 :
387 :
388 :
389 :
390 : // This rule is originally from Keast:
391 : // Patrick Keast,
392 : // Moderate Degree Tetrahedral Quadrature Formulas,
393 : // Computer Methods in Applied Mechanics and Engineering,
394 : // Volume 55, Number 3, May 1986, pages 339-348.
395 : //
396 : // It is accurate on 6th-degree polynomials and has 24 points
397 : // vs. 64 for the comparable conical product rule.
398 : //
399 : // Values copied 24th June 2008 from:
400 : // http://people.scs.fsu.edu/~burkardt/f_src/keast/keast.f90
401 639 : case SIXTH:
402 : {
403 639 : _points.resize (24);
404 639 : _weights.resize(24);
405 :
406 : // The raw data for the quadrature rule.
407 639 : const Real rule_data[4][4] = {
408 : {0.356191386222544953e+00_R , 0.214602871259151684e+00_R , 0._R, 0.00665379170969464506e+00_R}, // 4
409 : {0.877978124396165982e+00_R , 0.0406739585346113397e+00_R, 0._R, 0.00167953517588677620e+00_R}, // 4
410 : {0.0329863295731730594e+00_R, 0.322337890142275646e+00_R , 0._R, 0.00922619692394239843e+00_R}, // 4
411 : {0.0636610018750175299e+00_R, 0.269672331458315867e+00_R , 0.603005664791649076e+00_R, 0.00803571428571428248e+00_R} // 12
412 : };
413 :
414 :
415 : // Now call the keast routine to generate _points and _weights
416 639 : keast_rule(rule_data, 4);
417 :
418 18 : return;
419 : }
420 :
421 :
422 : // Keast's 31 point, 7th-order rule contains points on the reference
423 : // element boundary, so we've decided not to include it here.
424 : //
425 : // Keast's 8th-order rule has 45 points. and a negative
426 : // weight, so if you've explicitly disallowed such rules
427 : // you will fall through to the conical product rules
428 : // below.
429 36382 : case SEVENTH:
430 : case EIGHTH:
431 : {
432 36382 : if (allow_rules_with_negative_weights)
433 : {
434 36382 : _points.resize (45);
435 36382 : _weights.resize(45);
436 :
437 : // The raw data for the quadrature rule.
438 36382 : const Real rule_data[7][4] = {
439 : {0.250000000000000000e+00_R, 0._R, 0._R, -0.393270066412926145e-01_R}, // 1
440 : {0.617587190300082967e+00_R, 0.127470936566639015e+00_R, 0._R, 0.408131605934270525e-02_R}, // 4
441 : {0.903763508822103123e+00_R, 0.320788303926322960e-01_R, 0._R, 0.658086773304341943e-03_R}, // 4
442 : {0.497770956432810185e-01_R, 0._R, 0.450222904356718978e+00_R, 0.438425882512284693e-02_R}, // 6
443 : {0.183730447398549945e+00_R, 0._R, 0.316269552601450060e+00_R, 0.138300638425098166e-01_R}, // 6
444 : {0.231901089397150906e+00_R, 0.229177878448171174e-01_R, 0.513280033360881072e+00_R, 0.424043742468372453e-02_R}, // 12
445 : {0.379700484718286102e-01_R, 0.730313427807538396e+00_R, 0.193746475248804382e+00_R, 0.223873973961420164e-02_R} // 12
446 : };
447 :
448 :
449 : // Now call the keast routine to generate _points and _weights
450 36382 : keast_rule(rule_data, 7);
451 :
452 2784 : return;
453 0 : } // end if (allow_rules_with_negative_weights)
454 : // Note: if !allow_rules_with_negative_weights, fall through to next case.
455 : }
456 :
457 : libmesh_fallthrough();
458 :
459 :
460 : // Fall back on Grundmann-Moller or Conical Product rules at high orders.
461 : default:
462 : {
463 6110 : if ((allow_rules_with_negative_weights) && (get_order() < 34))
464 : {
465 : // The Grundmann-Moller rules are defined to arbitrary order and
466 : // can have significantly fewer evaluation points than conical product
467 : // rules. If you allow rules with negative weights, the GM rules
468 : // will be more efficient for degree up to 33 (for degree 34 and
469 : // higher, CP is more efficient!) but may be more susceptible
470 : // to round-off error. Safest is to disallow rules with negative
471 : // weights, but this decision should be made on a case-by-case basis.
472 260 : QGrundmann_Moller gm_rule(3, _order);
473 6110 : gm_rule.init(*this);
474 :
475 : // Swap points and weights with the about-to-be destroyed rule.
476 260 : _points.swap (gm_rule.get_points() );
477 260 : _weights.swap(gm_rule.get_weights());
478 :
479 260 : return;
480 : }
481 :
482 : else
483 : {
484 : // The following quadrature rules are generated as
485 : // conical products. These tend to be non-optimal
486 : // (use too many points, cluster points in certain
487 : // regions of the domain) but they are quite easy to
488 : // automatically generate using a 1D Gauss rule on
489 : // [0,1] and two 1D Jacobi-Gauss rules on [0,1].
490 0 : QConical conical_rule(3, _order);
491 0 : conical_rule.init(*this);
492 :
493 : // Swap points and weights with the about-to-be destroyed rule.
494 0 : _points.swap (conical_rule.get_points() );
495 0 : _weights.swap(conical_rule.get_weights());
496 :
497 0 : return;
498 : }
499 : }
500 : }
501 : } // end case TET
502 :
503 :
504 :
505 : //---------------------------------------------
506 : // Prism quadrature rules
507 26046 : case PRISM6:
508 : case PRISM15:
509 : case PRISM18:
510 : case PRISM20:
511 : case PRISM21:
512 : {
513 : // We compute the 3D quadrature rule as a tensor
514 : // product of the 1D quadrature rule and a 2D
515 : // triangle quadrature rule
516 :
517 29832 : QGauss q1D(1,get_order());
518 27939 : QGauss q2D(2,_order);
519 :
520 : // Initialize the 2D rule (1D is pre-initialized)
521 26046 : q2D.init(TRI3, _p_level, /*simple_type_only=*/true);
522 :
523 26046 : tensor_product_prism(q1D, q2D);
524 :
525 1893 : return;
526 : }
527 :
528 :
529 :
530 : //---------------------------------------------
531 : // Pyramid
532 18819 : case PYRAMID5:
533 : case PYRAMID13:
534 : case PYRAMID14:
535 : case PYRAMID18:
536 : {
537 : // We compute the Pyramid rule as a conical product of a
538 : // Jacobi rule with alpha==2 on the interval [0,1] two 1D
539 : // Gauss rules with suitably modified points. The idea comes
540 : // from: Stroud, A.H. "Approximate Calculation of Multiple
541 : // Integrals."
542 19837 : QConical conical_rule(3, _order);
543 18819 : conical_rule.init(*this);
544 :
545 : // Swap points and weights with the about-to-be destroyed rule.
546 1018 : _points.swap (conical_rule.get_points() );
547 1018 : _weights.swap(conical_rule.get_weights());
548 :
549 1018 : return;
550 :
551 : }
552 :
553 :
554 : //---------------------------------------------
555 : // Arbitrary polyhedron quadrature rules
556 587 : case C0POLYHEDRON:
557 : {
558 632 : QGauss tet_rule(3, _order);
559 587 : tet_rule.init(TET4, _p_level, true);
560 :
561 45 : std::vector<Point> & tetpoints = tet_rule.get_points();
562 45 : std::vector<Real> & tetweights = tet_rule.get_weights();
563 :
564 90 : std::size_t numtetpts = tetpoints.size();
565 :
566 : // C0Polyhedron requires the newer Quadrature API
567 587 : if (!_elem)
568 0 : libmesh_error();
569 :
570 45 : libmesh_assert(_elem->type() == C0POLYHEDRON);
571 :
572 45 : const C0Polyhedron & poly = *cast_ptr<const C0Polyhedron *>(_elem);
573 :
574 587 : std::size_t numtets = poly.n_subelements();
575 587 : _points.resize(numtetpts*numtets);
576 587 : _weights.resize(numtetpts*numtets);
577 :
578 3522 : for (std::size_t t = 0; t != numtets; ++t)
579 : {
580 2935 : auto master_points = poly.master_subelement(t);
581 :
582 225 : const Point v01 = master_points[1] - master_points[0];
583 225 : const Point v02 = master_points[2] - master_points[0];
584 225 : const Point v03 = master_points[3] - master_points[0];
585 :
586 : // The factor of one sixth from the tetweights cancels out
587 : // the factor of six here, so we don't need to do so
588 : // ourselves.
589 : const Real six_master_tet_vol =
590 2710 : triple_product(v01, v02, v03);
591 :
592 24630 : for (std::size_t i = 0; i != numtetpts; ++i)
593 : {
594 21695 : _points[numtetpts*t+i] =
595 1710 : master_points[0] +
596 5130 : v01 * tetpoints[i](0) +
597 5130 : v02 * tetpoints[i](1) +
598 6840 : v03 * tetpoints[i](2);
599 25115 : _weights[numtetpts*t+i] = tetweights[i] *
600 : six_master_tet_vol;
601 : }
602 : }
603 45 : return;
604 : }
605 :
606 : //---------------------------------------------
607 : // Unsupported type
608 0 : default:
609 0 : libmesh_error_msg("ERROR: Unsupported type: " << Utility::enum_to_string(_type));
610 : }
611 : #endif
612 : }
613 :
614 : } // namespace libMesh
|
