face_inf_quad4.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 
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 origin (this->origin());
120 
121  /*
122  * determine the minimal distance of the base from the origin
123  * use size_sq() instead of size(), it is slightly faster
124  */
125  const Real min_distance_sq = std::min((Point(this->point(0)-origin)).size_sq(),
126  (Point(this->point(1)-origin)).size_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-origin).size_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 
166 AutoPtr<Elem> InfQuad4::build_side (const unsigned int i,
167  bool proxy) const
168 {
169  // libmesh_assert_less (i, this->n_sides());
170 
171  if (proxy)
172  {
173  switch (i)
174  {
175  // base
176  case 0:
177  {
178  AutoPtr<Elem> ap(new Side<Edge2,InfQuad4>(this,i));
179  return ap;
180  }
181  // ifem edges
182  case 1:
183  case 2:
184  {
185  AutoPtr<Elem> ap(new Side<InfEdge2,InfQuad4>(this,i));
186  return ap;
187  }
188 
189  default:
190  libmesh_error();
191  }
192  }
193 
194  else
195  {
196  // Create NULL pointer to be initialized, returned later.
197  AutoPtr<Elem> edge(NULL);
198 
199  switch (i)
200  {
201  case 0:
202  {
203  edge.reset(new Edge2);
204 
205  edge->set_node(0) = this->get_node(0);
206  edge->set_node(1) = this->get_node(1);
207 
208  break;
209  }
210 
211  case 1:
212  {
213  // adjacent to another infinite element
214  edge.reset(new InfEdge2);
215 
216  edge->set_node(0) = this->get_node(1);
217  edge->set_node(1) = this->get_node(3);
218 
219  break;
220  }
221 
222  case 2:
223  {
224  // adjacent to another infinite element
225  edge.reset(new InfEdge2);
226 
227  edge->set_node(0) = this->get_node(0);
228  edge->set_node(1) = this->get_node(2);
229 
230  break;
231  }
232  default:
233  {
234  libmesh_error();
235  }
236  }
237 
238  edge->subdomain_id() = this->subdomain_id();
239  return edge;
240  }
241 
242  // How did we get here
243  libmesh_error();
244  AutoPtr<Elem> ap(NULL); return ap;
245 }
246 
247 
248 void InfQuad4::connectivity(const unsigned int libmesh_dbg_var(sf),
249  const IOPackage iop,
250  std::vector<dof_id_type>& conn) const
251 {
252  libmesh_assert_less (sf, this->n_sub_elem());
253  libmesh_assert_not_equal_to (iop, INVALID_IO_PACKAGE);
254 
255  conn.resize(4);
256 
257  switch (iop)
258  {
259  case TECPLOT:
260  {
261  conn[0] = this->node(0)+1;
262  conn[1] = this->node(1)+1;
263  conn[2] = this->node(3)+1;
264  conn[3] = this->node(2)+1;
265  return;
266  }
267  case VTK:
268  {
269  conn[0] = this->node(0);
270  conn[1] = this->node(1);
271  conn[2] = this->node(3);
272  conn[3] = this->node(2);
273  }
274  default:
275  libmesh_error();
276  }
277 
278  libmesh_error();
279 }
280 
281 } // namespace libMesh
282 
283 
284 #endif // ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS

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

Hosted By:
SourceForge.net Logo