Line data Source code
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 :
23 : namespace 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 :
32 0 : QuadratureType QMonomial::type() const
33 : {
34 0 : return QMONOMIAL;
35 : }
36 :
37 0 : std::unique_ptr<QBase> QMonomial::clone() const
38 : {
39 0 : return std::make_unique<QMonomial>(*this);
40 : }
41 :
42 28 : void QMonomial::wissmann_rule(const Real rule_data[][3],
43 : const unsigned int n_pts)
44 : {
45 175 : for (unsigned int i=0, c=0; i<n_pts; ++i)
46 : {
47 147 : _points[c] = Point( rule_data[i][0], rule_data[i][1] );
48 147 : _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 147 : if (rule_data[i][0] != static_cast<Real>(0.0))
53 : {
54 98 : _points[c] = Point( -rule_data[i][0], rule_data[i][1] );
55 126 : _weights[c++] = rule_data[i][2];
56 : }
57 : }
58 28 : }
59 :
60 :
61 :
62 70 : void QMonomial::stroud_rule(const Real rule_data[][3],
63 : const unsigned int * rule_symmetry,
64 : const unsigned int n_pts)
65 : {
66 553 : for (unsigned int i=0, c=0; i<n_pts; ++i)
67 : {
68 : const Real
69 483 : x=rule_data[i][0],
70 483 : y=rule_data[i][1],
71 483 : wt=rule_data[i][2];
72 :
73 483 : switch(rule_symmetry[i])
74 : {
75 8 : case 0: // Single point (no symmetry)
76 : {
77 28 : _points[c] = Point( x, y);
78 28 : _weights[c++] = wt;
79 :
80 28 : break;
81 : }
82 22 : case 1: // Fully-symmetric (x,y)
83 : {
84 77 : _points[c] = Point( x, y);
85 77 : _weights[c++] = wt;
86 :
87 77 : _points[c] = Point(-x, y);
88 77 : _weights[c++] = wt;
89 :
90 77 : _points[c] = Point( x,-y);
91 77 : _weights[c++] = wt;
92 :
93 77 : _points[c] = Point(-x,-y);
94 77 : _weights[c++] = wt;
95 :
96 77 : _points[c] = Point( y, x);
97 77 : _weights[c++] = wt;
98 :
99 77 : _points[c] = Point(-y, x);
100 77 : _weights[c++] = wt;
101 :
102 77 : _points[c] = Point( y,-x);
103 77 : _weights[c++] = wt;
104 :
105 77 : _points[c] = Point(-y,-x);
106 77 : _weights[c++] = wt;
107 :
108 77 : break;
109 : }
110 22 : case 2: // Fully-symmetric (x,x)
111 : {
112 77 : _points[c] = Point( x, x);
113 77 : _weights[c++] = wt;
114 :
115 77 : _points[c] = Point(-x, x);
116 77 : _weights[c++] = wt;
117 :
118 77 : _points[c] = Point( x,-x);
119 77 : _weights[c++] = wt;
120 :
121 77 : _points[c] = Point(-x,-x);
122 77 : _weights[c++] = wt;
123 :
124 77 : break;
125 : }
126 18 : case 3: // Fully-symmetric (x,0)
127 : {
128 18 : libmesh_assert_equal_to (y, 0.0);
129 :
130 63 : _points[c] = Point( x,0.);
131 63 : _weights[c++] = wt;
132 :
133 63 : _points[c] = Point(-x,0.);
134 63 : _weights[c++] = wt;
135 :
136 63 : _points[c] = Point(0., x);
137 63 : _weights[c++] = wt;
138 :
139 63 : _points[c] = Point(0.,-x);
140 63 : _weights[c++] = wt;
141 :
142 63 : break;
143 : }
144 64 : case 4: // Rotational invariant
145 : {
146 224 : _points[c] = Point( x, y);
147 224 : _weights[c++] = wt;
148 :
149 224 : _points[c] = Point(-x,-y);
150 224 : _weights[c++] = wt;
151 :
152 224 : _points[c] = Point(-y, x);
153 224 : _weights[c++] = wt;
154 :
155 224 : _points[c] = Point( y,-x);
156 224 : _weights[c++] = wt;
157 :
158 224 : break;
159 : }
160 0 : case 5: // Partial symmetry (Wissman's rules)
161 : {
162 0 : libmesh_assert_not_equal_to (x, 0.0);
163 :
164 0 : _points[c] = Point( x, y);
165 0 : _weights[c++] = wt;
166 :
167 0 : _points[c] = Point(-x, y);
168 0 : _weights[c++] = wt;
169 :
170 0 : break;
171 : }
172 2 : case 6: // Rectangular symmetry
173 : {
174 7 : _points[c] = Point( x, y);
175 7 : _weights[c++] = wt;
176 :
177 7 : _points[c] = Point(-x, y);
178 7 : _weights[c++] = wt;
179 :
180 7 : _points[c] = Point(-x,-y);
181 7 : _weights[c++] = wt;
182 :
183 7 : _points[c] = Point( x,-y);
184 7 : _weights[c++] = wt;
185 :
186 7 : break;
187 : }
188 2 : case 7: // Central symmetry
189 : {
190 2 : libmesh_assert_equal_to (x, 0.0);
191 2 : libmesh_assert_not_equal_to (y, 0.0);
192 :
193 7 : _points[c] = Point(0., y);
194 7 : _weights[c++] = wt;
195 :
196 7 : _points[c] = Point(0.,-y);
197 7 : _weights[c++] = wt;
198 :
199 7 : break;
200 : }
201 0 : default:
202 0 : libmesh_error_msg("Unknown symmetry!");
203 : } // end switch(rule_symmetry[i])
204 : }
205 70 : }
206 :
207 :
208 :
209 :
210 35 : void QMonomial::kim_rule(const Real rule_data[][4],
211 : const unsigned int * rule_id,
212 : const unsigned int n_pts)
213 : {
214 126 : for (unsigned int i=0, c=0; i<n_pts; ++i)
215 : {
216 : const Real
217 91 : x=rule_data[i][0],
218 91 : y=rule_data[i][1],
219 91 : z=rule_data[i][2],
220 91 : wt=rule_data[i][3];
221 :
222 91 : switch(rule_id[i])
223 : {
224 6 : case 0: // (0,0,0) 1 permutation
225 : {
226 27 : _points[c] = Point( x, y, z); _weights[c++] = wt;
227 :
228 21 : break;
229 : }
230 10 : case 1: // (x,0,0) 6 permutations
231 : {
232 45 : _points[c] = Point( x, 0., 0.); _weights[c++] = wt;
233 45 : _points[c] = Point(0., -x, 0.); _weights[c++] = wt;
234 45 : _points[c] = Point(-x, 0., 0.); _weights[c++] = wt;
235 45 : _points[c] = Point(0., x, 0.); _weights[c++] = wt;
236 45 : _points[c] = Point(0., 0., -x); _weights[c++] = wt;
237 45 : _points[c] = Point(0., 0., x); _weights[c++] = wt;
238 :
239 35 : break;
240 : }
241 0 : case 2: // (x,x,0) 12 permutations
242 : {
243 0 : _points[c] = Point( x, x, 0.); _weights[c++] = wt;
244 0 : _points[c] = Point( x, -x, 0.); _weights[c++] = wt;
245 0 : _points[c] = Point(-x, -x, 0.); _weights[c++] = wt;
246 0 : _points[c] = Point(-x, x, 0.); _weights[c++] = wt;
247 0 : _points[c] = Point( x, 0., -x); _weights[c++] = wt;
248 0 : _points[c] = Point( x, 0., x); _weights[c++] = wt;
249 0 : _points[c] = Point(0., x, -x); _weights[c++] = wt;
250 0 : _points[c] = Point(0., x, x); _weights[c++] = wt;
251 0 : _points[c] = Point(0., -x, -x); _weights[c++] = wt;
252 0 : _points[c] = Point(-x, 0., -x); _weights[c++] = wt;
253 0 : _points[c] = Point(0., -x, x); _weights[c++] = wt;
254 0 : _points[c] = Point(-x, 0., x); _weights[c++] = wt;
255 :
256 0 : break;
257 : }
258 0 : case 3: // (x,y,0) 24 permutations
259 : {
260 0 : _points[c] = Point( x, y, 0.); _weights[c++] = wt;
261 0 : _points[c] = Point( y, -x, 0.); _weights[c++] = wt;
262 0 : _points[c] = Point(-x, -y, 0.); _weights[c++] = wt;
263 0 : _points[c] = Point(-y, x, 0.); _weights[c++] = wt;
264 0 : _points[c] = Point( x, 0., -y); _weights[c++] = wt;
265 0 : _points[c] = Point( x, -y, 0.); _weights[c++] = wt;
266 0 : _points[c] = Point( x, 0., y); _weights[c++] = wt;
267 0 : _points[c] = Point(0., y, -x); _weights[c++] = wt;
268 0 : _points[c] = Point(-x, y, 0.); _weights[c++] = wt;
269 0 : _points[c] = Point(0., y, x); _weights[c++] = wt;
270 0 : _points[c] = Point( y, 0., -x); _weights[c++] = wt;
271 0 : _points[c] = Point(0., -y, -x); _weights[c++] = wt;
272 0 : _points[c] = Point(-y, 0., -x); _weights[c++] = wt;
273 0 : _points[c] = Point( y, x, 0.); _weights[c++] = wt;
274 0 : _points[c] = Point(-y, -x, 0.); _weights[c++] = wt;
275 0 : _points[c] = Point( y, 0., x); _weights[c++] = wt;
276 0 : _points[c] = Point(0., -y, x); _weights[c++] = wt;
277 0 : _points[c] = Point(-y, 0., x); _weights[c++] = wt;
278 0 : _points[c] = Point(-x, 0., y); _weights[c++] = wt;
279 0 : _points[c] = Point(0., -x, -y); _weights[c++] = wt;
280 0 : _points[c] = Point(0., -x, y); _weights[c++] = wt;
281 0 : _points[c] = Point(-x, 0., -y); _weights[c++] = wt;
282 0 : _points[c] = Point(0., x, y); _weights[c++] = wt;
283 0 : _points[c] = Point(0., x, -y); _weights[c++] = wt;
284 :
285 0 : break;
286 : }
287 4 : case 4: // (x,x,x) 8 permutations
288 : {
289 18 : _points[c] = Point( x, x, x); _weights[c++] = wt;
290 18 : _points[c] = Point( x, -x, x); _weights[c++] = wt;
291 18 : _points[c] = Point(-x, -x, x); _weights[c++] = wt;
292 18 : _points[c] = Point(-x, x, x); _weights[c++] = wt;
293 18 : _points[c] = Point( x, x, -x); _weights[c++] = wt;
294 18 : _points[c] = Point( x, -x, -x); _weights[c++] = wt;
295 18 : _points[c] = Point(-x, x, -x); _weights[c++] = wt;
296 18 : _points[c] = Point(-x, -x, -x); _weights[c++] = wt;
297 :
298 14 : break;
299 : }
300 0 : case 5: // (x,x,z) 24 permutations
301 : {
302 0 : _points[c] = Point( x, x, z); _weights[c++] = wt;
303 0 : _points[c] = Point( x, -x, z); _weights[c++] = wt;
304 0 : _points[c] = Point(-x, -x, z); _weights[c++] = wt;
305 0 : _points[c] = Point(-x, x, z); _weights[c++] = wt;
306 0 : _points[c] = Point( x, z, -x); _weights[c++] = wt;
307 0 : _points[c] = Point( x, -x, -z); _weights[c++] = wt;
308 0 : _points[c] = Point( x, -z, x); _weights[c++] = wt;
309 0 : _points[c] = Point( z, x, -x); _weights[c++] = wt;
310 0 : _points[c] = Point(-x, x, -z); _weights[c++] = wt;
311 0 : _points[c] = Point(-z, x, x); _weights[c++] = wt;
312 0 : _points[c] = Point( x, -z, -x); _weights[c++] = wt;
313 0 : _points[c] = Point(-z, -x, -x); _weights[c++] = wt;
314 0 : _points[c] = Point(-x, z, -x); _weights[c++] = wt;
315 0 : _points[c] = Point( x, x, -z); _weights[c++] = wt;
316 0 : _points[c] = Point(-x, -x, -z); _weights[c++] = wt;
317 0 : _points[c] = Point( x, z, x); _weights[c++] = wt;
318 0 : _points[c] = Point( z, -x, x); _weights[c++] = wt;
319 0 : _points[c] = Point(-x, -z, x); _weights[c++] = wt;
320 0 : _points[c] = Point(-x, z, x); _weights[c++] = wt;
321 0 : _points[c] = Point( z, -x, -x); _weights[c++] = wt;
322 0 : _points[c] = Point(-z, -x, x); _weights[c++] = wt;
323 0 : _points[c] = Point(-x, -z, -x); _weights[c++] = wt;
324 0 : _points[c] = Point( z, x, x); _weights[c++] = wt;
325 0 : _points[c] = Point(-z, x, -x); _weights[c++] = wt;
326 :
327 0 : break;
328 : }
329 6 : case 6: // (x,y,z) 24 permutations
330 : {
331 27 : _points[c] = Point( x, y, z); _weights[c++] = wt;
332 27 : _points[c] = Point( y, -x, z); _weights[c++] = wt;
333 27 : _points[c] = Point(-x, -y, z); _weights[c++] = wt;
334 27 : _points[c] = Point(-y, x, z); _weights[c++] = wt;
335 27 : _points[c] = Point( x, z, -y); _weights[c++] = wt;
336 27 : _points[c] = Point( x, -y, -z); _weights[c++] = wt;
337 27 : _points[c] = Point( x, -z, y); _weights[c++] = wt;
338 27 : _points[c] = Point( z, y, -x); _weights[c++] = wt;
339 27 : _points[c] = Point(-x, y, -z); _weights[c++] = wt;
340 27 : _points[c] = Point(-z, y, x); _weights[c++] = wt;
341 27 : _points[c] = Point( y, -z, -x); _weights[c++] = wt;
342 27 : _points[c] = Point(-z, -y, -x); _weights[c++] = wt;
343 27 : _points[c] = Point(-y, z, -x); _weights[c++] = wt;
344 27 : _points[c] = Point( y, x, -z); _weights[c++] = wt;
345 27 : _points[c] = Point(-y, -x, -z); _weights[c++] = wt;
346 27 : _points[c] = Point( y, z, x); _weights[c++] = wt;
347 27 : _points[c] = Point( z, -y, x); _weights[c++] = wt;
348 27 : _points[c] = Point(-y, -z, x); _weights[c++] = wt;
349 27 : _points[c] = Point(-x, z, y); _weights[c++] = wt;
350 27 : _points[c] = Point( z, -x, -y); _weights[c++] = wt;
351 27 : _points[c] = Point(-z, -x, y); _weights[c++] = wt;
352 27 : _points[c] = Point(-x, -z, -y); _weights[c++] = wt;
353 27 : _points[c] = Point( z, x, y); _weights[c++] = wt;
354 27 : _points[c] = Point(-z, x, -y); _weights[c++] = wt;
355 :
356 21 : break;
357 : }
358 0 : default:
359 0 : libmesh_error_msg("Unknown rule ID: " << rule_id[i] << "!");
360 : } // end switch(rule_id[i])
361 : }
362 35 : }
363 :
364 : } // namespace libMesh
|
