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