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_monomial.h"
22 : #include "libmesh/quadrature_gauss.h"
23 :
24 : namespace libMesh
25 : {
26 :
27 :
28 4828 : void QMonomial::init_3D()
29 : {
30 :
31 4828 : switch (_type)
32 : {
33 : //---------------------------------------------
34 : // Hex quadrature rules
35 1207 : case HEX8:
36 : case HEX20:
37 : case HEX27:
38 : {
39 68 : switch(get_order())
40 : {
41 :
42 : // The CONSTANT/FIRST rule is the 1-point Gauss "product" rule...we fall
43 : // through to the default case for this rule.
44 :
45 142 : case SECOND:
46 : case THIRD:
47 : {
48 : // A degree 3, 6-point, "rotationally-symmetric" rule by
49 : // Kim and Song, Comm. Korean Math. Soc vol. 13, no. 4, 1998, pp. 913-931.
50 : //
51 : // Warning: this rule contains points on the boundary of the reference
52 : // element, and therefore may be unsuitable for some problems. The alternative
53 : // would be a 2x2x2 Gauss product rule.
54 142 : const Real data[1][4] =
55 : {
56 : {1, 0, 0, Real(4)/3}
57 : };
58 :
59 142 : const unsigned int rule_id[1] = {
60 : 1 // (x,0,0) -> 6 permutations
61 : };
62 :
63 142 : _points.resize(6);
64 142 : _weights.resize(6);
65 :
66 142 : kim_rule(data, rule_id, 1);
67 4 : return;
68 : } // end case SECOND,THIRD
69 :
70 142 : case FOURTH:
71 : case FIFTH:
72 : {
73 : // A degree 5, 13-point rule by Stroud,
74 : // AH Stroud, "Some Fifth Degree Integration Formulas for Symmetric Regions II.",
75 : // Numerische Mathematik 9, pp. 460-468 (1967).
76 : //
77 : // This rule is provably minimal in the number of points. The equations given for
78 : // the n-cube on pg. 466 of the paper for mu/gamma and gamma are wrong, at least for
79 : // the n=3 case. The analytical values given here were computed by me [JWP] in Maple.
80 :
81 : // Convenient intermediate values.
82 4 : const Real sqrt19 = std::sqrt(Real(19)); // ~4.35889894354
83 4 : const Real tp = std::sqrt(71440 + 6802*sqrt19); // ~317.945326454
84 :
85 : // Point data for permutations.
86 4 : const Real eta = 0;
87 :
88 4 : const Real lambda = std::sqrt(Real(1919)/3285 - 148*sqrt19/3285 + 4*tp/3285);
89 : // 8.8030440669930978047737818209860e-01_R;
90 :
91 4 : const Real xi = -std::sqrt(Real(1121)/3285 + 74*sqrt19/3285 - 2*tp/3285);
92 : // -4.9584817142571115281421242364290e-01_R;
93 :
94 4 : const Real mu = std::sqrt(Real(1121)/3285 + 74*sqrt19/3285 + 2*tp/3285);
95 : // 7.9562142216409541542982482567580e-01_R;
96 :
97 4 : const Real gamma = std::sqrt(Real(1919)/3285 - 148*sqrt19/3285 - 4*tp/3285);
98 : // 2.5293711744842581347389255929324e-02_R;
99 :
100 : // Weights: the centroid weight is given analytically. Weight B (resp C) goes
101 : // with the {lambda,xi} (resp {gamma,mu}) permutation. The single-precision
102 : // results reported by Stroud are given for reference.
103 :
104 4 : const Real A = Real(32)/19; // ~1.684210560
105 : // Stroud: 0.21052632 * 8.0;
106 :
107 4 : const Real B = Real(1) / (Real(260072)/133225 - 1520*sqrt19/133225 + (133 - 37*sqrt19)*tp/133225);
108 : // 5.4498735127757671684690782180890e-01_R; // Stroud: 0.068123420 * 8.0 = 0.544987360;
109 :
110 4 : const Real C = Real(1) / (Real(260072)/133225 - 1520*sqrt19/133225 - (133 - 37*sqrt19)*tp/133225);
111 : // 5.0764422766979170420572375713840e-01_R; // Stroud: 0.063455527 * 8.0 = 0.507644216;
112 :
113 142 : _points.resize(13);
114 142 : _weights.resize(13);
115 :
116 4 : unsigned int c=0;
117 :
118 : // Point with weight A (origin)
119 142 : _points[c] = Point(eta, eta, eta);
120 142 : _weights[c++] = A;
121 :
122 : // Points with weight B
123 142 : _points[c] = Point(lambda, xi, xi);
124 142 : _weights[c++] = B;
125 142 : _points[c] = -_points[c-1];
126 142 : _weights[c++] = B;
127 :
128 142 : _points[c] = Point(xi, lambda, xi);
129 142 : _weights[c++] = B;
130 142 : _points[c] = -_points[c-1];
131 142 : _weights[c++] = B;
132 :
133 142 : _points[c] = Point(xi, xi, lambda);
134 142 : _weights[c++] = B;
135 142 : _points[c] = -_points[c-1];
136 142 : _weights[c++] = B;
137 :
138 : // Points with weight C
139 142 : _points[c] = Point(mu, mu, gamma);
140 142 : _weights[c++] = C;
141 142 : _points[c] = -_points[c-1];
142 142 : _weights[c++] = C;
143 :
144 142 : _points[c] = Point(mu, gamma, mu);
145 142 : _weights[c++] = C;
146 142 : _points[c] = -_points[c-1];
147 142 : _weights[c++] = C;
148 :
149 142 : _points[c] = Point(gamma, mu, mu);
150 142 : _weights[c++] = C;
151 142 : _points[c] = -_points[c-1];
152 142 : _weights[c++] = C;
153 :
154 142 : return;
155 :
156 :
157 : // // A degree 5, 14-point, "rotationally-symmetric" rule by
158 : // // Kim and Song, Comm. Korean Math. Soc vol. 13, no. 4, 1998, pp. 913-931.
159 : // // Was also reported in Stroud's 1971 book.
160 : // const Real data[2][4] =
161 : // {
162 : // {7.95822425754221463264548820476135e-01_R, 0.00000000000000000000000000000000e+00_R, 0.00000000000000000000000000000000e+00_R, 8.86426592797783933518005540166204e-01_R},
163 : // {7.58786910639328146269034278112267e-01_R, 7.58786910639328146269034278112267e-01_R, 7.58786910639328146269034278112267e-01_R, 3.35180055401662049861495844875346e-01_R}
164 : // };
165 :
166 : // const unsigned int rule_id[2] = {
167 : // 1, // (x,0,0) -> 6 permutations
168 : // 4 // (x,x,x) -> 8 permutations
169 : // };
170 :
171 : // _points.resize(14);
172 : // _weights.resize(14);
173 :
174 : // kim_rule(data, rule_id, 2);
175 : // return;
176 : } // end case FOURTH,FIFTH
177 :
178 :
179 142 : case SIXTH:
180 : case SEVENTH:
181 : {
182 142 : if (allow_rules_with_negative_weights)
183 : {
184 : // A degree 7, 31-point, "rotationally-symmetric" rule by
185 : // Kim and Song, Comm. Korean Math. Soc vol. 13, no. 4, 1998, pp. 913-931.
186 : // This rule contains a negative weight, so only use it if such type of
187 : // rules are allowed.
188 142 : const Real data[3][4] =
189 : {
190 : {0.00000000000000000000000000000000e+00_R, 0.00000000000000000000000000000000e+00_R, 0.00000000000000000000000000000000e+00_R, -1.27536231884057971014492753623188e+00_R},
191 : {5.85540043769119907612630781744060e-01_R, 0.00000000000000000000000000000000e+00_R, 0.00000000000000000000000000000000e+00_R, 8.71111111111111111111111111111111e-01_R},
192 : {6.94470135991704766602025803883310e-01_R, 9.37161638568208038511047377665396e-01_R, 4.15659267604065126239606672567031e-01_R, 1.68695652173913043478260869565217e-01_R}
193 : };
194 :
195 142 : const unsigned int rule_id[3] = {
196 : 0, // (0,0,0) -> 1 permutation
197 : 1, // (x,0,0) -> 6 permutations
198 : 6 // (x,y,z) -> 24 permutations
199 : };
200 :
201 142 : _points.resize(31);
202 142 : _weights.resize(31);
203 :
204 142 : kim_rule(data, rule_id, 3);
205 4 : return;
206 : } // end if (allow_rules_with_negative_weights)
207 :
208 :
209 : // A degree 7, 34-point, "fully-symmetric" rule, first published in
210 : // P.C. Hammer and A.H. Stroud, "Numerical Evaluation of Multiple Integrals II",
211 : // Mathematical Tables and Other Aids to Computation, vol 12., no 64, 1958, pp. 272-280
212 : //
213 : // This rule happens to fall under the same general
214 : // construction as the Kim rules, so we've re-used
215 : // that code here. Stroud gives 16 digits for his rule,
216 : // and this is the most accurate version I've found.
217 : //
218 : // For comparison, a SEVENTH-order Gauss product rule
219 : // (which integrates tri-7th order polynomials) would
220 : // have 4^3=64 points.
221 : const Real
222 0 : r = std::sqrt(Real(6)/7), // ~0.92582009977
223 0 : s = std::sqrt((Real(960) - 3*std::sqrt(Real(28798))) / 2726), // ~0.40670318642
224 0 : t = std::sqrt((Real(960) + 3*std::sqrt(Real(28798))) / 2726), // ~0.73411252875
225 0 : B1 = Real(8624)/29160, // ~0.29574759945
226 0 : B2 = Real(2744)/29160, // ~0.09410150891
227 0 : B3 = 8*(774*t*t - 230)/(9720*(t*t-s*s)), // ~0.41233386227
228 0 : B4 = 8*(230 - 774*s*s)/(9720*(t*t-s*s)); // ~0.22470317477
229 :
230 0 : const Real data[4][4] =
231 : {
232 : {r, 0, 0, B1},
233 : {r, r, 0, B2},
234 : {s, s, s, B3},
235 : {t, t, t, B4}
236 : };
237 :
238 0 : const unsigned int rule_id[4] = {
239 : 1, // (x,0,0) -> 6 permutations
240 : 2, // (x,x,0) -> 12 permutations
241 : 4, // (x,x,x) -> 8 permutations
242 : 4 // (x,x,x) -> 8 permutations
243 : };
244 :
245 0 : _points.resize(34);
246 0 : _weights.resize(34);
247 :
248 0 : kim_rule(data, rule_id, 4);
249 0 : return;
250 :
251 :
252 : // // A degree 7, 38-point, "rotationally-symmetric" rule by
253 : // // Kim and Song, Comm. Korean Math. Soc vol. 13, no. 4, 1998, pp. 913-931.
254 : // //
255 : // // This rule is obviously inferior to the 34-point rule above...
256 : // const Real data[3][4] =
257 : //{
258 : // {9.01687807821291289082811566285950e-01_R, 0.00000000000000000000000000000000e+00_R, 0.00000000000000000000000000000000e+00_R, 2.95189738262622903181631100062774e-01_R},
259 : // {4.08372221499474674069588900002128e-01_R, 4.08372221499474674069588900002128e-01_R, 4.08372221499474674069588900002128e-01_R, 4.04055417266200582425904380777126e-01_R},
260 : // {8.59523090201054193116477875786220e-01_R, 8.59523090201054193116477875786220e-01_R, 4.14735913727987720499709244748633e-01_R, 1.24850759678944080062624098058597e-01_R}
261 : //};
262 : //
263 : // const unsigned int rule_id[3] = {
264 : //1, // (x,0,0) -> 6 permutations
265 : //4, // (x,x,x) -> 8 permutations
266 : //5 // (x,x,z) -> 24 permutations
267 : // };
268 : //
269 : // _points.resize(38);
270 : // _weights.resize(38);
271 : //
272 : // kim_rule(data, rule_id, 3);
273 : // return;
274 : } // end case SIXTH,SEVENTH
275 :
276 71 : case EIGHTH:
277 : {
278 : // A degree 8, 47-point, "rotationally-symmetric" rule by
279 : // Kim and Song, Comm. Korean Math. Soc vol. 13, no. 4, 1998, pp. 913-931.
280 : //
281 : // A EIGHTH-order Gauss product rule (which integrates tri-8th order polynomials)
282 : // would have 5^3=125 points.
283 71 : const Real data[5][4] =
284 : {
285 : {0.00000000000000000000000000000000e+00_R, 0.00000000000000000000000000000000e+00_R, 0.00000000000000000000000000000000e+00_R, 4.51903714875199690490763818699555e-01_R},
286 : {7.82460796435951590652813975429717e-01_R, 0.00000000000000000000000000000000e+00_R, 0.00000000000000000000000000000000e+00_R, 2.99379177352338919703385618576171e-01_R},
287 : {4.88094669706366480526729301468686e-01_R, 4.88094669706366480526729301468686e-01_R, 4.88094669706366480526729301468686e-01_R, 3.00876159371240019939698689791164e-01_R},
288 : {8.62218927661481188856422891110042e-01_R, 8.62218927661481188856422891110042e-01_R, 8.62218927661481188856422891110042e-01_R, 4.94843255877038125738173175714853e-02_R},
289 : {2.81113909408341856058098281846420e-01_R, 9.44196578292008195318687494773744e-01_R, 6.97574833707236996779391729948984e-01_R, 1.22872389222467338799199767122592e-01_R}
290 : };
291 :
292 71 : const unsigned int rule_id[5] = {
293 : 0, // (0,0,0) -> 1 permutation
294 : 1, // (x,0,0) -> 6 permutations
295 : 4, // (x,x,x) -> 8 permutations
296 : 4, // (x,x,x) -> 8 permutations
297 : 6 // (x,y,z) -> 24 permutations
298 : };
299 :
300 71 : _points.resize(47);
301 71 : _weights.resize(47);
302 :
303 71 : kim_rule(data, rule_id, 5);
304 2 : return;
305 : } // end case EIGHTH
306 :
307 :
308 : // By default: construct and use a Gauss quadrature rule
309 20 : default:
310 : {
311 : // Break out and fall down into the default: case for the
312 : // outer switch statement.
313 20 : break;
314 : }
315 :
316 244 : } // end switch(_order + 2*p)
317 : } // end case HEX8/20/27
318 :
319 : libmesh_fallthrough();
320 :
321 : // By default: construct and use a Gauss quadrature rule
322 : default:
323 : {
324 4453 : QGauss gauss_rule(3, _order);
325 4331 : gauss_rule.init(*this);
326 :
327 : // Swap points and weights with the about-to-be destroyed rule.
328 122 : _points.swap (gauss_rule.get_points() );
329 122 : _weights.swap(gauss_rule.get_weights());
330 :
331 122 : return;
332 : }
333 : } // end switch (_type)
334 : }
335 :
336 : } // namespace libMesh
|
