quadrature_monomial_3D.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_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, static_cast<Real>(4.0L/3.0L)}
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.L + 6802.L*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 = 32.0L / 19.0L;
106  // Stroud: 0.21052632 * 8.0 = 1.684210560;
107 
108  const Real B = 1.L / ( 260072.L/133225.L - 1520*sqrt19/133225.L + (133.L - 37.L*sqrt19)*tp/133225.L );
109  // 5.4498735127757671684690782180890e-01L; // Stroud: 0.068123420 * 8.0 = 0.544987360;
110 
111  const Real C = 1.L / ( 260072.L/133225.L - 1520*sqrt19/133225.L - (133.L - 37.L*sqrt19)*tp/133225.L );
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  case SIXTH:
180  case SEVENTH:
181  {
183  {
184  // A degree 7, 31-point, "rotationally-symmetric" rule by
185  // Kim and Song, Comm. Korean Math. Soc vol. 13, no. 4, 1998, pp. 913-931.
186  // This rule contains a negative weight, so only use it if such type of
187  // rules are allowed.
188  const Real data[3][4] =
189  {
190  {0.00000000000000000000000000000000e+00L, 0.00000000000000000000000000000000e+00L, 0.00000000000000000000000000000000e+00L, -1.27536231884057971014492753623188e+00L},
191  {5.85540043769119907612630781744060e-01L, 0.00000000000000000000000000000000e+00L, 0.00000000000000000000000000000000e+00L, 8.71111111111111111111111111111111e-01L},
192  {6.94470135991704766602025803883310e-01L, 9.37161638568208038511047377665396e-01L, 4.15659267604065126239606672567031e-01L, 1.68695652173913043478260869565217e-01L}
193  };
194 
195  const unsigned int rule_id[3] = {
196  0, // (0,0,0) -> 1 permutation
197  1, // (x,0,0) -> 6 permutations
198  6 // (x,y,z) -> 24 permutations
199  };
200 
201  _points.resize(31);
202  _weights.resize(31);
203 
204  kim_rule(data, rule_id, 3);
205  return;
206  } // end if (allow_rules_with_negative_weights)
207 
208 
209  // A degree 7, 34-point, "fully-symmetric" rule, first published in
210  // P.C. Hammer and A.H. Stroud, "Numerical Evaluation of Multiple Integrals II",
211  // Mathmatical Tables and Other Aids to Computation, vol 12., no 64, 1958, pp. 272-280
212  //
213  // This rule happens to fall under the same general
214  // construction as the Kim rules, so we've re-used
215  // that code here. Stroud gives 16 digits for his rule,
216  // and this is the most accurate version I've found.
217  //
218  // For comparison, a SEVENTH-order Gauss product rule
219  // (which integrates tri-7th order polynomials) would
220  // have 4^3=64 points.
221  const Real
222  r = std::sqrt(6.L/7.L),
223  s = std::sqrt((960.L - 3.L*std::sqrt(28798.L)) / 2726.L),
224  t = std::sqrt((960.L + 3.L*std::sqrt(28798.L)) / 2726.L),
225  B1 = 8624.L / 29160.L,
226  B2 = 2744.L / 29160.L,
227  B3 = 8.L*(774.L*t*t - 230.L)/(9720.L*(t*t-s*s)),
228  B4 = 8.L*(230.L - 774.L*s*s)/(9720.L*(t*t-s*s));
229 
230  const Real data[4][4] =
231  {
232  {r, 0.L, 0.L, B1},
233  {r, r, 0.L, B2},
234  {s, s, s, B3},
235  {t, t, t, B4}
236  };
237 
238  const unsigned int rule_id[4] = {
239  1, // (x,0,0) -> 6 permutations
240  2, // (x,x,0) -> 12 permutations
241  4, // (x,x,x) -> 8 permutations
242  4 // (x,x,x) -> 8 permutations
243  };
244 
245  _points.resize(34);
246  _weights.resize(34);
247 
248  kim_rule(data, rule_id, 4);
249  return;
250 
251 
252 // // A degree 7, 38-point, "rotationally-symmetric" rule by
253 // // Kim and Song, Comm. Korean Math. Soc vol. 13, no. 4, 1998, pp. 913-931.
254 // //
255 // // This rule is obviously inferior to the 34-point rule above...
256 // const Real data[3][4] =
257 // {
258 // {9.01687807821291289082811566285950e-01L, 0.00000000000000000000000000000000e+00L, 0.00000000000000000000000000000000e+00L, 2.95189738262622903181631100062774e-01L},
259 // {4.08372221499474674069588900002128e-01L, 4.08372221499474674069588900002128e-01L, 4.08372221499474674069588900002128e-01L, 4.04055417266200582425904380777126e-01L},
260 // {8.59523090201054193116477875786220e-01L, 8.59523090201054193116477875786220e-01L, 4.14735913727987720499709244748633e-01L, 1.24850759678944080062624098058597e-01L}
261 // };
262 //
263 // const unsigned int rule_id[3] = {
264 // 1, // (x,0,0) -> 6 permutations
265 // 4, // (x,x,x) -> 8 permutations
266 // 5 // (x,x,z) -> 24 permutations
267 // };
268 //
269 // _points.resize(38);
270 // _weights.resize(38);
271 //
272 // kim_rule(data, rule_id, 3);
273 // return;
274  } // end case SIXTH,SEVENTH
275 
276  case EIGHTH:
277  {
278  // A degree 8, 47-point, "rotationally-symmetric" rule by
279  // Kim and Song, Comm. Korean Math. Soc vol. 13, no. 4, 1998, pp. 913-931.
280  //
281  // A EIGHTH-order Gauss product rule (which integrates tri-8th order polynomials)
282  // would have 5^3=125 points.
283  const Real data[5][4] =
284  {
285  {0.00000000000000000000000000000000e+00L, 0.00000000000000000000000000000000e+00L, 0.00000000000000000000000000000000e+00L, 4.51903714875199690490763818699555e-01L},
286  {7.82460796435951590652813975429717e-01L, 0.00000000000000000000000000000000e+00L, 0.00000000000000000000000000000000e+00L, 2.99379177352338919703385618576171e-01L},
287  {4.88094669706366480526729301468686e-01L, 4.88094669706366480526729301468686e-01L, 4.88094669706366480526729301468686e-01L, 3.00876159371240019939698689791164e-01L},
288  {8.62218927661481188856422891110042e-01L, 8.62218927661481188856422891110042e-01L, 8.62218927661481188856422891110042e-01L, 4.94843255877038125738173175714853e-02L},
289  {2.81113909408341856058098281846420e-01L, 9.44196578292008195318687494773744e-01L, 6.97574833707236996779391729948984e-01L, 1.22872389222467338799199767122592e-01L}
290  };
291 
292  const unsigned int rule_id[5] = {
293  0, // (0,0,0) -> 1 permutation
294  1, // (x,0,0) -> 6 permutations
295  4, // (x,x,x) -> 8 permutations
296  4, // (x,x,x) -> 8 permutations
297  6 // (x,y,z) -> 24 permutations
298  };
299 
300  _points.resize(47);
301  _weights.resize(47);
302 
303  kim_rule(data, rule_id, 5);
304  return;
305  } // end case EIGHTH
306 
307 
308  // By default: construct and use a Gauss quadrature rule
309  default:
310  {
311  // Break out and fall down into the default: case for the
312  // outer switch statement.
313  break;
314  }
315 
316  } // end switch(_order + 2*p)
317  } // end case HEX8/20/27
318 
319 
320  // By default: construct and use a Gauss quadrature rule
321  default:
322  {
323  QGauss gauss_rule(3, _order);
324  gauss_rule.init(type_in, p);
325 
326  // Swap points and weights with the about-to-be destroyed rule.
327  _points.swap (gauss_rule.get_points() );
328  _weights.swap(gauss_rule.get_weights());
329 
330  return;
331  }
332  } // end switch (type_in)
333 }
334 
335 } // namespace libMesh

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

Hosted By:
SourceForge.net Logo