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/enum_to_string.h"
24 :
25 : #include "libmesh/face_c0polygon.h"
26 :
27 : namespace libMesh
28 : {
29 :
30 :
31 677193 : void QGauss::init_2D()
32 : {
33 : #if LIBMESH_DIM > 1
34 :
35 : //-----------------------------------------------------------------------
36 : // 2D quadrature rules
37 677193 : switch (_type)
38 : {
39 :
40 :
41 : //---------------------------------------------
42 : // Quadrilateral quadrature rules
43 488975 : case QUAD4:
44 : case QUADSHELL4:
45 : case QUAD8:
46 : case QUADSHELL8:
47 : case QUAD9:
48 : case QUADSHELL9:
49 : {
50 : // We compute the 2D quadrature rule as a tensor
51 : // product of the 1D quadrature rule.
52 : //
53 : // For QUADs, a quadrature rule of order 'p' must be able to integrate
54 : // bilinear (p=1), biquadratic (p=2), bicubic (p=3), etc. polynomials of the form
55 : //
56 : // (x^p + x^{p-1} + ... + 1) * (y^p + y^{p-1} + ... + 1)
57 : //
58 : // These polynomials have terms *up to* degree 2p but they are *not* complete
59 : // polynomials of degree 2p. For example, when p=2 we have
60 : // 1
61 : // x y
62 : // x^2 xy y^2
63 : // yx^2 xy^2
64 : // x^2y^2
65 558622 : QGauss q1D(1,get_order());
66 488975 : tensor_product_quad( q1D );
67 34925 : return;
68 : }
69 :
70 :
71 : //---------------------------------------------
72 : // Triangle quadrature rules
73 187418 : case TRI3:
74 : case TRISHELL3:
75 : case TRI3SUBDIVISION:
76 : case TRI6:
77 : case TRI7:
78 : {
79 200810 : switch(get_order())
80 : {
81 10346 : case CONSTANT:
82 : case FIRST:
83 : {
84 : // Exact for linears
85 10346 : _points.resize(1);
86 10346 : _weights.resize(1);
87 :
88 10346 : _points[0](0) = Real(1)/3;
89 10346 : _points[0](1) = Real(1)/3;
90 :
91 10346 : _weights[0] = 0.5;
92 :
93 10346 : return;
94 : }
95 7233 : case SECOND:
96 : {
97 : // Exact for quadratics
98 7233 : _points.resize(3);
99 7233 : _weights.resize(3);
100 :
101 : // Alternate rule with points on ref. elt. boundaries.
102 : // Not ideal for problems with material coefficient discontinuities
103 : // aligned along element boundaries.
104 : // _points[0](0) = .5;
105 : // _points[0](1) = .5;
106 : // _points[1](0) = 0.;
107 : // _points[1](1) = .5;
108 : // _points[2](0) = .5;
109 : // _points[2](1) = .0;
110 :
111 7233 : _points[0](0) = Real(2)/3;
112 7233 : _points[0](1) = Real(1)/6;
113 :
114 7233 : _points[1](0) = Real(1)/6;
115 7233 : _points[1](1) = Real(2)/3;
116 :
117 7233 : _points[2](0) = Real(1)/6;
118 7233 : _points[2](1) = Real(1)/6;
119 :
120 :
121 7233 : _weights[0] = Real(1)/6;
122 7233 : _weights[1] = Real(1)/6;
123 7233 : _weights[2] = Real(1)/6;
124 :
125 7233 : return;
126 : }
127 85328 : case THIRD:
128 : {
129 : // Exact for cubics
130 85328 : _points.resize(4);
131 85328 : _weights.resize(4);
132 :
133 : // This rule is formed from a tensor product of
134 : // appropriately-scaled Gauss and Jacobi rules. (See
135 : // also the QConical quadrature class, this is a
136 : // hard-coded version of one of those rules.) For high
137 : // orders these rules generally have too many points,
138 : // but at extremely low order they are competitive and
139 : // have the additional benefit of having all positive
140 : // weights.
141 85328 : _points[0](0) = 1.5505102572168219018027159252941e-01_R;
142 85328 : _points[0](1) = 1.7855872826361642311703513337422e-01_R;
143 85328 : _points[1](0) = 6.4494897427831780981972840747059e-01_R;
144 85328 : _points[1](1) = 7.5031110222608118177475598324603e-02_R;
145 85328 : _points[2](0) = 1.5505102572168219018027159252941e-01_R;
146 85328 : _points[2](1) = 6.6639024601470138670269327409637e-01_R;
147 85328 : _points[3](0) = 6.4494897427831780981972840747059e-01_R;
148 85328 : _points[3](1) = 2.8001991549907407200279599420481e-01_R;
149 :
150 85328 : _weights[0] = 1.5902069087198858469718450103758e-01_R;
151 85328 : _weights[1] = 9.0979309128011415302815498962418e-02_R;
152 85328 : _weights[2] = 1.5902069087198858469718450103758e-01_R;
153 85328 : _weights[3] = 9.0979309128011415302815498962418e-02_R;
154 :
155 85328 : return;
156 :
157 :
158 : // The following third-order rule is quite commonly cited
159 : // in the literature and most likely works fine. However,
160 : // we generally prefer a rule with all positive weights
161 : // and an equal number of points, when available.
162 : //
163 : // (allow_rules_with_negative_weights)
164 : // {
165 : // // Exact for cubics
166 : // _points.resize(4);
167 : // _weights.resize(4);
168 : //
169 : // _points[0](0) = .33333333333333333333333333333333;
170 : // _points[0](1) = .33333333333333333333333333333333;
171 : //
172 : // _points[1](0) = .2;
173 : // _points[1](1) = .6;
174 : //
175 : // _points[2](0) = .2;
176 : // _points[2](1) = .2;
177 : //
178 : // _points[3](0) = .6;
179 : // _points[3](1) = .2;
180 : //
181 : //
182 : // _weights[0] = -27./96.;
183 : // _weights[1] = 25./96.;
184 : // _weights[2] = 25./96.;
185 : // _weights[3] = 25./96.;
186 : //
187 : // return;
188 : // } // end if (allow_rules_with_negative_weights)
189 : // Note: if !allow_rules_with_negative_weights, fall through to next case.
190 : }
191 :
192 :
193 :
194 : // A degree 4 rule with six points. This rule can be found in many places
195 : // including:
196 : //
197 : // J.N. Lyness and D. Jespersen, Moderate degree symmetric
198 : // quadrature rules for the triangle, J. Inst. Math. Appl. 15 (1975),
199 : // 19--32.
200 : //
201 : // We used the code in:
202 : // L. Zhang, T. Cui, and H. Liu. "A set of symmetric quadrature rules
203 : // on triangles and tetrahedra" Journal of Computational Mathematics,
204 : // v. 27, no. 1, 2009, pp. 89-96.
205 : // to generate additional precision.
206 2778 : case FOURTH:
207 : {
208 84 : const unsigned int n_wts = 2;
209 2778 : const Real wts[n_wts] =
210 : {
211 : 1.1169079483900573284750350421656140e-01_R,
212 : 5.4975871827660933819163162450105264e-02_R
213 : };
214 :
215 2778 : const Real a[n_wts] =
216 : {
217 : 4.4594849091596488631832925388305199e-01_R,
218 : 9.1576213509770743459571463402201508e-02_R
219 : };
220 :
221 2778 : const Real b[n_wts] = {0., 0.}; // not used
222 2778 : const unsigned int permutation_ids[n_wts] = {3, 3};
223 :
224 2778 : dunavant_rule2(wts, a, b, permutation_ids, n_wts); // 6 total points
225 :
226 84 : return;
227 : }
228 :
229 :
230 :
231 : // Exact for quintics
232 : // Can be found in "Quadrature on Simplices of Arbitrary
233 : // Dimension" by Walkington.
234 31514 : case FIFTH:
235 : {
236 2681 : const unsigned int n_wts = 3;
237 31514 : const Real wts[n_wts] =
238 : {
239 : Real(9)/80, // ~0.1125
240 : Real(31)/480 + std::sqrt(Real(15))/2400, // ~0.06619707639
241 : Real(31)/480 - std::sqrt(Real(15))/2400 // ~0.06296959027
242 : };
243 :
244 31514 : const Real a[n_wts] =
245 : {
246 : 0., // 'a' parameter not used for origin permutation
247 : Real(2)/7 + std::sqrt(Real(15))/21, // ~0.4701420641
248 : Real(2)/7 - std::sqrt(Real(15))/21 // ~0.1012865073
249 : };
250 :
251 31514 : const Real b[n_wts] = {0., 0., 0.}; // not used
252 31514 : const unsigned int permutation_ids[n_wts] = {1, 3, 3};
253 :
254 31514 : dunavant_rule2(wts, a, b, permutation_ids, n_wts); // 7 total points
255 :
256 2681 : return;
257 : }
258 :
259 :
260 :
261 : // A degree 6 rule with 12 points. This rule can be found in many places
262 : // including:
263 : //
264 : // J.N. Lyness and D. Jespersen, Moderate degree symmetric
265 : // quadrature rules for the triangle, J. Inst. Math. Appl. 15 (1975),
266 : // 19--32.
267 : //
268 : // We used the code in:
269 : // L. Zhang, T. Cui, and H. Liu. "A set of symmetric quadrature rules
270 : // on triangles and tetrahedra" Journal of Computational Mathematics,
271 : // v. 27, no. 1, 2009, pp. 89-96.
272 : // to generate additional precision.
273 : //
274 : // Note that the following 7th-order Ro3-invariant rule also has only 12 points,
275 : // which technically makes it the superior rule. This one is here for completeness.
276 1420 : case SIXTH:
277 : {
278 40 : const unsigned int n_wts = 3;
279 1420 : const Real wts[n_wts] =
280 : {
281 : 5.8393137863189683012644805692789721e-02_R,
282 : 2.5422453185103408460468404553434492e-02_R,
283 : 4.1425537809186787596776728210221227e-02_R
284 : };
285 :
286 1420 : const Real a[n_wts] =
287 : {
288 : 2.4928674517091042129163855310701908e-01_R,
289 : 6.3089014491502228340331602870819157e-02_R,
290 : 3.1035245103378440541660773395655215e-01_R
291 : };
292 :
293 1420 : const Real b[n_wts] =
294 : {
295 : 0.,
296 : 0.,
297 : 6.3650249912139864723014259441204970e-01_R
298 : };
299 :
300 1420 : const unsigned int permutation_ids[n_wts] = {3, 3, 6}; // 12 total points
301 :
302 1420 : dunavant_rule2(wts, a, b, permutation_ids, n_wts);
303 :
304 40 : return;
305 : }
306 :
307 :
308 : // A degree 7 rule with 12 points. This rule can be found in:
309 : //
310 : // K. Gatermann, The construction of symmetric cubature
311 : // formulas for the square and the triangle, Computing 40
312 : // (1988), 229--240.
313 : //
314 : // This rule, which is provably minimal in the number of
315 : // integration points, is said to be 'Ro3 invariant' which
316 : // means that a given set of barycentric coordinates
317 : // (z1,z2,z3) implies the quadrature points (z1,z2),
318 : // (z3,z1), (z2,z3) which are formed by taking the first
319 : // two entries in cyclic permutations of the barycentric
320 : // point. Barycentric coordinates are related in the
321 : // sense that: z3 = 1 - z1 - z2.
322 : //
323 : // The 12-point sixth-order rule for triangles given in
324 : // Flaherty's (http://www.cs.rpi.edu/~flaherje/FEM/fem6.ps)
325 : // lecture notes has been removed in favor of this rule
326 : // which is higher-order (for the same number of
327 : // quadrature points) and has a few more digits of
328 : // precision in the points and weights. Some 10-point
329 : // degree 6 rules exist for the triangle but they have
330 : // quadrature points outside the region of integration.
331 26582 : case SEVENTH:
332 : {
333 26582 : _points.resize (12);
334 26582 : _weights.resize(12);
335 :
336 1715 : const unsigned int nrows=4;
337 :
338 : // In each of the rows below, the first two entries are (z1, z2) which imply
339 : // z3. The third entry is the weight for each of the points in the cyclic permutation.
340 : // The original publication tabulated about 16 decimal digits for each point and weight
341 : // parameter. The additional digits shown here were obtained using a code in the
342 : // mp-quadrature library, https://github.com/jwpeterson/mp-quadrature
343 26582 : const Real rule_data[nrows][3] = {
344 : {6.2382265094402118173683000996350e-02_R, 6.7517867073916085442557131050869e-02_R, 2.6517028157436251428754180460739e-02_R}, // group A
345 : {5.5225456656926611737479190275645e-02_R, 3.2150249385198182266630784919920e-01_R, 4.3881408714446055036769903139288e-02_R}, // group B
346 : {3.4324302945097146469630642483938e-02_R, 6.6094919618673565761198031019780e-01_R, 2.8775042784981585738445496900219e-02_R}, // group C
347 : {5.1584233435359177925746338682643e-01_R, 2.7771616697639178256958187139372e-01_R, 6.7493187009802774462697086166421e-02_R} // group D
348 : };
349 :
350 132910 : for (unsigned int i=0, offset=0; i<nrows; ++i)
351 : {
352 106328 : _points[offset + 0] = Point(rule_data[i][0], rule_data[i][1]); // (z1,z2)
353 106328 : _points[offset + 1] = Point(1.-rule_data[i][0]-rule_data[i][1], rule_data[i][0]); // (z3,z1)
354 106328 : _points[offset + 2] = Point(rule_data[i][1], 1.-rule_data[i][0]-rule_data[i][1]); // (z2,z3)
355 :
356 : // All these points get the same weight
357 113188 : _weights[offset + 0] = rule_data[i][2];
358 106328 : _weights[offset + 1] = rule_data[i][2];
359 106328 : _weights[offset + 2] = rule_data[i][2];
360 :
361 : // Increment offset
362 106328 : offset += 3;
363 : }
364 :
365 1715 : return;
366 :
367 :
368 : // // The following is an inferior 7th-order Lyness-style rule with 15 points.
369 : // // It's here only for completeness and the Ro3-invariant rule above should
370 : // // be used instead!
371 : // const unsigned int n_wts = 3;
372 : // const Real wts[n_wts] =
373 : // {
374 : // 2.6538900895116205835977487499847719e-02_R,
375 : // 3.5426541846066783659206291623201826e-02_R,
376 : // 3.4637341039708446756138297960207647e-02_R
377 : // };
378 : //
379 : // const Real a[n_wts] =
380 : // {
381 : // 6.4930513159164863078379776030396538e-02_R,
382 : // 2.8457558424917033519741605734978046e-01_R,
383 : // 3.1355918438493150795585190219862865e-01_R
384 : // };
385 : //
386 : // const Real b[n_wts] =
387 : // {
388 : // 0.,
389 : // 1.9838447668150671917987659863332941e-01_R,
390 : // 4.3863471792372471511798695971295936e-02_R
391 : // };
392 : //
393 : // const unsigned int permutation_ids[n_wts] = {3, 6, 6}; // 15 total points
394 : //
395 : // dunavant_rule2(wts, a, b, permutation_ids, n_wts);
396 : //
397 : // return;
398 : }
399 :
400 :
401 :
402 :
403 : // Another Dunavant rule. This one has all positive weights. This rule has
404 : // 16 points while a comparable conical product rule would have 5*5=25.
405 : //
406 : // It was copied 23rd June 2008 from:
407 : // http://people.scs.fsu.edu/~burkardt/f_src/dunavant/dunavant.f90
408 : //
409 : // Additional precision obtained from the code in:
410 : // L. Zhang, T. Cui, and H. Liu. "A set of symmetric quadrature rules
411 : // on triangles and tetrahedra" Journal of Computational Mathematics,
412 : // v. 27, no. 1, 2009, pp. 89-96.
413 3383 : case EIGHTH:
414 : {
415 204 : const unsigned int n_wts = 5;
416 3383 : const Real wts[n_wts] =
417 : {
418 : 7.2157803838893584125545555244532310e-02_R,
419 : 4.7545817133642312396948052194292159e-02_R,
420 : 5.1608685267359125140895775146064515e-02_R,
421 : 1.6229248811599040155462964170890299e-02_R,
422 : 1.3615157087217497132422345036954462e-02_R
423 : };
424 :
425 3383 : const Real a[n_wts] =
426 : {
427 : 0.0, // 'a' parameter not used for origin permutation
428 : 4.5929258829272315602881551449416932e-01_R,
429 : 1.7056930775176020662229350149146450e-01_R,
430 : 5.0547228317030975458423550596598947e-02_R,
431 : 2.6311282963463811342178578628464359e-01_R,
432 : };
433 :
434 3383 : const Real b[n_wts] =
435 : {
436 : 0.,
437 : 0.,
438 : 0.,
439 : 0.,
440 : 7.2849239295540428124100037917606196e-01_R
441 : };
442 :
443 3383 : const unsigned int permutation_ids[n_wts] = {1, 3, 3, 3, 6}; // 16 total points
444 :
445 3383 : dunavant_rule2(wts, a, b, permutation_ids, n_wts);
446 :
447 204 : return;
448 : }
449 :
450 :
451 :
452 : // Another Dunavant rule. This one has all positive weights. This rule has 19
453 : // points. The comparable conical product rule would have 25.
454 : // It was copied 23rd June 2008 from:
455 : // http://people.scs.fsu.edu/~burkardt/f_src/dunavant/dunavant.f90
456 : //
457 : // Additional precision obtained from the code in:
458 : // L. Zhang, T. Cui, and H. Liu. "A set of symmetric quadrature rules
459 : // on triangles and tetrahedra" Journal of Computational Mathematics,
460 : // v. 27, no. 1, 2009, pp. 89-96.
461 14706 : case NINTH:
462 : {
463 1105 : const unsigned int n_wts = 6;
464 14706 : const Real wts[n_wts] =
465 : {
466 : 4.8567898141399416909620991253644315e-02_R,
467 : 1.5667350113569535268427415643604658e-02_R,
468 : 1.2788837829349015630839399279499912e-02_R,
469 : 3.8913770502387139658369678149701978e-02_R,
470 : 3.9823869463605126516445887132022637e-02_R,
471 : 2.1641769688644688644688644688644689e-02_R
472 : };
473 :
474 14706 : const Real a[n_wts] =
475 : {
476 : 0.0, // 'a' parameter not used for origin permutation
477 : 4.8968251919873762778370692483619280e-01_R,
478 : 4.4729513394452709865106589966276365e-02_R,
479 : 4.3708959149293663726993036443535497e-01_R,
480 : 1.8820353561903273024096128046733557e-01_R,
481 : 2.2196298916076569567510252769319107e-01_R
482 : };
483 :
484 14706 : const Real b[n_wts] =
485 : {
486 : 0.,
487 : 0.,
488 : 0.,
489 : 0.,
490 : 0.,
491 : 7.4119859878449802069007987352342383e-01_R
492 : };
493 :
494 14706 : const unsigned int permutation_ids[n_wts] = {1, 3, 3, 3, 3, 6}; // 19 total points
495 :
496 14706 : dunavant_rule2(wts, a, b, permutation_ids, n_wts);
497 :
498 1105 : return;
499 : }
500 :
501 :
502 : // Another Dunavant rule with all positive weights. This rule has 25
503 : // points. The comparable conical product rule would have 36.
504 : // It was copied 23rd June 2008 from:
505 : // http://people.scs.fsu.edu/~burkardt/f_src/dunavant/dunavant.f90
506 : //
507 : // Additional precision obtained from the code in:
508 : // L. Zhang, T. Cui, and H. Liu. "A set of symmetric quadrature rules
509 : // on triangles and tetrahedra" Journal of Computational Mathematics,
510 : // v. 27, no. 1, 2009, pp. 89-96.
511 426 : case TENTH:
512 : {
513 12 : const unsigned int n_wts = 6;
514 426 : const Real wts[n_wts] =
515 : {
516 : 4.5408995191376790047643297550014267e-02_R,
517 : 1.8362978878233352358503035945683300e-02_R,
518 : 2.2660529717763967391302822369298659e-02_R,
519 : 3.6378958422710054302157588309680344e-02_R,
520 : 1.4163621265528742418368530791049552e-02_R,
521 : 4.7108334818664117299637354834434138e-03_R
522 : };
523 :
524 426 : const Real a[n_wts] =
525 : {
526 : 0.0, // 'a' parameter not used for origin permutation
527 : 4.8557763338365737736750753220812615e-01_R,
528 : 1.0948157548503705479545863134052284e-01_R,
529 : 3.0793983876412095016515502293063162e-01_R,
530 : 2.4667256063990269391727646541117681e-01_R,
531 : 6.6803251012200265773540212762024737e-02_R
532 : };
533 :
534 426 : const Real b[n_wts] =
535 : {
536 : 0.,
537 : 0.,
538 : 0.,
539 : 5.5035294182099909507816172659300821e-01_R,
540 : 7.2832390459741092000873505358107866e-01_R,
541 : 9.2365593358750027664630697761508843e-01_R
542 : };
543 :
544 426 : const unsigned int permutation_ids[n_wts] = {1, 3, 3, 6, 6, 6}; // 25 total points
545 :
546 426 : dunavant_rule2(wts, a, b, permutation_ids, n_wts);
547 :
548 12 : return;
549 : }
550 :
551 :
552 : // Dunavant's 11th-order rule contains points outside the region of
553 : // integration, and is thus unacceptable for our FEM calculations.
554 : //
555 : // This 30-point, 11th-order rule was obtained by me [JWP] using the code in
556 : //
557 : // Additional precision obtained from the code in:
558 : // L. Zhang, T. Cui, and H. Liu. "A set of symmetric quadrature rules
559 : // on triangles and tetrahedra" Journal of Computational Mathematics,
560 : // v. 27, no. 1, 2009, pp. 89-96.
561 : //
562 : // Note: the 28-point 11th-order rule obtained by Zhang in the paper above
563 : // does not appear to be unique. It is a solution in the sense that it
564 : // minimizes the error in the least-squares minimization problem, but
565 : // it involves too many unknowns and the Jacobian is therefore singular
566 : // when attempting to improve the solution via Newton's method.
567 1714 : case ELEVENTH:
568 : {
569 138 : const unsigned int n_wts = 6;
570 1714 : const Real wts[n_wts] =
571 : {
572 : 3.6089021198604635216985338480426484e-02_R,
573 : 2.1607717807680420303346736867931050e-02_R,
574 : 3.1144524293927978774861144478241807e-03_R,
575 : 2.9086855161081509446654185084988077e-02_R,
576 : 8.4879241614917017182977532679947624e-03_R,
577 : 1.3795732078224796530729242858347546e-02_R
578 : };
579 :
580 1714 : const Real a[n_wts] =
581 : {
582 : 3.9355079629947969884346551941969960e-01_R,
583 : 4.7979065808897448654107733982929214e-01_R,
584 : 5.1003445645828061436081405648347852e-03_R,
585 : 2.6597620190330158952732822450744488e-01_R,
586 : 2.8536418538696461608233522814483715e-01_R,
587 : 1.3723536747817085036455583801851025e-01_R
588 : };
589 :
590 1714 : const Real b[n_wts] =
591 : {
592 : 0.,
593 : 0.,
594 : 5.6817155788572446538150614865768991e-02_R,
595 : 1.2539956353662088473247489775203396e-01_R,
596 : 1.2409970153698532116262152247041742e-02_R,
597 : 5.2792057988217708934207928630851643e-02_R
598 : };
599 :
600 1714 : const unsigned int permutation_ids[n_wts] = {3, 3, 6, 6, 6, 6}; // 30 total points
601 :
602 1714 : dunavant_rule2(wts, a, b, permutation_ids, n_wts);
603 :
604 138 : return;
605 : }
606 :
607 :
608 :
609 :
610 : // Another Dunavant rule with all positive weights. This rule has 33
611 : // points. The comparable conical product rule would have 36 (ELEVENTH) or 49 (TWELFTH).
612 : //
613 : // It was copied 23rd June 2008 from:
614 : // http://people.scs.fsu.edu/~burkardt/f_src/dunavant/dunavant.f90
615 : //
616 : // Additional precision obtained from the code in:
617 : // L. Zhang, T. Cui, and H. Liu. "A set of symmetric quadrature rules
618 : // on triangles and tetrahedra" Journal of Computational Mathematics,
619 : // v. 27, no. 1, 2009, pp. 89-96.
620 284 : case TWELFTH:
621 : {
622 8 : const unsigned int n_wts = 8;
623 284 : const Real wts[n_wts] =
624 : {
625 : 3.0831305257795086169332418926151771e-03_R,
626 : 3.1429112108942550177135256546441273e-02_R,
627 : 1.7398056465354471494664198647499687e-02_R,
628 : 2.1846272269019201067728631278737487e-02_R,
629 : 1.2865533220227667708895461535782215e-02_R,
630 : 1.1178386601151722855919538351159995e-02_R,
631 : 8.6581155543294461858210504055170332e-03_R,
632 : 2.0185778883190464758914349626118386e-02_R
633 : };
634 :
635 284 : const Real a[n_wts] =
636 : {
637 : 2.1317350453210370246856975515728246e-02_R,
638 : 2.7121038501211592234595134039689474e-01_R,
639 : 1.2757614554158592467389632515428357e-01_R,
640 : 4.3972439229446027297973662348436108e-01_R,
641 : 4.8821738977380488256466206525881104e-01_R,
642 : 2.8132558098993954824813069297455275e-01_R,
643 : 1.1625191590759714124135414784260182e-01_R,
644 : 2.7571326968551419397479634607976398e-01_R
645 : };
646 :
647 284 : const Real b[n_wts] =
648 : {
649 : 0.,
650 : 0.,
651 : 0.,
652 : 0.,
653 : 0.,
654 : 6.9583608678780342214163552323607254e-01_R,
655 : 8.5801403354407263059053661662617818e-01_R,
656 : 6.0894323577978780685619243776371007e-01_R
657 : };
658 :
659 284 : const unsigned int permutation_ids[n_wts] = {3, 3, 3, 3, 3, 6, 6, 6}; // 33 total points
660 :
661 284 : dunavant_rule2(wts, a, b, permutation_ids, n_wts);
662 :
663 8 : return;
664 : }
665 :
666 :
667 : // Another Dunavant rule with all positive weights. This rule has 37
668 : // points. The comparable conical product rule would have 49 points.
669 : //
670 : // It was copied 23rd June 2008 from:
671 : // http://people.scs.fsu.edu/~burkardt/f_src/dunavant/dunavant.f90
672 : //
673 : // A second rule with additional precision obtained from the code in:
674 : // L. Zhang, T. Cui, and H. Liu. "A set of symmetric quadrature rules
675 : // on triangles and tetrahedra" Journal of Computational Mathematics,
676 : // v. 27, no. 1, 2009, pp. 89-96.
677 284 : case THIRTEENTH:
678 : {
679 8 : const unsigned int n_wts = 9;
680 284 : const Real wts[n_wts] =
681 : {
682 : 3.3980018293415822140887212340442440e-02_R,
683 : 2.7800983765226664353628733005230734e-02_R,
684 : 2.9139242559599990702383541756669905e-02_R,
685 : 3.0261685517695859208964000161454122e-03_R,
686 : 1.1997200964447365386855399725479827e-02_R,
687 : 1.7320638070424185232993414255459110e-02_R,
688 : 7.4827005525828336316229285664517190e-03_R,
689 : 1.2089519905796909568722872786530380e-02_R,
690 : 4.7953405017716313612975450830554457e-03_R
691 : };
692 :
693 284 : const Real a[n_wts] =
694 : {
695 : 0., // 'a' parameter not used for origin permutation
696 : 4.2694141425980040602081253503137421e-01_R,
697 : 2.2137228629183290065481255470507908e-01_R,
698 : 2.1509681108843183869291313534052083e-02_R,
699 : 4.8907694645253934990068971909020439e-01_R,
700 : 3.0844176089211777465847185254124531e-01_R,
701 : 1.1092204280346339541286954522167452e-01_R,
702 : 1.6359740106785048023388790171095725e-01_R,
703 : 2.7251581777342966618005046435408685e-01_R
704 : };
705 :
706 284 : const Real b[n_wts] =
707 : {
708 : 0.,
709 : 0.,
710 : 0.,
711 : 0.,
712 : 0.,
713 : 6.2354599555367557081585435318623659e-01_R,
714 : 8.6470777029544277530254595089569318e-01_R,
715 : 7.4850711589995219517301859578870965e-01_R,
716 : 7.2235779312418796526062013230478405e-01_R
717 : };
718 :
719 284 : const unsigned int permutation_ids[n_wts] = {1, 3, 3, 3, 3, 6, 6, 6, 6}; // 37 total points
720 :
721 284 : dunavant_rule2(wts, a, b, permutation_ids, n_wts);
722 :
723 8 : return;
724 : }
725 :
726 :
727 : // Another Dunavant rule. This rule has 42 points, while
728 : // a comparable conical product rule would have 64.
729 : //
730 : // It was copied 23rd June 2008 from:
731 : // http://people.scs.fsu.edu/~burkardt/f_src/dunavant/dunavant.f90
732 : //
733 : // Additional precision obtained from the code in:
734 : // L. Zhang, T. Cui, and H. Liu. "A set of symmetric quadrature rules
735 : // on triangles and tetrahedra" Journal of Computational Mathematics,
736 : // v. 27, no. 1, 2009, pp. 89-96.
737 284 : case FOURTEENTH:
738 : {
739 8 : const unsigned int n_wts = 10;
740 284 : const Real wts[n_wts] =
741 : {
742 : 1.0941790684714445320422472981662986e-02_R,
743 : 1.6394176772062675320655489369312672e-02_R,
744 : 2.5887052253645793157392455083198201e-02_R,
745 : 2.1081294368496508769115218662093065e-02_R,
746 : 7.2168498348883338008549607403266583e-03_R,
747 : 2.4617018012000408409130117545210774e-03_R,
748 : 1.2332876606281836981437622591818114e-02_R,
749 : 1.9285755393530341614244513905205430e-02_R,
750 : 7.2181540567669202480443459995079017e-03_R,
751 : 2.5051144192503358849300465412445582e-03_R
752 : };
753 :
754 284 : const Real a[n_wts] =
755 : {
756 : 4.8896391036217863867737602045239024e-01_R,
757 : 4.1764471934045392250944082218564344e-01_R,
758 : 2.7347752830883865975494428326269856e-01_R,
759 : 1.7720553241254343695661069046505908e-01_R,
760 : 6.1799883090872601267478828436935788e-02_R,
761 : 1.9390961248701048178250095054529511e-02_R,
762 : 1.7226668782135557837528960161365733e-01_R,
763 : 3.3686145979634500174405519708892539e-01_R,
764 : 2.9837288213625775297083151805961273e-01_R,
765 : 1.1897449769695684539818196192990548e-01_R
766 : };
767 :
768 284 : const Real b[n_wts] =
769 : {
770 : 0.,
771 : 0.,
772 : 0.,
773 : 0.,
774 : 0.,
775 : 0.,
776 : 7.7060855477499648258903327416742796e-01_R,
777 : 5.7022229084668317349769621336235426e-01_R,
778 : 6.8698016780808783735862715402031306e-01_R,
779 : 8.7975717137017112951457163697460183e-01_R
780 : };
781 :
782 284 : const unsigned int permutation_ids[n_wts]
783 : = {3, 3, 3, 3, 3, 3, 6, 6, 6, 6}; // 42 total points
784 :
785 284 : dunavant_rule2(wts, a, b, permutation_ids, n_wts);
786 :
787 8 : return;
788 : }
789 :
790 :
791 : // This 49-point rule was found by me [JWP] using the code in:
792 : //
793 : // L. Zhang, T. Cui, and H. Liu. "A set of symmetric quadrature rules
794 : // on triangles and tetrahedra" Journal of Computational Mathematics,
795 : // v. 27, no. 1, 2009, pp. 89-96.
796 : //
797 : // A 54-point, 15th-order rule is reported by
798 : //
799 : // Stephen Wandzura, Hong Xiao,
800 : // Symmetric Quadrature Rules on a Triangle,
801 : // Computers and Mathematics with Applications,
802 : // Volume 45, Number 12, June 2003, pages 1829-1840.
803 : //
804 : // can be found here:
805 : // http://people.scs.fsu.edu/~burkardt/f_src/wandzura/wandzura.f90
806 : //
807 : // but this 49-point rule is superior.
808 284 : case FIFTEENTH:
809 : {
810 8 : const unsigned int n_wts = 11;
811 284 : const Real wts[n_wts] =
812 : {
813 : 2.4777380743035579804788826970198951e-02_R,
814 : 9.2433943023307730591540642828347660e-03_R,
815 : 2.2485768962175402793245929133296627e-03_R,
816 : 6.7052581900064143760518398833360903e-03_R,
817 : 1.9011381726930579256700190357527956e-02_R,
818 : 1.4605445387471889398286155981802858e-02_R,
819 : 1.5087322572773133722829435011138258e-02_R,
820 : 1.5630213780078803020711746273129099e-02_R,
821 : 6.1808086085778203192616856133701233e-03_R,
822 : 3.2209366452594664857296985751120513e-03_R,
823 : 5.8747373242569702667677969985668817e-03_R
824 : };
825 :
826 284 : const Real a[n_wts] =
827 : {
828 : 0.0, // 'a' parameter not used for origin
829 : 7.9031013655541635005816956762252155e-02_R,
830 : 1.8789501810770077611247984432284226e-02_R,
831 : 4.9250168823249670532514526605352905e-01_R,
832 : 4.0886316907744105975059040108092775e-01_R,
833 : 5.3877851064220142445952549348423733e-01_R,
834 : 2.0250549804829997692885033941362673e-01_R,
835 : 5.5349674918711643207148086558288110e-01_R,
836 : 7.8345022567320812359258882143250181e-01_R,
837 : 8.9514624528794883409864566727625002e-01_R,
838 : 3.2515745241110782862789881780746490e-01_R
839 : };
840 :
841 284 : const Real b[n_wts] =
842 : {
843 : 0.,
844 : 0.,
845 : 0.,
846 : 0.,
847 : 0.,
848 : 1.9412620368774630292701241080996842e-01_R,
849 : 9.8765911355712115933807754318089099e-02_R,
850 : 7.7663767064308164090246588765178087e-02_R,
851 : 2.1594628433980258573654682690950798e-02_R,
852 : 1.2563596287784997705599005477153617e-02_R,
853 : 1.5082654870922784345283124845552190e-02_R
854 : };
855 :
856 284 : const unsigned int permutation_ids[n_wts]
857 : = {1, 3, 3, 3, 3, 6, 6, 6, 6, 6, 6}; // 49 total points
858 :
859 284 : dunavant_rule2(wts, a, b, permutation_ids, n_wts);
860 :
861 8 : return;
862 : }
863 :
864 :
865 :
866 :
867 : // Dunavant's 16th-order rule contains points outside the region of
868 : // integration, and is thus unacceptable for our FEM calculations.
869 : //
870 : // This 55-point, 16th-order rule was obtained by me [JWP] using the code in
871 : //
872 : // Additional precision obtained from the code in:
873 : // L. Zhang, T. Cui, and H. Liu. "A set of symmetric quadrature rules
874 : // on triangles and tetrahedra" Journal of Computational Mathematics,
875 : // v. 27, no. 1, 2009, pp. 89-96.
876 : //
877 : // Note: the 55-point 16th-order rule obtained by Zhang in the paper above
878 : // does not appear to be unique. It is a solution in the sense that it
879 : // minimizes the error in the least-squares minimization problem, but
880 : // it involves too many unknowns and the Jacobian is therefore singular
881 : // when attempting to improve the solution via Newton's method.
882 284 : case SIXTEENTH:
883 : {
884 8 : const unsigned int n_wts = 12;
885 284 : const Real wts[n_wts] =
886 : {
887 : 2.2668082505910087151996321171534230e-02_R,
888 : 8.4043060714818596159798961899306135e-03_R,
889 : 1.0850949634049747713966288634484161e-03_R,
890 : 7.2252773375423638869298219383808751e-03_R,
891 : 1.2997715227338366024036316182572871e-02_R,
892 : 2.0054466616677715883228810959112227e-02_R,
893 : 9.7299841600417010281624372720122710e-03_R,
894 : 1.1651974438298104227427176444311766e-02_R,
895 : 9.1291185550484450744725847363097389e-03_R,
896 : 3.5568614040947150231712567900113671e-03_R,
897 : 5.8355861686234326181790822005304303e-03_R,
898 : 4.7411314396804228041879331486234396e-03_R
899 : };
900 :
901 284 : const Real a[n_wts] =
902 : {
903 : 0.0, // 'a' parameter not used for centroid weight
904 : 8.5402539407933203673769900926355911e-02_R,
905 : 1.2425572001444092841183633409631260e-02_R,
906 : 4.9174838341891594024701017768490960e-01_R,
907 : 4.5669426695387464162068900231444462e-01_R,
908 : 4.8506759880447437974189793537259677e-01_R,
909 : 2.0622099278664205707909858461264083e-01_R,
910 : 3.2374950270039093446805340265853956e-01_R,
911 : 7.3834330556606586255186213302750029e-01_R,
912 : 9.1210673061680792565673823935174611e-01_R,
913 : 6.6129919222598721544966837350891531e-01_R,
914 : 1.7807138906021476039088828811346122e-01_R
915 : };
916 :
917 284 : const Real b[n_wts] =
918 : {
919 : 0.0,
920 : 0.0,
921 : 0.0,
922 : 0.0,
923 : 0.0,
924 : 3.2315912848634384647700266402091638e-01_R,
925 : 1.5341553679414688425981898952416987e-01_R,
926 : 7.4295478991330687632977899141707872e-02_R,
927 : 7.1278762832147862035977841733532020e-02_R,
928 : 1.6623223223705792825395256602140459e-02_R,
929 : 1.4160772533794791868984026749196156e-02_R,
930 : 1.4539694958941854654807449467759690e-02_R
931 : };
932 :
933 284 : const unsigned int permutation_ids[n_wts]
934 : = {1, 3, 3, 3, 3, 6, 6, 6, 6, 6, 6, 6}; // 55 total points
935 :
936 284 : dunavant_rule2(wts, a, b, permutation_ids, n_wts);
937 :
938 8 : return;
939 : }
940 :
941 :
942 : // Dunavant's 17th-order rule has 61 points, while a
943 : // comparable conical product rule would have 81 (16th and 17th orders).
944 : //
945 : // It can be found here:
946 : // http://people.scs.fsu.edu/~burkardt/f_src/dunavant/dunavant.f90
947 : //
948 : // Zhang reports an identical rule in:
949 : // L. Zhang, T. Cui, and H. Liu. "A set of symmetric quadrature rules
950 : // on triangles and tetrahedra" Journal of Computational Mathematics,
951 : // v. 27, no. 1, 2009, pp. 89-96.
952 : //
953 : // Note: the 61-point 17th-order rule obtained by Dunavant and Zhang
954 : // does not appear to be unique. It is a solution in the sense that it
955 : // minimizes the error in the least-squares minimization problem, but
956 : // it involves too many unknowns and the Jacobian is therefore singular
957 : // when attempting to improve the solution via Newton's method.
958 : //
959 : // Therefore, we prefer the following 63-point rule which
960 : // I [JWP] found. It appears to be more accurate than the
961 : // rule reported by Dunavant and Zhang, even though it has
962 : // a few more points.
963 142 : case SEVENTEENTH:
964 : {
965 4 : const unsigned int n_wts = 12;
966 142 : const Real wts[n_wts] =
967 : {
968 : 1.7464603792572004485690588092246146e-02_R,
969 : 5.9429003555801725246549713984660076e-03_R,
970 : 1.2490753345169579649319736639588729e-02_R,
971 : 1.5386987188875607593083456905596468e-02_R,
972 : 1.1185807311917706362674684312990270e-02_R,
973 : 1.0301845740670206831327304917180007e-02_R,
974 : 1.1767783072977049696840016810370464e-02_R,
975 : 3.8045312849431209558329128678945240e-03_R,
976 : 4.5139302178876351271037137230354382e-03_R,
977 : 2.2178812517580586419412547665472893e-03_R,
978 : 5.2216271537483672304731416553063103e-03_R,
979 : 9.8381136389470256422419930926212114e-04_R
980 : };
981 :
982 142 : const Real a[n_wts] =
983 : {
984 : 2.8796825754667362165337965123570514e-01_R,
985 : 4.9216175986208465345536805750663939e-01_R,
986 : 4.6252866763171173685916780827044612e-01_R,
987 : 1.6730292951631792248498303276090273e-01_R,
988 : 1.5816335500814652972296428532213019e-01_R,
989 : 1.6352252138387564873002458959679529e-01_R,
990 : 6.2447680488959768233910286168417367e-01_R,
991 : 8.7317249935244454285263604347964179e-01_R,
992 : 3.4428164322282694677972239461699271e-01_R,
993 : 9.1584484467813674010523309855340209e-02_R,
994 : 2.0172088013378989086826623852040632e-01_R,
995 : 9.6538762758254643474731509845084691e-01_R
996 : };
997 :
998 142 : const Real b[n_wts] =
999 : {
1000 : 0.0,
1001 : 0.0,
1002 : 0.0,
1003 : 3.4429160695501713926320695771253348e-01_R,
1004 : 2.2541623431550639817203145525444726e-01_R,
1005 : 8.0670083153531811694942222940484991e-02_R,
1006 : 6.5967451375050925655738829747288190e-02_R,
1007 : 4.5677879890996762665044366994439565e-02_R,
1008 : 1.1528411723154215812386518751976084e-02_R,
1009 : 9.3057714323900610398389176844165892e-03_R,
1010 : 1.5916814107619812717966560404970160e-02_R,
1011 : 1.0734733163764032541125434215228937e-02_R
1012 : };
1013 :
1014 142 : const unsigned int permutation_ids[n_wts]
1015 : = {3, 3, 3, 6, 6, 6, 6, 6, 6, 6, 6, 6}; // 63 total points
1016 :
1017 142 : dunavant_rule2(wts, a, b, permutation_ids, n_wts);
1018 :
1019 4 : return;
1020 :
1021 : // _points.resize (61);
1022 : // _weights.resize(61);
1023 :
1024 : // // The raw data for the quadrature rule.
1025 : // const Real p[15][4] = {
1026 : // { 1./3., 0., 0., 0.033437199290803e+00 / 2.0}, // 1-perm
1027 : // {0.005658918886452e+00, 0.497170540556774e+00, 0., 0.005093415440507e+00 / 2.0}, // 3-perm
1028 : // {0.035647354750751e+00, 0.482176322624625e+00, 0., 0.014670864527638e+00 / 2.0}, // 3-perm
1029 : // {0.099520061958437e+00, 0.450239969020782e+00, 0., 0.024350878353672e+00 / 2.0}, // 3-perm
1030 : // {0.199467521245206e+00, 0.400266239377397e+00, 0., 0.031107550868969e+00 / 2.0}, // 3-perm
1031 : // {0.495717464058095e+00, 0.252141267970953e+00, 0., 0.031257111218620e+00 / 2.0}, // 3-perm
1032 : // {0.675905990683077e+00, 0.162047004658461e+00, 0., 0.024815654339665e+00 / 2.0}, // 3-perm
1033 : // {0.848248235478508e+00, 0.075875882260746e+00, 0., 0.014056073070557e+00 / 2.0}, // 3-perm
1034 : // {0.968690546064356e+00, 0.015654726967822e+00, 0., 0.003194676173779e+00 / 2.0}, // 3-perm
1035 : // {0.010186928826919e+00, 0.334319867363658e+00, 0.655493203809423e+00, 0.008119655318993e+00 / 2.0}, // 6-perm
1036 : // {0.135440871671036e+00, 0.292221537796944e+00, 0.572337590532020e+00, 0.026805742283163e+00 / 2.0}, // 6-perm
1037 : // {0.054423924290583e+00, 0.319574885423190e+00, 0.626001190286228e+00, 0.018459993210822e+00 / 2.0}, // 6-perm
1038 : // {0.012868560833637e+00, 0.190704224192292e+00, 0.796427214974071e+00, 0.008476868534328e+00 / 2.0}, // 6-perm
1039 : // {0.067165782413524e+00, 0.180483211648746e+00, 0.752351005937729e+00, 0.018292796770025e+00 / 2.0}, // 6-perm
1040 : // {0.014663182224828e+00, 0.080711313679564e+00, 0.904625504095608e+00, 0.006665632004165e+00 / 2.0} // 6-perm
1041 : // };
1042 :
1043 :
1044 : // // Now call the dunavant routine to generate _points and _weights
1045 : // dunavant_rule(p, 15);
1046 :
1047 : // return;
1048 : }
1049 :
1050 :
1051 :
1052 : // Dunavant's 18th-order rule contains points outside the region and is therefore unsuitable
1053 : // for our FEM calculations. His 19th-order rule has 73 points, compared with 100 points for
1054 : // a comparable-order conical product rule.
1055 : //
1056 : // It was copied 23rd June 2008 from:
1057 : // http://people.scs.fsu.edu/~burkardt/f_src/dunavant/dunavant.f90
1058 284 : case EIGHTTEENTH:
1059 : case NINETEENTH:
1060 : {
1061 284 : _points.resize (73);
1062 284 : _weights.resize(73);
1063 :
1064 : // The raw data for the quadrature rule.
1065 284 : const Real rule_data[17][4] = {
1066 : { 1./3., 0., 0., 0.032906331388919e+00 / 2.0}, // 1-perm
1067 : {0.020780025853987e+00, 0.489609987073006e+00, 0., 0.010330731891272e+00 / 2.0}, // 3-perm
1068 : {0.090926214604215e+00, 0.454536892697893e+00, 0., 0.022387247263016e+00 / 2.0}, // 3-perm
1069 : {0.197166638701138e+00, 0.401416680649431e+00, 0., 0.030266125869468e+00 / 2.0}, // 3-perm
1070 : {0.488896691193805e+00, 0.255551654403098e+00, 0., 0.030490967802198e+00 / 2.0}, // 3-perm
1071 : {0.645844115695741e+00, 0.177077942152130e+00, 0., 0.024159212741641e+00 / 2.0}, // 3-perm
1072 : {0.779877893544096e+00, 0.110061053227952e+00, 0., 0.016050803586801e+00 / 2.0}, // 3-perm
1073 : {0.888942751496321e+00, 0.055528624251840e+00, 0., 0.008084580261784e+00 / 2.0}, // 3-perm
1074 : {0.974756272445543e+00, 0.012621863777229e+00, 0., 0.002079362027485e+00 / 2.0}, // 3-perm
1075 : {0.003611417848412e+00, 0.395754787356943e+00, 0.600633794794645e+00, 0.003884876904981e+00 / 2.0}, // 6-perm
1076 : {0.134466754530780e+00, 0.307929983880436e+00, 0.557603261588784e+00, 0.025574160612022e+00 / 2.0}, // 6-perm
1077 : {0.014446025776115e+00, 0.264566948406520e+00, 0.720987025817365e+00, 0.008880903573338e+00 / 2.0}, // 6-perm
1078 : {0.046933578838178e+00, 0.358539352205951e+00, 0.594527068955871e+00, 0.016124546761731e+00 / 2.0}, // 6-perm
1079 : {0.002861120350567e+00, 0.157807405968595e+00, 0.839331473680839e+00, 0.002491941817491e+00 / 2.0}, // 6-perm
1080 : {0.223861424097916e+00, 0.075050596975911e+00, 0.701087978926173e+00, 0.018242840118951e+00 / 2.0}, // 6-perm
1081 : {0.034647074816760e+00, 0.142421601113383e+00, 0.822931324069857e+00, 0.010258563736199e+00 / 2.0}, // 6-perm
1082 : {0.010161119296278e+00, 0.065494628082938e+00, 0.924344252620784e+00, 0.003799928855302e+00 / 2.0} // 6-perm
1083 : };
1084 :
1085 :
1086 : // Now call the dunavant routine to generate _points and _weights
1087 284 : dunavant_rule(rule_data, 17);
1088 :
1089 8 : return;
1090 : }
1091 :
1092 :
1093 : // 20th-order rule by Wandzura.
1094 : //
1095 : // Stephen Wandzura, Hong Xiao,
1096 : // Symmetric Quadrature Rules on a Triangle,
1097 : // Computers and Mathematics with Applications,
1098 : // Volume 45, Number 12, June 2003, pages 1829-1840.
1099 : //
1100 : // Wandzura's work extends the work of Dunavant by providing degree
1101 : // 5,10,15,20,25, and 30 rules with positive weights for the triangle.
1102 : //
1103 : // Copied on 3rd July 2008 from:
1104 : // http://people.scs.fsu.edu/~burkardt/f_src/wandzura/wandzura.f90
1105 142 : case TWENTIETH:
1106 : {
1107 : // The equivalent conical product rule would have 121 points
1108 142 : _points.resize (85);
1109 142 : _weights.resize(85);
1110 :
1111 : // The raw data for the quadrature rule.
1112 142 : const Real rule_data[19][4] = {
1113 : {0.33333333333333e+00, 0.0, 0.0, 0.2761042699769952e-01 / 2.0}, // 1-perm
1114 : {0.00150064932443e+00, 0.49924967533779e+00, 0.0, 0.1779029547326740e-02 / 2.0}, // 3-perm
1115 : {0.09413975193895e+00, 0.45293012403052e+00, 0.0, 0.2011239811396117e-01 / 2.0}, // 3-perm
1116 : {0.20447212408953e+00, 0.39776393795524e+00, 0.0, 0.2681784725933157e-01 / 2.0}, // 3-perm
1117 : {0.47099959493443e+00, 0.26450020253279e+00, 0.0, 0.2452313380150201e-01 / 2.0}, // 3-perm
1118 : {0.57796207181585e+00, 0.21101896409208e+00, 0.0, 0.1639457841069539e-01 / 2.0}, // 3-perm
1119 : {0.78452878565746e+00, 0.10773560717127e+00, 0.0, 0.1479590739864960e-01 / 2.0}, // 3-perm
1120 : {0.92186182432439e+00, 0.03906908783780e+00, 0.0, 0.4579282277704251e-02 / 2.0}, // 3-perm
1121 : {0.97765124054134e+00, 0.01117437972933e+00, 0.0, 0.1651826515576217e-02 / 2.0}, // 3-perm
1122 : {0.00534961818734e+00, 0.06354966590835e+00, 0.93110071590431e+00, 0.2349170908575584e-02 / 2.0}, // 6-perm
1123 : {0.00795481706620e+00, 0.15710691894071e+00, 0.83493826399309e+00, 0.4465925754181793e-02 / 2.0}, // 6-perm
1124 : {0.01042239828126e+00, 0.39564211436437e+00, 0.59393548735436e+00, 0.6099566807907972e-02 / 2.0}, // 6-perm
1125 : {0.01096441479612e+00, 0.27316757071291e+00, 0.71586801449097e+00, 0.6891081327188203e-02 / 2.0}, // 6-perm
1126 : {0.03856671208546e+00, 0.10178538248502e+00, 0.85964790542952e+00, 0.7997475072478163e-02 / 2.0}, // 6-perm
1127 : {0.03558050781722e+00, 0.44665854917641e+00, 0.51776094300637e+00, 0.7386134285336024e-02 / 2.0}, // 6-perm
1128 : {0.04967081636276e+00, 0.19901079414950e+00, 0.75131838948773e+00, 0.1279933187864826e-01 / 2.0}, // 6-perm
1129 : {0.05851972508433e+00, 0.32426118369228e+00, 0.61721909122339e+00, 0.1725807117569655e-01 / 2.0}, // 6-perm
1130 : {0.12149778700439e+00, 0.20853136321013e+00, 0.66997084978547e+00, 0.1867294590293547e-01 / 2.0}, // 6-perm
1131 : {0.14071084494394e+00, 0.32317056653626e+00, 0.53611858851980e+00, 0.2281822405839526e-01 / 2.0} // 6-perm
1132 : };
1133 :
1134 :
1135 : // Now call the dunavant routine to generate _points and _weights
1136 142 : dunavant_rule(rule_data, 19);
1137 :
1138 4 : return;
1139 : }
1140 :
1141 :
1142 :
1143 : // 25th-order rule by Wandzura.
1144 : //
1145 : // Stephen Wandzura, Hong Xiao,
1146 : // Symmetric Quadrature Rules on a Triangle,
1147 : // Computers and Mathematics with Applications,
1148 : // Volume 45, Number 12, June 2003, pages 1829-1840.
1149 : //
1150 : // Wandzura's work extends the work of Dunavant by providing degree
1151 : // 5,10,15,20,25, and 30 rules with positive weights for the triangle.
1152 : //
1153 : // Copied on 3rd July 2008 from:
1154 : // http://people.scs.fsu.edu/~burkardt/f_src/wandzura/wandzura.f90
1155 : // case TWENTYFIRST: // fall through to 121 point conical product rule below
1156 0 : case TWENTYSECOND:
1157 : case TWENTYTHIRD:
1158 : case TWENTYFOURTH:
1159 : case TWENTYFIFTH:
1160 : {
1161 : // The equivalent conical product rule would have 169 points
1162 0 : _points.resize (126);
1163 0 : _weights.resize(126);
1164 :
1165 : // The raw data for the quadrature rule.
1166 0 : const Real rule_data[26][4] = {
1167 : {0.02794648307317e+00, 0.48602675846341e+00, 0.0, 0.8005581880020417e-02 / 2.0}, // 3-perm
1168 : {0.13117860132765e+00, 0.43441069933617e+00, 0.0, 0.1594707683239050e-01 / 2.0}, // 3-perm
1169 : {0.22022172951207e+00, 0.38988913524396e+00, 0.0, 0.1310914123079553e-01 / 2.0}, // 3-perm
1170 : {0.40311353196039e+00, 0.29844323401980e+00, 0.0, 0.1958300096563562e-01 / 2.0}, // 3-perm
1171 : {0.53191165532526e+00, 0.23404417233737e+00, 0.0, 0.1647088544153727e-01 / 2.0}, // 3-perm
1172 : {0.69706333078196e+00, 0.15146833460902e+00, 0.0, 0.8547279074092100e-02 / 2.0}, // 3-perm
1173 : {0.77453221290801e+00, 0.11273389354599e+00, 0.0, 0.8161885857226492e-02 / 2.0}, // 3-perm
1174 : {0.84456861581695e+00, 0.07771569209153e+00, 0.0, 0.6121146539983779e-02 / 2.0}, // 3-perm
1175 : {0.93021381277141e+00, 0.03489309361430e+00, 0.0, 0.2908498264936665e-02 / 2.0}, // 3-perm
1176 : {0.98548363075813e+00, 0.00725818462093e+00, 0.0, 0.6922752456619963e-03 / 2.0}, // 3-perm
1177 : {0.00129235270444e+00, 0.22721445215336e+00, 0.77149319514219e+00, 0.1248289199277397e-02 / 2.0}, // 6-perm
1178 : {0.00539970127212e+00, 0.43501055485357e+00, 0.55958974387431e+00, 0.3404752908803022e-02 / 2.0}, // 6-perm
1179 : {0.00638400303398e+00, 0.32030959927220e+00, 0.67330639769382e+00, 0.3359654326064051e-02 / 2.0}, // 6-perm
1180 : {0.00502821150199e+00, 0.09175032228001e+00, 0.90322146621800e+00, 0.1716156539496754e-02 / 2.0}, // 6-perm
1181 : {0.00682675862178e+00, 0.03801083585872e+00, 0.95516240551949e+00, 0.1480856316715606e-02 / 2.0}, // 6-perm
1182 : {0.01001619963993e+00, 0.15742521848531e+00, 0.83255858187476e+00, 0.3511312610728685e-02 / 2.0}, // 6-perm
1183 : {0.02575781317339e+00, 0.23988965977853e+00, 0.73435252704808e+00, 0.7393550149706484e-02 / 2.0}, // 6-perm
1184 : {0.03022789811992e+00, 0.36194311812606e+00, 0.60782898375402e+00, 0.7983087477376558e-02 / 2.0}, // 6-perm
1185 : {0.03050499010716e+00, 0.08355196095483e+00, 0.88594304893801e+00, 0.4355962613158041e-02 / 2.0}, // 6-perm
1186 : {0.04595654736257e+00, 0.14844322073242e+00, 0.80560023190501e+00, 0.7365056701417832e-02 / 2.0}, // 6-perm
1187 : {0.06744280054028e+00, 0.28373970872753e+00, 0.64881749073219e+00, 0.1096357284641955e-01 / 2.0}, // 6-perm
1188 : {0.07004509141591e+00, 0.40689937511879e+00, 0.52305553346530e+00, 0.1174996174354112e-01 / 2.0}, // 6-perm
1189 : {0.08391152464012e+00, 0.19411398702489e+00, 0.72197448833499e+00, 0.1001560071379857e-01 / 2.0}, // 6-perm
1190 : {0.12037553567715e+00, 0.32413434700070e+00, 0.55549011732214e+00, 0.1330964078762868e-01 / 2.0}, // 6-perm
1191 : {0.14806689915737e+00, 0.22927748355598e+00, 0.62265561728665e+00, 0.1415444650522614e-01 / 2.0}, // 6-perm
1192 : {0.19177186586733e+00, 0.32561812259598e+00, 0.48261001153669e+00, 0.1488137956116801e-01 / 2.0} // 6-perm
1193 : };
1194 :
1195 :
1196 : // Now call the dunavant routine to generate _points and _weights
1197 0 : dunavant_rule(rule_data, 26);
1198 :
1199 0 : return;
1200 : }
1201 :
1202 :
1203 :
1204 : // 30th-order rule by Wandzura.
1205 : //
1206 : // Stephen Wandzura, Hong Xiao,
1207 : // Symmetric Quadrature Rules on a Triangle,
1208 : // Computers and Mathematics with Applications,
1209 : // Volume 45, Number 12, June 2003, pages 1829-1840.
1210 : //
1211 : // Wandzura's work extends the work of Dunavant by providing degree
1212 : // 5,10,15,20,25, and 30 rules with positive weights for the triangle.
1213 : //
1214 : // Copied on 3rd July 2008 from:
1215 : // http://people.scs.fsu.edu/~burkardt/f_src/wandzura/wandzura.f90
1216 0 : case TWENTYSIXTH:
1217 : case TWENTYSEVENTH:
1218 : case TWENTYEIGHTH:
1219 : case TWENTYNINTH:
1220 : case THIRTIETH:
1221 : {
1222 : // The equivalent conical product rule would have 256 points
1223 0 : _points.resize (175);
1224 0 : _weights.resize(175);
1225 :
1226 : // The raw data for the quadrature rule.
1227 0 : const Real rule_data[36][4] = {
1228 : {0.33333333333333e+00, 0.0, 0.0, 0.1557996020289920e-01 / 2.0}, // 1-perm
1229 : {0.00733011643277e+00, 0.49633494178362e+00, 0.0, 0.3177233700534134e-02 / 2.0}, // 3-perm
1230 : {0.08299567580296e+00, 0.45850216209852e+00, 0.0, 0.1048342663573077e-01 / 2.0}, // 3-perm
1231 : {0.15098095612541e+00, 0.42450952193729e+00, 0.0, 0.1320945957774363e-01 / 2.0}, // 3-perm
1232 : {0.23590585989217e+00, 0.38204707005392e+00, 0.0, 0.1497500696627150e-01 / 2.0}, // 3-perm
1233 : {0.43802430840785e+00, 0.28098784579608e+00, 0.0, 0.1498790444338419e-01 / 2.0}, // 3-perm
1234 : {0.54530204829193e+00, 0.22734897585403e+00, 0.0, 0.1333886474102166e-01 / 2.0}, // 3-perm
1235 : {0.65088177698254e+00, 0.17455911150873e+00, 0.0, 0.1088917111390201e-01 / 2.0}, // 3-perm
1236 : {0.75348314559713e+00, 0.12325842720144e+00, 0.0, 0.8189440660893461e-02 / 2.0}, // 3-perm
1237 : {0.83983154221561e+00, 0.08008422889220e+00, 0.0, 0.5575387588607785e-02 / 2.0}, // 3-perm
1238 : {0.90445106518420e+00, 0.04777446740790e+00, 0.0, 0.3191216473411976e-02 / 2.0}, // 3-perm
1239 : {0.95655897063972e+00, 0.02172051468014e+00, 0.0, 0.1296715144327045e-02 / 2.0}, // 3-perm
1240 : {0.99047064476913e+00, 0.00476467761544e+00, 0.0, 0.2982628261349172e-03 / 2.0}, // 3-perm
1241 : {0.00092537119335e+00, 0.41529527091331e+00, 0.58377935789334e+00, 0.9989056850788964e-03 / 2.0}, // 6-perm
1242 : {0.00138592585556e+00, 0.06118990978535e+00, 0.93742416435909e+00, 0.4628508491732533e-03 / 2.0}, // 6-perm
1243 : {0.00368241545591e+00, 0.16490869013691e+00, 0.83140889440718e+00, 0.1234451336382413e-02 / 2.0}, // 6-perm
1244 : {0.00390322342416e+00, 0.02503506223200e+00, 0.97106171434384e+00, 0.5707198522432062e-03 / 2.0}, // 6-perm
1245 : {0.00323324815501e+00, 0.30606446515110e+00, 0.69070228669389e+00, 0.1126946125877624e-02 / 2.0}, // 6-perm
1246 : {0.00646743211224e+00, 0.10707328373022e+00, 0.88645928415754e+00, 0.1747866949407337e-02 / 2.0}, // 6-perm
1247 : {0.00324747549133e+00, 0.22995754934558e+00, 0.76679497516308e+00, 0.1182818815031657e-02 / 2.0}, // 6-perm
1248 : {0.00867509080675e+00, 0.33703663330578e+00, 0.65428827588746e+00, 0.1990839294675034e-02 / 2.0}, // 6-perm
1249 : {0.01559702646731e+00, 0.05625657618206e+00, 0.92814639735063e+00, 0.1900412795035980e-02 / 2.0}, // 6-perm
1250 : {0.01797672125369e+00, 0.40245137521240e+00, 0.57957190353391e+00, 0.4498365808817451e-02 / 2.0}, // 6-perm
1251 : {0.01712424535389e+00, 0.24365470201083e+00, 0.73922105263528e+00, 0.3478719460274719e-02 / 2.0}, // 6-perm
1252 : {0.02288340534658e+00, 0.16538958561453e+00, 0.81172700903888e+00, 0.4102399036723953e-02 / 2.0}, // 6-perm
1253 : {0.03273759728777e+00, 0.09930187449585e+00, 0.86796052821639e+00, 0.4021761549744162e-02 / 2.0}, // 6-perm
1254 : {0.03382101234234e+00, 0.30847833306905e+00, 0.65770065458860e+00, 0.6033164660795066e-02 / 2.0}, // 6-perm
1255 : {0.03554761446002e+00, 0.46066831859211e+00, 0.50378406694787e+00, 0.3946290302129598e-02 / 2.0}, // 6-perm
1256 : {0.05053979030687e+00, 0.21881529945393e+00, 0.73064491023920e+00, 0.6644044537680268e-02 / 2.0}, // 6-perm
1257 : {0.05701471491573e+00, 0.37920955156027e+00, 0.56377573352399e+00, 0.8254305856078458e-02 / 2.0}, // 6-perm
1258 : {0.06415280642120e+00, 0.14296081941819e+00, 0.79288637416061e+00, 0.6496056633406411e-02 / 2.0}, // 6-perm
1259 : {0.08050114828763e+00, 0.28373128210592e+00, 0.63576756960645e+00, 0.9252778144146602e-02 / 2.0}, // 6-perm
1260 : {0.10436706813453e+00, 0.19673744100444e+00, 0.69889549086103e+00, 0.9164920726294280e-02 / 2.0}, // 6-perm
1261 : {0.11384489442875e+00, 0.35588914121166e+00, 0.53026596435959e+00, 0.1156952462809767e-01 / 2.0}, // 6-perm
1262 : {0.14536348771552e+00, 0.25981868535191e+00, 0.59481782693256e+00, 0.1176111646760917e-01 / 2.0}, // 6-perm
1263 : {0.18994565282198e+00, 0.32192318123130e+00, 0.48813116594672e+00, 0.1382470218216540e-01 / 2.0} // 6-perm
1264 : };
1265 :
1266 :
1267 : // Now call the dunavant routine to generate _points and _weights
1268 0 : dunavant_rule(rule_data, 36);
1269 :
1270 0 : return;
1271 : }
1272 :
1273 :
1274 : // By default, we fall back on the conical product rules. If the user
1275 : // requests an order higher than what is currently available in the 1D
1276 : // rules, an error will be thrown from the respective 1D code.
1277 0 : default:
1278 : {
1279 : // The following quadrature rules are generated as
1280 : // conical products. These tend to be non-optimal
1281 : // (use too many points, cluster points in certain
1282 : // regions of the domain) but they are quite easy to
1283 : // automatically generate using a 1D Gauss rule on
1284 : // [0,1] and two 1D Jacobi-Gauss rules on [0,1].
1285 0 : QConical conical_rule(2, _order);
1286 0 : conical_rule.init(*this);
1287 :
1288 : // Swap points and weights with the about-to-be destroyed rule.
1289 0 : _points.swap (conical_rule.get_points() );
1290 0 : _weights.swap(conical_rule.get_weights());
1291 :
1292 0 : return;
1293 : }
1294 : }
1295 : }
1296 :
1297 :
1298 : //---------------------------------------------
1299 : // Arbitrary polygon quadrature rules
1300 800 : case C0POLYGON:
1301 : {
1302 851 : QGauss tri_rule(2, _order);
1303 800 : tri_rule.init(TRI3, _p_level, true);
1304 :
1305 51 : std::vector<Point> & tripoints = tri_rule.get_points();
1306 51 : std::vector<Real> & triweights = tri_rule.get_weights();
1307 :
1308 102 : std::size_t numtripts = tripoints.size();
1309 :
1310 : // C0Polygon requires the newer Quadrature API
1311 800 : if (!_elem)
1312 0 : libmesh_error();
1313 :
1314 51 : libmesh_assert(_elem->type() == C0POLYGON);
1315 :
1316 51 : const C0Polygon & poly = *cast_ptr<const C0Polygon *>(_elem);
1317 :
1318 800 : std::size_t numtris = poly.n_subtriangles();
1319 800 : _points.resize(numtripts*numtris);
1320 800 : _weights.resize(numtripts*numtris);
1321 3129 : for (std::size_t t = 0; t != numtris; ++t)
1322 : {
1323 2329 : auto master_points = poly.master_subtriangle(t);
1324 :
1325 : // The factor of one half from the triweights cancels out
1326 : // the factor of two here, so we don't need to do so
1327 : // ourselves.
1328 : const Real twice_master_tri_area =
1329 2329 : (- master_points[1](1) * master_points[2](0)
1330 2329 : - master_points[0](1) * master_points[1](0)
1331 2329 : + master_points[0](1) * master_points[2](0)
1332 2329 : + master_points[0](0) * master_points[1](1)
1333 2329 : - master_points[0](0) * master_points[2](1)
1334 2329 : + master_points[1](0) * master_points[2](1));
1335 :
1336 151 : const Point v01 = master_points[1] - master_points[0];
1337 151 : const Point v02 = master_points[2] - master_points[0];
1338 :
1339 13049 : for (std::size_t i = 0; i != numtripts; ++i)
1340 : {
1341 10720 : _points[numtripts*t+i] =
1342 721 : master_points[0] +
1343 2163 : v01 * tripoints[i](0) +
1344 2884 : v02 * tripoints[i](1);
1345 12162 : _weights[numtripts*t+i] = triweights[i] *
1346 : twice_master_tri_area;
1347 : }
1348 : }
1349 51 : return;
1350 : }
1351 :
1352 : //---------------------------------------------
1353 : // Unsupported type
1354 0 : default:
1355 0 : libmesh_error_msg("Element type not supported:" << Utility::enum_to_string(_type));
1356 : }
1357 : #endif
1358 : }
1359 :
1360 : } // namespace libMesh
|
