cell_hex.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2014 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 namespace libMesh
28 {
29 
30 
31 
32 
33 // ------------------------------------------------------------
34 // Hex class member functions
35 dof_id_type Hex::key (const unsigned int s) const
36 {
37  libmesh_assert_less (s, this->n_sides());
38 
39  // Think of a unit cube: (-1,1) x (-1,1)x (-1,1)
40  switch (s)
41  {
42  case 0: // the face at z = -1
43 
44  return
45  this->compute_key (this->node(0),
46  this->node(3),
47  this->node(2),
48  this->node(1));
49 
50  case 1: // the face at y = -1
51 
52  return
53  this->compute_key (this->node(0),
54  this->node(1),
55  this->node(5),
56  this->node(4));
57 
58  case 2: // the face at x = 1
59 
60  return
61  this->compute_key (this->node(1),
62  this->node(2),
63  this->node(6),
64  this->node(5));
65 
66  case 3: // the face at y = 1
67 
68  return
69  this->compute_key (this->node(2),
70  this->node(3),
71  this->node(7),
72  this->node(6));
73 
74  case 4: // the face at x = -1
75 
76  return
77  this->compute_key (this->node(3),
78  this->node(0),
79  this->node(4),
80  this->node(7));
81 
82  case 5: // the face at z = 1
83 
84  return
85  this->compute_key (this->node(4),
86  this->node(5),
87  this->node(6),
88  this->node(7));
89  }
90 
91  // We'll never get here.
92  libmesh_error();
93  return 0;
94 }
95 
96 
97 
98 AutoPtr<Elem> Hex::side (const unsigned int i) const
99 {
100  libmesh_assert_less (i, this->n_sides());
101 
102 
103 
104  Elem* face = new Quad4;
105 
106  // Think of a unit cube: (-1,1) x (-1,1)x (-1,1)
107  switch (i)
108  {
109  case 0: // the face at z = -1
110  {
111  face->set_node(0) = this->get_node(0);
112  face->set_node(1) = this->get_node(3);
113  face->set_node(2) = this->get_node(2);
114  face->set_node(3) = this->get_node(1);
115 
116  AutoPtr<Elem> ap(face);
117  return ap;
118  }
119  case 1: // the face at y = -1
120  {
121  face->set_node(0) = this->get_node(0);
122  face->set_node(1) = this->get_node(1);
123  face->set_node(2) = this->get_node(5);
124  face->set_node(3) = this->get_node(4);
125 
126  AutoPtr<Elem> ap(face);
127  return ap;
128  }
129  case 2: // the face at x = 1
130  {
131  face->set_node(0) = this->get_node(1);
132  face->set_node(1) = this->get_node(2);
133  face->set_node(2) = this->get_node(6);
134  face->set_node(3) = this->get_node(5);
135 
136  AutoPtr<Elem> ap(face);
137  return ap;
138  }
139  case 3: // the face at y = 1
140  {
141  face->set_node(0) = this->get_node(2);
142  face->set_node(1) = this->get_node(3);
143  face->set_node(2) = this->get_node(7);
144  face->set_node(3) = this->get_node(6);
145 
146  AutoPtr<Elem> ap(face);
147  return ap;
148  }
149  case 4: // the face at x = -1
150  {
151  face->set_node(0) = this->get_node(3);
152  face->set_node(1) = this->get_node(0);
153  face->set_node(2) = this->get_node(4);
154  face->set_node(3) = this->get_node(7);
155 
156  AutoPtr<Elem> ap(face);
157  return ap;
158  }
159  case 5: // the face at z = 1
160  {
161  face->set_node(0) = this->get_node(4);
162  face->set_node(1) = this->get_node(5);
163  face->set_node(2) = this->get_node(6);
164  face->set_node(3) = this->get_node(7);
165 
166  AutoPtr<Elem> ap(face);
167  return ap;
168  }
169  default:
170  {
171  libmesh_error();
172  AutoPtr<Elem> ap(face);
173  return ap;
174  }
175  }
176 
177  // We'll never get here.
178  libmesh_error();
179  AutoPtr<Elem> ap(face);
180  return ap;
181 }
182 
183 
184 
185 bool Hex::is_child_on_side(const unsigned int c,
186  const unsigned int s) const
187 {
188  libmesh_assert_less (c, this->n_children());
189  libmesh_assert_less (s, this->n_sides());
190 
191  // This array maps the Hex8 node numbering to the Hex8 child
192  // numbering. I.e.
193  // node 6 touches child 7, and
194  // node 7 touches child 6, etc.
195  const unsigned int node_child_map[8] = { 0, 1, 3, 2, 4, 5, 7, 6 };
196 
197  for (unsigned int i = 0; i != 4; ++i)
198  if (node_child_map[Hex8::side_nodes_map[s][i]] == c)
199  return true;
200 
201  return false;
202 }
203 
204 
205 
206 bool Hex::is_edge_on_side(const unsigned int e,
207  const unsigned int s) const
208 {
209  libmesh_assert_less (e, this->n_edges());
210  libmesh_assert_less (s, this->n_sides());
211 
212  return (is_node_on_side(Hex8::edge_nodes_map[e][0],s) &&
214 }
215 
216 
217 
218 unsigned int Hex::opposite_side(const unsigned int side_in) const
219 {
220  libmesh_assert_less (side_in, 6);
221  static const unsigned char hex_opposites[6] = {5, 3, 4, 1, 2, 0};
222  return hex_opposites[side_in];
223 }
224 
225 
226 
227 unsigned int Hex::opposite_node(const unsigned int node_in,
228  const unsigned int side_in) const
229 {
230  libmesh_assert_less (node_in, 26);
231  libmesh_assert_less (node_in, this->n_nodes());
232  libmesh_assert_less (side_in, this->n_sides());
233  libmesh_assert(this->is_node_on_side(node_in, side_in));
234 
235  static const unsigned char side05_nodes_map[] =
236  {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};
237  static const unsigned char side13_nodes_map[] =
238  {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};
239  static const unsigned char side24_nodes_map[] =
240  {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};
241 
242  switch (side_in)
243  {
244  case 0:
245  case 5:
246  return side05_nodes_map[node_in];
247  break;
248  case 1:
249  case 3:
250  return side13_nodes_map[node_in];
251  break;
252  case 2:
253  case 4:
254  return side24_nodes_map[node_in];
255  break;
256  }
257 
258  libmesh_error();
259  return 255;
260 }
261 
262 
263 
265 {
266  switch (q)
267  {
268 
273  case DIAGONAL:
274  {
275  // Diagonal between node 0 and node 6
276  const Real d06 = this->length(0,6);
277 
278  // Diagonal between node 3 and node 5
279  const Real d35 = this->length(3,5);
280 
281  // Diagonal between node 1 and node 7
282  const Real d17 = this->length(1,7);
283 
284  // Diagonal between node 2 and node 4
285  const Real d24 = this->length(2,4);
286 
287  // Find the biggest and smallest diagonals
288  const Real min = std::min(d06, std::min(d35, std::min(d17, d24)));
289  const Real max = std::max(d06, std::max(d35, std::max(d17, d24)));
290 
291  libmesh_assert_not_equal_to (max, 0.0);
292 
293  return min / max;
294 
295  break;
296  }
297 
302  case TAPER:
303  {
304 
308  const Real d01 = this->length(0,1);
309  const Real d12 = this->length(1,2);
310  const Real d23 = this->length(2,3);
311  const Real d03 = this->length(0,3);
312  const Real d45 = this->length(4,5);
313  const Real d56 = this->length(5,6);
314  const Real d67 = this->length(6,7);
315  const Real d47 = this->length(4,7);
316  const Real d04 = this->length(0,4);
317  const Real d15 = this->length(1,5);
318  const Real d37 = this->length(3,7);
319  const Real d26 = this->length(2,6);
320 
321  std::vector<Real> edge_ratios(12);
322  // Front
323  edge_ratios[0] = std::min(d01, d45) / std::max(d01, d45);
324  edge_ratios[1] = std::min(d04, d15) / std::max(d04, d15);
325 
326  // Right
327  edge_ratios[2] = std::min(d15, d26) / std::max(d15, d26);
328  edge_ratios[3] = std::min(d12, d56) / std::max(d12, d56);
329 
330  // Back
331  edge_ratios[4] = std::min(d67, d23) / std::max(d67, d23);
332  edge_ratios[5] = std::min(d26, d37) / std::max(d26, d37);
333 
334  // Left
335  edge_ratios[6] = std::min(d04, d37) / std::max(d04, d37);
336  edge_ratios[7] = std::min(d03, d47) / std::max(d03, d47);
337 
338  // Bottom
339  edge_ratios[8] = std::min(d01, d23) / std::max(d01, d23);
340  edge_ratios[9] = std::min(d03, d12) / std::max(d03, d12);
341 
342  // Top
343  edge_ratios[10] = std::min(d45, d67) / std::max(d45, d67);
344  edge_ratios[11] = std::min(d56, d47) / std::max(d56, d47);
345 
346  return *(std::min_element(edge_ratios.begin(), edge_ratios.end())) ;
347 
348  break;
349  }
350 
351 
356  case STRETCH:
357  {
358  const Real sqrt3 = 1.73205080756888;
359 
363  const Real d06 = this->length(0,6);
364  const Real d17 = this->length(1,7);
365  const Real d35 = this->length(3,5);
366  const Real d24 = this->length(2,4);
367  const Real max_diag = std::max(d06, std::max(d17, std::max(d35, d24)));
368 
369  libmesh_assert_not_equal_to ( max_diag, 0.0 );
370 
374  std::vector<Real> edges(12);
375  edges[0] = this->length(0,1);
376  edges[1] = this->length(1,2);
377  edges[2] = this->length(2,3);
378  edges[3] = this->length(0,3);
379  edges[4] = this->length(4,5);
380  edges[5] = this->length(5,6);
381  edges[6] = this->length(6,7);
382  edges[7] = this->length(4,7);
383  edges[8] = this->length(0,4);
384  edges[9] = this->length(1,5);
385  edges[10] = this->length(2,6);
386  edges[11] = this->length(3,7);
387 
388  const Real min_edge = *(std::min_element(edges.begin(), edges.end()));
389  return sqrt3 * min_edge / max_diag ;
390  break;
391  }
392 
393 
398  default:
399  {
400  return Elem::quality(q);
401  }
402  }
403 
404 
405  // Will never get here...
406  libmesh_error();
407  return 0.;
408 }
409 
410 
411 
412 std::pair<Real, Real> Hex::qual_bounds (const ElemQuality q) const
413 {
414  std::pair<Real, Real> bounds;
415 
416  switch (q)
417  {
418 
419  case ASPECT_RATIO:
420  bounds.first = 1.;
421  bounds.second = 4.;
422  break;
423 
424  case SKEW:
425  bounds.first = 0.;
426  bounds.second = 0.5;
427  break;
428 
429  case SHEAR:
430  case SHAPE:
431  bounds.first = 0.3;
432  bounds.second = 1.;
433  break;
434 
435  case CONDITION:
436  bounds.first = 1.;
437  bounds.second = 8.;
438  break;
439 
440  case JACOBIAN:
441  bounds.first = 0.5;
442  bounds.second = 1.;
443  break;
444 
445  case DISTORTION:
446  bounds.first = 0.6;
447  bounds.second = 1.;
448  break;
449 
450  case TAPER:
451  bounds.first = 0.;
452  bounds.second = 0.4;
453  break;
454 
455  case STRETCH:
456  bounds.first = 0.25;
457  bounds.second = 1.;
458  break;
459 
460  case DIAGONAL:
461  bounds.first = 0.65;
462  bounds.second = 1.;
463  break;
464 
465  case SIZE:
466  bounds.first = 0.5;
467  bounds.second = 1.;
468  break;
469 
470  default:
471  libMesh::out << "Warning: Invalid quality measure chosen." << std::endl;
472  bounds.first = -1;
473  bounds.second = -1;
474  }
475 
476  return bounds;
477 }
478 
479 
480 
481 const unsigned short int Hex::_second_order_vertex_child_number[27] =
482 {
483  99,99,99,99,99,99,99,99, // Vertices
484  0,1,2,0,0,1,2,3,4,5,6,5, // Edges
485  0,0,1,2,0,4, // Faces
486  0 // Interior
487 };
488 
489 
490 
491 const unsigned short int Hex::_second_order_vertex_child_index[27] =
492 {
493  99,99,99,99,99,99,99,99, // Vertices
494  1,2,3,3,4,5,6,7,5,6,7,7, // Edges
495  2,5,6,7,7,6, // Faces
496  6 // Interior
497 };
498 
499 
500 const unsigned short int Hex::_second_order_adjacent_vertices[12][2] =
501 {
502  { 0, 1}, // vertices adjacent to node 8
503  { 1, 2}, // vertices adjacent to node 9
504  { 2, 3}, // vertices adjacent to node 10
505  { 0, 3}, // vertices adjacent to node 11
506 
507  { 0, 4}, // vertices adjacent to node 12
508  { 1, 5}, // vertices adjacent to node 13
509  { 2, 6}, // vertices adjacent to node 14
510  { 3, 7}, // vertices adjacent to node 15
511 
512  { 4, 5}, // vertices adjacent to node 16
513  { 5, 6}, // vertices adjacent to node 17
514  { 6, 7}, // vertices adjacent to node 18
515  { 4, 7} // vertices adjacent to node 19
516 };
517 
518 } // namespace libMesh

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:04 UTC

Hosted By:
SourceForge.net Logo