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