libMesh
cell_inf_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 // Local includes
19 #include "libmesh/libmesh_config.h"
20 
21 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
22 
23 
24 // C++ includes
25 #include <algorithm> // for std::min, std::max
26 
27 // Local includes cont'd
28 #include "libmesh/cell_inf_hex.h"
29 #include "libmesh/cell_inf_hex8.h"
30 #include "libmesh/face_quad4.h"
31 #include "libmesh/face_inf_quad4.h"
32 #include "libmesh/fe_type.h"
33 #include "libmesh/fe_interface.h"
34 
35 namespace libMesh
36 {
37 
38 
39 
40 // ------------------------------------------------------------
41 // InfHex class static member initializations
42 
43 
44 // We need to require C++11...
45 const Real InfHex::_master_points[18][3] =
46  {
47  {-1, -1, 0},
48  {1, -1, 0},
49  {1, 1, 0},
50  {-1, 1, 0},
51  {-1, -1, 1},
52  {1, -1, 1},
53  {1, 1, 1},
54  {-1, 1, 1},
55  {0, -1, 0},
56  {1, 0, 0},
57  {0, 1, 0},
58  {-1, 0, 0},
59  {0, -1, 1},
60  {1, 0, 1},
61  {0, 1, 1},
62  {-1, 0, 1},
63  {0, 0, 0},
64  {0, 0, 1}
65  };
66 
67 
68 
69 
70 
71 
72 // ------------------------------------------------------------
73 // InfHex class member functions
74 dof_id_type InfHex::key (const unsigned int s) const
75 {
76  libmesh_assert_less (s, this->n_sides());
77 
78  // The order of the node ids does not matter, they are sorted by the
79  // compute_key() function.
80  return this->compute_key(this->node_id(InfHex8::side_nodes_map[s][0]),
81  this->node_id(InfHex8::side_nodes_map[s][1]),
82  this->node_id(InfHex8::side_nodes_map[s][2]),
83  this->node_id(InfHex8::side_nodes_map[s][3]));
84 }
85 
86 
87 
88 unsigned int InfHex::which_node_am_i(unsigned int side,
89  unsigned int side_node) const
90 {
91  libmesh_assert_less (side, this->n_sides());
92  libmesh_assert_less (side_node, 4);
93 
94  return InfHex8::side_nodes_map[side][side_node];
95 }
96 
97 
98 
99 UniquePtr<Elem> InfHex::side_ptr (const unsigned int i)
100 {
101  libmesh_assert_less (i, this->n_sides());
102 
103  // To be returned wrapped in a UniquePtr
104  Elem * face = libmesh_nullptr;
105 
106  // Think of a unit cube: (-1,1) x (-1,1) x (-1,1),
107  // with (in general) the normals pointing outwards
108  switch (i)
109  {
110  // the face at z = -1
111  // the base, where the infinite element couples to conventional
112  // elements
113  case 0:
114  {
115  // Oops, here we are, claiming the normal of the face
116  // elements point outwards -- and this is the exception:
117  // For the side built from the base face,
118  // the normal is pointing _into_ the element!
119  // Why is that? - In agreement with build_side_ptr(),
120  // which in turn _has_ to build the face in this
121  // way as to enable the cool way \p InfFE re-uses \p FE.
122  face = new Quad4;
123  break;
124  }
125 
126  // These faces connect to other infinite elements.
127  case 1: // the face at y = -1
128  case 2: // the face at x = 1
129  case 3: // the face at y = 1
130  case 4: // the face at x = -1
131  {
132  face = new InfQuad4;
133  break;
134  }
135 
136  default:
137  libmesh_error_msg("Invalid side i = " << i);
138  }
139 
140  // Set the nodes
141  for (unsigned n=0; n<face->n_nodes(); ++n)
142  face->set_node(n) = this->node_ptr(InfHex8::side_nodes_map[i][n]);
143 
144  return UniquePtr<Elem>(face);
145 }
146 
147 
148 
149 bool InfHex::is_child_on_side(const unsigned int c,
150  const unsigned int s) const
151 {
152  libmesh_assert_less (c, this->n_children());
153  libmesh_assert_less (s, this->n_sides());
154 
155  return (s == 0 || c+1 == s || c == s%4);
156 }
157 
158 
159 
160 bool InfHex::is_edge_on_side (const unsigned int e,
161  const unsigned int s) const
162 {
163  libmesh_assert_less (e, this->n_edges());
164  libmesh_assert_less (s, this->n_sides());
165 
166  return (is_node_on_side(InfHex8::edge_nodes_map[e][0],s) &&
167  is_node_on_side(InfHex8::edge_nodes_map[e][1],s));
168 }
169 
170 
171 
172 
174 {
175  switch (q)
176  {
177 
187  case DIAGONAL:
188  {
189  // Diagonal between node 0 and node 2
190  const Real d02 = this->length(0,2);
191 
192  // Diagonal between node 1 and node 3
193  const Real d13 = this->length(1,3);
194 
195  // Find the biggest and smallest diagonals
196  const Real min = std::min(d02, d13);
197  const Real max = std::max(d02, d13);
198 
199  libmesh_assert_not_equal_to (max, 0.0);
200 
201  return min / max;
202 
203  break;
204  }
205 
213  case TAPER:
214  {
215 
219  const Real d01 = this->length(0,1);
220  const Real d12 = this->length(1,2);
221  const Real d23 = this->length(2,3);
222  const Real d03 = this->length(0,3);
223 
224  std::vector<Real> edge_ratios(2);
225 
226  // Bottom
227  edge_ratios[8] = std::min(d01, d23) / std::max(d01, d23);
228  edge_ratios[9] = std::min(d03, d12) / std::max(d03, d12);
229 
230  return *(std::min_element(edge_ratios.begin(), edge_ratios.end())) ;
231 
232  break;
233  }
234 
235 
243  case STRETCH:
244  {
248  const Real sqrt3 = 1.73205080756888;
249 
253  const Real d02 = this->length(0,2);
254  const Real d13 = this->length(1,3);
255  const Real max_diag = std::max(d02, d13);
256 
257  libmesh_assert_not_equal_to ( max_diag, 0.0 );
258 
262  std::vector<Real> edges(4);
263  edges[0] = this->length(0,1);
264  edges[1] = this->length(1,2);
265  edges[2] = this->length(2,3);
266  edges[3] = this->length(0,3);
267 
268  const Real min_edge = *(std::min_element(edges.begin(), edges.end()));
269  return sqrt3 * min_edge / max_diag ;
270  }
271 
272 
277  default:
278  return Elem::quality(q);
279  }
280 
281  libmesh_error_msg("We'll never get here!");
282  return 0.;
283 }
284 
285 
286 
287 std::pair<Real, Real> InfHex::qual_bounds (const ElemQuality) const
288 {
289  libmesh_not_implemented();
290 
291  std::pair<Real, Real> bounds;
292 
293  /*
294  switch (q)
295  {
296 
297  case ASPECT_RATIO:
298  bounds.first = 1.;
299  bounds.second = 4.;
300  break;
301 
302  case SKEW:
303  bounds.first = 0.;
304  bounds.second = 0.5;
305  break;
306 
307  case SHEAR:
308  case SHAPE:
309  bounds.first = 0.3;
310  bounds.second = 1.;
311  break;
312 
313  case CONDITION:
314  bounds.first = 1.;
315  bounds.second = 8.;
316  break;
317 
318  case JACOBIAN:
319  bounds.first = 0.5;
320  bounds.second = 1.;
321  break;
322 
323  case DISTORTION:
324  bounds.first = 0.6;
325  bounds.second = 1.;
326  break;
327 
328  case TAPER:
329  bounds.first = 0.;
330  bounds.second = 0.4;
331  break;
332 
333  case STRETCH:
334  bounds.first = 0.25;
335  bounds.second = 1.;
336  break;
337 
338  case DIAGONAL:
339  bounds.first = 0.65;
340  bounds.second = 1.;
341  break;
342 
343  case SIZE:
344  bounds.first = 0.5;
345  bounds.second = 1.;
346  break;
347 
348  default:
349  libMesh::out << "Warning: Invalid quality measure chosen." << std::endl;
350  bounds.first = -1;
351  bounds.second = -1;
352  }
353  */
354 
355  return bounds;
356 }
357 
358 
359 
360 
361 
362 const unsigned short int InfHex::_second_order_adjacent_vertices[8][2] =
363  {
364  { 0, 1}, // vertices adjacent to node 8
365  { 1, 2}, // vertices adjacent to node 9
366  { 2, 3}, // vertices adjacent to node 10
367  { 0, 3}, // vertices adjacent to node 11
368 
369  { 4, 5}, // vertices adjacent to node 12
370  { 5, 6}, // vertices adjacent to node 13
371  { 6, 7}, // vertices adjacent to node 14
372  { 4, 7} // vertices adjacent to node 15
373  };
374 
375 
376 
377 const unsigned short int InfHex::_second_order_vertex_child_number[18] =
378  {
379  99,99,99,99,99,99,99,99, // Vertices
380  0,1,2,0, // Edges
381  0,1,2,0,0, // Faces
382  0 // Interior
383  };
384 
385 
386 
387 const unsigned short int InfHex::_second_order_vertex_child_index[18] =
388  {
389  99,99,99,99,99,99,99,99, // Vertices
390  1,2,3,3, // Edges
391  5,6,7,7,2, // Faces
392  6 // Interior
393  };
394 
395 
396 bool InfHex::contains_point (const Point & p, Real tol) const
397 {
398  // For infinite elements with linear base interpolation:
399  // make use of the fact that infinite elements do not
400  // live inside the envelope. Use a fast scheme to
401  // check whether point \p p is inside or outside
402  // our relevant part of the envelope. Note that
403  // this is not exclusive: only when the distance is less,
404  // we are safe. Otherwise, we cannot say anything. The
405  // envelope may be non-spherical, the physical point may lie
406  // inside the envelope, outside the envelope, or even inside
407  // this infinite element. Therefore if this fails,
408  // fall back to the FEInterface::inverse_map()
409  const Point my_origin (this->origin());
410 
411  // determine the minimal distance of the base from the origin
412  // Use norm_sq() instead of norm(), it is faster
413  Point pt0_o(this->point(0) - my_origin);
414  Point pt1_o(this->point(1) - my_origin);
415  Point pt2_o(this->point(2) - my_origin);
416  Point pt3_o(this->point(3) - my_origin);
417  const Real min_distance_sq = std::min(pt0_o.norm_sq(),
418  std::min(pt1_o.norm_sq(),
419  std::min(pt2_o.norm_sq(),
420  pt3_o.norm_sq())));
421 
422  // work with 1% allowable deviation. We can still fall
423  // back to the InfFE::inverse_map()
424  const Real conservative_p_dist_sq = 1.01 * (Point(p - my_origin).norm_sq());
425 
426 
427 
428  if (conservative_p_dist_sq < min_distance_sq)
429  {
430  // the physical point is definitely not contained in the element
431  return false;
432  }
433 
434  // this captures the case that the point is not (almost) in the direction of the element.:
435  // first, project the problem onto the unit sphere:
436  Point p_o(p - my_origin);
437  pt0_o /= pt0_o.norm();
438  pt1_o /= pt1_o.norm();
439  pt2_o /= pt2_o.norm();
440  pt3_o /= pt3_o.norm();
441  p_o /= p_o.norm();
442 
443 
444  // now, check if it is in the projected face; using that the diagonal contains
445  // the largest distance between points in it
446  Real max_h = std::max((pt0_o - pt2_o).norm_sq(),
447  (pt1_o - pt2_o).norm_sq())*1.01;
448 
449  if ((p_o - pt0_o).norm_sq() > max_h ||
450  (p_o - pt1_o).norm_sq() > max_h ||
451  (p_o - pt2_o).norm_sq() > max_h ||
452  (p_o - pt3_o).norm_sq() > max_h )
453  {
454  // the physical point is definitely not contained in the element
455  return false;
456  }
457 
458  // Declare a basic FEType. Will use default in the base,
459  // and something else (not important) in radial direction.
460  FEType fe_type(default_order());
461 
462  const Point mapped_point = FEInterface::inverse_map(dim(),
463  fe_type,
464  this,
465  p,
466  tol,
467  false);
468 
469  return FEInterface::on_reference_element(mapped_point, this->type(), tol);
470 }
471 
472 } // namespace libMesh
473 
474 
475 
476 
477 #endif // ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:178
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
virtual ElemType type() const =0
The QUAD4 is an element in 2D composed of 4 nodes.
Definition: face_quad4.h:49
virtual bool is_edge_on_side(const unsigned int e, const unsigned int s) const libmesh_override
Definition: cell_inf_hex.C:160
The INFQUAD4 is an infinite element in 2D composed of 4 nodes.
UniquePtr< Elem > side(const unsigned int i) const
Definition: elem.h:2105
This is the base class from which all geometric element types are derived.
Definition: elem.h:89
Real norm() const
Definition: type_vector.h:909
const class libmesh_nullptr_t libmesh_nullptr
virtual Point origin() const libmesh_override
Definition: cell_inf.h:71
The libMesh namespace provides an interface to certain functionality in the library.
long double max(long double a, double b)
static const unsigned short int _second_order_adjacent_vertices[8][2]
For higher-order elements, namely InfHex16 and InfHex18, the matrices for adjacent vertices of second...
Definition: cell_inf_hex.h:185
virtual unsigned int n_children() const libmesh_override
Definition: cell_inf_hex.h:108
virtual bool contains_point(const Point &p, Real tol=TOLERANCE) const libmesh_override
Definition: cell_inf_hex.C:396
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
virtual unsigned int n_nodes() const =0
static const Real _master_points[18][3]
Master element node locations.
Definition: cell_inf_hex.h:200
virtual dof_id_type key() const
Definition: elem.C:503
virtual unsigned int n_sides() const libmesh_override
Definition: cell_inf_hex.h:79
const Node * node_ptr(const unsigned int i) const
Definition: elem.h:1874
virtual unsigned int which_node_am_i(unsigned int side, unsigned int side_node) const libmesh_override
Definition: cell_inf_hex.C:88
virtual Real quality(const ElemQuality q) const
Definition: elem.C:1547
static Point inverse_map(const unsigned int dim, const FEType &fe_t, const Elem *elem, const Point &p, const Real tolerance=TOLERANCE, const bool secure=true)
Definition: fe_interface.C:569
virtual bool is_child_on_side(const unsigned int c, const unsigned int s) const libmesh_override
Definition: cell_inf_hex.C:149
virtual Order default_order() const =0
Real norm_sq() const
Definition: type_vector.h:940
static bool on_reference_element(const Point &p, const ElemType t, const Real eps=TOLERANCE)
Definition: fe_interface.C:631
virtual unsigned int dim() const libmesh_override
Definition: cell_inf.h:60
static const unsigned short int _second_order_vertex_child_index[18]
Vector that names the child vertex index for each second order node.
Definition: cell_inf_hex.h:195
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const Point & point(const unsigned int i) const
Definition: elem.h:1809
virtual UniquePtr< Elem > side_ptr(const unsigned int i) libmesh_override
Definition: cell_inf_hex.C:99
virtual Real quality(const ElemQuality q) const libmesh_override
Definition: cell_inf_hex.C:173
Real length(const unsigned int n1, const unsigned int n2) const
Definition: elem.C:492
virtual std::pair< Real, Real > qual_bounds(const ElemQuality q) const libmesh_override
Definition: cell_inf_hex.C:287
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
long double min(long double a, double b)
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
virtual unsigned int n_edges() const libmesh_override
Definition: cell_inf_hex.h:98
static const unsigned short int _second_order_vertex_child_number[18]
Vector that names a child sharing each second order node.
Definition: cell_inf_hex.h:190
uint8_t dof_id_type
Definition: id_types.h:64