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 : // C++ includes
20 : #include <algorithm> // for std::min, std::max
21 :
22 : // Local includes
23 : #include "libmesh/cell_hex.h"
24 : #include "libmesh/cell_hex8.h"
25 : #include "libmesh/face_quad4.h"
26 : #include "libmesh/enum_elem_quality.h"
27 : #include "libmesh/tensor_value.h"
28 :
29 : #include <array>
30 :
31 : namespace libMesh
32 : {
33 :
34 :
35 :
36 : // ------------------------------------------------------------
37 : // Hex class static member initializations
38 : const int Hex::num_sides;
39 : const int Hex::num_edges;
40 : const int Hex::num_children;
41 :
42 : const Real Hex::_master_points[27][3] =
43 : {
44 : {-1, -1, -1},
45 : {1, -1, -1},
46 : {1, 1, -1},
47 : {-1, 1, -1},
48 : {-1, -1, 1},
49 : {1, -1, 1},
50 : {1, 1, 1},
51 : {-1, 1, 1},
52 : {0, -1, -1},
53 : {1, 0, -1},
54 : {0, 1, -1},
55 : {-1, 0, -1},
56 : {-1, -1, 0},
57 : {1, -1, 0},
58 : {1, 1, 0},
59 : {-1, 1, 0},
60 : {0, -1, 1},
61 : {1, 0, 1},
62 : {0, 1, 1},
63 : {-1, 0, 1},
64 : {0, 0, -1},
65 : {0, -1, 0},
66 : {1, 0, 0},
67 : {0, 1, 0},
68 : {-1, 0, 0},
69 : {0, 0, 1}
70 : };
71 :
72 : const unsigned int Hex::edge_sides_map[12][2] =
73 : {
74 : {0, 1}, // Edge 0
75 : {0, 2}, // Edge 1
76 : {0, 3}, // Edge 2
77 : {0, 4}, // Edge 3
78 : {1, 4}, // Edge 4
79 : {1, 2}, // Edge 5
80 : {2, 3}, // Edge 6
81 : {3, 4}, // Edge 7
82 : {1, 5}, // Edge 8
83 : {2, 5}, // Edge 9
84 : {3, 5}, // Edge 10
85 : {4, 5} // Edge 11
86 : };
87 :
88 : const unsigned int Hex::adjacent_edges_map[/*num_vertices*/8][/*n_adjacent_edges*/3] =
89 : {
90 : {0, 3, 4}, // Edges adjacent to node 0
91 : {0, 1, 5}, // Edges adjacent to node 1
92 : {1, 2, 6}, // Edges adjacent to node 2
93 : {2, 3, 7}, // Edges adjacent to node 3
94 : {4, 8, 11}, // Edges adjacent to node 4
95 : {5, 8, 9}, // Edges adjacent to node 5
96 : {6, 9, 10}, // Edges adjacent to node 6
97 : {7, 10, 11} // Edges adjacent to node 7
98 : };
99 :
100 : // ------------------------------------------------------------
101 : // Hex class member functions
102 144576 : dof_id_type Hex::key (const unsigned int s) const
103 : {
104 12048 : libmesh_assert_less (s, this->n_sides());
105 :
106 156624 : return this->compute_key(this->node_id(Hex8::side_nodes_map[s][0]),
107 144576 : this->node_id(Hex8::side_nodes_map[s][1]),
108 144576 : this->node_id(Hex8::side_nodes_map[s][2]),
109 156624 : this->node_id(Hex8::side_nodes_map[s][3]));
110 : }
111 :
112 :
113 :
114 33988992 : dof_id_type Hex::low_order_key (const unsigned int s) const
115 : {
116 882916 : libmesh_assert_less (s, this->n_sides());
117 :
118 34871896 : return this->compute_key(this->node_id(Hex8::side_nodes_map[s][0]),
119 33988992 : this->node_id(Hex8::side_nodes_map[s][1]),
120 33988992 : this->node_id(Hex8::side_nodes_map[s][2]),
121 34871908 : this->node_id(Hex8::side_nodes_map[s][3]));
122 : }
123 :
124 :
125 :
126 256468423 : unsigned int Hex::local_side_node(unsigned int side,
127 : unsigned int side_node) const
128 : {
129 19011538 : libmesh_assert_less (side, this->n_sides());
130 19011538 : libmesh_assert_less (side_node, Hex8::nodes_per_side);
131 :
132 256468423 : return Hex8::side_nodes_map[side][side_node];
133 : }
134 :
135 :
136 :
137 162503440 : unsigned int Hex::local_edge_node(unsigned int edge,
138 : unsigned int edge_node) const
139 : {
140 12194510 : libmesh_assert_less (edge, this->n_edges());
141 12194510 : libmesh_assert_less (edge_node, Hex8::nodes_per_edge);
142 :
143 162503440 : return Hex8::edge_nodes_map[edge][edge_node];
144 : }
145 :
146 :
147 :
148 204467 : std::unique_ptr<Elem> Hex::side_ptr (const unsigned int i)
149 : {
150 11154 : libmesh_assert_less (i, this->n_sides());
151 :
152 204467 : std::unique_ptr<Elem> face = std::make_unique<Quad4>();
153 :
154 1022335 : for (auto n : face->node_index_range())
155 862484 : face->set_node(n, this->node_ptr(Hex8::side_nodes_map[i][n]));
156 :
157 204467 : return face;
158 0 : }
159 :
160 :
161 :
162 29032425 : void Hex::side_ptr (std::unique_ptr<Elem> & side,
163 : const unsigned int i)
164 : {
165 29032425 : this->simple_side_ptr<Hex,Hex8>(side, i, QUAD4);
166 29032425 : }
167 :
168 :
169 :
170 20604866 : bool Hex::is_child_on_side(const unsigned int c,
171 : const unsigned int s) const
172 : {
173 1299400 : libmesh_assert_less (c, this->n_children());
174 1299400 : libmesh_assert_less (s, this->n_sides());
175 :
176 : // This array maps the Hex8 node numbering to the Hex8 child
177 : // numbering. I.e.
178 : // node 6 touches child 7, and
179 : // node 7 touches child 6, etc.
180 20604866 : const unsigned int node_child_map[8] = { 0, 1, 3, 2, 4, 5, 7, 6 };
181 :
182 76527173 : for (unsigned int i = 0; i != 4; ++i)
183 66517451 : if (node_child_map[Hex8::side_nodes_map[s][i]] == c)
184 592456 : return true;
185 :
186 706944 : return false;
187 : }
188 :
189 :
190 :
191 18644356 : bool Hex::is_edge_on_side(const unsigned int e,
192 : const unsigned int s) const
193 : {
194 10970484 : libmesh_assert_less (e, this->n_edges());
195 10970484 : libmesh_assert_less (s, this->n_sides());
196 :
197 18644356 : return (edge_sides_map[e][0] == s || edge_sides_map[e][1] == s);
198 : }
199 :
200 :
201 :
202 186624 : std::vector<unsigned int> Hex::sides_on_edge(const unsigned int e) const
203 : {
204 15552 : libmesh_assert_less(e, this->n_edges());
205 :
206 186624 : return {edge_sides_map[e][0], edge_sides_map[e][1]};
207 : }
208 :
209 :
210 :
211 0 : unsigned int Hex::opposite_side(const unsigned int side_in) const
212 : {
213 0 : libmesh_assert_less (side_in, 6);
214 : static const unsigned char hex_opposites[6] = {5, 3, 4, 1, 2, 0};
215 0 : return hex_opposites[side_in];
216 : }
217 :
218 :
219 :
220 0 : unsigned int Hex::opposite_node(const unsigned int node_in,
221 : const unsigned int side_in) const
222 : {
223 0 : libmesh_assert_less (node_in, 26);
224 0 : libmesh_assert_less (node_in, this->n_nodes());
225 0 : libmesh_assert_less (side_in, this->n_sides());
226 0 : libmesh_assert(this->is_node_on_side(node_in, side_in));
227 :
228 : static const unsigned char side05_nodes_map[] =
229 : {4, 5, 6, 7, 0, 1, 2, 3, 16, 17, 18, 19, 255, 255, 255, 255, 8, 9, 10, 11, 25, 255, 255, 255, 255, 20};
230 : static const unsigned char side13_nodes_map[] =
231 : {3, 2, 1, 0, 7, 6, 5, 4, 10, 255, 8, 255, 15, 14, 13, 12, 18, 255, 16, 255, 255, 23, 255, 21, 255, 255};
232 : static const unsigned char side24_nodes_map[] =
233 : {1, 0, 3, 2, 5, 4, 7, 6, 255, 11, 255, 9, 13, 12, 15, 14, 255, 19, 255, 17, 255, 255, 24, 255, 22, 255};
234 :
235 0 : switch (side_in)
236 : {
237 0 : case 0:
238 : case 5:
239 0 : return side05_nodes_map[node_in];
240 0 : case 1:
241 : case 3:
242 0 : return side13_nodes_map[node_in];
243 0 : case 2:
244 : case 4:
245 0 : return side24_nodes_map[node_in];
246 0 : default:
247 0 : libmesh_error_msg("Unsupported side_in = " << side_in);
248 : }
249 : }
250 :
251 :
252 :
253 : bool
254 2760 : Hex::is_flipped() const
255 : {
256 2592 : return (triple_product(this->point(1)-this->point(0),
257 2592 : this->point(3)-this->point(0),
258 3096 : this->point(4)-this->point(0)) < 0);
259 : }
260 :
261 :
262 : std::vector<unsigned int>
263 837080 : Hex::edges_adjacent_to_node(const unsigned int n) const
264 : {
265 69600 : libmesh_assert_less(n, this->n_nodes());
266 :
267 : // For vertices, we use the Hex::adjacent_edges_map, otherwise each
268 : // of the mid-edge nodes is adjacent only to the edge it is on, and
269 : // face/internal nodes are not adjacent to any edge.
270 837080 : if (this->is_vertex(n))
271 397288 : return {std::begin(adjacent_edges_map[n]), std::end(adjacent_edges_map[n])};
272 470208 : else if (this->is_edge(n))
273 364032 : return {n - this->n_vertices()};
274 :
275 8848 : libmesh_assert(this->is_face(n) || this->is_internal(n));
276 97328 : return {};
277 : }
278 :
279 :
280 128440 : Real Hex::quality (const ElemQuality q) const
281 : {
282 128440 : switch (q)
283 : {
284 :
285 : #if LIBMESH_DIM >= 3
286 : /**
287 : * Compute the min/max diagonal ratio.
288 : * Source: CUBIT User's Manual.
289 : */
290 7488 : case DIAGONAL:
291 : {
292 : // Diagonal between node 0 and node 6
293 7488 : const Real d06 = this->length(0,6);
294 :
295 : // Diagonal between node 3 and node 5
296 7488 : const Real d35 = this->length(3,5);
297 :
298 : // Diagonal between node 1 and node 7
299 7488 : const Real d17 = this->length(1,7);
300 :
301 : // Diagonal between node 2 and node 4
302 7488 : const Real d24 = this->length(2,4);
303 :
304 : // Find the biggest and smallest diagonals
305 7488 : const Real min = std::min(d06, std::min(d35, std::min(d17, d24)));
306 7488 : const Real max = std::max(d06, std::max(d35, std::max(d17, d24)));
307 :
308 624 : libmesh_assert_not_equal_to (max, 0.0);
309 :
310 7488 : return min / max;
311 :
312 : break;
313 : }
314 :
315 : /**
316 : * Minimum ratio of lengths derived from opposite edges.
317 : * Source: CUBIT User's Manual.
318 : */
319 7488 : case TAPER:
320 : {
321 :
322 : /**
323 : * Compute the side lengths.
324 : */
325 7488 : const Real d01 = this->length(0,1);
326 7488 : const Real d12 = this->length(1,2);
327 7488 : const Real d23 = this->length(2,3);
328 7488 : const Real d03 = this->length(0,3);
329 7488 : const Real d45 = this->length(4,5);
330 7488 : const Real d56 = this->length(5,6);
331 7488 : const Real d67 = this->length(6,7);
332 7488 : const Real d47 = this->length(4,7);
333 7488 : const Real d04 = this->length(0,4);
334 7488 : const Real d15 = this->length(1,5);
335 7488 : const Real d37 = this->length(3,7);
336 7488 : const Real d26 = this->length(2,6);
337 :
338 7488 : std::vector<Real> edge_ratios(12);
339 : // Front
340 7488 : edge_ratios[0] = std::min(d01, d45) / std::max(d01, d45);
341 7488 : edge_ratios[1] = std::min(d04, d15) / std::max(d04, d15);
342 :
343 : // Right
344 7488 : edge_ratios[2] = std::min(d15, d26) / std::max(d15, d26);
345 7488 : edge_ratios[3] = std::min(d12, d56) / std::max(d12, d56);
346 :
347 : // Back
348 7488 : edge_ratios[4] = std::min(d67, d23) / std::max(d67, d23);
349 7488 : edge_ratios[5] = std::min(d26, d37) / std::max(d26, d37);
350 :
351 : // Left
352 7488 : edge_ratios[6] = std::min(d04, d37) / std::max(d04, d37);
353 7488 : edge_ratios[7] = std::min(d03, d47) / std::max(d03, d47);
354 :
355 : // Bottom
356 7488 : edge_ratios[8] = std::min(d01, d23) / std::max(d01, d23);
357 7488 : edge_ratios[9] = std::min(d03, d12) / std::max(d03, d12);
358 :
359 : // Top
360 7488 : edge_ratios[10] = std::min(d45, d67) / std::max(d45, d67);
361 7488 : edge_ratios[11] = std::min(d56, d47) / std::max(d56, d47);
362 :
363 7488 : return *(std::min_element(edge_ratios.begin(), edge_ratios.end())) ;
364 :
365 : break;
366 : }
367 :
368 :
369 : /**
370 : * Minimum edge length divided by max diagonal length.
371 : * Source: CUBIT User's Manual.
372 : */
373 7488 : case STRETCH:
374 : {
375 624 : const Real sqrt3 = 1.73205080756888;
376 :
377 : /**
378 : * Compute the maximum diagonal.
379 : */
380 7488 : const Real d06 = this->length(0,6);
381 7488 : const Real d17 = this->length(1,7);
382 7488 : const Real d35 = this->length(3,5);
383 7488 : const Real d24 = this->length(2,4);
384 7488 : const Real max_diag = std::max(d06, std::max(d17, std::max(d35, d24)));
385 :
386 624 : libmesh_assert_not_equal_to ( max_diag, 0.0 );
387 :
388 : /**
389 : * Compute the minimum edge length.
390 : */
391 7488 : std::vector<Real> edges(12);
392 7488 : edges[0] = this->length(0,1);
393 7488 : edges[1] = this->length(1,2);
394 7488 : edges[2] = this->length(2,3);
395 7488 : edges[3] = this->length(0,3);
396 7488 : edges[4] = this->length(4,5);
397 7488 : edges[5] = this->length(5,6);
398 7488 : edges[6] = this->length(6,7);
399 7488 : edges[7] = this->length(4,7);
400 7488 : edges[8] = this->length(0,4);
401 7488 : edges[9] = this->length(1,5);
402 7488 : edges[10] = this->length(2,6);
403 7488 : edges[11] = this->length(3,7);
404 :
405 7488 : const Real min_edge = *(std::min_element(edges.begin(), edges.end()));
406 7488 : return sqrt3 * min_edge / max_diag ;
407 : }
408 :
409 :
410 14976 : case SHAPE:
411 : case SKEW:
412 : {
413 : // From: P. Knupp, "Algebraic mesh quality metrics for
414 : // unstructured initial meshes," Finite Elements in Analysis
415 : // and Design 39, 2003, p. 217-241, Sections 6.2 and 6.3.
416 :
417 : // Make local copies of points, we will access these several
418 : // times below.
419 : const Point
420 19968 : x0 = point(0), x1 = point(1), x2 = point(2), x3 = point(3),
421 18720 : x4 = point(4), x5 = point(5), x6 = point(6), x7 = point(7);
422 :
423 : // The columns of the Jacobian matrices are:
424 : // \vec{x}_{\xi} = \vec{a1}*eta*zeta + \vec{b1}*eta + \vec{c1}*zeta + \vec{d1}
425 : // \vec{x}_{\eta} = \vec{a2}*xi*zeta + \vec{b2}*xi + \vec{c2}*zeta + \vec{d2}
426 : // \vec{x}_{\zeta} = \vec{a3}*xi*eta + \vec{b3}*xi + \vec{c3}*eta + \vec{d3}
427 : // where the ai, bi, ci, and di are constants defined below.
428 1248 : const Point a1 = -x0 + x1 - x2 + x3 + x4 - x5 + x6 - x7;
429 1248 : const Point b1 = x0 - x1 + x2 - x3 + x4 - x5 + x6 - x7;
430 1248 : const Point c1 = x0 - x1 - x2 + x3 - x4 + x5 + x6 - x7;
431 1248 : const Point d1 = -x0 + x1 + x2 - x3 - x4 + x5 + x6 - x7;
432 :
433 1248 : const Point a2 = a1;
434 1248 : const Point b2 = b1;
435 1248 : const Point c2 = x0 + x1 - x2 - x3 - x4 - x5 + x6 + x7;
436 1248 : const Point d2 = -x0 - x1 + x2 + x3 - x4 - x5 + x6 + x7;
437 :
438 1248 : const Point a3 = a1;
439 1248 : const Point b3 = c1;
440 1248 : const Point c3 = c2;
441 1248 : const Point d3 = -x0 - x1 - x2 - x3 + x4 + x5 + x6 + x7;
442 :
443 : // Form the nodal Jacobians. These were computed using a
444 : // Python script and the formulas above. Note that we are
445 : // actually computing the Jacobian _columns_ and passing them
446 : // to the RealTensor constructor which expects _rows_, but
447 : // it's OK because we are only interested in determinants and
448 : // products which are not affected by taking the transpose.
449 : std::array<RealTensor, 8> A =
450 : {{
451 : RealTensor(d1, d2, d3),
452 1248 : RealTensor(d1, b2 + d2, b3 + d3),
453 1248 : RealTensor(b1 + d1, b2 + d2, a3 + b3 + c3 + d3),
454 1248 : RealTensor(b1 + d1, d2, c3 + d3),
455 1248 : RealTensor(c1 + d1, c2 + d2, d3),
456 1248 : RealTensor(c1 + d1, a2 + b2 + c2 + d2, b3 + d3),
457 1248 : RealTensor(a1 + b1 + c1 + d1, a2 + b2 + c2 + d2, a3 + b3 + c3 + d3),
458 1248 : RealTensor(a1 + b1 + c1 + d1, c2 + d2, c3 + d3)
459 11232 : }};
460 :
461 : // Compute Nodal areas, alpha_k = det(A_k).
462 : // If any of these are zero or negative, we return zero
463 : // (lowest possible value) for the quality, since the formulas
464 : // below require positive nodal areas.
465 : std::array<Real, 8> alpha;
466 134784 : for (unsigned int k=0; k<alpha.size(); ++k)
467 : {
468 119808 : alpha[k] = A[k].det();
469 119808 : if (alpha[k] <= 0.)
470 0 : return 0.;
471 : }
472 :
473 : // Compute metric tensors, T_k = A_k^T * A_k.
474 14976 : std::array<RealTensor, 8> T;
475 134784 : for (unsigned int k=0; k<T.size(); ++k)
476 129792 : T[k] = A[k] * A[k].transpose();
477 :
478 : // Compute and return the shape metric. These only use the
479 : // diagonal entries of the T_k.
480 1248 : Real den = 0.;
481 14976 : if (q == SHAPE)
482 : {
483 67392 : for (unsigned int k=0; k<T.size(); ++k)
484 64896 : den += T[k].tr() / std::pow(alpha[k], 2./3.);
485 7488 : return (den == 0.) ? 0 : (24. / den);
486 : }
487 : else
488 : {
489 67392 : for (unsigned int k=0; k<T.size(); ++k)
490 59904 : den += std::pow(std::sqrt(T[k](0,0) * T[k](1,1) * T[k](2,2)) / alpha[k], 2./3.);
491 7488 : return (den == 0.) ? 0 : (8. / den);
492 : }
493 : }
494 :
495 : // Verdict/CUBIT "skew" metric: the maximum |cos A| over the
496 : // three pairs of principal axes, where A is the angle between a
497 : // pair of axes. Each principal axis is the sum of the vectors
498 : // connecting the midpoints of opposite faces along one logical
499 : // direction. A value of 0 indicates a perfectly orthogonal
500 : // (unskewed) element; larger values (up to 1) indicate
501 : // increasing skew. This differs from the SKEW metric above,
502 : // which is Knupp's algebraic skew (1 is ideal).
503 : // See: C. J. Stimpson et al., "The Verdict Geometric Quality
504 : // Library," Sandia report SAND2007-1751, 2007.
505 7701 : case SKEW_ANGLE:
506 : {
507 : const Point
508 10221 : x0 = point(0), x1 = point(1), x2 = point(2), x3 = point(3),
509 9591 : x4 = point(4), x5 = point(5), x6 = point(6), x7 = point(7);
510 :
511 : // Principal axes, one per logical (xi, eta, zeta) direction.
512 : const Point
513 630 : X1 = (x1 - x0) + (x2 - x3) + (x5 - x4) + (x6 - x7),
514 630 : X2 = (x3 - x0) + (x2 - x1) + (x7 - x4) + (x6 - x5),
515 630 : X3 = (x4 - x0) + (x5 - x1) + (x6 - x2) + (x7 - x3);
516 :
517 7071 : const Real n1 = X1.norm(), n2 = X2.norm(), n3 = X3.norm();
518 :
519 : // Degenerate element: return 0 (the Verdict convention) if any
520 : // principal axis has zero length.
521 7701 : if (n1 == 0. || n2 == 0. || n3 == 0.)
522 0 : return 0.;
523 :
524 : // Normalize, then take the largest |cos| among the three
525 : // pairs of principal axes.
526 630 : const Point X1h = X1 / n1, X2h = X2 / n2, X3h = X3 / n3;
527 10851 : return std::max({std::abs(X1h * X2h),
528 630 : std::abs(X1h * X3h),
529 1260 : std::abs(X2h * X3h)});
530 : }
531 : #endif // LIBMESH_DIM >= 3
532 :
533 : /**
534 : * I don't know what to do for this metric.
535 : * Maybe the base class knows...
536 : */
537 83299 : default:
538 83299 : return Elem::quality(q);
539 : }
540 : }
541 :
542 :
543 :
544 4896 : std::pair<Real, Real> Hex::qual_bounds (const ElemQuality q) const
545 : {
546 408 : std::pair<Real, Real> bounds;
547 :
548 408 : switch (q)
549 : {
550 :
551 24 : case ASPECT_RATIO:
552 24 : bounds.first = 1.;
553 24 : bounds.second = 4.;
554 24 : break;
555 :
556 24 : case SKEW_ANGLE:
557 24 : bounds.first = 0.;
558 24 : bounds.second = 0.5;
559 24 : break;
560 :
561 72 : case SKEW:
562 : case SHEAR:
563 : case SHAPE:
564 72 : bounds.first = 0.3;
565 72 : bounds.second = 1.;
566 72 : break;
567 :
568 24 : case CONDITION:
569 24 : bounds.first = 1.;
570 24 : bounds.second = 8.;
571 24 : break;
572 :
573 24 : case SCALED_JACOBIAN:
574 24 : bounds.first = 0.5;
575 24 : bounds.second = 1.;
576 24 : break;
577 :
578 24 : case DISTORTION:
579 24 : bounds.first = 0.6;
580 24 : bounds.second = 1.;
581 24 : break;
582 :
583 24 : case TAPER:
584 : // TAPER is 1 for an untapered element and decreases toward 0 with
585 : // increasing taper (see Hex::quality), so the good range runs up
586 : // to 1, not down from 0.
587 24 : bounds.first = 0.4;
588 24 : bounds.second = 1.;
589 24 : break;
590 :
591 24 : case STRETCH:
592 24 : bounds.first = 0.25;
593 24 : bounds.second = 1.;
594 24 : break;
595 :
596 24 : case DIAGONAL:
597 24 : bounds.first = 0.65;
598 24 : bounds.second = 1.;
599 24 : break;
600 :
601 24 : case SIZE:
602 24 : bounds.first = 0.5;
603 24 : bounds.second = 1.;
604 24 : break;
605 :
606 120 : default:
607 : libmesh_warning("Invalid quality measure chosen.");
608 120 : bounds.first = -1;
609 120 : bounds.second = -1;
610 : }
611 :
612 4896 : return bounds;
613 : }
614 :
615 :
616 :
617 13912996 : bool Hex::on_reference_element(const Point & p,
618 : const Real eps) const
619 : {
620 1518015 : const Real & xi = p(0);
621 1518015 : const Real & eta = p(1);
622 1518015 : const Real & zeta = p(2);
623 :
624 : // The reference hexahedral element is [-1,1]^3.
625 25930329 : return ((xi >= -1.-eps) &&
626 12017333 : (xi <= 1.+eps) &&
627 9281512 : (eta >= -1.-eps) &&
628 8052649 : (eta <= 1.+eps) &&
629 23368185 : (zeta >= -1.-eps) &&
630 15150120 : (zeta <= 1.+eps));
631 : }
632 :
633 :
634 :
635 : const unsigned short int Hex::_second_order_vertex_child_number[27] =
636 : {
637 : 99,99,99,99,99,99,99,99, // Vertices
638 : 0,1,2,0,0,1,2,3,4,5,6,5, // Edges
639 : 0,0,1,2,0,4, // Faces
640 : 0 // Interior
641 : };
642 :
643 :
644 :
645 : const unsigned short int Hex::_second_order_vertex_child_index[27] =
646 : {
647 : 99,99,99,99,99,99,99,99, // Vertices
648 : 1,2,3,3,4,5,6,7,5,6,7,7, // Edges
649 : 2,5,6,7,7,6, // Faces
650 : 6 // Interior
651 : };
652 :
653 :
654 : const unsigned short int Hex::_second_order_adjacent_vertices[12][2] =
655 : {
656 : { 0, 1}, // vertices adjacent to node 8
657 : { 1, 2}, // vertices adjacent to node 9
658 : { 2, 3}, // vertices adjacent to node 10
659 : { 0, 3}, // vertices adjacent to node 11
660 :
661 : { 0, 4}, // vertices adjacent to node 12
662 : { 1, 5}, // vertices adjacent to node 13
663 : { 2, 6}, // vertices adjacent to node 14
664 : { 3, 7}, // vertices adjacent to node 15
665 :
666 : { 4, 5}, // vertices adjacent to node 16
667 : { 5, 6}, // vertices adjacent to node 17
668 : { 6, 7}, // vertices adjacent to node 18
669 : { 4, 7} // vertices adjacent to node 19
670 : };
671 :
672 :
673 : #ifdef LIBMESH_ENABLE_AMR
674 :
675 : // We number 125 "possible node locations" for a 2x2x2 refinement of
676 : // hexes with up to 3x3x3 nodes each
677 : const int Hex::_child_node_lookup[8][27] =
678 : {
679 : // node lookup for child 0 (near node 0)
680 : { 0, 2, 12, 10, 50, 52, 62, 60, 1, 7, 11, 5, 25, 27, 37, 35,
681 : 51, 57, 61, 55, 6, 26, 32, 36, 30, 56, 31},
682 :
683 : // node lookup for child 1 (near node 1)
684 : { 2, 4, 14, 12, 52, 54, 64, 62, 3, 9, 13, 7, 27, 29, 39, 37,
685 : 53, 59, 63, 57, 8, 28, 34, 38, 32, 58, 33},
686 :
687 : // node lookup for child 2 (near node 3)
688 : { 10, 12, 22, 20, 60, 62, 72, 70, 11, 17, 21, 15, 35, 37, 47, 45,
689 : 61, 67, 71, 65, 16, 36, 42, 46, 40, 66, 41},
690 :
691 : // node lookup for child 3 (near node 2)
692 : { 12, 14, 24, 22, 62, 64, 74, 72, 13, 19, 23, 17, 37, 39, 49, 47,
693 : 63, 69, 73, 67, 18, 38, 44, 48, 42, 68, 43},
694 :
695 : // node lookup for child 4 (near node 4)
696 : { 50, 52, 62, 60, 100, 102, 112, 110, 51, 57, 61, 55, 75, 77, 87, 85,
697 : 101, 107, 111, 105, 56, 76, 82, 86, 80, 106, 81},
698 :
699 : // node lookup for child 5 (near node 5)
700 : { 52, 54, 64, 62, 102, 104, 114, 112, 53, 59, 63, 57, 77, 79, 89, 87,
701 : 103, 109, 113, 107, 58, 78, 84, 88, 82, 108, 93},
702 :
703 : // node lookup for child 6 (near node 7)
704 : { 60, 62, 72, 70, 110, 112, 122, 120, 61, 67, 71, 65, 85, 87, 97, 95,
705 : 111, 117, 121, 115, 66, 86, 92, 96, 90, 116, 91},
706 :
707 : // node lookup for child 7 (near node 6)
708 : { 62, 64, 74, 72, 112, 114, 124, 122, 63, 69, 73, 67, 87, 89, 99, 97,
709 : 113, 119, 123, 117, 68, 88, 94, 98, 92, 118, 103}
710 : };
711 :
712 : #endif // LIBMESH_ENABLE_AMR
713 :
714 :
715 : } // namespace libMesh
|
