libMesh
Loading...
Searching...
No Matches
quadrature_monomial.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// libMesh includes
20#include "libmesh/quadrature_monomial.h"
21#include "libmesh/enum_quadrature_type.h"
22
23namespace libMesh
24{
25
26
27// See the files:
28// quadrature_monomial_2D.C
29// quadrature_monomial_3D.C
30// for implementation of specific element types.
31
32QuadratureType QMonomial::type() const
33{
34 return QMONOMIAL;
35}
36
37std::unique_ptr<QBase> QMonomial::clone() const
38{
39 return std::make_unique<QMonomial>(*this);
40}
41
42void QMonomial::wissmann_rule(const Real rule_data[][3],
43 const unsigned int n_pts)
44{
45 for (unsigned int i=0, c=0; i<n_pts; ++i)
46 {
47 _points[c] = Point( rule_data[i][0], rule_data[i][1] );
48 _weights[c++] = rule_data[i][2];
49
50 // This may be an (x1,x2) -> (-x1,x2) point, in which case
51 // we will also generate the mirror point using the same weight.
52 if (rule_data[i][0] != static_cast<Real>(0.0))
53 {
54 _points[c] = Point( -rule_data[i][0], rule_data[i][1] );
55 _weights[c++] = rule_data[i][2];
56 }
57 }
58}
59
60
61
62void QMonomial::stroud_rule(const Real rule_data[][3],
63 const unsigned int * rule_symmetry,
64 const unsigned int n_pts)
65{
66 for (unsigned int i=0, c=0; i<n_pts; ++i)
67 {
68 const Real
69 x=rule_data[i][0],
70 y=rule_data[i][1],
71 wt=rule_data[i][2];
72
73 switch(rule_symmetry[i])
74 {
75 case 0: // Single point (no symmetry)
76 {
77 _points[c] = Point( x, y);
78 _weights[c++] = wt;
79
80 break;
81 }
82 case 1: // Fully-symmetric (x,y)
83 {
84 _points[c] = Point( x, y);
85 _weights[c++] = wt;
86
87 _points[c] = Point(-x, y);
88 _weights[c++] = wt;
89
90 _points[c] = Point( x,-y);
91 _weights[c++] = wt;
92
93 _points[c] = Point(-x,-y);
94 _weights[c++] = wt;
95
96 _points[c] = Point( y, x);
97 _weights[c++] = wt;
98
99 _points[c] = Point(-y, x);
100 _weights[c++] = wt;
101
102 _points[c] = Point( y,-x);
103 _weights[c++] = wt;
104
105 _points[c] = Point(-y,-x);
106 _weights[c++] = wt;
107
108 break;
109 }
110 case 2: // Fully-symmetric (x,x)
111 {
112 _points[c] = Point( x, x);
113 _weights[c++] = wt;
114
115 _points[c] = Point(-x, x);
116 _weights[c++] = wt;
117
118 _points[c] = Point( x,-x);
119 _weights[c++] = wt;
120
121 _points[c] = Point(-x,-x);
122 _weights[c++] = wt;
123
124 break;
125 }
126 case 3: // Fully-symmetric (x,0)
127 {
128 libmesh_assert_equal_to (y, 0.0);
129
130 _points[c] = Point( x,0.);
131 _weights[c++] = wt;
132
133 _points[c] = Point(-x,0.);
134 _weights[c++] = wt;
135
136 _points[c] = Point(0., x);
137 _weights[c++] = wt;
138
139 _points[c] = Point(0.,-x);
140 _weights[c++] = wt;
141
142 break;
143 }
144 case 4: // Rotational invariant
145 {
146 _points[c] = Point( x, y);
147 _weights[c++] = wt;
148
149 _points[c] = Point(-x,-y);
150 _weights[c++] = wt;
151
152 _points[c] = Point(-y, x);
153 _weights[c++] = wt;
154
155 _points[c] = Point( y,-x);
156 _weights[c++] = wt;
157
158 break;
159 }
160 case 5: // Partial symmetry (Wissman's rules)
161 {
162 libmesh_assert_not_equal_to (x, 0.0);
163
164 _points[c] = Point( x, y);
165 _weights[c++] = wt;
166
167 _points[c] = Point(-x, y);
168 _weights[c++] = wt;
169
170 break;
171 }
172 case 6: // Rectangular symmetry
173 {
174 _points[c] = Point( x, y);
175 _weights[c++] = wt;
176
177 _points[c] = Point(-x, y);
178 _weights[c++] = wt;
179
180 _points[c] = Point(-x,-y);
181 _weights[c++] = wt;
182
183 _points[c] = Point( x,-y);
184 _weights[c++] = wt;
185
186 break;
187 }
188 case 7: // Central symmetry
189 {
190 libmesh_assert_equal_to (x, 0.0);
191 libmesh_assert_not_equal_to (y, 0.0);
192
193 _points[c] = Point(0., y);
194 _weights[c++] = wt;
195
196 _points[c] = Point(0.,-y);
197 _weights[c++] = wt;
198
199 break;
200 }
201 default:
202 libmesh_error_msg("Unknown symmetry!");
203 } // end switch(rule_symmetry[i])
204 }
205}
206
207
208
209
210void QMonomial::kim_rule(const Real rule_data[][4],
211 const unsigned int * rule_id,
212 const unsigned int n_pts)
213{
214 for (unsigned int i=0, c=0; i<n_pts; ++i)
215 {
216 const Real
217 x=rule_data[i][0],
218 y=rule_data[i][1],
219 z=rule_data[i][2],
220 wt=rule_data[i][3];
221
222 switch(rule_id[i])
223 {
224 case 0: // (0,0,0) 1 permutation
225 {
226 _points[c] = Point( x, y, z); _weights[c++] = wt;
227
228 break;
229 }
230 case 1: // (x,0,0) 6 permutations
231 {
232 _points[c] = Point( x, 0., 0.); _weights[c++] = wt;
233 _points[c] = Point(0., -x, 0.); _weights[c++] = wt;
234 _points[c] = Point(-x, 0., 0.); _weights[c++] = wt;
235 _points[c] = Point(0., x, 0.); _weights[c++] = wt;
236 _points[c] = Point(0., 0., -x); _weights[c++] = wt;
237 _points[c] = Point(0., 0., x); _weights[c++] = wt;
238
239 break;
240 }
241 case 2: // (x,x,0) 12 permutations
242 {
243 _points[c] = Point( x, x, 0.); _weights[c++] = wt;
244 _points[c] = Point( x, -x, 0.); _weights[c++] = wt;
245 _points[c] = Point(-x, -x, 0.); _weights[c++] = wt;
246 _points[c] = Point(-x, x, 0.); _weights[c++] = wt;
247 _points[c] = Point( x, 0., -x); _weights[c++] = wt;
248 _points[c] = Point( x, 0., x); _weights[c++] = wt;
249 _points[c] = Point(0., x, -x); _weights[c++] = wt;
250 _points[c] = Point(0., x, x); _weights[c++] = wt;
251 _points[c] = Point(0., -x, -x); _weights[c++] = wt;
252 _points[c] = Point(-x, 0., -x); _weights[c++] = wt;
253 _points[c] = Point(0., -x, x); _weights[c++] = wt;
254 _points[c] = Point(-x, 0., x); _weights[c++] = wt;
255
256 break;
257 }
258 case 3: // (x,y,0) 24 permutations
259 {
260 _points[c] = Point( x, y, 0.); _weights[c++] = wt;
261 _points[c] = Point( y, -x, 0.); _weights[c++] = wt;
262 _points[c] = Point(-x, -y, 0.); _weights[c++] = wt;
263 _points[c] = Point(-y, x, 0.); _weights[c++] = wt;
264 _points[c] = Point( x, 0., -y); _weights[c++] = wt;
265 _points[c] = Point( x, -y, 0.); _weights[c++] = wt;
266 _points[c] = Point( x, 0., y); _weights[c++] = wt;
267 _points[c] = Point(0., y, -x); _weights[c++] = wt;
268 _points[c] = Point(-x, y, 0.); _weights[c++] = wt;
269 _points[c] = Point(0., y, x); _weights[c++] = wt;
270 _points[c] = Point( y, 0., -x); _weights[c++] = wt;
271 _points[c] = Point(0., -y, -x); _weights[c++] = wt;
272 _points[c] = Point(-y, 0., -x); _weights[c++] = wt;
273 _points[c] = Point( y, x, 0.); _weights[c++] = wt;
274 _points[c] = Point(-y, -x, 0.); _weights[c++] = wt;
275 _points[c] = Point( y, 0., x); _weights[c++] = wt;
276 _points[c] = Point(0., -y, x); _weights[c++] = wt;
277 _points[c] = Point(-y, 0., x); _weights[c++] = wt;
278 _points[c] = Point(-x, 0., y); _weights[c++] = wt;
279 _points[c] = Point(0., -x, -y); _weights[c++] = wt;
280 _points[c] = Point(0., -x, y); _weights[c++] = wt;
281 _points[c] = Point(-x, 0., -y); _weights[c++] = wt;
282 _points[c] = Point(0., x, y); _weights[c++] = wt;
283 _points[c] = Point(0., x, -y); _weights[c++] = wt;
284
285 break;
286 }
287 case 4: // (x,x,x) 8 permutations
288 {
289 _points[c] = Point( x, x, x); _weights[c++] = wt;
290 _points[c] = Point( x, -x, x); _weights[c++] = wt;
291 _points[c] = Point(-x, -x, x); _weights[c++] = wt;
292 _points[c] = Point(-x, x, x); _weights[c++] = wt;
293 _points[c] = Point( x, x, -x); _weights[c++] = wt;
294 _points[c] = Point( x, -x, -x); _weights[c++] = wt;
295 _points[c] = Point(-x, x, -x); _weights[c++] = wt;
296 _points[c] = Point(-x, -x, -x); _weights[c++] = wt;
297
298 break;
299 }
300 case 5: // (x,x,z) 24 permutations
301 {
302 _points[c] = Point( x, x, z); _weights[c++] = wt;
303 _points[c] = Point( x, -x, z); _weights[c++] = wt;
304 _points[c] = Point(-x, -x, z); _weights[c++] = wt;
305 _points[c] = Point(-x, x, z); _weights[c++] = wt;
306 _points[c] = Point( x, z, -x); _weights[c++] = wt;
307 _points[c] = Point( x, -x, -z); _weights[c++] = wt;
308 _points[c] = Point( x, -z, x); _weights[c++] = wt;
309 _points[c] = Point( z, x, -x); _weights[c++] = wt;
310 _points[c] = Point(-x, x, -z); _weights[c++] = wt;
311 _points[c] = Point(-z, x, x); _weights[c++] = wt;
312 _points[c] = Point( x, -z, -x); _weights[c++] = wt;
313 _points[c] = Point(-z, -x, -x); _weights[c++] = wt;
314 _points[c] = Point(-x, z, -x); _weights[c++] = wt;
315 _points[c] = Point( x, x, -z); _weights[c++] = wt;
316 _points[c] = Point(-x, -x, -z); _weights[c++] = wt;
317 _points[c] = Point( x, z, x); _weights[c++] = wt;
318 _points[c] = Point( z, -x, x); _weights[c++] = wt;
319 _points[c] = Point(-x, -z, x); _weights[c++] = wt;
320 _points[c] = Point(-x, z, x); _weights[c++] = wt;
321 _points[c] = Point( z, -x, -x); _weights[c++] = wt;
322 _points[c] = Point(-z, -x, x); _weights[c++] = wt;
323 _points[c] = Point(-x, -z, -x); _weights[c++] = wt;
324 _points[c] = Point( z, x, x); _weights[c++] = wt;
325 _points[c] = Point(-z, x, -x); _weights[c++] = wt;
326
327 break;
328 }
329 case 6: // (x,y,z) 24 permutations
330 {
331 _points[c] = Point( x, y, z); _weights[c++] = wt;
332 _points[c] = Point( y, -x, z); _weights[c++] = wt;
333 _points[c] = Point(-x, -y, z); _weights[c++] = wt;
334 _points[c] = Point(-y, x, z); _weights[c++] = wt;
335 _points[c] = Point( x, z, -y); _weights[c++] = wt;
336 _points[c] = Point( x, -y, -z); _weights[c++] = wt;
337 _points[c] = Point( x, -z, y); _weights[c++] = wt;
338 _points[c] = Point( z, y, -x); _weights[c++] = wt;
339 _points[c] = Point(-x, y, -z); _weights[c++] = wt;
340 _points[c] = Point(-z, y, x); _weights[c++] = wt;
341 _points[c] = Point( y, -z, -x); _weights[c++] = wt;
342 _points[c] = Point(-z, -y, -x); _weights[c++] = wt;
343 _points[c] = Point(-y, z, -x); _weights[c++] = wt;
344 _points[c] = Point( y, x, -z); _weights[c++] = wt;
345 _points[c] = Point(-y, -x, -z); _weights[c++] = wt;
346 _points[c] = Point( y, z, x); _weights[c++] = wt;
347 _points[c] = Point( z, -y, x); _weights[c++] = wt;
348 _points[c] = Point(-y, -z, x); _weights[c++] = wt;
349 _points[c] = Point(-x, z, y); _weights[c++] = wt;
350 _points[c] = Point( z, -x, -y); _weights[c++] = wt;
351 _points[c] = Point(-z, -x, y); _weights[c++] = wt;
352 _points[c] = Point(-x, -z, -y); _weights[c++] = wt;
353 _points[c] = Point( z, x, y); _weights[c++] = wt;
354 _points[c] = Point(-z, x, -y); _weights[c++] = wt;
355
356 break;
357 }
358 default:
359 libmesh_error_msg("Unknown rule ID: " << rule_id[i] << "!");
360 } // end switch(rule_id[i])
361 }
362}
363
364} // namespace libMesh
libMesh::Point
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition point.h:40
libMesh::QBase::_points
std::vector< Point > _points
The locations of the quadrature points in reference element space.
Definition quadrature.h:409
libMesh::QBase::_weights
std::vector< Real > _weights
The quadrature weights.
Definition quadrature.h:415
libMesh::QMonomial::kim_rule
void kim_rule(const Real rule_data[][4], const unsigned int *rule_id, const unsigned int n_pts)
Rules from Kim and Song, Comm.
Definition quadrature_monomial.C:210
libMesh::QMonomial::clone
virtual std::unique_ptr< QBase > clone() const override
Definition quadrature_monomial.C:37
libMesh::QMonomial::stroud_rule
void stroud_rule(const Real rule_data[][3], const unsigned int *rule_symmetry, const unsigned int n_pts)
Stroud's rules for quads and hexes can have one of several different types of symmetry.
Definition quadrature_monomial.C:62
libMesh::QMonomial::wissmann_rule
void wissmann_rule(const Real rule_data[][3], const unsigned int n_pts)
Wissmann published three interesting "partially symmetric" rules for integrating degree 4,...
Definition quadrature_monomial.C:42
libMesh::QMonomial::type
virtual QuadratureType type() const override
Definition quadrature_monomial.C:32
libMesh
The libMesh namespace provides an interface to certain functionality in the library.
libMesh::QuadratureType
QuadratureType
Defines an enum for currently available quadrature rules.
Definition enum_quadrature_type.h:33
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition libmesh_common.h:144
Generated on Fri Aug 21 2026 13:32:30 for libMesh by doxygen 1.9.8