libMesh
cell_inf_prism.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2024 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 #include "libmesh/libmesh_config.h"
19 
20 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
21 
22 // Local includes
23 #include "libmesh/cell_inf_prism.h"
24 #include "libmesh/cell_inf_prism6.h"
25 #include "libmesh/face_tri3.h"
26 #include "libmesh/face_inf_quad4.h"
27 #include "libmesh/fe_type.h"
28 #include "libmesh/fe_interface.h"
29 #include "libmesh/inf_fe_map.h"
30 #include "libmesh/enum_order.h"
31 
32 namespace libMesh
33 {
34 
35 
36 // ------------------------------------------------------------
37 // InfPrism class static member initializations
38 
39 
40 // We need to require C++11...
41 const Real InfPrism::_master_points[12][3] =
42  {
43  {0, 0, 0},
44  {1, 0, 0},
45  {0, 1, 0},
46  {0, 0, 1},
47  {1, 0, 1},
48  {0, 1, 1},
49  {.5, 0, 0},
50  {.5, .5, 0},
51  {0, .5, 0},
52  {.5, 0, 1},
53  {.5, .5, 1},
54  {0, .5, 1}
55  };
56 
57 const unsigned int InfPrism::edge_sides_map[6][2] =
58  {
59  {0, 1}, // Edge 0
60  {0, 2}, // Edge 1
61  {0, 3}, // Edge 2
62  {1, 3}, // Edge 3
63  {1, 2}, // Edge 4
64  {2, 3} // Edge 5
65  };
66 
67 // ------------------------------------------------------------
68 // InfPrism class member functions
69 dof_id_type InfPrism::key (const unsigned int s) const
70 {
71  libmesh_assert_less (s, this->n_sides());
72 
73  switch (s)
74  {
75  case 0: // the triangular face at z=-1, base face
76  return this->compute_key (this->node_id(InfPrism6::side_nodes_map[s][0]),
78  this->node_id(InfPrism6::side_nodes_map[s][2]));
79 
80  case 1: // the quad face at y=0
81  case 2: // the other quad face
82  case 3: // the quad face at x=0
83  return this->compute_key (this->node_id(InfPrism6::side_nodes_map[s][0]),
86  this->node_id(InfPrism6::side_nodes_map[s][3]));
87 
88  default:
89  libmesh_error_msg("Invalid side s = " << s);
90  }
91 }
92 
93 
94 
95 dof_id_type InfPrism::low_order_key (const unsigned int s) const
96 {
97  libmesh_assert_less (s, this->n_sides());
98 
99  switch (s)
100  {
101  case 0: // the triangular face at z=-1, base face
102  return this->compute_key (this->node_id(InfPrism6::side_nodes_map[s][0]),
103  this->node_id(InfPrism6::side_nodes_map[s][1]),
104  this->node_id(InfPrism6::side_nodes_map[s][2]));
105 
106  case 1: // the quad face at y=0
107  case 2: // the other quad face
108  case 3: // the quad face at x=0
109  return this->compute_key (this->node_id(InfPrism6::side_nodes_map[s][0]),
110  this->node_id(InfPrism6::side_nodes_map[s][1]),
111  this->node_id(InfPrism6::side_nodes_map[s][2]),
112  this->node_id(InfPrism6::side_nodes_map[s][3]));
113 
114  default:
115  libmesh_error_msg("Invalid side s = " << s);
116  }
117 }
118 
119 
120 
121 unsigned int InfPrism::local_side_node(unsigned int side,
122  unsigned int side_node) const
123 {
124  libmesh_assert_less (side, this->n_sides());
125 
126  // Never more than 4 nodes per side.
127  libmesh_assert_less(side_node, InfPrism6::nodes_per_side);
128 
129  // Some sides have 3 nodes.
130  libmesh_assert(side != 0 || side_node < 3);
131 
132  return InfPrism6::side_nodes_map[side][side_node];
133 }
134 
135 
136 
137 unsigned int InfPrism::local_edge_node(unsigned int edge,
138  unsigned int edge_node) const
139 {
140  libmesh_assert_less(edge, this->n_edges());
141  libmesh_assert_less(edge_node, InfPrism6::nodes_per_edge);
142 
143  return InfPrism6::edge_nodes_map[edge][edge_node];
144 }
145 
146 
147 
148 std::unique_ptr<Elem> InfPrism::side_ptr (const unsigned int i)
149 {
150  libmesh_assert_less (i, this->n_sides());
151 
152  std::unique_ptr<Elem> face;
153 
154  switch (i)
155  {
156  case 0: // the triangular face at z=-1, base face
157  {
158  face = std::make_unique<Tri3>();
159  break;
160  }
161 
162  case 1: // the quad face at y=0
163  case 2: // the other quad face
164  case 3: // the quad face at x=0
165  {
166  face = std::make_unique<InfQuad4>();
167  break;
168  }
169 
170  default:
171  libmesh_error_msg("Invalid side i = " << i);
172  }
173 
174  // Set the nodes
175  for (auto n : face->node_index_range())
176  face->set_node(n) = this->node_ptr(InfPrism6::side_nodes_map[i][n]);
177 
178  return face;
179 }
180 
181 
182 
183 void InfPrism::side_ptr (std::unique_ptr<Elem> & side,
184  const unsigned int i)
185 {
186  libmesh_assert_less (i, this->n_sides());
187 
188  switch (i)
189  {
190  // the base face
191  case 0:
192  {
193  if (!side.get() || side->type() != TRI3)
194  {
195  side = this->side_ptr(i);
196  return;
197  }
198  break;
199  }
200 
201  // connecting to another infinite element
202  case 1:
203  case 2:
204  case 3:
205  case 4:
206  {
207  if (!side.get() || side->type() != INFQUAD4)
208  {
209  side = this->side_ptr(i);
210  return;
211  }
212  break;
213  }
214 
215  default:
216  libmesh_error_msg("Invalid side i = " << i);
217  }
218 
219  side->subdomain_id() = this->subdomain_id();
220 
221  // Set the nodes
222  for (auto n : side->node_index_range())
223  side->set_node(n) = this->node_ptr(InfPrism6::side_nodes_map[i][n]);
224 }
225 
226 
227 
228 bool InfPrism::is_child_on_side(const unsigned int c,
229  const unsigned int s) const
230 {
231  libmesh_assert_less (c, this->n_children());
232  libmesh_assert_less (s, this->n_sides());
233 
234  return (s == 0 || c+1 == s || c == s%3);
235 }
236 
237 
238 
239 bool InfPrism::is_edge_on_side (const unsigned int e,
240  const unsigned int s) const
241 {
242  libmesh_assert_less (e, this->n_edges());
243  libmesh_assert_less (s, this->n_sides());
244 
245  return (edge_sides_map[e][0] == s || edge_sides_map[e][1] == s);
246 }
247 
248 bool InfPrism::contains_point (const Point & p, Real tol) const
249 {
250  // For infinite elements with linear base interpolation:
251  // make use of the fact that infinite elements do not
252  // live inside the envelope. Use a fast scheme to
253  // check whether point \p p is inside or outside
254  // our relevant part of the envelope. Note that
255  // this is not exclusive: only when the distance is less,
256  // we are safe. Otherwise, we cannot say anything. The
257  // envelope may be non-spherical, the physical point may lie
258  // inside the envelope, outside the envelope, or even inside
259  // this infinite element. Therefore if this fails,
260  // fall back to the inverse_map()
261  const Point my_origin (this->origin());
262 
263  // determine the minimal distance of the base from the origin
264  // use norm_sq(), it is faster than norm() and produces
265  // the same behavior. Shift my_origin to center first:
266  Point pt0_o(this->point(0) - my_origin);
267  Point pt1_o(this->point(1) - my_origin);
268  Point pt2_o(this->point(2) - my_origin);
269  const Real tmp_min_distance_sq = std::min(pt0_o.norm_sq(),
270  std::min(pt1_o.norm_sq(),
271  pt2_o.norm_sq()));
272 
273  // For a coarse grid, it is important to account for the fact
274  // that the sides are not spherical, thus the closest point
275  // can be closer than all edges.
276  // This is an estimator using Pythagoras:
277  const Real min_distance_sq = tmp_min_distance_sq
278  - .5*std::max((point(0)-point(1)).norm_sq(),
279  std::max((point(0)-point(2)).norm_sq(),
280  (point(1)-point(2)).norm_sq()));
281 
282  // work with 1% allowable deviation. We can still fall
283  // back to the InfFE::inverse_map()
284  const Real conservative_p_dist_sq = 1.01 * (Point(p - my_origin).norm_sq());
285 
286  if (conservative_p_dist_sq < min_distance_sq)
287  {
288  // the physical point is definitely not contained in the element:
289  return false;
290  }
291 
292  // this captures the case that the point is not (almost) in the direction of the element.:
293  // first, project the problem onto the unit sphere:
294  Point p_o(p - my_origin);
295  pt0_o /= pt0_o.norm();
296  pt1_o /= pt1_o.norm();
297  pt2_o /= pt2_o.norm();
298  p_o /= p_o.norm();
299 
300  // now, check if it is in the projected face; by comparing the distance of
301  // any point in the element to \p p with the largest distance between this point
302  // to any other point in the element.
303  if ((p_o - pt0_o).norm_sq() > std::max((pt0_o - pt1_o).norm_sq(), (pt0_o - pt2_o).norm_sq()) ||
304  (p_o - pt1_o).norm_sq() > std::max((pt1_o - pt2_o).norm_sq(), (pt1_o - pt0_o).norm_sq()) ||
305  (p_o - pt2_o).norm_sq() > std::max((pt2_o - pt0_o).norm_sq(), (pt2_o - pt1_o).norm_sq()) )
306  {
307  // the physical point is definitely not contained in the element
308  return false;
309  }
310 
311 
312  // It seems that \p p is at least close to this element.
313  //
314  // Declare a basic FEType. Will use default in the base,
315  // and something else (not important) in radial direction.
316  FEType fe_type(default_order());
317 
318  const Point mapped_point = InfFEMap::inverse_map(dim(), this, p,
319  tol, false);
320 
321  return FEInterface::on_reference_element(mapped_point, this->type(), tol);
322 }
323 
324 std::vector<unsigned>
325 InfPrism::sides_on_edge(const unsigned int e) const
326 {
327  libmesh_assert_less(e, n_edges());
328  return {std::begin(edge_sides_map[e]), std::end(edge_sides_map[e])};
329 }
330 
331 
332 bool
334 {
335  return (triple_product(this->point(1)-this->point(0),
336  this->point(2)-this->point(0),
337  this->point(3)-this->point(0)) < 0);
338 }
339 
340 
341 } // namespace libMesh
342 
343 
344 
345 #endif // ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:182
auto norm() const -> decltype(std::norm(T()))
Definition: type_vector.h:929
static const Real _master_points[12][3]
Master element node locations.
virtual dof_id_type key() const
Definition: elem.C:563
virtual unsigned int n_edges() const override final
The libMesh namespace provides an interface to certain functionality in the library.
virtual bool is_flipped() const override final
static Point inverse_map(const unsigned int dim, const Elem *elem, const Point &p, const Real tolerance=TOLERANCE, const bool secure=true)
Definition: inf_fe_map.C:96
virtual unsigned int n_children() const override final
std::vector< unsigned int > sides_on_edge(const unsigned int e) const override final
virtual unsigned int local_side_node(unsigned int side, unsigned int side_node) const override
T triple_product(const TypeVector< T > &a, const TypeVector< T > &b, const TypeVector< T > &c)
Definition: type_vector.h:1068
auto norm_sq() const -> decltype(std::norm(T()))
Definition: type_vector.h:948
static const unsigned int side_nodes_map[num_sides][nodes_per_side]
This maps the node of the side to element node numbers.
virtual unsigned int local_edge_node(unsigned int edge, unsigned int edge_node) const override
virtual std::unique_ptr< Elem > side_ptr(const unsigned int i) override final
libmesh_assert(ctx)
virtual bool is_child_on_side(const unsigned int c, const unsigned int s) const override final
static bool on_reference_element(const Point &p, const ElemType t, const Real eps=TOLERANCE)
Definition: fe_interface.C:990
virtual bool is_edge_on_side(const unsigned int e, const unsigned int s) const override final
static const unsigned int edge_nodes_map[num_edges][nodes_per_edge]
This maps the node of the edge to element node numbers.
virtual Point origin() const override
Definition: cell_inf.h:77
static const unsigned int edge_sides_map[6][2]
This maps each edge to the sides that contain said edge.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
subdomain_id_type subdomain_id() const
Definition: elem.h:2391
static const int nodes_per_side
const Node * node_ptr(const unsigned int i) const
Definition: elem.h:2331
auto norm_sq(const T &a) -> decltype(std::norm(a))
Definition: tensor_tools.h:104
virtual unsigned short dim() const override
Definition: cell_inf.h:66
virtual bool contains_point(const Point &p, Real tol=TOLERANCE) const override
static const int nodes_per_edge
static dof_id_type compute_key(dof_id_type n0)
Definition: elem.h:3131
virtual Order default_order() const =0
virtual unsigned int n_sides() const override final
virtual ElemType type() const =0
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
dof_id_type node_id(const unsigned int i) const
Definition: elem.h:2299
const Point & point(const unsigned int i) const
Definition: elem.h:2277
virtual dof_id_type low_order_key(const unsigned int s) const override
uint8_t dof_id_type
Definition: id_types.h:67