libMesh
face_inf_quad4.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 
20 // Local includes
21 #include "libmesh/libmesh_config.h"
22 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
23 
24 
25 // Local includes cont'd
26 #include "libmesh/face_inf_quad4.h"
27 #include "libmesh/fe_interface.h"
28 #include "libmesh/fe_type.h"
29 #include "libmesh/side.h"
30 #include "libmesh/edge_edge2.h"
31 #include "libmesh/edge_inf_edge2.h"
32 
33 namespace libMesh
34 {
35 
36 
37 
38 // ------------------------------------------------------------
39 // InfQuad4 class static member initializations
40 const unsigned int InfQuad4::side_nodes_map[3][2] =
41  {
42  {0, 1}, // Side 0
43  {1, 3}, // Side 1
44  {0, 2} // Side 2
45  };
46 
47 
48 
49 #ifdef LIBMESH_ENABLE_AMR
50 
51 const float InfQuad4::_embedding_matrix[2][4][4] =
52  {
53  // embedding matrix for child 0
54  {
55  // 0 1 2 3rd parent node
56  {1.0, 0.0, 0.0, 0.0}, // 0th child node
57  {0.5, 0.5, 0.0, 0.0}, // 1
58  {0.0, 0.0, 1.0, 0.0}, // 2
59  {0.0, 0.0, 0.5, 0.5} // 3
60  },
61 
62  // embedding matrix for child 1
63  {
64  // 0 1 2 3
65  {0.5, 0.5, 0.0, 0.0}, // 0
66  {0.0, 1.0, 0.0, 0.0}, // 1
67  {0.0, 0.0, 0.5, 0.5}, // 2
68  {0.0, 0.0, 0.0, 1.0} // 3
69  }
70  };
71 
72 #endif
73 
74 
75 // ------------------------------------------------------------
76 // InfQuad4 class member functions
77 
78 bool InfQuad4::is_vertex(const unsigned int i) const
79 {
80  if (i < 2)
81  return true;
82  return false;
83 }
84 
85 bool InfQuad4::is_edge(const unsigned int i) const
86 {
87  if (i < 2)
88  return false;
89  return true;
90 }
91 
92 bool InfQuad4::is_face(const unsigned int) const
93 {
94  return false;
95 }
96 
97 bool InfQuad4::is_node_on_side(const unsigned int n,
98  const unsigned int s) const
99 {
100  libmesh_assert_less (s, n_sides());
101  for (unsigned int i = 0; i != 2; ++i)
102  if (side_nodes_map[s][i] == n)
103  return true;
104  return false;
105 }
106 
107 bool InfQuad4::contains_point (const Point & p, Real tol) const
108 {
109  /*
110  * make use of the fact that infinite elements do not
111  * live inside the envelope. Use a fast scheme to
112  * check whether point \p p is inside or outside
113  * our relevant part of the envelope. Note that
114  * this is not exclusive: the point may be outside
115  * the envelope, but contained in another infinite element.
116  * Therefore, if the distance is greater, do fall back
117  * to the scheme of using FEInterface::inverse_map().
118  */
119  const Point my_origin (this->origin());
120 
121  /*
122  * determine the minimal distance of the base from the origin
123  * use norm_sq() instead of norm(), it is slightly faster
124  */
125  const Real min_distance_sq = std::min((Point(this->point(0)-my_origin)).norm_sq(),
126  (Point(this->point(1)-my_origin)).norm_sq());
127 
128  /*
129  * work with 1% allowable deviation. Can still fall
130  * back to the InfFE::inverse_map()
131  */
132  const Real conservative_p_dist_sq = 1.01 * (Point(p-my_origin).norm_sq());
133 
134  if (conservative_p_dist_sq < min_distance_sq)
135  {
136  /*
137  * the physical point is definitely not contained
138  * in the element, return false.
139  */
140  return false;
141  }
142  else
143  {
144  /*
145  * cannot say anything, fall back to the FEInterface::inverse_map()
146  *
147  * Declare a basic FEType. Will use default in the base,
148  * and something else (not important) in radial direction.
149  */
150  FEType fe_type(default_order());
151 
152  const Point mapped_point = FEInterface::inverse_map(dim(),
153  fe_type,
154  this,
155  p,
156  tol,
157  false);
158 
159  return FEInterface::on_reference_element(mapped_point, this->type(), tol);
160  }
161 }
162 
163 
164 
165 
167  bool proxy)
168 {
169  // libmesh_assert_less (i, this->n_sides());
170 
171  if (proxy)
172  {
173  switch (i)
174  {
175  // base
176  case 0:
177  return UniquePtr<Elem>(new Side<Edge2,InfQuad4>(this,i));
178 
179  // ifem edges
180  case 1:
181  case 2:
182  return UniquePtr<Elem>(new Side<InfEdge2,InfQuad4>(this,i));
183 
184  default:
185  libmesh_error_msg("Invalid side i = " << i);
186  }
187  }
188 
189  else
190  {
191  // Create NULL pointer to be initialized, returned later.
192  Elem * edge = libmesh_nullptr;
193 
194  switch (i)
195  {
196  case 0:
197  {
198  edge = new Edge2;
199  break;
200  }
201 
202  // adjacent to another infinite element
203  case 1:
204  case 2:
205  {
206  edge = new InfEdge2;
207  break;
208  }
209 
210  default:
211  libmesh_error_msg("Invalid side i = " << i);
212  }
213 
214  edge->subdomain_id() = this->subdomain_id();
215 
216  // Set the nodes
217  for (unsigned n=0; n<edge->n_nodes(); ++n)
218  edge->set_node(n) = this->node_ptr(InfQuad4::side_nodes_map[i][n]);
219 
220  return UniquePtr<Elem>(edge);
221  }
222 
223  libmesh_error_msg("We'll never get here!");
224  return UniquePtr<Elem>();
225 }
226 
227 
228 void InfQuad4::connectivity(const unsigned int libmesh_dbg_var(sf),
229  const IOPackage iop,
230  std::vector<dof_id_type> & conn) const
231 {
232  libmesh_assert_less (sf, this->n_sub_elem());
233  libmesh_assert_not_equal_to (iop, INVALID_IO_PACKAGE);
234 
235  conn.resize(4);
236 
237  switch (iop)
238  {
239  case TECPLOT:
240  {
241  conn[0] = this->node_id(0)+1;
242  conn[1] = this->node_id(1)+1;
243  conn[2] = this->node_id(3)+1;
244  conn[3] = this->node_id(2)+1;
245  return;
246  }
247  case VTK:
248  {
249  conn[0] = this->node_id(0);
250  conn[1] = this->node_id(1);
251  conn[2] = this->node_id(3);
252  conn[3] = this->node_id(2);
253  }
254  default:
255  libmesh_error_msg("Unsupported IO package " << iop);
256  }
257 }
258 
259 } // namespace libMesh
260 
261 
262 #endif // ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:178
virtual Node *& set_node(const unsigned int i)
Definition: elem.h:1941
virtual unsigned int dim() const libmesh_override
Definition: face_inf_quad.h:89
virtual bool is_face(const unsigned int i) const libmesh_override
virtual bool is_vertex(const unsigned int i) const libmesh_override
virtual bool is_node_on_side(const unsigned int n, const unsigned int s) const libmesh_override
static const float _embedding_matrix[2][4][4]
Matrix that computes new nodal locations/solution values from current nodes/solution.
virtual Point origin() const libmesh_override
This is the base class from which all geometric element types are derived.
Definition: elem.h:89
const class libmesh_nullptr_t libmesh_nullptr
The InfEdge2 is an infinite element in 1D composed of 2 nodes.
virtual bool contains_point(const Point &p, Real tol=TOLERANCE) const libmesh_override
static const unsigned int side_nodes_map[3][2]
This maps the node of the side to element node numbers.
The libMesh namespace provides an interface to certain functionality in the library.
virtual ElemType type() const libmesh_override
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
virtual unsigned int n_nodes() const =0
virtual unsigned int n_sub_elem() const libmesh_override
const Node * node_ptr(const unsigned int i) const
Definition: elem.h:1874
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
subdomain_id_type subdomain_id() const
Definition: elem.h:1951
PetscErrorCode Vec Mat libmesh_dbg_var(j)
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
This defines the Side class.
Definition: side.h:48
virtual Order default_order() const libmesh_override
virtual bool is_edge(const unsigned int i) const libmesh_override
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const Point & point(const unsigned int i) const
Definition: elem.h:1809
The Edge2 is an element in 1D composed of 2 nodes.
Definition: edge_edge2.h:43
dof_id_type node_id(const unsigned int i) const
Definition: elem.h:1831
IOPackage
libMesh interfaces with several different software packages for the purposes of creating, reading, and writing mesh files.
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_sides() const libmesh_override
Definition: face_inf_quad.h:96
virtual UniquePtr< Elem > build_side_ptr(const unsigned int i, bool proxy) libmesh_override
virtual void connectivity(const unsigned int sf, const IOPackage iop, std::vector< dof_id_type > &conn) const libmesh_override