libMesh
cell_inf_prism.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 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
21 
22 
23 // C++ includes
24 // include <algorithm>
25 
26 // Local includes cont'd
27 #include "libmesh/cell_inf_prism.h"
28 #include "libmesh/cell_inf_prism6.h"
29 #include "libmesh/face_tri3.h"
30 #include "libmesh/face_inf_quad4.h"
31 #include "libmesh/fe_type.h"
32 #include "libmesh/fe_interface.h"
33 
34 namespace libMesh
35 {
36 
37 
38 // ------------------------------------------------------------
39 // InfPrism class static member initializations
40 
41 
42 // We need to require C++11...
43 const Real InfPrism::_master_points[12][3] =
44  {
45  {0, 0, 0},
46  {1, 0, 0},
47  {0, 1, 0},
48  {0, 0, 1},
49  {1, 0, 1},
50  {0, 1, 1},
51  {.5, 0, 0},
52  {.5, .5, 0},
53  {0, .5, 0},
54  {.5, 0, 1},
55  {.5, .5, 1},
56  {0, .5, 1}
57  };
58 
59 
60 
61 // ------------------------------------------------------------
62 // InfPrism class member functions
63 dof_id_type InfPrism::key (const unsigned int s) const
64 {
65  libmesh_assert_less (s, this->n_sides());
66 
67  switch (s)
68  {
69  case 0: // the triangular face at z=-1, base face
70  return this->compute_key (this->node_id(InfPrism6::side_nodes_map[s][0]),
71  this->node_id(InfPrism6::side_nodes_map[s][1]),
72  this->node_id(InfPrism6::side_nodes_map[s][2]));
73 
74  case 1: // the quad face at y=0
75  case 2: // the other quad face
76  case 3: // the quad face at x=0
77  return this->compute_key (this->node_id(InfPrism6::side_nodes_map[s][0]),
78  this->node_id(InfPrism6::side_nodes_map[s][1]),
79  this->node_id(InfPrism6::side_nodes_map[s][2]),
80  this->node_id(InfPrism6::side_nodes_map[s][3]));
81 
82  default:
83  libmesh_error_msg("Invalid side s = " << s);
84  }
85 
86  libmesh_error_msg("We'll never get here!");
87  return 0;
88 }
89 
90 
91 
92 unsigned int InfPrism::which_node_am_i(unsigned int side,
93  unsigned int side_node) const
94 {
95  libmesh_assert_less (side, this->n_sides());
96 
97  // Never more than 4 nodes per side.
98  libmesh_assert_less(side_node, 4);
99 
100  // Some sides have 3 nodes.
101  libmesh_assert(side != 0 || side_node < 3);
102 
103  return InfPrism6::side_nodes_map[side][side_node];
104 }
105 
106 
107 
108 UniquePtr<Elem> InfPrism::side_ptr (const unsigned int i)
109 {
110  libmesh_assert_less (i, this->n_sides());
111 
112  Elem * face = libmesh_nullptr;
113 
114  switch (i)
115  {
116  case 0: // the triangular face at z=-1, base face
117  {
118  face = new Tri3;
119  break;
120  }
121 
122  case 1: // the quad face at y=0
123  case 2: // the other quad face
124  case 3: // the quad face at x=0
125  {
126  face = new InfQuad4;
127  break;
128  }
129 
130  default:
131  libmesh_error_msg("Invalid side i = " << i);
132  }
133 
134  // Set the nodes
135  for (unsigned n=0; n<face->n_nodes(); ++n)
136  face->set_node(n) = this->node_ptr(InfPrism6::side_nodes_map[i][n]);
137 
138  return UniquePtr<Elem>(face);
139 }
140 
141 
142 
143 bool InfPrism::is_child_on_side(const unsigned int c,
144  const unsigned int s) const
145 {
146  libmesh_assert_less (c, this->n_children());
147  libmesh_assert_less (s, this->n_sides());
148 
149  return (s == 0 || c+1 == s || c == s%3);
150 }
151 
152 
153 
154 bool InfPrism::is_edge_on_side (const unsigned int e,
155  const unsigned int s) const
156 {
157  libmesh_assert_less (e, this->n_edges());
158  libmesh_assert_less (s, this->n_sides());
159 
160  return (is_node_on_side(InfPrism6::edge_nodes_map[e][0],s) &&
161  is_node_on_side(InfPrism6::edge_nodes_map[e][1],s));
162 }
163 
164 bool InfPrism::contains_point (const Point & p, Real tol) const
165 {
166  // For infinite elements with linear base interpolation:
167  // make use of the fact that infinite elements do not
168  // live inside the envelope. Use a fast scheme to
169  // check whether point \p p is inside or outside
170  // our relevant part of the envelope. Note that
171  // this is not exclusive: only when the distance is less,
172  // we are safe. Otherwise, we cannot say anything. The
173  // envelope may be non-spherical, the physical point may lie
174  // inside the envelope, outside the envelope, or even inside
175  // this infinite element. Therefore if this fails,
176  // fall back to the FEInterface::inverse_map()
177  const Point my_origin (this->origin());
178 
179  // determine the minimal distance of the base from the origin
180  // use norm_sq(), it is faster than norm() and produces
181  // the same behavior. Shift my_origin to center first:
182  Point pt0_o(this->point(0) - my_origin);
183  Point pt1_o(this->point(1) - my_origin);
184  Point pt2_o(this->point(2) - my_origin);
185  const Real min_distance_sq = std::min(pt0_o.norm_sq(),
186  std::min(pt1_o.norm_sq(),
187  pt2_o.norm_sq()));
188 
189  // work with 1% allowable deviation. We can still fall
190  // back to the InfFE::inverse_map()
191  const Real conservative_p_dist_sq = 1.01 * (Point(p - my_origin).norm_sq());
192 
193  if (conservative_p_dist_sq < min_distance_sq)
194  {
195  // the physical point is definitely not contained in the element:
196  return false;
197  }
198 
199  // this captures the case that the point is not (almost) in the direction of the element.:
200  // first, project the problem onto the unit sphere:
201  Point p_o(p - my_origin);
202  pt0_o /= pt0_o.norm();
203  pt1_o /= pt1_o.norm();
204  pt2_o /= pt2_o.norm();
205  p_o /= p_o.norm();
206 
207  // now, check if it is in the projected face; by comparing the distance of
208  // any point in the element to \p p with the largest distance between this point
209  // to any other point in the element.
210  if ((p_o - pt0_o).norm_sq() > std::max((pt0_o - pt1_o).norm_sq(), (pt0_o - pt2_o).norm_sq()) ||
211  (p_o - pt1_o).norm_sq() > std::max((pt1_o - pt2_o).norm_sq(), (pt1_o - pt0_o).norm_sq()) ||
212  (p_o - pt2_o).norm_sq() > std::max((pt2_o - pt0_o).norm_sq(), (pt2_o - pt1_o).norm_sq()) )
213  {
214  // the physical point is definitely not contained in the element
215  return false;
216  }
217 
218 
219  // It seems that \p p is at least close to this element.
220  //
221  // Declare a basic FEType. Will use default in the base,
222  // and something else (not important) in radial direction.
223  FEType fe_type(default_order());
224 
225  const Point mapped_point = FEInterface::inverse_map(dim(),
226  fe_type,
227  this,
228  p,
229  tol,
230  false);
231 
232  return FEInterface::on_reference_element(mapped_point, this->type(), tol);
233 }
234 
235 
236 } // namespace libMesh
237 
238 
239 
240 #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
The Tri3 is an element in 2D composed of 3 nodes.
Definition: face_tri3.h:54
virtual Node *& set_node(const unsigned int i)
Definition: elem.h:1941
virtual unsigned int n_sides() const libmesh_override
virtual ElemType type() const =0
static const Real _master_points[12][3]
Master element node locations.
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
virtual bool is_edge_on_side(const unsigned int e, const unsigned int s) const libmesh_override
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)
virtual unsigned int which_node_am_i(unsigned int side, unsigned int side_node) const libmesh_override
virtual UniquePtr< Elem > side_ptr(const unsigned int i) libmesh_override
libmesh_assert(j)
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
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 bool is_child_on_side(const unsigned int c, const unsigned int s) const libmesh_override
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 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 n_edges() const libmesh_override
virtual unsigned int dim() const libmesh_override
Definition: cell_inf.h:60
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const Point & point(const unsigned int i) const
Definition: elem.h:1809
virtual unsigned int n_children() const libmesh_override
dof_id_type node_id(const unsigned int i) const
Definition: elem.h:1831
static dof_id_type compute_key(dof_id_type n0)
Definition: elem.h:2621
virtual bool contains_point(const Point &p, Real tol=TOLERANCE) const libmesh_override
long double min(long double a, double b)
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
uint8_t dof_id_type
Definition: id_types.h:64