cell_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 
19 // C++ includes
20 
21 // Local includes
22 #include "libmesh/side.h"
23 #include "libmesh/cell_prism6.h"
24 #include "libmesh/edge_edge2.h"
25 #include "libmesh/face_quad4.h"
26 #include "libmesh/face_tri3.h"
27 
28 namespace libMesh
29 {
30 
31 
32 
33 // ------------------------------------------------------------
34 // Prism6 class static member initializations
35 const unsigned int Prism6::side_nodes_map[5][4] =
36 {
37  {0, 2, 1, 99}, // Side 0
38  {0, 1, 4, 3}, // Side 1
39  {1, 2, 5, 4}, // Side 2
40  {2, 0, 3, 5}, // Side 3
41  {3, 4, 5, 99} // Side 4
42 };
43 
44 const unsigned int Prism6::side_elems_map[5][4] =
45 {
46  {0, 1, 2, 3}, // Side 0
47  {0, 1, 4, 5}, // Side 1
48  {1, 2, 5, 6}, // Side 2
49  {0, 2, 4, 6}, // Side 3
50  {4, 5, 6, 7} // Side 4
51 };
52 
53 const unsigned int Prism6::edge_nodes_map[9][2] =
54 {
55  {0, 1}, // Side 0
56  {1, 2}, // Side 1
57  {0, 2}, // Side 2
58  {0, 3}, // Side 3
59  {1, 4}, // Side 4
60  {2, 5}, // Side 5
61  {3, 4}, // Side 6
62  {4, 5}, // Side 7
63  {3, 5} // Side 8
64 };
65 
66 
67 // ------------------------------------------------------------
68 // Prism6 class member functions
69 
70 bool Prism6::is_vertex(const unsigned int) const
71 {
72  return true;
73 }
74 
75 bool Prism6::is_edge(const unsigned int) const
76 {
77  return false;
78 }
79 
80 bool Prism6::is_face(const unsigned int) const
81 {
82  return false;
83 }
84 
85 bool Prism6::is_node_on_side(const unsigned int n,
86  const unsigned int s) const
87 {
88  libmesh_assert_less (s, n_sides());
89  for (unsigned int i = 0; i != 4; ++i)
90  if (side_nodes_map[s][i] == n)
91  return true;
92  return false;
93 }
94 
95 bool Prism6::is_node_on_edge(const unsigned int n,
96  const unsigned int e) const
97 {
98  libmesh_assert_less (e, n_edges());
99  for (unsigned int i = 0; i != 2; ++i)
100  if (edge_nodes_map[e][i] == n)
101  return true;
102  return false;
103 }
104 
105 
106 
108 {
109  // Make sure z edges are affine
110  Point v = this->point(3) - this->point(0);
111  if (!v.relative_fuzzy_equals(this->point(4) - this->point(1)) ||
112  !v.relative_fuzzy_equals(this->point(5) - this->point(2)))
113  return false;
114  return true;
115 }
116 
117 
118 
119 AutoPtr<Elem> Prism6::build_side (const unsigned int i,
120  bool proxy) const
121 {
122  libmesh_assert_less (i, this->n_sides());
123 
124  if (proxy)
125  {
126  switch(i)
127  {
128  case 0:
129  case 4:
130  {
131  AutoPtr<Elem> face(new Side<Tri3,Prism6>(this,i));
132  return face;
133  }
134 
135  case 1:
136  case 2:
137  case 3:
138  {
139  AutoPtr<Elem> face(new Side<Quad4,Prism6>(this,i));
140  return face;
141  }
142 
143  default:
144  {
145  libmesh_error();
146  }
147  }
148  }
149 
150  else
151  {
152  // Create NULL pointer to be initialized, returned later.
153  AutoPtr<Elem> face(NULL);
154 
155  switch (i)
156  {
157  case 0: // the triangular face at z=-1
158  {
159  face.reset(new Tri3);
160 
161  face->set_node(0) = this->get_node(0);
162  face->set_node(1) = this->get_node(2);
163  face->set_node(2) = this->get_node(1);
164 
165  break;
166  }
167  case 1: // the quad face at y=0
168  {
169  face.reset(new Quad4);
170 
171  face->set_node(0) = this->get_node(0);
172  face->set_node(1) = this->get_node(1);
173  face->set_node(2) = this->get_node(4);
174  face->set_node(3) = this->get_node(3);
175 
176  break;
177  }
178  case 2: // the other quad face
179  {
180  face.reset(new Quad4);
181 
182  face->set_node(0) = this->get_node(1);
183  face->set_node(1) = this->get_node(2);
184  face->set_node(2) = this->get_node(5);
185  face->set_node(3) = this->get_node(4);
186 
187  break;
188  }
189  case 3: // the quad face at x=0
190  {
191  face.reset(new Quad4);
192 
193  face->set_node(0) = this->get_node(2);
194  face->set_node(1) = this->get_node(0);
195  face->set_node(2) = this->get_node(3);
196  face->set_node(3) = this->get_node(5);
197 
198  break;
199  }
200  case 4: // the triangular face at z=1
201  {
202  face.reset(new Tri3);
203 
204  face->set_node(0) = this->get_node(3);
205  face->set_node(1) = this->get_node(4);
206  face->set_node(2) = this->get_node(5);
207 
208  break;
209  }
210  default:
211  {
212  libmesh_error();
213  }
214  }
215 
216  face->subdomain_id() = this->subdomain_id();
217  return face;
218  }
219 
220  // We'll never get here.
221  libmesh_error();
222  AutoPtr<Elem> ap(NULL); return ap;
223 }
224 
225 
226 
227 AutoPtr<Elem> Prism6::build_edge (const unsigned int i) const
228 {
229  libmesh_assert_less (i, this->n_edges());
230 
231  return AutoPtr<Elem>(new SideEdge<Edge2,Prism6>(this,i));
232 }
233 
234 
235 
236 void Prism6::connectivity(const unsigned int libmesh_dbg_var(sc),
237  const IOPackage iop,
238  std::vector<dof_id_type>& conn) const
239 {
241  libmesh_assert_less (sc, this->n_sub_elem());
242  libmesh_assert_not_equal_to (iop, INVALID_IO_PACKAGE);
243 
244  switch (iop)
245  {
246  case TECPLOT:
247  {
248  conn.resize(8);
249  conn[0] = this->node(0)+1;
250  conn[1] = this->node(1)+1;
251  conn[2] = this->node(2)+1;
252  conn[3] = this->node(2)+1;
253  conn[4] = this->node(3)+1;
254  conn[5] = this->node(4)+1;
255  conn[6] = this->node(5)+1;
256  conn[7] = this->node(5)+1;
257  return;
258  }
259 
260  case VTK:
261  {
262  conn.resize(6);
263  conn[0] = this->node(0);
264  conn[1] = this->node(2);
265  conn[2] = this->node(1);
266  conn[3] = this->node(3);
267  conn[4] = this->node(5);
268  conn[5] = this->node(4);
269  return;
270  }
271 
272  default:
273  libmesh_error();
274  }
275 
276  libmesh_error();
277 }
278 
279 
280 
281 #ifdef LIBMESH_ENABLE_AMR
282 
283 const float Prism6::_embedding_matrix[8][6][6] =
284 {
285  // embedding matrix for child 0
286  {
287  // 0 1 2 3 4 5
288  { 1.0, 0.0, 0.0, 0.0, 0.0, 0.0}, // 0
289  { 0.5, 0.5, 0.0, 0.0, 0.0, 0.0}, // 1
290  { 0.5, 0.0, 0.5, 0.0, 0.0, 0.0}, // 2
291  { 0.5, 0.0, 0.0, 0.5, 0.0, 0.0}, // 3
292  { .25, .25, 0.0, .25, .25, 0.0}, // 4
293  { .25, 0.0, .25, .25, 0.0, .25} // 5
294  },
295 
296  // embedding matrix for child 1
297  {
298  // 0 1 2 3 4 5
299  { 0.5, 0.5, 0.0, 0.0, 0.0, 0.0}, // 0
300  { 0.0, 1.0, 0.0, 0.0, 0.0, 0.0}, // 1
301  { 0.0, 0.5, 0.5, 0.0, 0.0, 0.0}, // 2
302  { .25, .25, 0.0, .25, .25, 0.0}, // 3
303  { 0.0, 0.5, 0.0, 0.0, 0.5, 0.0}, // 4
304  { 0.0, .25, .25, 0.0, .25, .25} // 5
305  },
306 
307  // embedding matrix for child 2
308  {
309  // 0 1 2 3 4 5
310  { 0.5, 0.0, 0.5, 0.0, 0.0, 0.0}, // 0
311  { 0.0, 0.5, 0.5, 0.0, 0.0, 0.0}, // 1
312  { 0.0, 0.0, 1.0, 0.0, 0.0, 0.0}, // 2
313  { .25, 0.0, .25, .25, 0.0, .25}, // 3
314  { 0.0, .25, .25, 0.0, .25, .25}, // 4
315  { 0.0, 0.0, 0.5, 0.0, 0.0, 0.5} // 5
316  },
317 
318  // embedding matrix for child 3
319  {
320  // 0 1 2 3 4 5
321  { 0.5, 0.5, 0.0, 0.0, 0.0, 0.0}, // 0
322  { 0.0, 0.5, 0.5, 0.0, 0.0, 0.0}, // 1
323  { 0.5, 0.0, 0.5, 0.0, 0.0, 0.0}, // 2
324  { .25, .25, 0.0, .25, .25, 0.0}, // 3
325  { 0.0, .25, .25, 0.0, .25, .25}, // 4
326  { .25, 0.0, .25, .25, 0.0, .25} // 5
327  },
328 
329  // embedding matrix for child 4
330  {
331  // 0 1 2 3 4 5
332  { 0.5, 0.0, 0.0, 0.5, 0.0, 0.0}, // 0
333  { .25, .25, 0.0, .25, .25, 0.0}, // 1
334  { .25, 0.0, .25, .25, 0.0, .25}, // 2
335  { 0.0, 0.0, 0.0, 1.0, 0.0, 0.0}, // 3
336  { 0.0, 0.0, 0.0, 0.5, 0.5, 0.0}, // 4
337  { 0.0, 0.0, 0.0, 0.5, 0.0, 0.5} // 5
338  },
339 
340  // embedding matrix for child 5
341  {
342  // 0 1 2 3 4 5
343  { .25, .25, 0.0, .25, .25, 0.0}, // 0
344  { 0.0, 0.5, 0.0, 0.0, 0.5, 0.0}, // 1
345  { 0.0, .25, .25, 0.0, .25, .25}, // 2
346  { 0.0, 0.0, 0.0, 0.5, 0.5, 0.0}, // 3
347  { 0.0, 0.0, 0.0, 0.0, 1.0, 0.0}, // 4
348  { 0.0, 0.0, 0.0, 0.0, 0.5, 0.5} // 5
349  },
350 
351  // embedding matrix for child 6
352  {
353  // 0 1 2 3 4 5
354  { .25, 0.0, .25, .25, 0.0, .25}, // 0
355  { 0.0, .25, .25, 0.0, .25, .25}, // 1
356  { 0.0, 0.0, 0.5, 0.0, 0.0, 0.5}, // 2
357  { 0.0, 0.0, 0.0, 0.5, 0.0, 0.5}, // 3
358  { 0.0, 0.0, 0.0, 0.0, 0.5, 0.5}, // 4
359  { 0.0, 0.0, 0.0, 0.0, 0.0, 1.0} // 5
360  },
361 
362  // embedding matrix for child 7
363  {
364  // 0 1 2 3 4 5
365  { .25, .25, 0.0, .25, .25, 0.0}, // 0
366  { 0.0, .25, .25, 0.0, .25, .25}, // 1
367  { .25, 0.0, .25, .25, 0.0, .25}, // 2
368  { 0.0, 0.0, 0.0, 0.5, 0.5, 0.0}, // 3
369  { 0.0, 0.0, 0.0, 0.0, 0.5, 0.5}, // 4
370  { 0.0, 0.0, 0.0, 0.5, 0.0, 0.5} // 5
371  }
372 };
373 
374 #endif
375 
376 
377 
379 {
380  // The volume of the prism is computed by splitting
381  // it into 2 tetrahedra and 3 pyramids with bilinear bases.
382  // Then the volume formulae for the tetrahedron and pyramid
383  // are applied and summed to obtain the prism's volume.
384 
385  static const unsigned char sub_pyr[3][4] =
386  {
387  {0, 1, 4, 3},
388  {1, 2, 5, 4},
389  {0, 3, 5, 2}
390  };
391 
392  static const unsigned char sub_tet[2][3] =
393  {
394  {0, 1, 2},
395  {5, 4, 3}
396  };
397 
398  // The centroid is a convenient point to use
399  // for the apex of all the pyramids.
400  const Point R = this->centroid();
401 
402  // temporary storage for Nodes which form the base of the
403  // subelements
404  Node* base[4];
405 
406  // volume accumulation variable
407  Real vol=0.;
408 
409  // Add up the sub-pyramid volumes
410  for (unsigned int n=0; n<3; ++n)
411  {
412  // Set the nodes of the pyramid base
413  for (unsigned int i=0; i<4; ++i)
414  base[i] = this->_nodes[sub_pyr[n][i]];
415 
416  // Compute diff vectors
417  Point a ( *base[0] - R );
418  Point b ( *base[1] - *base[3] );
419  Point c ( *base[2] - *base[0] );
420  Point d ( *base[3] - *base[0] );
421  Point e ( *base[1] - *base[0] );
422 
423  // Compute pyramid volume
424  Real sub_vol = (1./6.)*(a*(b.cross(c))) + (1./12.)*(c*(d.cross(e)));
425 
426  libmesh_assert (sub_vol>0.);
427 
428  vol += sub_vol;
429  }
430 
431 
432  // Add up the sub-tet volumes
433  for (unsigned int n=0; n<2; ++n)
434  {
435  // Set the nodes of the pyramid base
436  for (unsigned int i=0; i<3; ++i)
437  base[i] = this->_nodes[sub_tet[n][i]];
438 
439  // The volume of a tetrahedron is 1/6 the box product formed
440  // by its base and apex vectors
441  Point a ( R - *base[0] );
442 
443  // b is the vector pointing from 0 to 1
444  Point b ( *base[1] - *base[0] );
445 
446  // c is the vector pointing from 0 to 2
447  Point c ( *base[2] - *base[0] );
448 
449  Real sub_vol = (1.0 / 6.0) * (a * (b.cross(c)));
450 
451  libmesh_assert (sub_vol>0.);
452 
453  vol += sub_vol;
454  }
455 
456 
457  // Done with all sub-volumes, so return
458  return vol;
459 }
460 
461 } // namespace libMesh

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

Hosted By:
SourceForge.net Logo