quadrature_gauss_3D.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 #include "libmesh/quadrature_gm.h"
24 
25 namespace libMesh
26 {
27 
28 void QGauss::init_3D(const ElemType type_in,
29  unsigned int p)
30 {
31 #if LIBMESH_DIM == 3
32 
33  //-----------------------------------------------------------------------
34  // 3D quadrature rules
35  switch (type_in)
36  {
37  //---------------------------------------------
38  // Hex quadrature rules
39  case HEX8:
40  case HEX20:
41  case HEX27:
42  {
43  // We compute the 3D quadrature rule as a tensor
44  // product of the 1D quadrature rule.
45  QGauss q1D(1,_order);
46  q1D.init(EDGE2,p);
47 
48  tensor_product_hex( q1D );
49 
50  return;
51  }
52 
53 
54 
55  //---------------------------------------------
56  // Tetrahedral quadrature rules
57  case TET4:
58  case TET10:
59  {
60  switch(_order + 2*p)
61  {
62  // Taken from pg. 222 of "The finite element method," vol. 1
63  // ed. 5 by Zienkiewicz & Taylor
64  case CONSTANT:
65  case FIRST:
66  {
67  // Exact for linears
68  _points.resize(1);
69  _weights.resize(1);
70 
71 
72  _points[0](0) = .25;
73  _points[0](1) = .25;
74  _points[0](2) = .25;
75 
76  _weights[0] = Real(1)/6;
77 
78  return;
79  }
80  case SECOND:
81  {
82  // Exact for quadratics
83  _points.resize(4);
84  _weights.resize(4);
85 
86 
87  const Real a = .585410196624969;
88  const Real b = .138196601125011;
89 
90  _points[0](0) = a;
91  _points[0](1) = b;
92  _points[0](2) = b;
93 
94  _points[1](0) = b;
95  _points[1](1) = a;
96  _points[1](2) = b;
97 
98  _points[2](0) = b;
99  _points[2](1) = b;
100  _points[2](2) = a;
101 
102  _points[3](0) = b;
103  _points[3](1) = b;
104  _points[3](2) = b;
105 
106 
107 
108  _weights[0] = Real(1)/24;
109  _weights[1] = _weights[0];
110  _weights[2] = _weights[0];
111  _weights[3] = _weights[0];
112 
113  return;
114  }
115 
116 
117 
118  // Can be found in the class notes
119  // http://www.cs.rpi.edu/~flaherje/FEM/fem6.ps
120  // by Flaherty.
121  //
122  // Caution: this rule has a negative weight and may be
123  // unsuitable for some problems.
124  // Exact for cubics.
125  //
126  // Note: Keast (see ref. elsewhere in this file) also gives
127  // a third-order rule with positive weights, but it contains points
128  // on the ref. elt. boundary, making it less suitable for FEM calculations.
129  case THIRD:
130  {
132  {
133  _points.resize(5);
134  _weights.resize(5);
135 
136 
137  _points[0](0) = .25;
138  _points[0](1) = .25;
139  _points[0](2) = .25;
140 
141  _points[1](0) = .5;
142  _points[1](1) = Real(1)/6;
143  _points[1](2) = Real(1)/6;
144 
145  _points[2](0) = Real(1)/6;
146  _points[2](1) = .5;
147  _points[2](2) = Real(1)/6;
148 
149  _points[3](0) = Real(1)/6;
150  _points[3](1) = Real(1)/6;
151  _points[3](2) = .5;
152 
153  _points[4](0) = Real(1)/6;
154  _points[4](1) = Real(1)/6;
155  _points[4](2) = Real(1)/6;
156 
157 
158  _weights[0] = Real(-2)/15;
159  _weights[1] = .075;
160  _weights[2] = _weights[1];
161  _weights[3] = _weights[1];
162  _weights[4] = _weights[1];
163 
164  return;
165  } // end if (allow_rules_with_negative_weights)
166  else
167  {
168  // If a rule with positive weights is required, the 2x2x2 conical
169  // product rule is third-order accurate and has less points than
170  // the next-available positive-weight rule at FIFTH order.
171  QConical conical_rule(3, _order);
172  conical_rule.init(type_in, p);
173 
174  // Swap points and weights with the about-to-be destroyed rule.
175  _points.swap (conical_rule.get_points() );
176  _weights.swap(conical_rule.get_weights());
177 
178  return;
179  }
180  // Note: if !allow_rules_with_negative_weights, fall through to next case.
181  }
182 
183 
184 
185  // Originally a Keast rule,
186  // Patrick Keast,
187  // Moderate Degree Tetrahedral Quadrature Formulas,
188  // Computer Methods in Applied Mechanics and Engineering,
189  // Volume 55, Number 3, May 1986, pages 339-348.
190  //
191  // Can also be found the class notes
192  // http://www.cs.rpi.edu/~flaherje/FEM/fem6.ps
193  // by Flaherty.
194  //
195  // Caution: this rule has a negative weight and may be
196  // unsuitable for some problems.
197  case FOURTH:
198  {
200  {
201  _points.resize(11);
202  _weights.resize(11);
203 
204  // The raw data for the quadrature rule.
205  const Real rule_data[3][4] = {
206  {0.250000000000000000e+00, 0., 0., -0.131555555555555556e-01}, // 1
207  {0.785714285714285714e+00, 0.714285714285714285e-01, 0., 0.762222222222222222e-02}, // 4
208  {0.399403576166799219e+00, 0., 0.100596423833200785e+00, 0.248888888888888889e-01} // 6
209  };
210 
211 
212  // Now call the keast routine to generate _points and _weights
213  keast_rule(rule_data, 3);
214 
215  return;
216  } // end if (allow_rules_with_negative_weights)
217  // Note: if !allow_rules_with_negative_weights, fall through to next case.
218  }
219 
220  libmesh_fallthrough();
221 
222 
223  // Walkington's fifth-order 14-point rule from
224  // "Quadrature on Simplices of Arbitrary Dimension"
225  //
226  // We originally had a Keast rule here, but this rule had
227  // more points than an equivalent rule by Walkington and
228  // also contained points on the boundary of the ref. elt,
229  // making it less suitable for FEM calculations.
230  case FIFTH:
231  {
232  _points.resize(14);
233  _weights.resize(14);
234 
235  // permutations of these points and suitably-modified versions of
236  // these points are the quadrature point locations
237  const Real a[3] = {0.31088591926330060980, // a1 from the paper
238  0.092735250310891226402, // a2 from the paper
239  0.045503704125649649492}; // a3 from the paper
240 
241  // weights. a[] and wt[] are the only floating-point inputs required
242  // for this rule.
243  const Real wt[3] = {0.018781320953002641800, // w1 from the paper
244  0.012248840519393658257, // w2 from the paper
245  0.0070910034628469110730}; // w3 from the paper
246 
247  // The first two sets of 4 points are formed in a similar manner
248  for (unsigned int i=0; i<2; ++i)
249  {
250  // Where we will insert values into _points and _weights
251  const unsigned int offset=4*i;
252 
253  // Stuff points and weights values into their arrays
254  const Real b = 1. - 3.*a[i];
255 
256  // Here are the permutations. Order of these is not important,
257  // all have the same weight
258  _points[offset + 0] = Point(a[i], a[i], a[i]);
259  _points[offset + 1] = Point(a[i], b, a[i]);
260  _points[offset + 2] = Point( b, a[i], a[i]);
261  _points[offset + 3] = Point(a[i], a[i], b);
262 
263  // These 4 points all have the same weights
264  for (unsigned int j=0; j<4; ++j)
265  _weights[offset + j] = wt[i];
266  } // end for
267 
268 
269  {
270  // The third set contains 6 points and is formed a little differently
271  const unsigned int offset = 8;
272  const Real b = 0.5*(1. - 2.*a[2]);
273 
274  // Here are the permutations. Order of these is not important,
275  // all have the same weight
276  _points[offset + 0] = Point(b , b, a[2]);
277  _points[offset + 1] = Point(b , a[2], a[2]);
278  _points[offset + 2] = Point(a[2], a[2], b);
279  _points[offset + 3] = Point(a[2], b, a[2]);
280  _points[offset + 4] = Point( b, a[2], b);
281  _points[offset + 5] = Point(a[2], b, b);
282 
283  // These 6 points all have the same weights
284  for (unsigned int j=0; j<6; ++j)
285  _weights[offset + j] = wt[2];
286  }
287 
288 
289  // Original rule by Keast, unsuitable because it has points on the
290  // reference element boundary!
291  // _points.resize(15);
292  // _weights.resize(15);
293 
294  // _points[0](0) = 0.25;
295  // _points[0](1) = 0.25;
296  // _points[0](2) = 0.25;
297 
298  // {
299  // const Real a = 0.;
300  // const Real b = Real(1)/3;
301 
302  // _points[1](0) = a;
303  // _points[1](1) = b;
304  // _points[1](2) = b;
305 
306  // _points[2](0) = b;
307  // _points[2](1) = a;
308  // _points[2](2) = b;
309 
310  // _points[3](0) = b;
311  // _points[3](1) = b;
312  // _points[3](2) = a;
313 
314  // _points[4](0) = b;
315  // _points[4](1) = b;
316  // _points[4](2) = b;
317  // }
318  // {
319  // const Real a = Real(8)/11;
320  // const Real b = Real(1)/11;
321 
322  // _points[5](0) = a;
323  // _points[5](1) = b;
324  // _points[5](2) = b;
325 
326  // _points[6](0) = b;
327  // _points[6](1) = a;
328  // _points[6](2) = b;
329 
330  // _points[7](0) = b;
331  // _points[7](1) = b;
332  // _points[7](2) = a;
333 
334  // _points[8](0) = b;
335  // _points[8](1) = b;
336  // _points[8](2) = b;
337  // }
338  // {
339  // const Real a = 0.066550153573664;
340  // const Real b = 0.433449846426336;
341 
342  // _points[9](0) = b;
343  // _points[9](1) = a;
344  // _points[9](2) = a;
345 
346  // _points[10](0) = a;
347  // _points[10](1) = a;
348  // _points[10](2) = b;
349 
350  // _points[11](0) = a;
351  // _points[11](1) = b;
352  // _points[11](2) = b;
353 
354  // _points[12](0) = b;
355  // _points[12](1) = a;
356  // _points[12](2) = b;
357 
358  // _points[13](0) = b;
359  // _points[13](1) = b;
360  // _points[13](2) = a;
361 
362  // _points[14](0) = a;
363  // _points[14](1) = b;
364  // _points[14](2) = a;
365  // }
366 
367  // _weights[0] = 0.030283678097089;
368  // _weights[1] = 0.006026785714286;
369  // _weights[2] = _weights[1];
370  // _weights[3] = _weights[1];
371  // _weights[4] = _weights[1];
372  // _weights[5] = 0.011645249086029;
373  // _weights[6] = _weights[5];
374  // _weights[7] = _weights[5];
375  // _weights[8] = _weights[5];
376  // _weights[9] = 0.010949141561386;
377  // _weights[10] = _weights[9];
378  // _weights[11] = _weights[9];
379  // _weights[12] = _weights[9];
380  // _weights[13] = _weights[9];
381  // _weights[14] = _weights[9];
382 
383  return;
384  }
385 
386 
387 
388 
389  // This rule is originally from Keast:
390  // Patrick Keast,
391  // Moderate Degree Tetrahedral Quadrature Formulas,
392  // Computer Methods in Applied Mechanics and Engineering,
393  // Volume 55, Number 3, May 1986, pages 339-348.
394  //
395  // It is accurate on 6th-degree polynomials and has 24 points
396  // vs. 64 for the comparable conical product rule.
397  //
398  // Values copied 24th June 2008 from:
399  // http://people.scs.fsu.edu/~burkardt/f_src/keast/keast.f90
400  case SIXTH:
401  {
402  _points.resize (24);
403  _weights.resize(24);
404 
405  // The raw data for the quadrature rule.
406  const Real rule_data[4][4] = {
407  {0.356191386222544953e+00 , 0.214602871259151684e+00 , 0., 0.00665379170969464506e+00}, // 4
408  {0.877978124396165982e+00 , 0.0406739585346113397e+00, 0., 0.00167953517588677620e+00}, // 4
409  {0.0329863295731730594e+00, 0.322337890142275646e+00 , 0., 0.00922619692394239843e+00}, // 4
410  {0.0636610018750175299e+00, 0.269672331458315867e+00 , 0.603005664791649076e+00, 0.00803571428571428248e+00} // 12
411  };
412 
413 
414  // Now call the keast routine to generate _points and _weights
415  keast_rule(rule_data, 4);
416 
417  return;
418  }
419 
420 
421  // Keast's 31 point, 7th-order rule contains points on the reference
422  // element boundary, so we've decided not to include it here.
423  //
424  // Keast's 8th-order rule has 45 points. and a negative
425  // weight, so if you've explicitly disallowed such rules
426  // you will fall through to the conical product rules
427  // below.
428  case SEVENTH:
429  case EIGHTH:
430  {
432  {
433  _points.resize (45);
434  _weights.resize(45);
435 
436  // The raw data for the quadrature rule.
437  const Real rule_data[7][4] = {
438  {0.250000000000000000e+00, 0., 0., -0.393270066412926145e-01}, // 1
439  {0.617587190300082967e+00, 0.127470936566639015e+00, 0., 0.408131605934270525e-02}, // 4
440  {0.903763508822103123e+00, 0.320788303926322960e-01, 0., 0.658086773304341943e-03}, // 4
441  {0.497770956432810185e-01, 0., 0.450222904356718978e+00, 0.438425882512284693e-02}, // 6
442  {0.183730447398549945e+00, 0., 0.316269552601450060e+00, 0.138300638425098166e-01}, // 6
443  {0.231901089397150906e+00, 0.229177878448171174e-01, 0.513280033360881072e+00, 0.424043742468372453e-02}, // 12
444  {0.379700484718286102e-01, 0.730313427807538396e+00, 0.193746475248804382e+00, 0.223873973961420164e-02} // 12
445  };
446 
447 
448  // Now call the keast routine to generate _points and _weights
449  keast_rule(rule_data, 7);
450 
451  return;
452  } // end if (allow_rules_with_negative_weights)
453  // Note: if !allow_rules_with_negative_weights, fall through to next case.
454  }
455 
456  libmesh_fallthrough();
457 
458 
459  // Fall back on Grundmann-Moller or Conical Product rules at high orders.
460  default:
461  {
462  if ((allow_rules_with_negative_weights) && (_order + 2*p < 34))
463  {
464  // The Grundmann-Moller rules are defined to arbitrary order and
465  // can have significantly fewer evaluation points than conical product
466  // rules. If you allow rules with negative weights, the GM rules
467  // will be more efficient for degree up to 33 (for degree 34 and
468  // higher, CP is more efficient!) but may be more susceptible
469  // to round-off error. Safest is to disallow rules with negative
470  // weights, but this decision should be made on a case-by-case basis.
471  QGrundmann_Moller gm_rule(3, _order);
472  gm_rule.init(type_in, p);
473 
474  // Swap points and weights with the about-to-be destroyed rule.
475  _points.swap (gm_rule.get_points() );
476  _weights.swap(gm_rule.get_weights());
477 
478  return;
479  }
480 
481  else
482  {
483  // The following quadrature rules are generated as
484  // conical products. These tend to be non-optimal
485  // (use too many points, cluster points in certain
486  // regions of the domain) but they are quite easy to
487  // automatically generate using a 1D Gauss rule on
488  // [0,1] and two 1D Jacobi-Gauss rules on [0,1].
489  QConical conical_rule(3, _order);
490  conical_rule.init(type_in, p);
491 
492  // Swap points and weights with the about-to-be destroyed rule.
493  _points.swap (conical_rule.get_points() );
494  _weights.swap(conical_rule.get_weights());
495 
496  return;
497  }
498  }
499  }
500  } // end case TET4,TET10
501 
502 
503 
504  //---------------------------------------------
505  // Prism quadrature rules
506  case PRISM6:
507  case PRISM15:
508  case PRISM18:
509  {
510  // We compute the 3D quadrature rule as a tensor
511  // product of the 1D quadrature rule and a 2D
512  // triangle quadrature rule
513 
514  QGauss q1D(1,_order);
515  QGauss q2D(2,_order);
516 
517  // Initialize
518  q1D.init(EDGE2,p);
519  q2D.init(TRI3,p);
520 
521  tensor_product_prism(q1D, q2D);
522 
523  return;
524  }
525 
526 
527 
528  //---------------------------------------------
529  // Pyramid
530  case PYRAMID5:
531  case PYRAMID13:
532  case PYRAMID14:
533  {
534  // We compute the Pyramid rule as a conical product of a
535  // Jacobi rule with alpha==2 on the interval [0,1] two 1D
536  // Gauss rules with suitably modified points. The idea comes
537  // from: Stroud, A.H. "Approximate Calculation of Multiple
538  // Integrals."
539  QConical conical_rule(3, _order);
540  conical_rule.init(type_in, p);
541 
542  // Swap points and weights with the about-to-be destroyed rule.
543  _points.swap (conical_rule.get_points() );
544  _weights.swap(conical_rule.get_weights());
545 
546  return;
547 
548  }
549 
550 
551 
552  //---------------------------------------------
553  // Unsupported type
554  default:
555  libmesh_error_msg("ERROR: Unsupported type: " << type_in);
556  }
557 #endif
558 }
559 
560 } // namespace libMesh
virtual void init(const ElemType type=INVALID_ELEM, unsigned int p_level=0)
Definition: quadrature.C:28
bool allow_rules_with_negative_weights
Definition: quadrature.h:246
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 tensor_product_prism(const QBase &q1D, const QBase &q2D)
Definition: quadrature.C:181
virtual void init_3D(const ElemType _type=INVALID_ELEM, unsigned int p_level=0) override
void tensor_product_hex(const QBase &q1D)
Definition: quadrature.C:154
Conical product quadrature rules for Tri and Tet elements.
Implements the quadrature rules of Grundmann and Moller in 2D and 3D.
Definition: quadrature_gm.h:96
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
void keast_rule(const Real rule_data[][4], const unsigned int n_pts)
const std::vector< Point > & get_points() const
Definition: quadrature.h:142
Implements 1, 2, and 3D "Gaussian" quadrature rules.
A geometric point in (x,y,z) space.
Definition: point.h:38