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 : // Local includes
20 : #include "libmesh/cell_hex8.h"
21 : #include "libmesh/edge_edge2.h"
22 : #include "libmesh/face_quad4.h"
23 : #include "libmesh/enum_io_package.h"
24 : #include "libmesh/enum_order.h"
25 : #include "libmesh/fe_lagrange_shape_1D.h"
26 :
27 : // C++ includes
28 : #include <array>
29 :
30 : namespace libMesh
31 : {
32 :
33 :
34 :
35 :
36 : // ------------------------------------------------------------
37 : // Hex8 class static member initializations
38 : const int Hex8::num_nodes;
39 : const int Hex8::nodes_per_side;
40 : const int Hex8::nodes_per_edge;
41 :
42 : // ------------------------------------------------------------
43 : // Hex8 class member functions
44 :
45 18430801 : bool Hex8::is_vertex(const unsigned int libmesh_dbg_var(n)) const
46 : {
47 1344568 : libmesh_assert_not_equal_to (n, invalid_uint);
48 18430801 : return true;
49 : }
50 :
51 0 : bool Hex8::is_edge(const unsigned int) const
52 : {
53 0 : return false;
54 : }
55 :
56 0 : bool Hex8::is_face(const unsigned int) const
57 : {
58 0 : return false;
59 : }
60 :
61 930114 : bool Hex8::is_node_on_side(const unsigned int n,
62 : const unsigned int s) const
63 : {
64 77245 : libmesh_assert_less (s, n_sides());
65 77245 : return std::find(std::begin(side_nodes_map[s]),
66 930114 : std::end(side_nodes_map[s]),
67 930114 : n) != std::end(side_nodes_map[s]);
68 : }
69 :
70 : std::vector<unsigned>
71 63750105 : Hex8::nodes_on_side(const unsigned int s) const
72 : {
73 4708962 : libmesh_assert_less(s, n_sides());
74 68459067 : return {std::begin(side_nodes_map[s]), std::end(side_nodes_map[s])};
75 : }
76 :
77 : std::vector<unsigned>
78 7764 : Hex8::nodes_on_edge(const unsigned int e) const
79 : {
80 600 : libmesh_assert_less(e, n_edges());
81 8364 : return {std::begin(edge_nodes_map[e]), std::end(edge_nodes_map[e])};
82 : }
83 :
84 1625184 : bool Hex8::is_node_on_edge(const unsigned int n,
85 : const unsigned int e) const
86 : {
87 62824 : libmesh_assert_less (e, n_edges());
88 62824 : return std::find(std::begin(edge_nodes_map[e]),
89 1625184 : std::end(edge_nodes_map[e]),
90 1625184 : n) != std::end(edge_nodes_map[e]);
91 : }
92 :
93 :
94 :
95 5907517 : bool Hex8::has_affine_map() const
96 : {
97 : // Make sure x-edge endpoints are affine
98 985236 : Point v = this->point(1) - this->point(0);
99 5825203 : if (!v.relative_fuzzy_equals(this->point(2) - this->point(3), affine_tol) ||
100 11719340 : !v.relative_fuzzy_equals(this->point(5) - this->point(4), affine_tol) ||
101 6303733 : !v.relative_fuzzy_equals(this->point(6) - this->point(7), affine_tol))
102 96506 : return false;
103 : // Make sure xz-faces are identical parallelograms
104 6297953 : v = this->point(4) - this->point(0);
105 5811011 : if (!v.relative_fuzzy_equals(this->point(7) - this->point(3), affine_tol))
106 24 : return false;
107 : // If all the above checks out, the map is affine
108 486930 : return true;
109 : }
110 :
111 :
112 :
113 137451576 : Order Hex8::default_order() const
114 : {
115 137451576 : return FIRST;
116 : }
117 :
118 :
119 :
120 1291068 : std::unique_ptr<Elem> Hex8::build_side_ptr (const unsigned int i)
121 : {
122 1291068 : return this->simple_build_side_ptr<Quad4, Hex8>(i);
123 : }
124 :
125 :
126 :
127 546488 : void Hex8::build_side_ptr (std::unique_ptr<Elem> & side,
128 : const unsigned int i)
129 : {
130 546488 : this->simple_build_side_ptr<Hex8>(side, i, QUAD4);
131 546488 : }
132 :
133 :
134 :
135 6379110 : std::unique_ptr<Elem> Hex8::build_edge_ptr (const unsigned int i)
136 : {
137 6379110 : return this->simple_build_edge_ptr<Edge2,Hex8>(i);
138 : }
139 :
140 :
141 :
142 135440 : void Hex8::build_edge_ptr (std::unique_ptr<Elem> & edge, const unsigned int i)
143 : {
144 135440 : this->simple_build_edge_ptr<Hex8>(edge, i, EDGE2);
145 135440 : }
146 :
147 :
148 :
149 234943 : void Hex8::connectivity(const unsigned int libmesh_dbg_var(sc),
150 : const IOPackage iop,
151 : std::vector<dof_id_type> & conn) const
152 : {
153 20253 : libmesh_assert(_nodes);
154 20253 : libmesh_assert_less (sc, this->n_sub_elem());
155 20253 : libmesh_assert_not_equal_to (iop, INVALID_IO_PACKAGE);
156 :
157 234943 : conn.resize(8);
158 :
159 234943 : switch (iop)
160 : {
161 234943 : case TECPLOT:
162 : {
163 255196 : conn[0] = this->node_id(0)+1;
164 234943 : conn[1] = this->node_id(1)+1;
165 234943 : conn[2] = this->node_id(2)+1;
166 234943 : conn[3] = this->node_id(3)+1;
167 234943 : conn[4] = this->node_id(4)+1;
168 234943 : conn[5] = this->node_id(5)+1;
169 234943 : conn[6] = this->node_id(6)+1;
170 234943 : conn[7] = this->node_id(7)+1;
171 234943 : return;
172 : }
173 :
174 0 : case VTK:
175 : {
176 0 : for (auto i : index_range(conn))
177 0 : conn[i] = this->node_id(i);
178 0 : return;
179 : }
180 :
181 0 : default:
182 0 : libmesh_error_msg("Unsupported IO package " << iop);
183 : }
184 : }
185 :
186 :
187 :
188 : #ifdef LIBMESH_ENABLE_AMR
189 :
190 : const Real Hex8::_embedding_matrix[Hex8::num_children][Hex8::num_nodes][Hex8::num_nodes] =
191 : {
192 : // The 8 children of the Hex-type elements can be thought of as being
193 : // associated with the 8 vertices of the Hex. Some of the children are
194 : // numbered the same as their corresponding vertex, while some are
195 : // not. The children which are numbered differently have been marked
196 : // with ** in the comments below.
197 :
198 : // embedding matrix for child 0 (child 0 is associated with vertex 0)
199 : {
200 : // 0 1 2 3 4 5 6 7
201 : { 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, // 0
202 : { 0.5, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, // 1
203 : { .25, .25, .25, .25, 0.0, 0.0, 0.0, 0.0}, // 2
204 : { 0.5, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0}, // 3
205 : { 0.5, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0}, // 4
206 : { .25, .25, 0.0, 0.0, .25, .25, 0.0, 0.0}, // 5
207 : {.125, .125, .125, .125, .125, .125, .125, .125}, // 6
208 : { .25, 0.0, 0.0, .25, .25, 0.0, 0.0, .25} // 7
209 : },
210 :
211 : // embedding matrix for child 1 (child 1 is associated with vertex 1)
212 : {
213 : // 0 1 2 3 4 5 6 7
214 : { 0.5, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, // 0
215 : { 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, // 1
216 : { 0.0, 0.5, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0}, // 2
217 : { .25, .25, .25, .25, 0.0, 0.0, 0.0, 0.0}, // 3
218 : { .25, .25, 0.0, 0.0, .25, .25, 0.0, 0.0}, // 4
219 : { 0.0, 0.5, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0}, // 5
220 : { 0.0, .25, .25, 0.0, 0.0, .25, .25, 0.0}, // 6
221 : {.125, .125, .125, .125, .125, .125, .125, .125} // 7
222 : },
223 :
224 : // embedding matrix for child 2 (child 2 is associated with vertex 3**)
225 : {
226 : // 0 1 2 3 4 5 6 7
227 : { 0.5, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0}, // 0
228 : { .25, .25, .25, .25, 0.0, 0.0, 0.0, 0.0}, // 1
229 : { 0.0, 0.0, 0.5, 0.5, 0.0, 0.0, 0.0, 0.0}, // 2
230 : { 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0}, // 3
231 : { .25, 0.0, 0.0, .25, .25, 0.0, 0.0, .25}, // 4
232 : {.125, .125, .125, .125, .125, .125, .125, .125}, // 5
233 : { 0.0, 0.0, .25, .25, 0.0, 0.0, .25, .25}, // 6
234 : { 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.5} // 7
235 : },
236 :
237 : // embedding matrix for child 3 (child 3 is associated with vertex 2**)
238 : {
239 : // 0 1 2 3 4 5 6 7
240 : { .25, .25, .25, .25, 0.0, 0.0, 0.0, 0.0}, // 0
241 : { 0.0, 0.5, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0}, // 1
242 : { 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0}, // 2
243 : { 0.0, 0.0, 0.5, 0.5, 0.0, 0.0, 0.0, 0.0}, // 3
244 : {.125, .125, .125, .125, .125, .125, .125, .125}, // 4
245 : { 0.0, .25, .25, 0.0, 0.0, .25, .25, 0.0}, // 5
246 : { 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.5, 0.0}, // 6
247 : { 0.0, 0.0, .25, .25, 0.0, 0.0, .25, .25} // 7
248 : },
249 :
250 : // embedding matrix for child 4 (child 4 is associated with vertex 4)
251 : {
252 : // 0 1 2 3 4 5 6 7
253 : { 0.5, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0}, // 0
254 : { .25, .25, 0.0, 0.0, .25, .25, 0.0, 0.0}, // 1
255 : {.125, .125, .125, .125, .125, .125, .125, .125}, // 2
256 : { .25, 0.0, 0.0, .25, .25, 0.0, 0.0, .25}, // 3
257 : { 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0}, // 4
258 : { 0.0, 0.0, 0.0, 0.0, 0.5, 0.5, 0.0, 0.0}, // 5
259 : { 0.0, 0.0, 0.0, 0.0, .25, .25, .25, .25}, // 6
260 : { 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.5} // 7
261 : },
262 :
263 : // embedding matrix for child 5 (child 5 is associated with vertex 5)
264 : {
265 : // 0 1 2 3 4 5 6 7
266 : { .25, .25, 0.0, 0.0, .25, .25, 0.0, 0.0}, // 0
267 : { 0.0, 0.5, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0}, // 1
268 : { 0.0, .25, .25, 0.0, 0.0, .25, .25, 0.0}, // 2
269 : {.125, .125, .125, .125, .125, .125, .125, .125}, // 3
270 : { 0.0, 0.0, 0.0, 0.0, 0.5, 0.5, 0.0, 0.0}, // 4
271 : { 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0}, // 5
272 : { 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.5, 0.0}, // 6
273 : { 0.0, 0.0, 0.0, 0.0, .25, .25, .25, .25} // 7
274 : },
275 :
276 : // embedding matrix for child 6 (child 6 is associated with vertex 7**)
277 : {
278 : // 0 1 2 3 4 5 6 7
279 : { .25, 0.0, 0.0, .25, .25, 0.0, 0.0, .25}, // 0
280 : {.125, .125, .125, .125, .125, .125, .125, .125}, // 1
281 : { 0.0, 0.0, .25, .25, 0.0, 0.0, .25, .25}, // 2
282 : { 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.5}, // 3
283 : { 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.5}, // 4
284 : { 0.0, 0.0, 0.0, 0.0, .25, .25, .25, .25}, // 5
285 : { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.5}, // 6
286 : { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0} // 7
287 : },
288 :
289 : // embedding matrix for child 7 (child 7 is associated with vertex 6**)
290 : {
291 : // 0 1 2 3 4 5 6 7
292 : {.125, .125, .125, .125, .125, .125, .125, .125}, // 0
293 : { 0.0, .25, .25, 0.0, 0.0, .25, .25, 0.0}, // 1
294 : { 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.5, 0.0}, // 2
295 : { 0.0, 0.0, .25, .25, 0.0, 0.0, .25, .25}, // 3
296 : { 0.0, 0.0, 0.0, 0.0, .25, .25, .25, .25}, // 4
297 : { 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.5, 0.0}, // 5
298 : { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0}, // 6
299 : { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.5} // 7
300 : }
301 : };
302 :
303 :
304 :
305 :
306 : #endif
307 :
308 :
309 :
310 9327 : Point Hex8::centroid_from_points(
311 : const Point & x0, const Point & x1, const Point & x2, const Point & x3,
312 : const Point & x4, const Point & x5, const Point & x6, const Point & x7)
313 : {
314 : // The Jacobian is dx/d(xi) dot (dx/d(eta) cross dx/d(zeta)), where
315 : // dx/d(xi) = a1*eta*zeta + b1*eta + c1*zeta + d1
316 : // dx/d(eta) = a2*xi*zeta + b2*xi + c2*zeta + d2
317 : // dx/d(zeta) = a3*xi*eta + b3*xi + c3*eta + d3
318 :
319 : // Notes:
320 : // 1.) Several of these coefficient vectors are equal, as noted below.
321 : // 2.) These are all off by a factor of 8, but this cancels when we
322 : // divide by the volume, which will also be off by the same
323 : // factor.
324 : Point
325 554 : a1 = -x0 + x1 - x2 + x3 + x4 - x5 + x6 - x7,
326 554 : a2 = a1,
327 554 : a3 = a1;
328 :
329 : Point
330 554 : b1 = x0 - x1 + x2 - x3 + x4 - x5 + x6 - x7,
331 554 : b2 = b1,
332 554 : b3 = x0 - x1 - x2 + x3 - x4 + x5 + x6 - x7;
333 :
334 : Point
335 554 : c1 = b3,
336 554 : c2 = x0 + x1 - x2 - x3 - x4 - x5 + x6 + x7,
337 554 : c3 = c2;
338 :
339 : Point
340 554 : d1 = -x0 + x1 + x2 - x3 - x4 + x5 + x6 - x7,
341 554 : d2 = -x0 - x1 + x2 + x3 - x4 - x5 + x6 + x7,
342 554 : d3 = -x0 - x1 - x2 - x3 + x4 + x5 + x6 + x7;
343 :
344 : // Use 2x2x2 quadrature to compute the integral of each basis
345 : // function (as defined on the [-1,1]^3 reference domain). We use
346 : // a quadrature rule which is exact for tri-cubics. The weights for
347 : // this rule are all equal to 1.
348 : static const Real q[2] = {-std::sqrt(Real(3))/3, std::sqrt(Real(3))/3.};
349 :
350 : // Indices for computing tensor product basis functions. This is
351 : // copied from fe_lagrange_shape_3D.C
352 : static const unsigned int i0[] = {0, 1, 1, 0, 0, 1, 1, 0};
353 : static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 1, 1};
354 : static const unsigned int i2[] = {0, 0, 0, 0, 1, 1, 1, 1};
355 :
356 : // Compute nodal volumes
357 9327 : std::array<Real, Hex8::num_nodes> V{};
358 :
359 27981 : for (const auto & xi : q)
360 55962 : for (const auto & eta : q)
361 111924 : for (const auto & zeta : q)
362 : {
363 70184 : Real jxw = triple_product(a1*eta*zeta + b1*eta + c1*zeta + d1,
364 74616 : a2*xi*zeta + b2*xi + c2*zeta + d2,
365 79048 : a3*xi*eta + b3*xi + c3*eta + d3);
366 :
367 671544 : for (int i=0; i<Hex8::num_nodes; ++i)
368 1158400 : V[i] += jxw *
369 596928 : fe_lagrange_1D_linear_shape(i0[i], xi) *
370 1158400 : fe_lagrange_1D_linear_shape(i1[i], eta) *
371 596928 : fe_lagrange_1D_linear_shape(i2[i], zeta);
372 : }
373 :
374 : // Compute centroid
375 : return
376 554 : (x0*V[0] + x1*V[1] + x2*V[2] + x3*V[3] + x4*V[4] + x5*V[5] + x6*V[6] + x7*V[7]) /
377 10435 : (V[0] + V[1] + V[2] + V[3] + V[4] + V[5] + V[6] + V[7]);
378 : }
379 :
380 :
381 2968 : Point Hex8::true_centroid () const
382 : {
383 : return Hex8::centroid_from_points
384 2536 : (point(0), point(1), point(2), point(3),
385 3184 : point(4), point(5), point(6), point(7));
386 : }
387 :
388 :
389 248807 : Real Hex8::volume () const
390 : {
391 : // This specialization is good for Lagrange mappings only in general
392 248807 : if (this->mapping_type() != LAGRANGE_MAP)
393 44 : return this->Elem::volume();
394 :
395 : // Make copies of our points. It makes the subsequent calculations a bit
396 : // shorter and avoids dereferencing the same pointer multiple times.
397 : Point
398 331219 : x0 = point(0), x1 = point(1), x2 = point(2), x3 = point(3),
399 310605 : x4 = point(4), x5 = point(5), x6 = point(6), x7 = point(7);
400 :
401 : // Construct constant data vectors. The notation is:
402 : // \vec{x}_{\xi} = \vec{a1}*eta*zeta + \vec{b1}*eta + \vec{c1}*zeta + \vec{d1}
403 : // \vec{x}_{\eta} = \vec{a2}*xi*zeta + \vec{b2}*xi + \vec{c2}*zeta + \vec{d2}
404 : // \vec{x}_{\zeta} = \vec{a3}*xi*eta + \vec{b3}*xi + \vec{c3}*eta + \vec{d3}
405 : // but it turns out that a1, a2, and a3 are not needed for the volume calculation.
406 :
407 : // Build up the 6 unique vectors which make up dx/dxi, dx/deta, and dx/dzeta.
408 : Point q[6] =
409 : {
410 21043 : /*b1*/ x0 - x1 + x2 - x3 + x4 - x5 + x6 - x7, /*=b2*/
411 21043 : /*c1*/ x0 - x1 - x2 + x3 - x4 + x5 + x6 - x7, /*=b3*/
412 21043 : /*d1*/ -x0 + x1 + x2 - x3 - x4 + x5 + x6 - x7,
413 21043 : /*c2*/ x0 + x1 - x2 - x3 - x4 - x5 + x6 + x7, /*=c3*/
414 21043 : /*d2*/ -x0 - x1 + x2 + x3 - x4 - x5 + x6 + x7,
415 21043 : /*d3*/ -x0 - x1 - x2 - x3 + x4 + x5 + x6 + x7
416 147301 : };
417 :
418 : // We could check for a linear element, but it's probably faster to
419 : // just compute the result...
420 : return
421 248763 : (triple_product(q[0], q[4], q[3]) +
422 248763 : triple_product(q[2], q[0], q[1]) +
423 248763 : triple_product(q[1], q[3], q[5])) / 192. +
424 248763 : triple_product(q[2], q[4], q[5]) / 64.;
425 : }
426 :
427 : BoundingBox
428 10295000 : Hex8::loose_bounding_box () const
429 : {
430 10295000 : return Elem::loose_bounding_box();
431 : }
432 :
433 :
434 : void
435 8532 : Hex8::permute(unsigned int perm_num)
436 : {
437 792 : libmesh_assert_less (perm_num, 24);
438 8532 : const unsigned int side = perm_num % 6;
439 8532 : const unsigned int rotate = perm_num / 6;
440 :
441 23760 : for (unsigned int i = 0; i != rotate; ++i)
442 : {
443 15228 : swap4nodes(0,1,2,3);
444 15228 : swap4nodes(4,5,6,7);
445 13716 : swap4neighbors(1,2,3,4);
446 : }
447 :
448 8532 : switch (side) {
449 96 : case 0:
450 96 : break;
451 1152 : case 1:
452 1152 : swap4nodes(3,7,4,0);
453 1152 : swap4nodes(2,6,5,1);
454 1056 : swap4neighbors(0,3,5,1);
455 1056 : break;
456 1152 : case 2:
457 1152 : swap4nodes(0,4,5,1);
458 1152 : swap4nodes(3,7,6,2);
459 1056 : swap4neighbors(0,4,5,2);
460 1056 : break;
461 1152 : case 3:
462 1152 : swap4nodes(0,4,7,3);
463 1152 : swap4nodes(1,5,6,2);
464 1056 : swap4neighbors(0,1,5,3);
465 1056 : break;
466 1152 : case 4:
467 1152 : swap4nodes(1,5,4,0);
468 1152 : swap4nodes(2,6,7,3);
469 1056 : swap4neighbors(0,2,5,4);
470 1056 : break;
471 2772 : case 5:
472 2772 : swap2nodes(0,7);
473 2772 : swap2nodes(1,6);
474 2772 : swap2nodes(2,5);
475 2772 : swap2nodes(3,4);
476 312 : swap2neighbors(0,5);
477 312 : swap2neighbors(1,3);
478 312 : break;
479 0 : default:
480 0 : libmesh_error();
481 : }
482 8532 : }
483 :
484 :
485 : void
486 288 : Hex8::flip(BoundaryInfo * boundary_info)
487 : {
488 24 : libmesh_assert(boundary_info);
489 :
490 288 : swap2nodes(0,1);
491 288 : swap2nodes(2,3);
492 288 : swap2nodes(4,5);
493 288 : swap2nodes(6,7);
494 24 : swap2neighbors(2,4);
495 288 : swap2boundarysides(2,4,boundary_info);
496 288 : swap2boundaryedges(1,3,boundary_info);
497 288 : swap2boundaryedges(4,5,boundary_info);
498 288 : swap2boundaryedges(6,7,boundary_info);
499 288 : swap2boundaryedges(9,11,boundary_info);
500 288 : }
501 :
502 :
503 : ElemType
504 126704454 : Hex8::side_type (const unsigned int libmesh_dbg_var(s)) const
505 : {
506 9378651 : libmesh_assert_less (s, 6);
507 126704454 : return QUAD4;
508 : }
509 :
510 :
511 : Point
512 852 : Hex8::side_vertex_average_normal(const unsigned int s) const
513 : {
514 24 : libmesh_assert_less (s, 6);
515 24 : libmesh_assert_equal_to(this->mapping_type(), LAGRANGE_MAP);
516 :
517 : // At the side vertex average, things simplify a bit
518 : // We get the side "plane" normal at all vertices, then average them
519 24 : Point normal;
520 900 : Point current_edge = this->point(side_nodes_map[s][1]) -
521 852 : this->point(side_nodes_map[s][0]);
522 4260 : for (auto i : make_range(4))
523 : {
524 3600 : const Point next_edge = this->point(side_nodes_map[s][(i + 2) % 4]) -
525 3408 : this->point(side_nodes_map[s][(i + 1) % 4]);
526 3312 : const Point normal_at_vertex = current_edge.cross(next_edge);
527 : // Note: unweighted and unnormalized sum matches the FE computation of normals
528 96 : normal += normal_at_vertex;
529 3312 : current_edge = next_edge;
530 : }
531 24 : normal /= 4;
532 852 : return normal.unit();
533 : }
534 :
535 :
536 : } // namespace libMesh
|
