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

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:07 UTC

Hosted By:
SourceForge.net Logo