1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2020 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
init_2D(const ElemType,unsigned int)28 void QMonomial::init_2D(const ElemType, unsigned int)
29 {
30
31 switch (_type)
32 {
33 //---------------------------------------------
34 // Quadrilateral quadrature rules
35 case QUAD4:
36 case QUADSHELL4:
37 case QUAD8:
38 case QUADSHELL8:
39 case QUAD9:
40 {
41 switch(get_order())
42 {
43 case SECOND:
44 {
45 // A degree=2 rule for the QUAD with 3 points.
46 // A tensor product degree-2 Gauss would have 4 points.
47 // This rule (or a variation on it) is probably available in
48 //
49 // A.H. Stroud, Approximate calculation of multiple integrals,
50 // Prentice-Hall, Englewood Cliffs, N.J., 1971.
51 //
52 // though I have never actually seen a reference for it.
53 // Luckily it's fairly easy to derive, which is what I've done
54 // here [JWP].
55 const Real
56 s=std::sqrt(Real(1)/3),
57 t=std::sqrt(Real(2)/3);
58
59 const Real data[2][3] =
60 {
61 {0.0, s, 2.0},
62 { t, -s, 1.0}
63 };
64
65 _points.resize(3);
66 _weights.resize(3);
67
68 wissmann_rule(data, 2);
69
70 return;
71 } // end case SECOND
72
73
74
75 // For third-order, fall through to default case, use 2x2 Gauss product rule.
76 // case THIRD:
77 // {
78 // } // end case THIRD
79
80 // Tabulated-in-double-precision rules aren't accurate enough for
81 // higher precision, so fall back on Gauss
82 #if !defined(LIBMESH_DEFAULT_TRIPLE_PRECISION) && !defined(LIBMESH_DEFAULT_QUADRUPLE_PRECISION)
83 case FOURTH:
84 {
85 // A pair of degree=4 rules for the QUAD "C2" due to
86 // Wissmann and Becker. These rules both have six points.
87 // A tensor product degree-4 Gauss would have 9 points.
88 //
89 // J. W. Wissmann and T. Becker, Partially symmetric cubature
90 // formulas for even degrees of exactness, SIAM J. Numer. Anal. 23
91 // (1986), 676--685.
92 const Real data[4][3] =
93 {
94 // First of 2 degree-4 rules given by Wissmann
95 {Real(0.0000000000000000e+00), Real(0.0000000000000000e+00), Real(1.1428571428571428e+00)},
96 {Real(0.0000000000000000e+00), Real(9.6609178307929590e-01), Real(4.3956043956043956e-01)},
97 {Real(8.5191465330460049e-01), Real(4.5560372783619284e-01), Real(5.6607220700753210e-01)},
98 {Real(6.3091278897675402e-01), Real(-7.3162995157313452e-01), Real(6.4271900178367668e-01)}
99 //
100 // Second of 2 degree-4 rules given by Wissmann. These both
101 // yield 4th-order accurate rules, I just chose the one that
102 // happened to contain the origin.
103 // {0.000000000000000, -0.356822089773090, 1.286412084888852},
104 // {0.000000000000000, 0.934172358962716, 0.491365692888926},
105 // {0.774596669241483, 0.390885162530071, 0.761883709085613},
106 // {0.774596669241483, -0.852765377881771, 0.349227402025498}
107 };
108
109 _points.resize(6);
110 _weights.resize(6);
111
112 wissmann_rule(data, 4);
113
114 return;
115 } // end case FOURTH
116 #endif
117
118
119
120
121 case FIFTH:
122 {
123 // A degree 5, 7-point rule due to Stroud.
124 //
125 // A.H. Stroud, Approximate calculation of multiple integrals,
126 // Prentice-Hall, Englewood Cliffs, N.J., 1971.
127 //
128 // This rule is provably minimal in the number of points.
129 // A tensor-product rule accurate for "bi-quintic" polynomials would have 9 points.
130 const Real data[3][3] =
131 {
132 { 0.L, 0.L, Real(8)/7 }, // 1
133 { 0.L, std::sqrt(Real(14)/15), Real(20)/63}, // 2
134 {std::sqrt(Real(3)/5), std::sqrt(Real(1)/3), Real(20)/36} // 4
135 };
136
137 const unsigned int symmetry[3] = {
138 0, // Origin
139 7, // Central Symmetry
140 6 // Rectangular
141 };
142
143 _points.resize (7);
144 _weights.resize(7);
145
146 stroud_rule(data, symmetry, 3);
147
148 return;
149 } // end case FIFTH
150
151
152
153
154 // Tabulated-in-double-precision rules aren't accurate enough for
155 // higher precision, so fall back on Gauss
156 #if !defined(LIBMESH_DEFAULT_TRIPLE_PRECISION) && !defined(LIBMESH_DEFAULT_QUADRUPLE_PRECISION)
157 case SIXTH:
158 {
159 // A pair of degree=6 rules for the QUAD "C2" due to
160 // Wissmann and Becker. These rules both have 10 points.
161 // A tensor product degree-6 Gauss would have 16 points.
162 //
163 // J. W. Wissmann and T. Becker, Partially symmetric cubature
164 // formulas for even degrees of exactness, SIAM J. Numer. Anal. 23
165 // (1986), 676--685.
166 const Real data[6][3] =
167 {
168 // First of 2 degree-6, 10 point rules given by Wissmann
169 // {0.000000000000000, 0.836405633697626, 0.455343245714174},
170 // {0.000000000000000, -0.357460165391307, 0.827395973202966},
171 // {0.888764014654765, 0.872101531193131, 0.144000884599645},
172 // {0.604857639464685, 0.305985162155427, 0.668259104262665},
173 // {0.955447506641064, -0.410270899466658, 0.225474004890679},
174 // {0.565459993438754, -0.872869311156879, 0.320896396788441}
175 //
176 // Second of 2 degree-6, 10 point rules given by Wissmann.
177 // Either of these will work, I just chose the one with points
178 // slightly further into the element interior.
179 {Real(0.0000000000000000e+00), Real(8.6983337525005900e-01), Real(3.9275059096434794e-01)},
180 {Real(0.0000000000000000e+00), Real(-4.7940635161211124e-01), Real(7.5476288124261053e-01)},
181 {Real(8.6374282634615388e-01), Real(8.0283751620765670e-01), Real(2.0616605058827902e-01)},
182 {Real(5.1869052139258234e-01), Real(2.6214366550805818e-01), Real(6.8999213848986375e-01)},
183 {Real(9.3397254497284950e-01), Real(-3.6309658314806653e-01), Real(2.6051748873231697e-01)},
184 {Real(6.0897753601635630e-01), Real(-8.9660863276245265e-01), Real(2.6956758608606100e-01)}
185 };
186
187 _points.resize(10);
188 _weights.resize(10);
189
190 wissmann_rule(data, 6);
191
192 return;
193 } // end case SIXTH
194 #endif
195
196
197
198
199 case SEVENTH:
200 {
201 // A degree 7, 12-point rule due to Tyler, can be found in Stroud's book
202 //
203 // A.H. Stroud, Approximate calculation of multiple integrals,
204 // Prentice-Hall, Englewood Cliffs, N.J., 1971.
205 //
206 // This rule is fully-symmetric and provably minimal in the number of points.
207 // A tensor-product rule accurate for "bi-septic" polynomials would have 16 points.
208 const Real
209 r = std::sqrt(Real(6)/7),
210 s = std::sqrt( (114 - 3*std::sqrt(Real(583))) / 287 ),
211 t = std::sqrt( (114 + 3*std::sqrt(Real(583))) / 287 ),
212 B1 = Real(196)/810,
213 B2 = 4 * (178981 + 2769*std::sqrt(Real(583))) / 1888920,
214 B3 = 4 * (178981 - 2769*std::sqrt(Real(583))) / 1888920;
215
216 const Real data[3][3] =
217 {
218 {r, 0.0, B1}, // 4
219 {s, 0.0, B2}, // 4
220 {t, 0.0, B3} // 4
221 };
222
223 const unsigned int symmetry[3] = {
224 3, // Full Symmetry, (x,0)
225 2, // Full Symmetry, (x,x)
226 2 // Full Symmetry, (x,x)
227 };
228
229 _points.resize (12);
230 _weights.resize(12);
231
232 stroud_rule(data, symmetry, 3);
233
234 return;
235 } // end case SEVENTH
236
237
238
239
240 // Tabulated-in-double-precision rules aren't accurate enough for
241 // higher precision, so fall back on Gauss
242 #if !defined(LIBMESH_DEFAULT_TRIPLE_PRECISION) && !defined(LIBMESH_DEFAULT_QUADRUPLE_PRECISION)
243 case EIGHTH:
244 {
245 // A pair of degree=8 rules for the QUAD "C2" due to
246 // Wissmann and Becker. These rules both have 16 points.
247 // A tensor product degree-6 Gauss would have 25 points.
248 //
249 // J. W. Wissmann and T. Becker, Partially symmetric cubature
250 // formulas for even degrees of exactness, SIAM J. Numer. Anal. 23
251 // (1986), 676--685.
252 const Real data[10][3] =
253 {
254 // First of 2 degree-8, 16 point rules given by Wissmann
255 // {0.000000000000000, 0.000000000000000, 0.055364705621440},
256 // {0.000000000000000, 0.757629177660505, 0.404389368726076},
257 // {0.000000000000000, -0.236871842255702, 0.533546604952635},
258 // {0.000000000000000, -0.989717929044527, 0.117054188786739},
259 // {0.639091304900370, 0.950520955645667, 0.125614417613747},
260 // {0.937069076924990, 0.663882736885633, 0.136544584733588},
261 // {0.537083530541494, 0.304210681724104, 0.483408479211257},
262 // {0.887188506449625, -0.236496718536120, 0.252528506429544},
263 // {0.494698820670197, -0.698953476086564, 0.361262323882172},
264 // {0.897495818279768, -0.900390774211580, 0.085464254086247}
265 //
266 // Second of 2 degree-8, 16 point rules given by Wissmann.
267 // Either of these will work, I just chose the one with points
268 // further into the element interior.
269 {Real(0.0000000000000000e+00), Real(6.5956013196034176e-01), Real(4.5027677630559029e-01)},
270 {Real(0.0000000000000000e+00), Real(-9.4914292304312538e-01), Real(1.6657042677781274e-01)},
271 {Real(9.5250946607156228e-01), Real(7.6505181955768362e-01), Real(9.8869459933431422e-02)},
272 {Real(5.3232745407420624e-01), Real(9.3697598108841598e-01), Real(1.5369674714081197e-01)},
273 {Real(6.8473629795173504e-01), Real(3.3365671773574759e-01), Real(3.9668697607290278e-01)},
274 {Real(2.3314324080140552e-01), Real(-7.9583272377396852e-02), Real(3.5201436794569501e-01)},
275 {Real(9.2768331930611748e-01), Real(-2.7224008061253425e-01), Real(1.8958905457779799e-01)},
276 {Real(4.5312068740374942e-01), Real(-6.1373535339802760e-01), Real(3.7510100114758727e-01)},
277 {Real(8.3750364042281223e-01), Real(-8.8847765053597136e-01), Real(1.2561879164007201e-01)}
278 };
279
280 _points.resize(16);
281 _weights.resize(16);
282
283 wissmann_rule(data, /*10*/ 9);
284
285 return;
286 } // end case EIGHTH
287
288
289
290
291 case NINTH:
292 {
293 // A degree 9, 17-point rule due to Moller.
294 //
295 // H.M. Moller, Kubaturformeln mit minimaler Knotenzahl,
296 // Numer. Math. 25 (1976), 185--200.
297 //
298 // This rule is provably minimal in the number of points.
299 // A tensor-product rule accurate for "bi-ninth" degree polynomials would have 25 points.
300 const Real data[5][3] =
301 {
302 {Real(0.0000000000000000e+00), Real(0.0000000000000000e+00), Real(5.2674897119341563e-01)}, // 1
303 {Real(6.3068011973166885e-01), Real(9.6884996636197772e-01), Real(8.8879378170198706e-02)}, // 4
304 {Real(9.2796164595956966e-01), Real(7.5027709997890053e-01), Real(1.1209960212959648e-01)}, // 4
305 {Real(4.5333982113564719e-01), Real(5.2373582021442933e-01), Real(3.9828243926207009e-01)}, // 4
306 {Real(8.5261572933366230e-01), Real(7.6208328192617173e-02), Real(2.6905133763978080e-01)} // 4
307 };
308
309 const unsigned int symmetry[5] = {
310 0, // Single point
311 4, // Rotational Invariant
312 4, // Rotational Invariant
313 4, // Rotational Invariant
314 4 // Rotational Invariant
315 };
316
317 _points.resize (17);
318 _weights.resize(17);
319
320 stroud_rule(data, symmetry, 5);
321
322 return;
323 } // end case NINTH
324
325
326
327
328 case TENTH:
329 case ELEVENTH:
330 {
331 // A degree 11, 24-point rule due to Cools and Haegemans.
332 //
333 // R. Cools and A. Haegemans, Another step forward in searching for
334 // cubature formulae with a minimal number of knots for the square,
335 // Computing 40 (1988), 139--146.
336 //
337 // P. Verlinden and R. Cools, The algebraic construction of a minimal
338 // cubature formula of degree 11 for the square, Cubature Formulas
339 // and their Applications (Russian) (Krasnoyarsk) (M.V. Noskov, ed.),
340 // 1994, pp. 13--23.
341 //
342 // This rule is provably minimal in the number of points.
343 // A tensor-product rule accurate for "bi-tenth" or "bi-eleventh" degree polynomials would have 36 points.
344 const Real data[6][3] =
345 {
346 {Real(6.9807610454956756e-01), Real(9.8263922354085547e-01), Real(4.8020763350723814e-02)}, // 4
347 {Real(9.3948638281673690e-01), Real(8.2577583590296393e-01), Real(6.6071329164550595e-02)}, // 4
348 {Real(9.5353952820153201e-01), Real(1.8858613871864195e-01), Real(9.7386777358668164e-02)}, // 4
349 {Real(3.1562343291525419e-01), Real(8.1252054830481310e-01), Real(2.1173634999894860e-01)}, // 4
350 {Real(7.1200191307533630e-01), Real(5.2532025036454776e-01), Real(2.2562606172886338e-01)}, // 4
351 {Real(4.2484724884866925e-01), Real(4.1658071912022368e-02), Real(3.5115871839824543e-01)} // 4
352 };
353
354 const unsigned int symmetry[6] = {
355 4, // Rotational Invariant
356 4, // Rotational Invariant
357 4, // Rotational Invariant
358 4, // Rotational Invariant
359 4, // Rotational Invariant
360 4 // Rotational Invariant
361 };
362
363 _points.resize (24);
364 _weights.resize(24);
365
366 stroud_rule(data, symmetry, 6);
367
368 return;
369 } // end case TENTH,ELEVENTH
370
371
372
373
374 case TWELFTH:
375 case THIRTEENTH:
376 {
377 // A degree 13, 33-point rule due to Cools and Haegemans.
378 //
379 // R. Cools and A. Haegemans, Another step forward in searching for
380 // cubature formulae with a minimal number of knots for the square,
381 // Computing 40 (1988), 139--146.
382 //
383 // A tensor-product rule accurate for "bi-12" or "bi-13" degree polynomials would have 49 points.
384 const Real data[9][3] =
385 {
386 {Real(0.0000000000000000e+00), Real(0.0000000000000000e+00), Real(3.0038211543122536e-01)}, // 1
387 {Real(9.8348668243987226e-01), Real(7.7880971155441942e-01), Real(2.9991838864499131e-02)}, // 4
388 {Real(8.5955600564163892e-01), Real(9.5729769978630736e-01), Real(3.8174421317083669e-02)}, // 4
389 {Real(9.5892517028753485e-01), Real(1.3818345986246535e-01), Real(6.0424923817749980e-02)}, // 4
390 {Real(3.9073621612946100e-01), Real(9.4132722587292523e-01), Real(7.7492738533105339e-02)}, // 4
391 {Real(8.5007667369974857e-01), Real(4.7580862521827590e-01), Real(1.1884466730059560e-01)}, // 4
392 {Real(6.4782163718701073e-01), Real(7.5580535657208143e-01), Real(1.2976355037000271e-01)}, // 4
393 {Real(7.0741508996444936e-02), Real(6.9625007849174941e-01), Real(2.1334158145718938e-01)}, // 4
394 {Real(4.0930456169403884e-01), Real(3.4271655604040678e-01), Real(2.5687074948196783e-01)} // 4
395 };
396
397 const unsigned int symmetry[9] = {
398 0, // Single point
399 4, // Rotational Invariant
400 4, // Rotational Invariant
401 4, // Rotational Invariant
402 4, // Rotational Invariant
403 4, // Rotational Invariant
404 4, // Rotational Invariant
405 4, // Rotational Invariant
406 4 // Rotational Invariant
407 };
408
409 _points.resize (33);
410 _weights.resize(33);
411
412 stroud_rule(data, symmetry, 9);
413
414 return;
415 } // end case TWELFTH,THIRTEENTH
416
417
418
419
420 case FOURTEENTH:
421 case FIFTEENTH:
422 {
423 // A degree-15, 48 point rule originally due to Rabinowitz and Richter,
424 // can be found in Cools' 1971 book.
425 //
426 // A.H. Stroud, Approximate calculation of multiple integrals,
427 // Prentice-Hall, Englewood Cliffs, N.J., 1971.
428 //
429 // The product Gauss rule for this order has 8^2=64 points.
430 const Real data[9][3] =
431 {
432 {Real(9.915377816777667e-01L), Real(0.0000000000000000e+00), Real(3.01245207981210e-02L)}, // 4
433 {Real(8.020163879230440e-01L), Real(0.0000000000000000e+00), Real(8.71146840209092e-02L)}, // 4
434 {Real(5.648674875232742e-01L), Real(0.0000000000000000e+00), Real(1.250080294351494e-01L)}, // 4
435 {Real(9.354392392539896e-01L), Real(0.0000000000000000e+00), Real(2.67651407861666e-02L)}, // 4
436 {Real(7.624563338825799e-01L), Real(0.0000000000000000e+00), Real(9.59651863624437e-02L)}, // 4
437 {Real(2.156164241427213e-01L), Real(0.0000000000000000e+00), Real(1.750832998343375e-01L)}, // 4
438 {Real(9.769662659711761e-01L), Real(6.684480048977932e-01L), Real(2.83136372033274e-02L)}, // 4
439 {Real(8.937128379503403e-01L), Real(3.735205277617582e-01L), Real(8.66414716025093e-02L)}, // 4
440 {Real(6.122485619312083e-01L), Real(4.078983303613935e-01L), Real(1.150144605755996e-01L)} // 4
441 };
442
443 const unsigned int symmetry[9] = {
444 3, // Full Symmetry, (x,0)
445 3, // Full Symmetry, (x,0)
446 3, // Full Symmetry, (x,0)
447 2, // Full Symmetry, (x,x)
448 2, // Full Symmetry, (x,x)
449 2, // Full Symmetry, (x,x)
450 1, // Full Symmetry, (x,y)
451 1, // Full Symmetry, (x,y)
452 1, // Full Symmetry, (x,y)
453 };
454
455 _points.resize (48);
456 _weights.resize(48);
457
458 stroud_rule(data, symmetry, 9);
459
460 return;
461 } // case FOURTEENTH, FIFTEENTH:
462
463
464
465
466 case SIXTEENTH:
467 case SEVENTEENTH:
468 {
469 // A degree 17, 60-point rule due to Cools and Haegemans.
470 //
471 // R. Cools and A. Haegemans, Another step forward in searching for
472 // cubature formulae with a minimal number of knots for the square,
473 // Computing 40 (1988), 139--146.
474 //
475 // A tensor-product rule accurate for "bi-14" or "bi-15" degree polynomials would have 64 points.
476 // A tensor-product rule accurate for "bi-16" or "bi-17" degree polynomials would have 81 points.
477 const Real data[10][3] =
478 {
479 {Real(9.8935307451260049e-01), Real(0.0000000000000000e+00), Real(2.0614915919990959e-02)}, // 4
480 {Real(3.7628520715797329e-01), Real(0.0000000000000000e+00), Real(1.2802571617990983e-01)}, // 4
481 {Real(9.7884827926223311e-01), Real(0.0000000000000000e+00), Real(5.5117395340318905e-03)}, // 4
482 {Real(8.8579472916411612e-01), Real(0.0000000000000000e+00), Real(3.9207712457141880e-02)}, // 4
483 {Real(1.7175612383834817e-01), Real(0.0000000000000000e+00), Real(7.6396945079863302e-02)}, // 4
484 {Real(5.9049927380600241e-01), Real(3.1950503663457394e-01), Real(1.4151372994997245e-01)}, // 8
485 {Real(7.9907913191686325e-01), Real(5.9797245192945738e-01), Real(8.3903279363797602e-02)}, // 8
486 {Real(8.0374396295874471e-01), Real(5.8344481776550529e-02), Real(6.0394163649684546e-02)}, // 8
487 {Real(9.3650627612749478e-01), Real(3.4738631616620267e-01), Real(5.7387752969212695e-02)}, // 8
488 {Real(9.8132117980545229e-01), Real(7.0600028779864611e-01), Real(2.1922559481863763e-02)}, // 8
489 };
490
491 const unsigned int symmetry[10] = {
492 3, // Fully symmetric (x,0)
493 3, // Fully symmetric (x,0)
494 2, // Fully symmetric (x,x)
495 2, // Fully symmetric (x,x)
496 2, // Fully symmetric (x,x)
497 1, // Fully symmetric (x,y)
498 1, // Fully symmetric (x,y)
499 1, // Fully symmetric (x,y)
500 1, // Fully symmetric (x,y)
501 1 // Fully symmetric (x,y)
502 };
503
504 _points.resize (60);
505 _weights.resize(60);
506
507 stroud_rule(data, symmetry, 10);
508
509 return;
510 } // end case FOURTEENTH through SEVENTEENTH
511 #endif
512
513
514
515 // By default: construct and use a Gauss quadrature rule
516 default:
517 {
518 // Break out and fall down into the default: case for the
519 // outer switch statement.
520 break;
521 }
522
523 } // end switch(_order + 2*p)
524 } // end case QUAD4/8/9
525
526 libmesh_fallthrough();
527
528 // By default: construct and use a Gauss quadrature rule
529 default:
530 {
531 QGauss gauss_rule(2, _order);
532 gauss_rule.init(_type, _p_level);
533
534 // Swap points and weights with the about-to-be destroyed rule.
535 _points.swap (gauss_rule.get_points() );
536 _weights.swap(gauss_rule.get_weights());
537
538 return;
539 }
540 } // end switch (_type)
541 }
542
543 } // namespace libMesh
544