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 :
20 : // Local includes
21 : #include "libmesh/quadrature_nodal.h"
22 : #include "libmesh/quadrature_trap.h"
23 : #include "libmesh/quadrature_simpson.h"
24 : #include "libmesh/enum_to_string.h"
25 : #include "libmesh/cell_c0polyhedron.h"
26 :
27 : namespace libMesh
28 : {
29 :
30 78545 : void QNodal::init_3D()
31 : {
32 : #if LIBMESH_DIM == 3
33 :
34 78545 : switch (_type)
35 : {
36 6044 : case TET4:
37 : case PRISM6:
38 : case HEX8:
39 : case PYRAMID5:
40 : {
41 : // Construct QTrap rule that matches our own nodal pyramid quadrature permissions
42 6532 : QTrap rule(/*dim=*/3, /*ignored*/_order);
43 6044 : rule.allow_nodal_pyramid_quadrature = this->allow_nodal_pyramid_quadrature;
44 6044 : rule.init(*this);
45 488 : _points.swap (rule.get_points());
46 488 : _weights.swap(rule.get_weights());
47 488 : return;
48 : }
49 :
50 100 : case PRISM15:
51 : {
52 : // A rule with 15 points which is exact for linears, and
53 : // naturally produces a lumped approximation to the mass
54 : // matrix. The quadrature points are numbered the same way as
55 : // the reference element nodes.
56 2388 : _points =
57 : {
58 : Point(0.,0.,-1), Point(+1,0.,-1), Point(0.,+1,-1),
59 : Point(0.,0.,+1), Point(+1,0.,+1), Point(0.,+1,+1),
60 : Point(.5,0.,-1), Point(.5,.5,-1), Point(0.,.5,-1),
61 : Point(0.,0.,0.), Point(+1,0.,0.), Point(0.,+1,0.),
62 : Point(.5,0.,+1), Point(.5,.5,+1), Point(0.,.5,+1),
63 1294 : };
64 :
65 : // vertex (wv), tri edge (wt), and quad edge (wq) weights are
66 : // obtained using the same approach that was used for the Quad8,
67 : // see quadrature_nodal_2D.C for details.
68 100 : Real wv = Real(1) / 34;
69 100 : Real wt = Real(4) / 51;
70 100 : Real wq = Real(2) / 17;
71 :
72 2388 : _weights = {wv, wv, wv, wv, wv, wv,
73 : wt, wt, wt,
74 : wq, wq, wq,
75 1294 : wt, wt, wt};
76 :
77 1294 : return;
78 : }
79 :
80 196 : case HEX20:
81 : {
82 : // A rule with 20 points which is exact for linears, and
83 : // naturally produces a lumped approximation to the mass
84 : // matrix. The quadrature points are numbered the same way as
85 : // the reference element nodes.
86 4500 : _points =
87 : {
88 : Point(-1,-1,-1), Point(+1,-1,-1), Point(+1,+1,-1), Point(-1,+1,-1),
89 : Point(-1,-1,+1), Point(+1,-1,+1), Point(+1,+1,+1), Point(-1,+1,+1),
90 : Point(0.,-1,-1), Point(+1,0.,-1), Point(0.,+1,-1), Point(-1,0.,-1),
91 : Point(-1,-1,0.), Point(+1,-1,0.), Point(+1,+1,0.), Point(-1,+1,0.),
92 : Point(0.,-1,+1), Point(+1,0.,+1), Point(0.,+1,+1), Point(-1,0.,+1)
93 2446 : };
94 :
95 : // vertex (wv), and edge (we) weights are obtained using the
96 : // same approach that was used for the Quad8, see
97 : // quadrature_nodal_2D.C for details.
98 196 : Real wv = Real(7) / 31;
99 196 : Real we = Real(16) / 31;
100 :
101 4500 : _weights = {wv, wv, wv, wv, wv, wv, wv, wv,
102 2446 : we, we, we, we, we, we, we, we, we, we, we, we};
103 :
104 2446 : return;
105 : }
106 :
107 36209 : case TET10:
108 : case PRISM18:
109 : case HEX27:
110 : case PYRAMID13:
111 : case PYRAMID14:
112 : {
113 : // Construct QSimpson rule that matches our own nodal pyramid quadrature permissions
114 39199 : QSimpson rule(/*dim=*/3, /*ignored*/_order);
115 36209 : rule.allow_nodal_pyramid_quadrature = this->allow_nodal_pyramid_quadrature;
116 36209 : rule.init(*this);
117 36209 : _points.swap (rule.get_points());
118 36209 : _weights.swap(rule.get_weights());
119 :
120 : // We can't do a proper Simpson rule for pyramids regardless
121 36209 : if (_type == PYRAMID13)
122 : {
123 2375 : _points.resize(13);
124 2375 : _weights.resize(13);
125 : }
126 :
127 2990 : return;
128 : }
129 :
130 98 : case PRISM20:
131 : {
132 2250 : _points =
133 : {
134 : Point(0.,0.,-1), Point(+1,0.,-1), Point(0.,+1,-1),
135 : Point(0.,0.,+1), Point(+1,0.,+1), Point(0.,+1,+1),
136 : Point(.5,0.,-1), Point(.5,.5,-1), Point(0.,.5,-1),
137 : Point(0.,0.,0.), Point(+1,0.,0.), Point(0.,+1,0.),
138 : Point(.5,0.,+1), Point(.5,.5,+1), Point(0.,.5,+1),
139 : Point(.5,0.,0.), Point(.5,.5,0.), Point(0.,.5,0.),
140 : Point(1/Real(3),1/Real(3),-1), Point(1/Real(3),1/Real(3),+1)
141 1223 : };
142 :
143 : // Symmetry gives us identical weights on vertices, triangle
144 : // edges, square edges, triangle faces, square faces; then we
145 : // have the midnode. Solving for weights which exactly
146 : // integrate cubics and xi^2*zeta^2 would give a unique answer
147 : // ... with wse=0 and wte<0. See Quad8 in
148 : // quadrature_nodal_2D.C for discussion of this problem.
149 : //
150 : // Dropping the xi^2 zeta^2 constraint gives us a rank-1 null
151 : // space to play with ... but adding anything from that null
152 : // space to push wte closer to positive just pushes wse
153 : // negative.
154 : //
155 : // Dropping exact integration of xi^3 type terms gives us a
156 : // rank-2 null space. I just found the max(min(weight)) from
157 : // that solution space using Octave. Someone who cares more
158 : // than me might want to repeat this exercise with better than
159 : // double precision...
160 :
161 98 : constexpr Real wv = 8.546754839782711e-03;
162 98 : constexpr Real wte = 8.548403936750599e-03;
163 98 : constexpr Real wse = wv; // here's our min(weight) constraint...
164 98 : constexpr Real wtf = 1.153811903370667e-01;
165 98 : constexpr Real wsf = 2.136754673824396e-01;
166 :
167 2250 : _weights = {wv, wv, wv, wv, wv, wv,
168 : wte, wte, wte, wse, wse, wse, wte, wte, wte,
169 1223 : wsf, wsf, wsf, wtf, wtf};
170 :
171 1223 : return;
172 : }
173 :
174 98 : case PRISM21:
175 : {
176 2250 : _points =
177 : {
178 : Point(0.,0.,-1), Point(+1,0.,-1), Point(0.,+1,-1),
179 : Point(0.,0.,+1), Point(+1,0.,+1), Point(0.,+1,+1),
180 : Point(.5,0.,-1), Point(.5,.5,-1), Point(0.,.5,-1),
181 : Point(0.,0.,0.), Point(+1,0.,0.), Point(0.,+1,0.),
182 : Point(.5,0.,+1), Point(.5,.5,+1), Point(0.,.5,+1),
183 : Point(.5,0.,0.), Point(.5,.5,0.), Point(0.,.5,0.),
184 : Point(1/Real(3),1/Real(3),-1), Point(1/Real(3),1/Real(3),+1),
185 : Point(1/Real(3),1/Real(3),0.)
186 1223 : };
187 :
188 : // Symmetry gives us identical weights on vertices, triangle
189 : // edges, square edges, triangle faces, square faces; then we
190 : // have the midnode. Solving for weights which exactly
191 : // integrate cubics, xi^2*zeta^2, and xi^3*zeta^2, gives a
192 : // unique answer:
193 98 : constexpr Real wv = Real(1)/120;
194 98 : constexpr Real wte = Real(1)/45;
195 98 : constexpr Real wse = Real(1)/30;
196 98 : constexpr Real wtf = Real(3)/40;
197 98 : constexpr Real wsf = Real(4)/45;
198 :
199 2250 : _weights = {wv, wv, wv, wv, wv, wv,
200 : wte, wte, wte, wse, wse, wse, wte, wte, wte,
201 1223 : wsf, wsf, wsf, wtf, wtf, Real(3)/10};
202 :
203 1223 : return;
204 : }
205 :
206 2375 : case PYRAMID18:
207 : {
208 2375 : libmesh_error_msg_if(!allow_nodal_pyramid_quadrature,
209 : "Nodal quadrature on Pyramid elements is not allowed by default since\n"
210 : "the Jacobian of the inverse element map is not well-defined at the Pyramid apex.\n"
211 : "Set the QBase::allow_nodal_pyramid_quadrature flag to true to ignore skip this check.");
212 :
213 4362 : _points =
214 : {
215 : Point(-1,-1,0.), Point(+1,-1,0.), Point(+1,+1,0.),
216 : Point(-1,+1,0.), Point(0.,0.,+1), Point(0.,-1,0.),
217 : Point(+1,0.,0.), Point(0.,+1,0.), Point(-1,0.,0.),
218 : Point(-.5,-.5,.5), Point(-.5,.5,.5), Point(.5,.5,.5),
219 : Point(-.5,.5,.5), Point(0.,0.,0.), Point(0.,-2/Real(3),1/Real(3)),
220 : Point(2/Real(3),0.,1/Real(3)), Point(0.,2/Real(3),1/Real(3)),
221 : Point(-2/Real(3),0.,1/Real(3))
222 2375 : };
223 :
224 : // Even with triangle faces to play with, I can't seem to get
225 : // exact integrals of any higher order functions without
226 : // negative weights. So I punt and just use QTrap weights.
227 4362 : _weights = {1/Real(4), 1/Real(4), 1/Real(4), 1/Real(4), 1/Real(3),
228 2375 : 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
229 :
230 2375 : return;
231 : }
232 :
233 27719 : case TET14:
234 : {
235 27719 : _points.resize(14);
236 27719 : _weights.resize(14);
237 :
238 27719 : _points[0](0) = 0.; _points[5](0) = .5;
239 27719 : _points[0](1) = 0.; _points[5](1) = .5;
240 27719 : _points[0](2) = 0.; _points[5](2) = 0.;
241 :
242 27719 : _points[1](0) = 1.; _points[6](0) = 0.;
243 27719 : _points[1](1) = 0.; _points[6](1) = .5;
244 27719 : _points[1](2) = 0.; _points[6](2) = 0.;
245 :
246 27719 : _points[2](0) = 0.; _points[7](0) = 0.;
247 27719 : _points[2](1) = 1.; _points[7](1) = 0.;
248 27719 : _points[2](2) = 0.; _points[7](2) = .5;
249 :
250 27719 : _points[3](0) = 0.; _points[8](0) = .5;
251 27719 : _points[3](1) = 0.; _points[8](1) = 0.;
252 27719 : _points[3](2) = 1.; _points[8](2) = .5;
253 :
254 27719 : _points[4](0) = .5; _points[9](0) = 0.;
255 27719 : _points[4](1) = 0.; _points[9](1) = .5;
256 27719 : _points[4](2) = 0.; _points[9](2) = .5;
257 :
258 :
259 27719 : _points[10](0) = 1/Real(3); _points[11](0) = 1/Real(3);
260 27719 : _points[10](1) = 1/Real(3); _points[11](1) = 0.;
261 27719 : _points[10](2) = 0.; _points[11](2) = 1/Real(3);
262 :
263 27719 : _points[12](0) = 1/Real(3); _points[13](0) = 0.;
264 27719 : _points[12](1) = 1/Real(3); _points[13](1) = 1/Real(3);
265 27719 : _points[12](2) = 1/Real(3); _points[13](2) = 1/Real(3);
266 :
267 : // RHS:
268 : //
269 : // The ``optimal'' nodal quadrature here, the one that
270 : // integrates every Lagrange polynomial on these nodes
271 : // exactly, produces an indefinite mass matrix ...
272 : //
273 : // const Real wv = 1/Real(240);
274 : // const Real we = 0;
275 : // const Real wf = 3/Real(80);
276 :
277 : // We could average that with our Tet10 rule and get:
278 : //
279 : // const Real wv = (1/Real(240)+1/Real(192))/2;
280 : // const Real we = Real(14)/576/2;
281 : // const Real wf = 3/Real(80)/2;
282 : //
283 : // Except our Tet10 rule won't actually exactly integrate
284 : // quadratics! (exactly integrating quadratics wouldn't even
285 : // have given a positive semidefinite mass matrix there...)
286 : //
287 : // John derived equations for wv and we based on symmetry and
288 : // the requirement to exactly integrate quadratics; within
289 : // those constraints we might pick the wf that maximizes the
290 : // minimum nodal weight:
291 : // const Real wf = Real(15)/440;
292 : // const Real wv = -1/Real(120) + wf/3;
293 : // const Real we = 1/Real(30) - wf*(Real(8)/9);
294 : //
295 : // But John also did the same clever optimization trick that
296 : // quadrature_nodal_2D.C discusses in the context of Quad8 and
297 : // Hex20 outputs, and for Tet14 that gives us:
298 2306 : const Real wv = Real(87)/47120;
299 2306 : const Real we = Real(164)/26505;
300 2306 : const Real wf = Real(1439)/47120;
301 :
302 27719 : _weights = {wv, wv, wv, wv, we, we, we, we, we, we, wf, wf, wf, wf};
303 27719 : return;
304 : }
305 :
306 : //---------------------------------------------
307 : // Arbitrary polygon quadrature rules
308 12 : case C0POLYHEDRON:
309 : {
310 : // Polyhedra require the newer Quadrature API
311 12 : if (!_elem)
312 0 : libmesh_error();
313 :
314 1 : libmesh_assert(_elem->type() == C0POLYHEDRON);
315 :
316 1 : const C0Polyhedron & poly = *cast_ptr<const C0Polyhedron *>(_elem);
317 :
318 12 : const Real master_poly_volume = _elem->volume();
319 :
320 : // For polyhedra the master element is the physical element.
321 : //
322 : // Using even weights is not the most effective nodal rule for
323 : // that case, but users who want an effective rule probably
324 : // want a non-nodal rule anyway.
325 12 : const unsigned int nn = poly.n_nodes();
326 12 : const Real weight = master_poly_volume / nn;
327 12 : _weights.resize(nn, weight);
328 :
329 12 : _points.resize(poly.n_nodes());
330 108 : for (auto n : make_range(nn))
331 96 : _points[n] = poly.master_point(n);
332 :
333 1 : return;
334 : }
335 :
336 :
337 0 : default:
338 0 : libmesh_error_msg("ERROR: Unsupported type: " << Utility::enum_to_string(_type));
339 : }
340 : #endif
341 : }
342 :
343 : } // namespace libMesh
|
