libMesh
quadrature_gm_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_gm.h"
22 
23 namespace libMesh
24 {
25 
26 
27 
28 void QGrundmann_Moller::init_3D(const ElemType type_in,
29  unsigned int p)
30 {
31  // Nearly all GM rules contain negative weights, so if you are not
32  // allowing rules with negative weights, we cannot continue!
33  if (!allow_rules_with_negative_weights)
34  libmesh_error_msg("You requested a Grundmann-Moller rule but\n" \
35  << "are not allowing rules with negative weights!\n" \
36  << "Either select a different quadrature class or\n" \
37  << "set allow_rules_with_negative_weights==true.");
38 
39  switch (type_in)
40  {
41  case TET4:
42  case TET10:
43  {
44  switch(_order + 2*p)
45  {
46  // We hard-code the first few orders based on output from
47  // the mp-quadrature library:
48  //
49  // https://code.google.com/p/mp-quadrature
50  //
51  // The points are ordered in such a way that the amount of
52  // round-off error in the quadrature calculations is
53  // (hopefully) minimized. These orderings were obtained
54  // via a simple permutation optimization strategy designed
55  // to produce the smallest overall average error while
56  // integrating all polynomials of degree <= d.
57  case CONSTANT:
58  case FIRST:
59  {
60  _points.resize(1);
61  _weights.resize(1);
62 
63  _points[0](0) = 1.L/4.L;
64  _points[0](1) = 1.L/4.L;
65  _points[0](2) = 1.L/4.L;
66 
67  _weights[0] = 1.L/6.L;
68  return;
69  }
70 
71  case SECOND:
72  case THIRD:
73  {
74  _points.resize(5);
75  _weights.resize(5);
76 
77  _points[0](0) = 1.L/6.L; _points[0](1) = 1.L/6.L; _points[0](2) = 1.L/2.L; _weights[0] = 3.L/40.L;
78  _points[1](0) = 1.L/6.L; _points[1](1) = 1.L/6.L; _points[1](2) = 1.L/6.L; _weights[1] = 3.L/40.L;
79  _points[2](0) = 1.L/4.L; _points[2](1) = 1.L/4.L; _points[2](2) = 1.L/4.L; _weights[2] = -2.L/15.L;
80  _points[3](0) = 1.L/2.L; _points[3](1) = 1.L/6.L; _points[3](2) = 1.L/6.L; _weights[3] = 3.L/40.L;
81  _points[4](0) = 1.L/6.L; _points[4](1) = 1.L/2.L; _points[4](2) = 1.L/6.L; _weights[4] = 3.L/40.L;
82  return;
83  }
84 
85  case FOURTH:
86  case FIFTH:
87  {
88  _points.resize(15);
89  _weights.resize(15);
90 
91  _points[ 0](0) = 1.L/8.L; _points[ 0](1) = 3.L/8.L; _points[ 0](2) = 1.L/8.L; _weights[ 0] = 16.L/315.L;
92  _points[ 1](0) = 1.L/8.L; _points[ 1](1) = 1.L/8.L; _points[ 1](2) = 5.L/8.L; _weights[ 1] = 16.L/315.L;
93  _points[ 2](0) = 3.L/8.L; _points[ 2](1) = 1.L/8.L; _points[ 2](2) = 1.L/8.L; _weights[ 2] = 16.L/315.L;
94  _points[ 3](0) = 1.L/6.L; _points[ 3](1) = 1.L/2.L; _points[ 3](2) = 1.L/6.L; _weights[ 3] = -27.L/280.L;
95  _points[ 4](0) = 3.L/8.L; _points[ 4](1) = 1.L/8.L; _points[ 4](2) = 3.L/8.L; _weights[ 4] = 16.L/315.L;
96  _points[ 5](0) = 1.L/8.L; _points[ 5](1) = 3.L/8.L; _points[ 5](2) = 3.L/8.L; _weights[ 5] = 16.L/315.L;
97  _points[ 6](0) = 1.L/6.L; _points[ 6](1) = 1.L/6.L; _points[ 6](2) = 1.L/6.L; _weights[ 6] = -27.L/280.L;
98  _points[ 7](0) = 1.L/6.L; _points[ 7](1) = 1.L/6.L; _points[ 7](2) = 1.L/2.L; _weights[ 7] = -27.L/280.L;
99  _points[ 8](0) = 1.L/8.L; _points[ 8](1) = 1.L/8.L; _points[ 8](2) = 1.L/8.L; _weights[ 8] = 16.L/315.L;
100  _points[ 9](0) = 1.L/4.L; _points[ 9](1) = 1.L/4.L; _points[ 9](2) = 1.L/4.L; _weights[ 9] = 2.L/45.L;
101  _points[10](0) = 1.L/8.L; _points[10](1) = 5.L/8.L; _points[10](2) = 1.L/8.L; _weights[10] = 16.L/315.L;
102  _points[11](0) = 1.L/2.L; _points[11](1) = 1.L/6.L; _points[11](2) = 1.L/6.L; _weights[11] = -27.L/280.L;
103  _points[12](0) = 1.L/8.L; _points[12](1) = 1.L/8.L; _points[12](2) = 3.L/8.L; _weights[12] = 16.L/315.L;
104  _points[13](0) = 5.L/8.L; _points[13](1) = 1.L/8.L; _points[13](2) = 1.L/8.L; _weights[13] = 16.L/315.L;
105  _points[14](0) = 3.L/8.L; _points[14](1) = 3.L/8.L; _points[14](2) = 1.L/8.L; _weights[14] = 16.L/315.L;
106  return;
107  }
108 
109  case SIXTH:
110  case SEVENTH:
111  {
112  _points.resize(35);
113  _weights.resize(35);
114 
115  _points[ 0](0) = 3.L/10.L; _points[ 0](1) = 1.L/10.L; _points[ 0](2) = 3.L/10.L; _weights[ 0] = 3125.L/72576.L;
116  _points[ 1](0) = 1.L/6.L; _points[ 1](1) = 1.L/2.L; _points[ 1](2) = 1.L/6.L; _weights[ 1] = 243.L/4480.L;
117  _points[ 2](0) = 1.L/6.L; _points[ 2](1) = 1.L/6.L; _points[ 2](2) = 1.L/2.L; _weights[ 2] = 243.L/4480.L;
118  _points[ 3](0) = 1.L/2.L; _points[ 3](1) = 1.L/10.L; _points[ 3](2) = 1.L/10.L; _weights[ 3] = 3125.L/72576.L;
119  _points[ 4](0) = 3.L/10.L; _points[ 4](1) = 1.L/10.L; _points[ 4](2) = 1.L/10.L; _weights[ 4] = 3125.L/72576.L;
120  _points[ 5](0) = 3.L/10.L; _points[ 5](1) = 3.L/10.L; _points[ 5](2) = 1.L/10.L; _weights[ 5] = 3125.L/72576.L;
121  _points[ 6](0) = 1.L/2.L; _points[ 6](1) = 1.L/6.L; _points[ 6](2) = 1.L/6.L; _weights[ 6] = 243.L/4480.L;
122  _points[ 7](0) = 3.L/10.L; _points[ 7](1) = 1.L/10.L; _points[ 7](2) = 1.L/2.L; _weights[ 7] = 3125.L/72576.L;
123  _points[ 8](0) = 7.L/10.L; _points[ 8](1) = 1.L/10.L; _points[ 8](2) = 1.L/10.L; _weights[ 8] = 3125.L/72576.L;
124  _points[ 9](0) = 3.L/8.L; _points[ 9](1) = 1.L/8.L; _points[ 9](2) = 1.L/8.L; _weights[ 9] = -256.L/2835.L;
125  _points[10](0) = 3.L/10.L; _points[10](1) = 3.L/10.L; _points[10](2) = 3.L/10.L; _weights[10] = 3125.L/72576.L;
126  _points[11](0) = 1.L/10.L; _points[11](1) = 1.L/2.L; _points[11](2) = 3.L/10.L; _weights[11] = 3125.L/72576.L;
127  _points[12](0) = 1.L/10.L; _points[12](1) = 1.L/10.L; _points[12](2) = 7.L/10.L; _weights[12] = 3125.L/72576.L;
128  _points[13](0) = 1.L/2.L; _points[13](1) = 1.L/10.L; _points[13](2) = 3.L/10.L; _weights[13] = 3125.L/72576.L;
129  _points[14](0) = 1.L/10.L; _points[14](1) = 7.L/10.L; _points[14](2) = 1.L/10.L; _weights[14] = 3125.L/72576.L;
130  _points[15](0) = 1.L/10.L; _points[15](1) = 1.L/2.L; _points[15](2) = 1.L/10.L; _weights[15] = 3125.L/72576.L;
131  _points[16](0) = 1.L/6.L; _points[16](1) = 1.L/6.L; _points[16](2) = 1.L/6.L; _weights[16] = 243.L/4480.L;
132  _points[17](0) = 3.L/8.L; _points[17](1) = 1.L/8.L; _points[17](2) = 3.L/8.L; _weights[17] = -256.L/2835.L;
133  _points[18](0) = 1.L/8.L; _points[18](1) = 1.L/8.L; _points[18](2) = 5.L/8.L; _weights[18] = -256.L/2835.L;
134  _points[19](0) = 1.L/10.L; _points[19](1) = 1.L/10.L; _points[19](2) = 3.L/10.L; _weights[19] = 3125.L/72576.L;
135  _points[20](0) = 1.L/8.L; _points[20](1) = 3.L/8.L; _points[20](2) = 3.L/8.L; _weights[20] = -256.L/2835.L;
136  _points[21](0) = 5.L/8.L; _points[21](1) = 1.L/8.L; _points[21](2) = 1.L/8.L; _weights[21] = -256.L/2835.L;
137  _points[22](0) = 1.L/8.L; _points[22](1) = 5.L/8.L; _points[22](2) = 1.L/8.L; _weights[22] = -256.L/2835.L;
138  _points[23](0) = 1.L/10.L; _points[23](1) = 3.L/10.L; _points[23](2) = 1.L/10.L; _weights[23] = 3125.L/72576.L;
139  _points[24](0) = 1.L/4.L; _points[24](1) = 1.L/4.L; _points[24](2) = 1.L/4.L; _weights[24] = -8.L/945.L;
140  _points[25](0) = 1.L/8.L; _points[25](1) = 1.L/8.L; _points[25](2) = 3.L/8.L; _weights[25] = -256.L/2835.L;
141  _points[26](0) = 3.L/8.L; _points[26](1) = 3.L/8.L; _points[26](2) = 1.L/8.L; _weights[26] = -256.L/2835.L;
142  _points[27](0) = 1.L/8.L; _points[27](1) = 3.L/8.L; _points[27](2) = 1.L/8.L; _weights[27] = -256.L/2835.L;
143  _points[28](0) = 1.L/10.L; _points[28](1) = 3.L/10.L; _points[28](2) = 1.L/2.L; _weights[28] = 3125.L/72576.L;
144  _points[29](0) = 3.L/10.L; _points[29](1) = 1.L/2.L; _points[29](2) = 1.L/10.L; _weights[29] = 3125.L/72576.L;
145  _points[30](0) = 1.L/10.L; _points[30](1) = 1.L/10.L; _points[30](2) = 1.L/2.L; _weights[30] = 3125.L/72576.L;
146  _points[31](0) = 1.L/2.L; _points[31](1) = 3.L/10.L; _points[31](2) = 1.L/10.L; _weights[31] = 3125.L/72576.L;
147  _points[32](0) = 1.L/8.L; _points[32](1) = 1.L/8.L; _points[32](2) = 1.L/8.L; _weights[32] = -256.L/2835.L;
148  _points[33](0) = 1.L/10.L; _points[33](1) = 3.L/10.L; _points[33](2) = 3.L/10.L; _weights[33] = 3125.L/72576.L;
149  _points[34](0) = 1.L/10.L; _points[34](1) = 1.L/10.L; _points[34](2) = 1.L/10.L; _weights[34] = 3125.L/72576.L;
150  return;
151  }
152 
153  default:
154  {
155  // Untested above _order=23 but should work...
156  gm_rule((_order + 2*p)/2, /*dim=*/3);
157  return;
158  }
159  } // end switch (order)
160  } // end case TET4, TET10
161 
162  default:
163  libmesh_error_msg("ERROR: Unsupported element type: " << type_in);
164  } // end switch (type_in)
165 }
166 
167 } // namespace libMesh
ElemType
Defines an enum for geometric element types.
The libMesh namespace provides an interface to certain functionality in the library.