libMesh
Loading...
Searching...
No Matches
quadrature_gm_3D.C
Go to the documentation of this file.
1// The libMesh Finite Element Library.
2// Copyright (C) 2002-2026 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#include "libmesh/enum_to_string.h"
23
24namespace libMesh
25{
26
27void QGrundmann_Moller::init_3D()
28{
29 // Nearly all GM rules contain negative weights, so if you are not
30 // allowing rules with negative weights, we cannot continue!
31 libmesh_error_msg_if(!allow_rules_with_negative_weights,
32 "You requested a Grundmann-Moller rule but\n"
33 "are not allowing rules with negative weights!\n"
34 "Either select a different quadrature class or\n"
35 "set allow_rules_with_negative_weights==true.");
36
37 switch (_type)
38 {
39 case TET4:
40 case TET10:
41 case TET14:
42 {
43 switch(get_order())
44 {
45 // We hard-code the first few orders based on output from
46 // the mp-quadrature library:
47 //
48 // https://code.google.com/p/mp-quadrature
49 //
50 // The points are ordered in such a way that the amount of
51 // round-off error in the quadrature calculations is
52 // (hopefully) minimized. These orderings were obtained
53 // via a simple permutation optimization strategy designed
54 // to produce the smallest overall average error while
55 // integrating all polynomials of degree <= d.
56 case CONSTANT:
57 case FIRST:
58 {
59 _points.resize(1);
60 _weights.resize(1);
61
62 _points[0](0) = Real(1)/4;
63 _points[0](1) = Real(1)/4;
64 _points[0](2) = Real(1)/4;
65
66 _weights[0] = Real(1)/6;
67 return;
68 }
69
70 case SECOND:
71 case THIRD:
72 {
73 _points.resize(5);
74 _weights.resize(5);
75
76 _points[0](0) = Real(1)/6; _points[0](1) = Real(1)/6; _points[0](2) = Real(1)/2; _weights[0] = Real(3)/40;
77 _points[1](0) = Real(1)/6; _points[1](1) = Real(1)/6; _points[1](2) = Real(1)/6; _weights[1] = Real(3)/40;
78 _points[2](0) = Real(1)/4; _points[2](1) = Real(1)/4; _points[2](2) = Real(1)/4; _weights[2] = -Real(2)/15;
79 _points[3](0) = Real(1)/2; _points[3](1) = Real(1)/6; _points[3](2) = Real(1)/6; _weights[3] = Real(3)/40;
80 _points[4](0) = Real(1)/6; _points[4](1) = Real(1)/2; _points[4](2) = Real(1)/6; _weights[4] = Real(3)/40;
81 return;
82 }
83
84 case FOURTH:
85 case FIFTH:
86 {
87 _points.resize(15);
88 _weights.resize(15);
89
90 _points[ 0](0) = Real(1)/8; _points[ 0](1) = Real(3)/8; _points[ 0](2) = Real(1)/8; _weights[ 0] = Real(16)/315;
91 _points[ 1](0) = Real(1)/8; _points[ 1](1) = Real(1)/8; _points[ 1](2) = Real(5)/8; _weights[ 1] = Real(16)/315;
92 _points[ 2](0) = Real(3)/8; _points[ 2](1) = Real(1)/8; _points[ 2](2) = Real(1)/8; _weights[ 2] = Real(16)/315;
93 _points[ 3](0) = Real(1)/6; _points[ 3](1) = Real(1)/2; _points[ 3](2) = Real(1)/6; _weights[ 3] = -Real(27)/280;
94 _points[ 4](0) = Real(3)/8; _points[ 4](1) = Real(1)/8; _points[ 4](2) = Real(3)/8; _weights[ 4] = Real(16)/315;
95 _points[ 5](0) = Real(1)/8; _points[ 5](1) = Real(3)/8; _points[ 5](2) = Real(3)/8; _weights[ 5] = Real(16)/315;
96 _points[ 6](0) = Real(1)/6; _points[ 6](1) = Real(1)/6; _points[ 6](2) = Real(1)/6; _weights[ 6] = -Real(27)/280;
97 _points[ 7](0) = Real(1)/6; _points[ 7](1) = Real(1)/6; _points[ 7](2) = Real(1)/2; _weights[ 7] = -Real(27)/280;
98 _points[ 8](0) = Real(1)/8; _points[ 8](1) = Real(1)/8; _points[ 8](2) = Real(1)/8; _weights[ 8] = Real(16)/315;
99 _points[ 9](0) = Real(1)/4; _points[ 9](1) = Real(1)/4; _points[ 9](2) = Real(1)/4; _weights[ 9] = Real(2)/45;
100 _points[10](0) = Real(1)/8; _points[10](1) = Real(5)/8; _points[10](2) = Real(1)/8; _weights[10] = Real(16)/315;
101 _points[11](0) = Real(1)/2; _points[11](1) = Real(1)/6; _points[11](2) = Real(1)/6; _weights[11] = -Real(27)/280;
102 _points[12](0) = Real(1)/8; _points[12](1) = Real(1)/8; _points[12](2) = Real(3)/8; _weights[12] = Real(16)/315;
103 _points[13](0) = Real(5)/8; _points[13](1) = Real(1)/8; _points[13](2) = Real(1)/8; _weights[13] = Real(16)/315;
104 _points[14](0) = Real(3)/8; _points[14](1) = Real(3)/8; _points[14](2) = Real(1)/8; _weights[14] = Real(16)/315;
105 return;
106 }
107
108 case SIXTH:
109 case SEVENTH:
110 {
111 _points.resize(35);
112 _weights.resize(35);
113
114 _points[ 0](0) = Real(3)/10; _points[ 0](1) = Real(1)/10; _points[ 0](2) = Real(3)/10; _weights[ 0] = Real(3125)/72576;
115 _points[ 1](0) = Real(1)/6; _points[ 1](1) = Real(1)/2; _points[ 1](2) = Real(1)/6; _weights[ 1] = Real(243)/4480;
116 _points[ 2](0) = Real(1)/6; _points[ 2](1) = Real(1)/6; _points[ 2](2) = Real(1)/2; _weights[ 2] = Real(243)/4480;
117 _points[ 3](0) = Real(1)/2; _points[ 3](1) = Real(1)/10; _points[ 3](2) = Real(1)/10; _weights[ 3] = Real(3125)/72576;
118 _points[ 4](0) = Real(3)/10; _points[ 4](1) = Real(1)/10; _points[ 4](2) = Real(1)/10; _weights[ 4] = Real(3125)/72576;
119 _points[ 5](0) = Real(3)/10; _points[ 5](1) = Real(3)/10; _points[ 5](2) = Real(1)/10; _weights[ 5] = Real(3125)/72576;
120 _points[ 6](0) = Real(1)/2; _points[ 6](1) = Real(1)/6; _points[ 6](2) = Real(1)/6; _weights[ 6] = Real(243)/4480;
121 _points[ 7](0) = Real(3)/10; _points[ 7](1) = Real(1)/10; _points[ 7](2) = Real(1)/2; _weights[ 7] = Real(3125)/72576;
122 _points[ 8](0) = Real(7)/10; _points[ 8](1) = Real(1)/10; _points[ 8](2) = Real(1)/10; _weights[ 8] = Real(3125)/72576;
123 _points[ 9](0) = Real(3)/8; _points[ 9](1) = Real(1)/8; _points[ 9](2) = Real(1)/8; _weights[ 9] = -Real(256)/2835;
124 _points[10](0) = Real(3)/10; _points[10](1) = Real(3)/10; _points[10](2) = Real(3)/10; _weights[10] = Real(3125)/72576;
125 _points[11](0) = Real(1)/10; _points[11](1) = Real(1)/2; _points[11](2) = Real(3)/10; _weights[11] = Real(3125)/72576;
126 _points[12](0) = Real(1)/10; _points[12](1) = Real(1)/10; _points[12](2) = Real(7)/10; _weights[12] = Real(3125)/72576;
127 _points[13](0) = Real(1)/2; _points[13](1) = Real(1)/10; _points[13](2) = Real(3)/10; _weights[13] = Real(3125)/72576;
128 _points[14](0) = Real(1)/10; _points[14](1) = Real(7)/10; _points[14](2) = Real(1)/10; _weights[14] = Real(3125)/72576;
129 _points[15](0) = Real(1)/10; _points[15](1) = Real(1)/2; _points[15](2) = Real(1)/10; _weights[15] = Real(3125)/72576;
130 _points[16](0) = Real(1)/6; _points[16](1) = Real(1)/6; _points[16](2) = Real(1)/6; _weights[16] = Real(243)/4480;
131 _points[17](0) = Real(3)/8; _points[17](1) = Real(1)/8; _points[17](2) = Real(3)/8; _weights[17] = -Real(256)/2835;
132 _points[18](0) = Real(1)/8; _points[18](1) = Real(1)/8; _points[18](2) = Real(5)/8; _weights[18] = -Real(256)/2835;
133 _points[19](0) = Real(1)/10; _points[19](1) = Real(1)/10; _points[19](2) = Real(3)/10; _weights[19] = Real(3125)/72576;
134 _points[20](0) = Real(1)/8; _points[20](1) = Real(3)/8; _points[20](2) = Real(3)/8; _weights[20] = -Real(256)/2835;
135 _points[21](0) = Real(5)/8; _points[21](1) = Real(1)/8; _points[21](2) = Real(1)/8; _weights[21] = -Real(256)/2835;
136 _points[22](0) = Real(1)/8; _points[22](1) = Real(5)/8; _points[22](2) = Real(1)/8; _weights[22] = -Real(256)/2835;
137 _points[23](0) = Real(1)/10; _points[23](1) = Real(3)/10; _points[23](2) = Real(1)/10; _weights[23] = Real(3125)/72576;
138 _points[24](0) = Real(1)/4; _points[24](1) = Real(1)/4; _points[24](2) = Real(1)/4; _weights[24] = -Real(8)/945;
139 _points[25](0) = Real(1)/8; _points[25](1) = Real(1)/8; _points[25](2) = Real(3)/8; _weights[25] = -Real(256)/2835;
140 _points[26](0) = Real(3)/8; _points[26](1) = Real(3)/8; _points[26](2) = Real(1)/8; _weights[26] = -Real(256)/2835;
141 _points[27](0) = Real(1)/8; _points[27](1) = Real(3)/8; _points[27](2) = Real(1)/8; _weights[27] = -Real(256)/2835;
142 _points[28](0) = Real(1)/10; _points[28](1) = Real(3)/10; _points[28](2) = Real(1)/2; _weights[28] = Real(3125)/72576;
143 _points[29](0) = Real(3)/10; _points[29](1) = Real(1)/2; _points[29](2) = Real(1)/10; _weights[29] = Real(3125)/72576;
144 _points[30](0) = Real(1)/10; _points[30](1) = Real(1)/10; _points[30](2) = Real(1)/2; _weights[30] = Real(3125)/72576;
145 _points[31](0) = Real(1)/2; _points[31](1) = Real(3)/10; _points[31](2) = Real(1)/10; _weights[31] = Real(3125)/72576;
146 _points[32](0) = Real(1)/8; _points[32](1) = Real(1)/8; _points[32](2) = Real(1)/8; _weights[32] = -Real(256)/2835;
147 _points[33](0) = Real(1)/10; _points[33](1) = Real(3)/10; _points[33](2) = Real(3)/10; _weights[33] = Real(3125)/72576;
148 _points[34](0) = Real(1)/10; _points[34](1) = Real(1)/10; _points[34](2) = Real(1)/10; _weights[34] = Real(3125)/72576;
149 return;
150 }
151
152 default:
153 {
154 // Untested above _order=23 but should work...
155 gm_rule(get_order()/2, /*dim=*/3);
156 return;
157 }
158 } // end switch (order)
159 } // end case TET
160
161 default:
162 libmesh_error_msg("ERROR: Unsupported element type: " << Utility::enum_to_string(_type));
163 } // end switch (_type)
164}
165
166} // namespace libMesh
libMesh::QBase::_points
std::vector< Point > _points
The locations of the quadrature points in reference element space.
Definition quadrature.h:409
libMesh::QBase::_type
ElemType _type
The type of element for which the current values have been computed.
Definition quadrature.h:391
libMesh::QBase::get_order
Order get_order() const
Definition quadrature.h:249
libMesh::QBase::allow_rules_with_negative_weights
bool allow_rules_with_negative_weights
Flag (default true) controlling the use of quadrature rules with negative weights.
Definition quadrature.h:301
libMesh::QBase::_weights
std::vector< Real > _weights
The quadrature weights.
Definition quadrature.h:415
libMesh::QGrundmann_Moller::init_3D
virtual void init_3D() override
Initializes the 3D quadrature rule by filling the points and weights vectors with the appropriate val...
Definition quadrature_gm_3D.C:27
libMesh::QGrundmann_Moller::gm_rule
void gm_rule(unsigned int s, unsigned int dim)
This routine is called from init_2D() and init_3D().
Definition quadrature_gm.C:57
libMesh::Utility::enum_to_string
std::string enum_to_string(const T e)
libMesh
The libMesh namespace provides an interface to certain functionality in the library.
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition libmesh_common.h:144
Generated on Mon Aug 24 2026 19:00:37 for libMesh by doxygen 1.9.8