libMesh
cell_hex.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 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 
27 
28 
29 namespace libMesh
30 {
31 
32 
33 
34 // ------------------------------------------------------------
35 // Hex class static member initializations
36 
37 
38 // We need to require C++11...
39 const Real Hex::_master_points[27][3] =
40  {
41  {-1, -1, -1},
42  {1, -1, -1},
43  {1, 1, -1},
44  {-1, 1, -1},
45  {-1, -1, 1},
46  {1, -1, 1},
47  {1, 1, 1},
48  {-1, 1, 1},
49  {0, -1, -1},
50  {1, 0, -1},
51  {0, 1, -1},
52  {-1, 0, -1},
53  {-1, -1, 0},
54  {1, -1, 0},
55  {1, 1, 0},
56  {-1, 1, 0},
57  {0, -1, 1},
58  {1, 0, 1},
59  {0, 1, 1},
60  {-1, 0, 1},
61  {0, 0, -1},
62  {0, -1, 0},
63  {1, 0, 0},
64  {0, 1, 0},
65  {-1, 0, 0},
66  {0, 0, 1}
67  };
68 
69 
70 
71 
72 // ------------------------------------------------------------
73 // Hex class member functions
74 dof_id_type Hex::key (const unsigned int s) const
75 {
76  libmesh_assert_less (s, this->n_sides());
77 
78  return this->compute_key(this->node_id(Hex8::side_nodes_map[s][0]),
79  this->node_id(Hex8::side_nodes_map[s][1]),
80  this->node_id(Hex8::side_nodes_map[s][2]),
81  this->node_id(Hex8::side_nodes_map[s][3]));
82 }
83 
84 
85 
86 unsigned int Hex::which_node_am_i(unsigned int side,
87  unsigned int side_node) const
88 {
89  libmesh_assert_less (side, this->n_sides());
90  libmesh_assert_less (side_node, 4);
91 
92  return Hex8::side_nodes_map[side][side_node];
93 }
94 
95 
96 
97 UniquePtr<Elem> Hex::side_ptr (const unsigned int i)
98 {
99  libmesh_assert_less (i, this->n_sides());
100 
101  Elem * face = new Quad4;
102 
103  for (unsigned n=0; n<face->n_nodes(); ++n)
104  face->set_node(n) = this->node_ptr(Hex8::side_nodes_map[i][n]);
105 
106  return UniquePtr<Elem>(face);
107 }
108 
109 
110 
111 bool Hex::is_child_on_side(const unsigned int c,
112  const unsigned int s) const
113 {
114  libmesh_assert_less (c, this->n_children());
115  libmesh_assert_less (s, this->n_sides());
116 
117  // This array maps the Hex8 node numbering to the Hex8 child
118  // numbering. I.e.
119  // node 6 touches child 7, and
120  // node 7 touches child 6, etc.
121  const unsigned int node_child_map[8] = { 0, 1, 3, 2, 4, 5, 7, 6 };
122 
123  for (unsigned int i = 0; i != 4; ++i)
124  if (node_child_map[Hex8::side_nodes_map[s][i]] == c)
125  return true;
126 
127  return false;
128 }
129 
130 
131 
132 bool Hex::is_edge_on_side(const unsigned int e,
133  const unsigned int s) const
134 {
135  libmesh_assert_less (e, this->n_edges());
136  libmesh_assert_less (s, this->n_sides());
137 
138  return (is_node_on_side(Hex8::edge_nodes_map[e][0],s) &&
139  is_node_on_side(Hex8::edge_nodes_map[e][1],s));
140 }
141 
142 
143 
144 unsigned int Hex::opposite_side(const unsigned int side_in) const
145 {
146  libmesh_assert_less (side_in, 6);
147  static const unsigned char hex_opposites[6] = {5, 3, 4, 1, 2, 0};
148  return hex_opposites[side_in];
149 }
150 
151 
152 
153 unsigned int Hex::opposite_node(const unsigned int node_in,
154  const unsigned int side_in) const
155 {
156  libmesh_assert_less (node_in, 26);
157  libmesh_assert_less (node_in, this->n_nodes());
158  libmesh_assert_less (side_in, this->n_sides());
159  libmesh_assert(this->is_node_on_side(node_in, side_in));
160 
161  static const unsigned char side05_nodes_map[] =
162  {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};
163  static const unsigned char side13_nodes_map[] =
164  {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};
165  static const unsigned char side24_nodes_map[] =
166  {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};
167 
168  switch (side_in)
169  {
170  case 0:
171  case 5:
172  return side05_nodes_map[node_in];
173  case 1:
174  case 3:
175  return side13_nodes_map[node_in];
176  case 2:
177  case 4:
178  return side24_nodes_map[node_in];
179  default:
180  libmesh_error_msg("Unsupported side_in = " << side_in);
181  }
182 
183  libmesh_error_msg("We'll never get here!");
184  return 255;
185 }
186 
187 
188 
190 {
191  switch (q)
192  {
193 
198  case DIAGONAL:
199  {
200  // Diagonal between node 0 and node 6
201  const Real d06 = this->length(0,6);
202 
203  // Diagonal between node 3 and node 5
204  const Real d35 = this->length(3,5);
205 
206  // Diagonal between node 1 and node 7
207  const Real d17 = this->length(1,7);
208 
209  // Diagonal between node 2 and node 4
210  const Real d24 = this->length(2,4);
211 
212  // Find the biggest and smallest diagonals
213  const Real min = std::min(d06, std::min(d35, std::min(d17, d24)));
214  const Real max = std::max(d06, std::max(d35, std::max(d17, d24)));
215 
216  libmesh_assert_not_equal_to (max, 0.0);
217 
218  return min / max;
219 
220  break;
221  }
222 
227  case TAPER:
228  {
229 
233  const Real d01 = this->length(0,1);
234  const Real d12 = this->length(1,2);
235  const Real d23 = this->length(2,3);
236  const Real d03 = this->length(0,3);
237  const Real d45 = this->length(4,5);
238  const Real d56 = this->length(5,6);
239  const Real d67 = this->length(6,7);
240  const Real d47 = this->length(4,7);
241  const Real d04 = this->length(0,4);
242  const Real d15 = this->length(1,5);
243  const Real d37 = this->length(3,7);
244  const Real d26 = this->length(2,6);
245 
246  std::vector<Real> edge_ratios(12);
247  // Front
248  edge_ratios[0] = std::min(d01, d45) / std::max(d01, d45);
249  edge_ratios[1] = std::min(d04, d15) / std::max(d04, d15);
250 
251  // Right
252  edge_ratios[2] = std::min(d15, d26) / std::max(d15, d26);
253  edge_ratios[3] = std::min(d12, d56) / std::max(d12, d56);
254 
255  // Back
256  edge_ratios[4] = std::min(d67, d23) / std::max(d67, d23);
257  edge_ratios[5] = std::min(d26, d37) / std::max(d26, d37);
258 
259  // Left
260  edge_ratios[6] = std::min(d04, d37) / std::max(d04, d37);
261  edge_ratios[7] = std::min(d03, d47) / std::max(d03, d47);
262 
263  // Bottom
264  edge_ratios[8] = std::min(d01, d23) / std::max(d01, d23);
265  edge_ratios[9] = std::min(d03, d12) / std::max(d03, d12);
266 
267  // Top
268  edge_ratios[10] = std::min(d45, d67) / std::max(d45, d67);
269  edge_ratios[11] = std::min(d56, d47) / std::max(d56, d47);
270 
271  return *(std::min_element(edge_ratios.begin(), edge_ratios.end())) ;
272 
273  break;
274  }
275 
276 
281  case STRETCH:
282  {
283  const Real sqrt3 = 1.73205080756888;
284 
288  const Real d06 = this->length(0,6);
289  const Real d17 = this->length(1,7);
290  const Real d35 = this->length(3,5);
291  const Real d24 = this->length(2,4);
292  const Real max_diag = std::max(d06, std::max(d17, std::max(d35, d24)));
293 
294  libmesh_assert_not_equal_to ( max_diag, 0.0 );
295 
299  std::vector<Real> edges(12);
300  edges[0] = this->length(0,1);
301  edges[1] = this->length(1,2);
302  edges[2] = this->length(2,3);
303  edges[3] = this->length(0,3);
304  edges[4] = this->length(4,5);
305  edges[5] = this->length(5,6);
306  edges[6] = this->length(6,7);
307  edges[7] = this->length(4,7);
308  edges[8] = this->length(0,4);
309  edges[9] = this->length(1,5);
310  edges[10] = this->length(2,6);
311  edges[11] = this->length(3,7);
312 
313  const Real min_edge = *(std::min_element(edges.begin(), edges.end()));
314  return sqrt3 * min_edge / max_diag ;
315  }
316 
317 
322  default:
323  return Elem::quality(q);
324  }
325 
326  libmesh_error_msg("We'll never get here!");
327  return 0.;
328 }
329 
330 
331 
332 std::pair<Real, Real> Hex::qual_bounds (const ElemQuality q) const
333 {
334  std::pair<Real, Real> bounds;
335 
336  switch (q)
337  {
338 
339  case ASPECT_RATIO:
340  bounds.first = 1.;
341  bounds.second = 4.;
342  break;
343 
344  case SKEW:
345  bounds.first = 0.;
346  bounds.second = 0.5;
347  break;
348 
349  case SHEAR:
350  case SHAPE:
351  bounds.first = 0.3;
352  bounds.second = 1.;
353  break;
354 
355  case CONDITION:
356  bounds.first = 1.;
357  bounds.second = 8.;
358  break;
359 
360  case JACOBIAN:
361  bounds.first = 0.5;
362  bounds.second = 1.;
363  break;
364 
365  case DISTORTION:
366  bounds.first = 0.6;
367  bounds.second = 1.;
368  break;
369 
370  case TAPER:
371  bounds.first = 0.;
372  bounds.second = 0.4;
373  break;
374 
375  case STRETCH:
376  bounds.first = 0.25;
377  bounds.second = 1.;
378  break;
379 
380  case DIAGONAL:
381  bounds.first = 0.65;
382  bounds.second = 1.;
383  break;
384 
385  case SIZE:
386  bounds.first = 0.5;
387  bounds.second = 1.;
388  break;
389 
390  default:
391  libMesh::out << "Warning: Invalid quality measure chosen." << std::endl;
392  bounds.first = -1;
393  bounds.second = -1;
394  }
395 
396  return bounds;
397 }
398 
399 
400 
401 const unsigned short int Hex::_second_order_vertex_child_number[27] =
402  {
403  99,99,99,99,99,99,99,99, // Vertices
404  0,1,2,0,0,1,2,3,4,5,6,5, // Edges
405  0,0,1,2,0,4, // Faces
406  0 // Interior
407  };
408 
409 
410 
411 const unsigned short int Hex::_second_order_vertex_child_index[27] =
412  {
413  99,99,99,99,99,99,99,99, // Vertices
414  1,2,3,3,4,5,6,7,5,6,7,7, // Edges
415  2,5,6,7,7,6, // Faces
416  6 // Interior
417  };
418 
419 
420 const unsigned short int Hex::_second_order_adjacent_vertices[12][2] =
421  {
422  { 0, 1}, // vertices adjacent to node 8
423  { 1, 2}, // vertices adjacent to node 9
424  { 2, 3}, // vertices adjacent to node 10
425  { 0, 3}, // vertices adjacent to node 11
426 
427  { 0, 4}, // vertices adjacent to node 12
428  { 1, 5}, // vertices adjacent to node 13
429  { 2, 6}, // vertices adjacent to node 14
430  { 3, 7}, // vertices adjacent to node 15
431 
432  { 4, 5}, // vertices adjacent to node 16
433  { 5, 6}, // vertices adjacent to node 17
434  { 6, 7}, // vertices adjacent to node 18
435  { 4, 7} // vertices adjacent to node 19
436  };
437 
438 
439 #ifdef LIBMESH_ENABLE_AMR
440 
441 // We number 125 "possible node locations" for a 2x2x2 refinement of
442 // hexes with up to 3x3x3 nodes each
443 const int Hex::_child_node_lookup[8][27] =
444  {
445  // node lookup for child 0 (near node 0)
446  { 0, 2, 12, 10, 50, 52, 62, 60, 1, 7, 11, 5, 25, 27, 37, 35,
447  51, 57, 61, 55, 6, 26, 32, 36, 30, 56, 31},
448 
449  // node lookup for child 1 (near node 1)
450  { 2, 4, 14, 12, 52, 54, 64, 62, 3, 9, 13, 7, 27, 29, 39, 37,
451  53, 59, 63, 57, 8, 28, 34, 38, 32, 58, 33},
452 
453  // node lookup for child 2 (near node 3)
454  { 10, 12, 22, 20, 60, 62, 72, 70, 11, 17, 21, 15, 35, 37, 47, 45,
455  61, 67, 71, 65, 16, 36, 42, 46, 40, 66, 41},
456 
457  // node lookup for child 3 (near node 2)
458  { 12, 14, 24, 22, 62, 64, 74, 72, 13, 19, 23, 17, 37, 39, 49, 47,
459  63, 69, 73, 67, 18, 38, 44, 48, 42, 68, 43},
460 
461  // node lookup for child 4 (near node 4)
462  { 50, 52, 62, 60, 100, 102, 112, 110, 51, 57, 61, 55, 75, 77, 87, 85,
463  101, 107, 111, 105, 56, 76, 82, 86, 80, 106, 81},
464 
465  // node lookup for child 5 (near node 5)
466  { 52, 54, 64, 62, 102, 104, 114, 112, 53, 59, 63, 57, 77, 79, 89, 87,
467  103, 109, 113, 107, 58, 78, 84, 88, 82, 108, 93},
468 
469  // node lookup for child 6 (near node 7)
470  { 60, 62, 72, 70, 110, 112, 122, 120, 61, 67, 71, 65, 85, 87, 97, 95,
471  111, 117, 121, 115, 66, 86, 92, 96, 90, 116, 91},
472 
473  // node lookup for child 7 (near node 6)
474  { 62, 64, 74, 72, 112, 114, 124, 122, 63, 69, 73, 67, 87, 89, 99, 97,
475  113, 119, 123, 117, 68, 88, 94, 98, 92, 118, 103}
476  };
477 
478 #endif // LIBMESH_ENABLE_AMR
479 
480 
481 } // namespace libMesh
virtual bool is_node_on_side(const unsigned int n, const unsigned int s) const =0
virtual Node *& set_node(const unsigned int i)
Definition: elem.h:1941
static const Real _master_points[27][3]
Master element node locations.
Definition: cell_hex.h:182
The QUAD4 is an element in 2D composed of 4 nodes.
Definition: face_quad4.h:49
UniquePtr< Elem > side(const unsigned int i) const
Definition: elem.h:2105
static const unsigned short int _second_order_adjacent_vertices[12][2]
Matrix that tells which vertices define the location of mid-side (or second-order) nodes...
Definition: cell_hex.h:167
This is the base class from which all geometric element types are derived.
Definition: elem.h:89
The libMesh namespace provides an interface to certain functionality in the library.
virtual unsigned int n_children() const libmesh_override
Definition: cell_hex.h:88
long double max(long double a, double b)
virtual std::pair< Real, Real > qual_bounds(const ElemQuality q) const libmesh_override
Definition: cell_hex.C:332
virtual unsigned int opposite_node(const unsigned int n, const unsigned int s) const libmesh_override
Definition: cell_hex.C:153
libmesh_assert(j)
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
virtual bool is_child_on_side(const unsigned int c, const unsigned int s) const libmesh_override
Definition: cell_hex.C:111
virtual unsigned int n_nodes() const =0
virtual dof_id_type key() const
Definition: elem.C:503
const Node * node_ptr(const unsigned int i) const
Definition: elem.h:1874
virtual Real quality(const ElemQuality q) const
Definition: elem.C:1547
static const int _child_node_lookup[8][27]
Lookup table from child id, child node id to "possible node location" (a simple dictionary-index in a...
Definition: cell_hex.h:188
virtual unsigned int which_node_am_i(unsigned int side, unsigned int side_node) const libmesh_override
Definition: cell_hex.C:86
virtual bool is_edge_on_side(const unsigned int e, const unsigned int s) const libmesh_override
Definition: cell_hex.C:132
virtual UniquePtr< Elem > side_ptr(const unsigned int i) libmesh_override
Definition: cell_hex.C:97
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual Real quality(const ElemQuality q) const libmesh_override
Definition: cell_hex.C:189
OStreamProxy out
static const unsigned short int _second_order_vertex_child_number[27]
Vector that names a child sharing each second order node.
Definition: cell_hex.h:172
Real length(const unsigned int n1, const unsigned int n2) const
Definition: elem.C:492
static const unsigned short int _second_order_vertex_child_index[27]
Vector that names the child vertex index for each second order node.
Definition: cell_hex.h:177
dof_id_type node_id(const unsigned int i) const
Definition: elem.h:1831
ElemQuality
Defines an enum for element quality metrics.
static dof_id_type compute_key(dof_id_type n0)
Definition: elem.h:2621
virtual unsigned int n_edges() const libmesh_override
Definition: cell_hex.h:78
long double min(long double a, double b)
virtual unsigned int n_sides() const libmesh_override
Definition: cell_hex.h:68
virtual unsigned int opposite_side(const unsigned int s) const libmesh_override
Definition: cell_hex.C:144
uint8_t dof_id_type
Definition: id_types.h:64