cell_inf_prism6.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 // Local includes
19 #include "libmesh/libmesh_config.h"
20 
21 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
22 
23 // C++ includes
24 
25 // Local includes cont'd
27 #include "libmesh/edge_edge2.h"
28 #include "libmesh/edge_inf_edge2.h"
29 #include "libmesh/fe_interface.h"
30 #include "libmesh/fe_type.h"
31 #include "libmesh/side.h"
32 #include "libmesh/face_inf_quad4.h"
33 #include "libmesh/face_tri3.h"
34 
35 namespace libMesh
36 {
37 
38 
39 // ------------------------------------------------------------
40 // InfPrism6 class static member initializations
41 const unsigned int InfPrism6::side_nodes_map[4][4] =
42 {
43  { 0, 1, 2, 99}, // Side 0
44  { 0, 1, 3, 4}, // Side 1
45  { 1, 2, 4, 5}, // Side 2
46  { 2, 0, 5, 3} // Side 3
47 };
48 
49 const unsigned int InfPrism6::edge_nodes_map[6][2] =
50 {
51  { 0, 1}, // Side 0
52  { 1, 2}, // Side 1
53  { 0, 2}, // Side 2
54  { 0, 3}, // Side 3
55  { 1, 4}, // Side 4
56  { 2, 5} // Side 5
57 };
58 
59 
60 // ------------------------------------------------------------
61 // InfPrism6 class member functions
62 
63 bool InfPrism6::is_vertex(const unsigned int i) const
64 {
65  if (i < 3)
66  return true;
67  return false;
68 }
69 
70 bool InfPrism6::is_edge(const unsigned int i) const
71 {
72  if (i < 3)
73  return false;
74  return true;
75 }
76 
77 bool InfPrism6::is_face(const unsigned int) const
78 {
79  return false;
80 }
81 
82 bool InfPrism6::is_node_on_side(const unsigned int n,
83  const unsigned int s) const
84 {
85  libmesh_assert_less (s, n_sides());
86  for (unsigned int i = 0; i != 4; ++i)
87  if (side_nodes_map[s][i] == n)
88  return true;
89  return false;
90 }
91 
92 bool InfPrism6::is_node_on_edge(const unsigned int n,
93  const unsigned int e) const
94 {
95  libmesh_assert_less (e, n_edges());
96  for (unsigned int i = 0; i != 2; ++i)
97  if (edge_nodes_map[e][i] == n)
98  return true;
99  return false;
100 }
101 
102 
103 AutoPtr<Elem> InfPrism6::build_side (const unsigned int i,
104  bool proxy) const
105 {
106  libmesh_assert_less (i, this->n_sides());
107 
108  if (proxy)
109  {
110  switch (i)
111  {
112  // base
113  case 0:
114  {
115  AutoPtr<Elem> ap(new Side<Tri3,InfPrism6>(this,i));
116  return ap;
117  }
118  // ifem sides
119  case 1:
120  case 2:
121  case 3:
122  {
123  AutoPtr<Elem> ap(new Side<InfQuad4,InfPrism6>(this,i));
124  return ap;
125  }
126  default:
127  libmesh_error();
128  }
129  }
130 
131  else
132  {
133  // Create NULL pointer to be initialized, returned later.
134  AutoPtr<Elem> face(NULL);
135 
136  switch (i)
137  {
138  case 0: // the triangular face at z=-1, base face
139  {
140  face.reset(new Tri3);
141 
142  // Note that for this face element, the normal points inward
143  face->set_node(0) = this->get_node(0);
144  face->set_node(1) = this->get_node(1);
145  face->set_node(2) = this->get_node(2);
146 
147  break;
148  }
149 
150  case 1: // the quad face at y=0
151  {
152  face.reset(new InfQuad4);
153 
154  face->set_node(0) = this->get_node(0);
155  face->set_node(1) = this->get_node(1);
156  face->set_node(2) = this->get_node(3);
157  face->set_node(3) = this->get_node(4);
158 
159  break;
160  }
161 
162  case 2: // the other quad face
163  {
164  face.reset(new InfQuad4);
165 
166  face->set_node(0) = this->get_node(1);
167  face->set_node(1) = this->get_node(2);
168  face->set_node(2) = this->get_node(4);
169  face->set_node(3) = this->get_node(5);
170 
171  break;
172  }
173 
174  case 3: // the quad face at x=0
175  {
176  face.reset(new InfQuad4);
177 
178  face->set_node(0) = this->get_node(2);
179  face->set_node(1) = this->get_node(0);
180  face->set_node(2) = this->get_node(5);
181  face->set_node(3) = this->get_node(3);
182 
183  break;
184  }
185 
186  default:
187  {
188  libmesh_error();
189  }
190  }
191 
192  face->subdomain_id() = this->subdomain_id();
193  return face;
194  }
195 
196 
197  // We will never get here... Look at the code above.
198  libmesh_error();
199  AutoPtr<Elem> ap(NULL); return ap;
200 }
201 
202 
203 AutoPtr<Elem> InfPrism6::build_edge (const unsigned int i) const
204 {
205  libmesh_assert_less (i, n_edges());
206 
207  if (i < 3)
208  return AutoPtr<Elem>(new SideEdge<Edge2,InfPrism6>(this,i));
209  return AutoPtr<Elem>(new SideEdge<InfEdge2,InfPrism6>(this,i));
210 }
211 
212 
213 bool InfPrism6::contains_point (const Point& p, Real tol) const
214 {
215  /*
216  * For infinite elements with linear base interpolation:
217  *
218  * make use of the fact that infinite elements do not
219  * live inside the envelope. Use a fast scheme to
220  * check whether point \p p is inside or outside
221  * our relevant part of the envelope. Note that
222  * this is not exclusive: only when the distance is less,
223  * we are safe. Otherwise, we cannot say anything. The
224  * envelope may be non-spherical, the physical point may lie
225  * inside the envelope, outside the envelope, or even inside
226  * this infinite element. Therefore if this fails,
227  * fall back to the FEInterface::inverse_map()
228  */
229  const Point origin (this->origin());
230 
231  /*
232  * determine the minimal distance of the base from the origin
233  * use size_sq(), it is faster than size() and produces
234  * the same behavior
235  */
236  const Real min_distance_sq = std::min((Point(this->point(0)-origin)).size_sq(),
237  std::min((Point(this->point(1)-origin)).size_sq(),
238  (Point(this->point(2)-origin)).size_sq()));
239 
240  /*
241  * work with 1% allowable deviation. We can still fall
242  * back to the InfFE::inverse_map()
243  */
244  const Real conservative_p_dist_sq = 1.01 * (Point(p-origin).size_sq());
245 
246 
247  if (conservative_p_dist_sq < min_distance_sq)
248  {
249  /*
250  * the physical point is definitely not contained in the element
251  */
252  return false;
253  }
254  else
255  {
256  /*
257  * Declare a basic FEType. Will use default in the base,
258  * and something else (not important) in radial direction.
259  */
260  FEType fe_type(default_order());
261 
262  const Point mapped_point = FEInterface::inverse_map(dim(),
263  fe_type,
264  this,
265  p,
266  tol,
267  false);
268 
269  return FEInterface::on_reference_element(mapped_point, this->type(), tol);
270  }
271 }
272 
273 
274 
275 
276 void InfPrism6::connectivity(const unsigned int libmesh_dbg_var(sc),
277  const IOPackage iop,
278  std::vector<dof_id_type>& conn) const
279 {
281  libmesh_assert_less (sc, this->n_sub_elem());
282  libmesh_assert_not_equal_to (iop, INVALID_IO_PACKAGE);
283 
284  switch (iop)
285  {
286  case TECPLOT:
287  {
288  conn.resize(8);
289  conn[0] = this->node(0)+1;
290  conn[1] = this->node(1)+1;
291  conn[2] = this->node(2)+1;
292  conn[3] = this->node(2)+1;
293  conn[4] = this->node(3)+1;
294  conn[5] = this->node(4)+1;
295  conn[6] = this->node(5)+1;
296  conn[7] = this->node(5)+1;
297  return;
298  }
299 
300  default:
301  libmesh_error();
302  }
303 
304  libmesh_error();
305 }
306 
307 
308 
309 
310 
311 #ifdef LIBMESH_ENABLE_AMR
312 
313 const float InfPrism6::_embedding_matrix[4][6][6] =
314 {
315  // embedding matrix for child 0
316  {
317  // 0 1 2 3 4 5 th parent Node
318  { 1.0, 0.0, 0.0, 0.0, 0.0, 0.0}, // 0th child N.
319  { 0.5, 0.5, 0.0, 0.0, 0.0, 0.0}, // 1
320  { 0.5, 0.0, 0.5, 0.0, 0.0, 0.0}, // 2
321  { 0.0, 0.0, 0.0, 1.0, 0.0, 0.0}, // 3
322  { 0.0, 0.0, 0.0, 0.5, 0.5, 0.0}, // 4
323  { 0.0, 0.0, 0.0, 0.5, 0.0, 0.5} // 5
324  },
325 
326  // embedding matrix for child 1
327  {
328  // 0 1 2 3 4 5 th parent Node
329  { 0.5, 0.5, 0.0, 0.0, 0.0, 0.0}, // 0th child N.
330  { 0.0, 1.0, 0.0, 0.0, 0.0, 0.0}, // 1
331  { 0.0, 0.5, 0.5, 0.0, 0.0, 0.0}, // 2
332  { 0.0, 0.0, 0.0, 0.5, 0.5, 0.0}, // 3
333  { 0.0, 0.0, 0.0, 0.0, 1.0, 0.0}, // 4
334  { 0.0, 0.0, 0.0, 0.0, 0.5, 0.5} // 5
335  },
336 
337  // embedding matrix for child 2
338  {
339  // 0 1 2 3 4 5 th parent Node
340  { 0.5, 0.0, 0.5, 0.0, 0.0, 0.0}, // 0th child N.
341  { 0.0, 0.5, 0.5, 0.0, 0.0, 0.0}, // 1
342  { 0.0, 0.0, 1.0, 0.0, 0.0, 0.0}, // 2
343  { 0.0, 0.0, 0.0, 0.5, 0.0, 0.5}, // 3
344  { 0.0, 0.0, 0.0, 0.0, 0.5, 0.5}, // 4
345  { 0.0, 0.0, 0.0, 0.0, 0.0, 1.0} // 5
346  },
347 
348  // embedding matrix for child 3
349  {
350  // 0 1 2 3 4 5 th parent Node
351  { 0.5, 0.5, 0.0, 0.0, 0.0, 0.0}, // 0th child N.
352  { 0.0, 0.5, 0.5, 0.0, 0.0, 0.0}, // 1
353  { 0.5, 0.0, 0.5, 0.0, 0.0, 0.0}, // 2
354  { 0.0, 0.0, 0.0, 0.5, 0.5, 0.0}, // 3
355  { 0.0, 0.0, 0.0, 0.0, 0.5, 0.5}, // 4
356  { 0.0, 0.0, 0.0, 0.5, 0.0, 0.5} // 5
357  }
358 };
359 
360 
361 
362 #endif
363 
364 } // namespace libMesh
365 
366 #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