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 :
20 : // Local includes
21 : #include "libmesh/quadrature_grid.h"
22 : #include "libmesh/enum_to_string.h"
23 :
24 : namespace libMesh
25 : {
26 :
27 :
28 1917 : void QGrid::init_3D()
29 : {
30 : #if LIBMESH_DIM == 3
31 :
32 1917 : switch (_type)
33 : {
34 : //---------------------------------------------
35 : // Hex quadrature rules
36 639 : case HEX8:
37 : case HEX20:
38 : case HEX27:
39 : {
40 : // We compute the 3D quadrature rule as a tensor
41 : // product of the 1D quadrature rule.
42 : //
43 : // We ignore p_level in QGrid, since the "order" here is a
44 : // user-requested number of grid points, nothing to do with
45 : // the order of exactly-integrated polynomial spaces
46 657 : QGrid q1D(1,_order);
47 :
48 639 : tensor_product_hex( q1D );
49 :
50 18 : return;
51 : }
52 :
53 :
54 :
55 : //---------------------------------------------
56 : // Tetrahedral quadrature rules
57 639 : case TET4:
58 : case TET10:
59 : case TET14:
60 : {
61 : // This is the composite centroid rule on a uniform refinement
62 : // of the reference tet with order+1 intervals per edge: one
63 : // point at the centroid of each of the (order+1)^3 subtets.
64 : // The refinement consists of upright tets, inverted tets,
65 : // and octahedra, and we split each octahedron into four tets
66 : // about the same diagonal. Every subtet has the same volume
67 : // regardless of which diagonal is chosen, so the weights are
68 : // equal. The rule is exact for linears and has O(dx^2) error
69 : // on smooth functions, like the midpoint grid on hexes.
70 639 : const unsigned int np = (_order+1)*(_order+1)*(_order+1);
71 : // Master tet has 1x1 triangle base, height 1, so volume = 1/6
72 639 : const Real weight = Real(1)/Real(6)/np;
73 639 : const Real dx = Real(1)/static_cast<int>(_order+1);
74 639 : _points.resize(np);
75 639 : _weights.resize(np);
76 :
77 639 : unsigned int pt = 0;
78 250992 : auto add_point = [this, &pt, weight, dx](Real x, Real y, Real z)
79 : {
80 213276 : _points[pt] = Point(x*dx, y*dx, z*dx);
81 210084 : _weights[pt] = weight;
82 210306 : pt++;
83 209277 : };
84 :
85 4491 : for (int i = 0; i != _order + 1; ++i)
86 : {
87 19383 : for (int j = 0; j != _order + 1 - i; ++j)
88 : {
89 66243 : for (int k = 0; k != _order + 1 - i - j; ++k)
90 : {
91 : // Centroid of the upright tet with vertices
92 : // (i,j,k), (i+1,j,k), (i,j+1,k), (i,j,k+1)
93 50694 : add_point(i+Real(1)/4, j+Real(1)/4, k+Real(1)/4);
94 :
95 : // The octahedron with vertices (i+1,j,k),
96 : // (i,j+1,k), (i,j,k+1), (i+1,j+1,k),
97 : // (i+1,j,k+1), (i,j+1,k+1) lies inside the
98 : // reference tet only if i+j+k < order. We split
99 : // it into four tets about the diagonal from
100 : // (i+1,j,k) to (i,j+1,k+1); their centroids are
101 : // offset from the octahedron center by
102 : // (-1/4,0,0), (1/4,0,0), (0,1/4,-1/4), and
103 : // (0,-1/4,1/4).
104 50694 : if (i + j + k < _order)
105 : {
106 35145 : const Real ci = i+Real(1)/2;
107 35145 : const Real cj = j+Real(1)/2;
108 35145 : const Real ck = k+Real(1)/2;
109 35145 : add_point(ci-Real(1)/4, cj, ck);
110 35145 : add_point(ci+Real(1)/4, cj, ck);
111 35145 : add_point(ci, cj+Real(1)/4, ck-Real(1)/4);
112 35145 : add_point(ci, cj-Real(1)/4, ck+Real(1)/4);
113 : }
114 :
115 : // Centroid of the inverted tet with vertices
116 : // (i+1,j+1,k), (i+1,j,k+1), (i,j+1,k+1),
117 : // (i+1,j+1,k+1), which lies inside the reference
118 : // tet only if i+j+k < order-1
119 50694 : if (i + j + k + 1 < _order)
120 23430 : add_point(i+Real(3)/4, j+Real(3)/4, k+Real(3)/4);
121 : }
122 : }
123 : }
124 18 : libmesh_assert_equal_to(pt, np);
125 18 : return;
126 : }
127 :
128 :
129 : // Prism quadrature rules
130 0 : case PRISM6:
131 : case PRISM15:
132 : case PRISM18:
133 : case PRISM20:
134 : case PRISM21:
135 : {
136 : // We compute the 3D quadrature rule as a tensor
137 : // product of the 1D quadrature rule and a 2D
138 : // triangle quadrature rule
139 : //
140 : // We ignore p_level in QGrid, since the "order" here is a
141 : // user-requested number of grid points, nothing to do with
142 : // the order of exactly-integrated polynomial spaces
143 0 : QGrid q1D(1,_order);
144 0 : QGrid q2D(2,_order);
145 :
146 : // Initialize the 2D rule (1D is pre-initialized)
147 0 : q2D.init(TRI3, 0, /*simple_type_only=*/true);
148 :
149 0 : tensor_product_prism(q1D, q2D);
150 :
151 0 : return;
152 : }
153 :
154 :
155 :
156 : //---------------------------------------------
157 : // Pyramid
158 639 : case PYRAMID5:
159 : case PYRAMID13:
160 : case PYRAMID14:
161 : case PYRAMID18:
162 : {
163 : // This is the composite centroid rule on a uniform refinement
164 : // of the reference pyramid with order+1 intervals per edge.
165 : // Layer k of the refinement, between z = k*dz and z =
166 : // (k+1)*dz, has an m x m grid of cells on its bottom face and
167 : // an (m-1) x (m-1) grid on its top face, where m = order+1-k.
168 : // It consists of
169 : // * m^2 upright pyramids, one on each bottom cell,
170 : // * (m-1)^2 inverted pyramids, one under each top cell, and
171 : // * 2m(m-1) tets, one per interior edge of the bottom grid,
172 : // between the two upright pyramids sharing that edge.
173 : // Each pyramid has twice the volume of each tet, so we split
174 : // every pyramid into two tets along a base diagonal, giving
175 : // 2(order+1)^3 tets of equal volume. We place one point at
176 : // the centroid of each tet. The rule is exact for linears
177 : // and has O(dz^2) error on smooth functions.
178 639 : const unsigned int np = 2*(_order+1)*(_order+1)*(_order+1);
179 639 : _points.resize(np);
180 639 : _weights.resize(np);
181 : // Master pyramid has 2x2 base, height 1, so volume = 4/3
182 639 : const Real weight = Real(4)/Real(3)/np;
183 639 : const Real dz = Real(1)/static_cast<int>(_order+1);
184 :
185 : // Lattice point (i,j) on the bottom of layer k is at
186 : // ((2i+k)*dz - 1, (2j+k)*dz - 1, k*dz). This map is affine,
187 : // so each centroid below is computed as the average of its
188 : // tet's vertices in (i,j,k) lattice coordinates.
189 639 : unsigned int pt = 0;
190 538272 : auto add_point = [this, &pt, weight, dz](Real i, Real j, Real k)
191 : {
192 453600 : _points[pt] = Point((2*i+k)*dz - 1, (2*j+k)*dz - 1, k*dz);
193 441504 : _weights[pt] = weight;
194 429408 : pt++;
195 417951 : };
196 :
197 4491 : for (int k = 0; k != _order + 1; ++k)
198 : {
199 108 : const int m = _order + 1 - k;
200 19383 : for (int i = 0; i != m; ++i)
201 : {
202 101388 : for (int j = 0; j != m; ++j)
203 : {
204 : // The upright pyramid with base vertices
205 : // (i,j,k), (i+1,j,k), (i+1,j+1,k), (i,j+1,k) and
206 : // apex (i,j,k+1), split along the diagonal from
207 : // (i,j,k) to (i+1,j+1,k)
208 85839 : add_point(i+Real(1)/2, j+Real(1)/4, k+Real(1)/4);
209 85839 : add_point(i+Real(1)/4, j+Real(1)/2, k+Real(1)/4);
210 :
211 : // The inverted pyramid with base vertices
212 : // (i,j,k+1), (i+1,j,k+1), (i+1,j+1,k+1),
213 : // (i,j+1,k+1) and apex (i+1,j+1,k), split along
214 : // the diagonal from (i,j,k+1) to (i+1,j+1,k+1)
215 85839 : if (i + 1 < m && j + 1 < m)
216 : {
217 58575 : add_point(i+Real(3)/4, j+Real(1)/2, k+Real(3)/4);
218 58575 : add_point(i+Real(1)/2, j+Real(3)/4, k+Real(3)/4);
219 : }
220 :
221 : // The tet with vertices (i+1,j,k), (i+1,j+1,k),
222 : // (i,j,k+1), (i+1,j,k+1), between the upright
223 : // pyramids on cells (i,j) and (i+1,j)
224 85839 : if (i + 1 < m)
225 70290 : add_point(i+Real(3)/4, j+Real(1)/4, k+Real(1)/2);
226 :
227 : // The tet with vertices (i,j+1,k), (i+1,j+1,k),
228 : // (i,j,k+1), (i,j+1,k+1), between the upright
229 : // pyramids on cells (i,j) and (i,j+1)
230 85839 : if (j + 1 < m)
231 70290 : add_point(i+Real(1)/4, j+Real(3)/4, k+Real(1)/2);
232 : }
233 : }
234 : }
235 18 : libmesh_assert_equal_to(pt, np);
236 18 : return;
237 : }
238 :
239 :
240 :
241 : //---------------------------------------------
242 : // Unsupported type
243 0 : default:
244 0 : libmesh_error_msg("ERROR: Unsupported type: " << Utility::enum_to_string(_type));
245 : }
246 : #endif
247 : }
248 :
249 : } // namespace libMesh
|
