Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2025 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 284 : void QMonomial::wissmann_rule(const Real rule_data[][3],
43 : const unsigned int n_pts)
44 : {
45 1775 : for (unsigned int i=0, c=0; i<n_pts; ++i)
46 : {
47 1491 : _points[c] = Point( rule_data[i][0], rule_data[i][1] );
48 1491 : _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 1491 : if (rule_data[i][0] != static_cast<Real>(0.0))
53 : {
54 994 : _points[c] = Point( -rule_data[i][0], rule_data[i][1] );
55 1022 : _weights[c++] = rule_data[i][2];
56 : }
57 : }
58 284 : }
59 :
60 :
61 :
62 710 : void QMonomial::stroud_rule(const Real rule_data[][3],
63 : const unsigned int * rule_symmetry,
64 : const unsigned int n_pts)
65 : {
66 5609 : for (unsigned int i=0, c=0; i<n_pts; ++i)
67 : {
68 : const Real
69 4899 : x=rule_data[i][0],
70 4899 : y=rule_data[i][1],
71 4899 : wt=rule_data[i][2];
72 :
73 4899 : switch(rule_symmetry[i])
74 : {
75 8 : case 0: // Single point (no symmetry)
76 : {
77 284 : _points[c] = Point( x, y);
78 284 : _weights[c++] = wt;
79 :
80 284 : break;
81 : }
82 22 : case 1: // Fully-symmetric (x,y)
83 : {
84 781 : _points[c] = Point( x, y);
85 781 : _weights[c++] = wt;
86 :
87 781 : _points[c] = Point(-x, y);
88 781 : _weights[c++] = wt;
89 :
90 781 : _points[c] = Point( x,-y);
91 781 : _weights[c++] = wt;
92 :
93 781 : _points[c] = Point(-x,-y);
94 781 : _weights[c++] = wt;
95 :
96 781 : _points[c] = Point( y, x);
97 781 : _weights[c++] = wt;
98 :
99 781 : _points[c] = Point(-y, x);
100 781 : _weights[c++] = wt;
101 :
102 781 : _points[c] = Point( y,-x);
103 781 : _weights[c++] = wt;
104 :
105 781 : _points[c] = Point(-y,-x);
106 781 : _weights[c++] = wt;
107 :
108 781 : break;
109 : }
110 22 : case 2: // Fully-symmetric (x,x)
111 : {
112 781 : _points[c] = Point( x, x);
113 781 : _weights[c++] = wt;
114 :
115 781 : _points[c] = Point(-x, x);
116 781 : _weights[c++] = wt;
117 :
118 781 : _points[c] = Point( x,-x);
119 781 : _weights[c++] = wt;
120 :
121 781 : _points[c] = Point(-x,-x);
122 781 : _weights[c++] = wt;
123 :
124 781 : break;
125 : }
126 18 : case 3: // Fully-symmetric (x,0)
127 : {
128 18 : libmesh_assert_equal_to (y, 0.0);
129 :
130 639 : _points[c] = Point( x,0.);
131 639 : _weights[c++] = wt;
132 :
133 639 : _points[c] = Point(-x,0.);
134 639 : _weights[c++] = wt;
135 :
136 639 : _points[c] = Point(0., x);
137 639 : _weights[c++] = wt;
138 :
139 639 : _points[c] = Point(0.,-x);
140 639 : _weights[c++] = wt;
141 :
142 639 : break;
143 : }
144 64 : case 4: // Rotational invariant
145 : {
146 2272 : _points[c] = Point( x, y);
147 2272 : _weights[c++] = wt;
148 :
149 2272 : _points[c] = Point(-x,-y);
150 2272 : _weights[c++] = wt;
151 :
152 2272 : _points[c] = Point(-y, x);
153 2272 : _weights[c++] = wt;
154 :
155 2272 : _points[c] = Point( y,-x);
156 2272 : _weights[c++] = wt;
157 :
158 2272 : 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 71 : _points[c] = Point( x, y);
175 71 : _weights[c++] = wt;
176 :
177 71 : _points[c] = Point(-x, y);
178 71 : _weights[c++] = wt;
179 :
180 71 : _points[c] = Point(-x,-y);
181 71 : _weights[c++] = wt;
182 :
183 71 : _points[c] = Point( x,-y);
184 71 : _weights[c++] = wt;
185 :
186 71 : 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 71 : _points[c] = Point(0., y);
194 71 : _weights[c++] = wt;
195 :
196 71 : _points[c] = Point(0.,-y);
197 71 : _weights[c++] = wt;
198 :
199 71 : break;
200 : }
201 0 : default:
202 0 : libmesh_error_msg("Unknown symmetry!");
203 : } // end switch(rule_symmetry[i])
204 : }
205 710 : }
206 :
207 :
208 :
209 :
210 355 : void QMonomial::kim_rule(const Real rule_data[][4],
211 : const unsigned int * rule_id,
212 : const unsigned int n_pts)
213 : {
214 1278 : for (unsigned int i=0, c=0; i<n_pts; ++i)
215 : {
216 : const Real
217 923 : x=rule_data[i][0],
218 923 : y=rule_data[i][1],
219 923 : z=rule_data[i][2],
220 923 : wt=rule_data[i][3];
221 :
222 923 : switch(rule_id[i])
223 : {
224 6 : case 0: // (0,0,0) 1 permutation
225 : {
226 219 : _points[c] = Point( x, y, z); _weights[c++] = wt;
227 :
228 213 : break;
229 : }
230 10 : case 1: // (x,0,0) 6 permutations
231 : {
232 365 : _points[c] = Point( x, 0., 0.); _weights[c++] = wt;
233 365 : _points[c] = Point(0., -x, 0.); _weights[c++] = wt;
234 365 : _points[c] = Point(-x, 0., 0.); _weights[c++] = wt;
235 365 : _points[c] = Point(0., x, 0.); _weights[c++] = wt;
236 365 : _points[c] = Point(0., 0., -x); _weights[c++] = wt;
237 365 : _points[c] = Point(0., 0., x); _weights[c++] = wt;
238 :
239 355 : 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 146 : _points[c] = Point( x, x, x); _weights[c++] = wt;
290 146 : _points[c] = Point( x, -x, x); _weights[c++] = wt;
291 146 : _points[c] = Point(-x, -x, x); _weights[c++] = wt;
292 146 : _points[c] = Point(-x, x, x); _weights[c++] = wt;
293 146 : _points[c] = Point( x, x, -x); _weights[c++] = wt;
294 146 : _points[c] = Point( x, -x, -x); _weights[c++] = wt;
295 146 : _points[c] = Point(-x, x, -x); _weights[c++] = wt;
296 146 : _points[c] = Point(-x, -x, -x); _weights[c++] = wt;
297 :
298 142 : 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 219 : _points[c] = Point( x, y, z); _weights[c++] = wt;
332 219 : _points[c] = Point( y, -x, z); _weights[c++] = wt;
333 219 : _points[c] = Point(-x, -y, z); _weights[c++] = wt;
334 219 : _points[c] = Point(-y, x, z); _weights[c++] = wt;
335 219 : _points[c] = Point( x, z, -y); _weights[c++] = wt;
336 219 : _points[c] = Point( x, -y, -z); _weights[c++] = wt;
337 219 : _points[c] = Point( x, -z, y); _weights[c++] = wt;
338 219 : _points[c] = Point( z, y, -x); _weights[c++] = wt;
339 219 : _points[c] = Point(-x, y, -z); _weights[c++] = wt;
340 219 : _points[c] = Point(-z, y, x); _weights[c++] = wt;
341 219 : _points[c] = Point( y, -z, -x); _weights[c++] = wt;
342 219 : _points[c] = Point(-z, -y, -x); _weights[c++] = wt;
343 219 : _points[c] = Point(-y, z, -x); _weights[c++] = wt;
344 219 : _points[c] = Point( y, x, -z); _weights[c++] = wt;
345 219 : _points[c] = Point(-y, -x, -z); _weights[c++] = wt;
346 219 : _points[c] = Point( y, z, x); _weights[c++] = wt;
347 219 : _points[c] = Point( z, -y, x); _weights[c++] = wt;
348 219 : _points[c] = Point(-y, -z, x); _weights[c++] = wt;
349 219 : _points[c] = Point(-x, z, y); _weights[c++] = wt;
350 219 : _points[c] = Point( z, -x, -y); _weights[c++] = wt;
351 219 : _points[c] = Point(-z, -x, y); _weights[c++] = wt;
352 219 : _points[c] = Point(-x, -z, -y); _weights[c++] = wt;
353 219 : _points[c] = Point( z, x, y); _weights[c++] = wt;
354 219 : _points[c] = Point(-z, x, -y); _weights[c++] = wt;
355 :
356 213 : break;
357 : }
358 0 : default:
359 0 : libmesh_error_msg("Unknown rule ID: " << rule_id[i] << "!");
360 : } // end switch(rule_id[i])
361 : }
362 355 : }
363 :
364 : } // namespace libMesh
|
