libMesh
quadrature_monomial_3D.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 
20 // Local includes
21 #include "libmesh/quadrature_monomial.h"
22 #include "libmesh/quadrature_gauss.h"
23 
24 namespace libMesh
25 {
26 
27 
28 void QMonomial::init_3D(const ElemType type_in,
29  unsigned int p)
30 {
31 
32  switch (type_in)
33  {
34  //---------------------------------------------
35  // Hex quadrature rules
36  case HEX8:
37  case HEX20:
38  case HEX27:
39  {
40  switch(_order + 2*p)
41  {
42 
43  // The CONSTANT/FIRST rule is the 1-point Gauss "product" rule...we fall
44  // through to the default case for this rule.
45 
46  case SECOND:
47  case THIRD:
48  {
49  // A degree 3, 6-point, "rotationally-symmetric" rule by
50  // Kim and Song, Comm. Korean Math. Soc vol. 13, no. 4, 1998, pp. 913-931.
51  //
52  // Warning: this rule contains points on the boundary of the reference
53  // element, and therefore may be unsuitable for some problems. The alternative
54  // would be a 2x2x2 Gauss product rule.
55  const Real data[1][4] =
56  {
57  {1.0L, 0.0L, 0.0L, Real(4)/3}
58  };
59 
60  const unsigned int rule_id[1] = {
61  1 // (x,0,0) -> 6 permutations
62  };
63 
64  _points.resize(6);
65  _weights.resize(6);
66 
67  kim_rule(data, rule_id, 1);
68  return;
69  } // end case SECOND,THIRD
70 
71  case FOURTH:
72  case FIFTH:
73  {
74  // A degree 5, 13-point rule by Stroud,
75  // AH Stroud, "Some Fifth Degree Integration Formulas for Symmetric Regions II.",
76  // Numerische Mathematik 9, pp. 460-468 (1967).
77  //
78  // This rule is provably minimal in the number of points. The equations given for
79  // the n-cube on pg. 466 of the paper for mu/gamma and gamma are wrong, at least for
80  // the n=3 case. The analytical values given here were computed by me [JWP] in Maple.
81 
82  // Convenient intermediate values.
83  const Real sqrt19 = std::sqrt(19.L);
84  const Real tp = std::sqrt(71440 + 6802*sqrt19);
85 
86  // Point data for permutations.
87  const Real eta = 0.00000000000000000000000000000000e+00L;
88 
89  const Real lambda = std::sqrt(1919.L/3285.L - 148.L*sqrt19/3285.L + 4.L*tp/3285.L);
90  // 8.8030440669930978047737818209860e-01L;
91 
92  const Real xi = -std::sqrt(1121.L/3285.L + 74.L*sqrt19/3285.L - 2.L*tp/3285.L);
93  // -4.9584817142571115281421242364290e-01L;
94 
95  const Real mu = std::sqrt(1121.L/3285.L + 74.L*sqrt19/3285.L + 2.L*tp/3285.L);
96  // 7.9562142216409541542982482567580e-01L;
97 
98  const Real gamma = std::sqrt(1919.L/3285.L - 148.L*sqrt19/3285.L - 4.L*tp/3285.L);
99  // 2.5293711744842581347389255929324e-02L;
100 
101  // Weights: the centroid weight is given analytically. Weight B (resp C) goes
102  // with the {lambda,xi} (resp {gamma,mu}) permutation. The single-precision
103  // results reported by Stroud are given for reference.
104 
105  const Real A = Real(32)/19;
106  // Stroud: 0.21052632 * 8.0 = 1.684210560;
107 
108  const Real B = Real(1) / (Real(260072)/133225 - 1520*sqrt19/133225 + (133 - 37*sqrt19)*tp/133225);
109  // 5.4498735127757671684690782180890e-01L; // Stroud: 0.068123420 * 8.0 = 0.544987360;
110 
111  const Real C = Real(1) / (Real(260072)/133225 - 1520*sqrt19/133225 - (133 - 37*sqrt19)*tp/133225);
112  // 5.0764422766979170420572375713840e-01L; // Stroud: 0.063455527 * 8.0 = 0.507644216;
113 
114  _points.resize(13);
115  _weights.resize(13);
116 
117  unsigned int c=0;
118 
119  // Point with weight A (origin)
120  _points[c] = Point(eta, eta, eta);
121  _weights[c++] = A;
122 
123  // Points with weight B
124  _points[c] = Point(lambda, xi, xi);
125  _weights[c++] = B;
126  _points[c] = -_points[c-1];
127  _weights[c++] = B;
128 
129  _points[c] = Point(xi, lambda, xi);
130  _weights[c++] = B;
131  _points[c] = -_points[c-1];
132  _weights[c++] = B;
133 
134  _points[c] = Point(xi, xi, lambda);
135  _weights[c++] = B;
136  _points[c] = -_points[c-1];
137  _weights[c++] = B;
138 
139  // Points with weight C
140  _points[c] = Point(mu, mu, gamma);
141  _weights[c++] = C;
142  _points[c] = -_points[c-1];
143  _weights[c++] = C;
144 
145  _points[c] = Point(mu, gamma, mu);
146  _weights[c++] = C;
147  _points[c] = -_points[c-1];
148  _weights[c++] = C;
149 
150  _points[c] = Point(gamma, mu, mu);
151  _weights[c++] = C;
152  _points[c] = -_points[c-1];
153  _weights[c++] = C;
154 
155  return;
156 
157 
158  // // A degree 5, 14-point, "rotationally-symmetric" rule by
159  // // Kim and Song, Comm. Korean Math. Soc vol. 13, no. 4, 1998, pp. 913-931.
160  // // Was also reported in Stroud's 1971 book.
161  // const Real data[2][4] =
162  // {
163  // {7.95822425754221463264548820476135e-01L, 0.00000000000000000000000000000000e+00L, 0.00000000000000000000000000000000e+00L, 8.86426592797783933518005540166204e-01L},
164  // {7.58786910639328146269034278112267e-01L, 7.58786910639328146269034278112267e-01L, 7.58786910639328146269034278112267e-01L, 3.35180055401662049861495844875346e-01L}
165  // };
166 
167  // const unsigned int rule_id[2] = {
168  // 1, // (x,0,0) -> 6 permutations
169  // 4 // (x,x,x) -> 8 permutations
170  // };
171 
172  // _points.resize(14);
173  // _weights.resize(14);
174 
175  // kim_rule(data, rule_id, 2);
176  // return;
177  } // end case FOURTH,FIFTH
178 
179 
180  case SIXTH:
181  case SEVENTH:
182  {
183  if (allow_rules_with_negative_weights)
184  {
185  // A degree 7, 31-point, "rotationally-symmetric" rule by
186  // Kim and Song, Comm. Korean Math. Soc vol. 13, no. 4, 1998, pp. 913-931.
187  // This rule contains a negative weight, so only use it if such type of
188  // rules are allowed.
189  const Real data[3][4] =
190  {
191  {0.00000000000000000000000000000000e+00L, 0.00000000000000000000000000000000e+00L, 0.00000000000000000000000000000000e+00L, -1.27536231884057971014492753623188e+00L},
192  {5.85540043769119907612630781744060e-01L, 0.00000000000000000000000000000000e+00L, 0.00000000000000000000000000000000e+00L, 8.71111111111111111111111111111111e-01L},
193  {6.94470135991704766602025803883310e-01L, 9.37161638568208038511047377665396e-01L, 4.15659267604065126239606672567031e-01L, 1.68695652173913043478260869565217e-01L}
194  };
195 
196  const unsigned int rule_id[3] = {
197  0, // (0,0,0) -> 1 permutation
198  1, // (x,0,0) -> 6 permutations
199  6 // (x,y,z) -> 24 permutations
200  };
201 
202  _points.resize(31);
203  _weights.resize(31);
204 
205  kim_rule(data, rule_id, 3);
206  return;
207  } // end if (allow_rules_with_negative_weights)
208 
209 
210  // A degree 7, 34-point, "fully-symmetric" rule, first published in
211  // P.C. Hammer and A.H. Stroud, "Numerical Evaluation of Multiple Integrals II",
212  // Mathematical Tables and Other Aids to Computation, vol 12., no 64, 1958, pp. 272-280
213  //
214  // This rule happens to fall under the same general
215  // construction as the Kim rules, so we've re-used
216  // that code here. Stroud gives 16 digits for his rule,
217  // and this is the most accurate version I've found.
218  //
219  // For comparison, a SEVENTH-order Gauss product rule
220  // (which integrates tri-7th order polynomials) would
221  // have 4^3=64 points.
222  const Real
223  r = std::sqrt(6.L/7.L),
224  s = std::sqrt((960.L - 3.L*std::sqrt(28798.L)) / 2726.L),
225  t = std::sqrt((960.L + 3.L*std::sqrt(28798.L)) / 2726.L),
226  B1 = Real(8624)/29160,
227  B2 = Real(2744)/29160,
228  B3 = 8*(774*t*t - 230)/(9720*(t*t-s*s)),
229  B4 = 8*(230 - 774*s*s)/(9720*(t*t-s*s));
230 
231  const Real data[4][4] =
232  {
233  {r, 0, 0, B1},
234  {r, r, 0, B2},
235  {s, s, s, B3},
236  {t, t, t, B4}
237  };
238 
239  const unsigned int rule_id[4] = {
240  1, // (x,0,0) -> 6 permutations
241  2, // (x,x,0) -> 12 permutations
242  4, // (x,x,x) -> 8 permutations
243  4 // (x,x,x) -> 8 permutations
244  };
245 
246  _points.resize(34);
247  _weights.resize(34);
248 
249  kim_rule(data, rule_id, 4);
250  return;
251 
252 
253  // // A degree 7, 38-point, "rotationally-symmetric" rule by
254  // // Kim and Song, Comm. Korean Math. Soc vol. 13, no. 4, 1998, pp. 913-931.
255  // //
256  // // This rule is obviously inferior to the 34-point rule above...
257  // const Real data[3][4] =
258  //{
259  // {9.01687807821291289082811566285950e-01L, 0.00000000000000000000000000000000e+00L, 0.00000000000000000000000000000000e+00L, 2.95189738262622903181631100062774e-01L},
260  // {4.08372221499474674069588900002128e-01L, 4.08372221499474674069588900002128e-01L, 4.08372221499474674069588900002128e-01L, 4.04055417266200582425904380777126e-01L},
261  // {8.59523090201054193116477875786220e-01L, 8.59523090201054193116477875786220e-01L, 4.14735913727987720499709244748633e-01L, 1.24850759678944080062624098058597e-01L}
262  //};
263  //
264  // const unsigned int rule_id[3] = {
265  //1, // (x,0,0) -> 6 permutations
266  //4, // (x,x,x) -> 8 permutations
267  //5 // (x,x,z) -> 24 permutations
268  // };
269  //
270  // _points.resize(38);
271  // _weights.resize(38);
272  //
273  // kim_rule(data, rule_id, 3);
274  // return;
275  } // end case SIXTH,SEVENTH
276 
277  case EIGHTH:
278  {
279  // A degree 8, 47-point, "rotationally-symmetric" rule by
280  // Kim and Song, Comm. Korean Math. Soc vol. 13, no. 4, 1998, pp. 913-931.
281  //
282  // A EIGHTH-order Gauss product rule (which integrates tri-8th order polynomials)
283  // would have 5^3=125 points.
284  const Real data[5][4] =
285  {
286  {0.00000000000000000000000000000000e+00L, 0.00000000000000000000000000000000e+00L, 0.00000000000000000000000000000000e+00L, 4.51903714875199690490763818699555e-01L},
287  {7.82460796435951590652813975429717e-01L, 0.00000000000000000000000000000000e+00L, 0.00000000000000000000000000000000e+00L, 2.99379177352338919703385618576171e-01L},
288  {4.88094669706366480526729301468686e-01L, 4.88094669706366480526729301468686e-01L, 4.88094669706366480526729301468686e-01L, 3.00876159371240019939698689791164e-01L},
289  {8.62218927661481188856422891110042e-01L, 8.62218927661481188856422891110042e-01L, 8.62218927661481188856422891110042e-01L, 4.94843255877038125738173175714853e-02L},
290  {2.81113909408341856058098281846420e-01L, 9.44196578292008195318687494773744e-01L, 6.97574833707236996779391729948984e-01L, 1.22872389222467338799199767122592e-01L}
291  };
292 
293  const unsigned int rule_id[5] = {
294  0, // (0,0,0) -> 1 permutation
295  1, // (x,0,0) -> 6 permutations
296  4, // (x,x,x) -> 8 permutations
297  4, // (x,x,x) -> 8 permutations
298  6 // (x,y,z) -> 24 permutations
299  };
300 
301  _points.resize(47);
302  _weights.resize(47);
303 
304  kim_rule(data, rule_id, 5);
305  return;
306  } // end case EIGHTH
307 
308 
309  // By default: construct and use a Gauss quadrature rule
310  default:
311  {
312  // Break out and fall down into the default: case for the
313  // outer switch statement.
314  break;
315  }
316 
317  } // end switch(_order + 2*p)
318  } // end case HEX8/20/27
319 
320  libmesh_fallthrough();
321 
322  // By default: construct and use a Gauss quadrature rule
323  default:
324  {
325  QGauss gauss_rule(3, _order);
326  gauss_rule.init(type_in, p);
327 
328  // Swap points and weights with the about-to-be destroyed rule.
329  _points.swap (gauss_rule.get_points() );
330  _weights.swap(gauss_rule.get_weights());
331 
332  return;
333  }
334  } // end switch (type_in)
335 }
336 
337 } // namespace libMesh
ElemType
Defines an enum for geometric element types.
The libMesh namespace provides an interface to certain functionality in the library.
Definition: assembly.h:38
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
static PetscErrorCode Mat * A
IterBase * data
Ideally this private member data should have protected access.