cell_prism15.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_prism15.h"
24 #include "libmesh/edge_edge3.h"
25 #include "libmesh/face_quad8.h"
26 #include "libmesh/face_tri6.h"
27 
28 namespace libMesh
29 {
30 
31 
32 
33 // ------------------------------------------------------------
34 // Prism15 class static member initializations
35 const unsigned int Prism15::side_nodes_map[5][8] =
36 {
37  {0, 2, 1, 8, 7, 6, 99, 99}, // Side 0
38  {0, 1, 4, 3, 6, 10, 12, 9}, // Side 1
39  {1, 2, 5, 4, 7, 11, 13, 10}, // Side 2
40  {2, 0, 3, 5, 8, 9, 14, 11}, // Side 3
41  {3, 4, 5, 12, 13, 14, 99, 99} // Side 4
42 };
43 
44 const unsigned int Prism15::edge_nodes_map[9][3] =
45 {
46  {0, 1, 6}, // Side 0
47  {1, 2, 7}, // Side 1
48  {0, 2, 8}, // Side 2
49  {0, 3, 9}, // Side 3
50  {1, 4, 10}, // Side 4
51  {2, 5, 11}, // Side 5
52  {3, 4, 12}, // Side 6
53  {4, 5, 13}, // Side 7
54  {3, 5, 14} // Side 8
55 };
56 
57 
58 // ------------------------------------------------------------
59 // Prism15 class member functions
60 
61 bool Prism15::is_vertex(const unsigned int i) const
62 {
63  if (i < 6)
64  return true;
65  return false;
66 }
67 
68 bool Prism15::is_edge(const unsigned int i) const
69 {
70  if (i < 6)
71  return false;
72  return true;
73 }
74 
75 bool Prism15::is_face(const unsigned int) const
76 {
77  return false;
78 }
79 
80 bool Prism15::is_node_on_side(const unsigned int n,
81  const unsigned int s) const
82 {
83  libmesh_assert_less (s, n_sides());
84  for (unsigned int i = 0; i != 8; ++i)
85  if (side_nodes_map[s][i] == n)
86  return true;
87  return false;
88 }
89 
90 bool Prism15::is_node_on_edge(const unsigned int n,
91  const unsigned int e) const
92 {
93  libmesh_assert_less (e, n_edges());
94  for (unsigned int i = 0; i != 3; ++i)
95  if (edge_nodes_map[e][i] == n)
96  return true;
97  return false;
98 }
99 
100 
101 
103 {
104  // Make sure z edges are affine
105  Point v = this->point(3) - this->point(0);
106  if (!v.relative_fuzzy_equals(this->point(4) - this->point(1)) ||
107  !v.relative_fuzzy_equals(this->point(5) - this->point(2)))
108  return false;
109  // Make sure edges are straight
110  v /= 2;
111  if (!v.relative_fuzzy_equals(this->point(9) - this->point(0)) ||
112  !v.relative_fuzzy_equals(this->point(10) - this->point(1)) ||
113  !v.relative_fuzzy_equals(this->point(11) - this->point(2)))
114  return false;
115  v = (this->point(1) - this->point(0))/2;
116  if (!v.relative_fuzzy_equals(this->point(6) - this->point(0)) ||
117  !v.relative_fuzzy_equals(this->point(12) - this->point(3)))
118  return false;
119  v = (this->point(2) - this->point(0))/2;
120  if (!v.relative_fuzzy_equals(this->point(8) - this->point(0)) ||
121  !v.relative_fuzzy_equals(this->point(14) - this->point(3)))
122  return false;
123  v = (this->point(2) - this->point(1))/2;
124  if (!v.relative_fuzzy_equals(this->point(7) - this->point(1)) ||
125  !v.relative_fuzzy_equals(this->point(13) - this->point(4)))
126  return false;
127  return true;
128 }
129 
130 
131 
132 AutoPtr<Elem> Prism15::build_side (const unsigned int i,
133  bool proxy) const
134 {
135  libmesh_assert_less (i, this->n_sides());
136 
137  if (proxy)
138  {
139  switch (i)
140  {
141  case 0: // the triangular face at z=-1
142  case 4:
143  {
144  AutoPtr<Elem> face(new Side<Tri6,Prism15>(this,i));
145  return face;
146  }
147 
148  case 1:
149  case 2:
150  case 3:
151  {
152  AutoPtr<Elem> face(new Side<Quad8,Prism15>(this,i));
153  return face;
154  }
155 
156  default:
157  {
158  libmesh_error();
159  }
160  }
161  }
162 
163  else
164  {
165  // Create NULL pointer to be initialized, returned later.
166  AutoPtr<Elem> face(NULL);
167 
168  switch (i)
169  {
170  case 0: // the triangular face at z=-1
171  {
172  face.reset(new Tri6);
173 
174  face->set_node(0) = this->get_node(0);
175  face->set_node(1) = this->get_node(2);
176  face->set_node(2) = this->get_node(1);
177  face->set_node(3) = this->get_node(8);
178  face->set_node(4) = this->get_node(7);
179  face->set_node(5) = this->get_node(6);
180 
181  break;
182  }
183  case 1: // the quad face at y=0
184  {
185  face.reset(new Quad8);
186 
187  face->set_node(0) = this->get_node(0);
188  face->set_node(1) = this->get_node(1);
189  face->set_node(2) = this->get_node(4);
190  face->set_node(3) = this->get_node(3);
191  face->set_node(4) = this->get_node(6);
192  face->set_node(5) = this->get_node(10);
193  face->set_node(6) = this->get_node(12);
194  face->set_node(7) = this->get_node(9);
195 
196  break;
197  }
198  case 2: // the other quad face
199  {
200  face.reset(new Quad8);
201 
202  face->set_node(0) = this->get_node(1);
203  face->set_node(1) = this->get_node(2);
204  face->set_node(2) = this->get_node(5);
205  face->set_node(3) = this->get_node(4);
206  face->set_node(4) = this->get_node(7);
207  face->set_node(5) = this->get_node(11);
208  face->set_node(6) = this->get_node(13);
209  face->set_node(7) = this->get_node(10);
210 
211  break;
212  }
213  case 3: // the quad face at x=0
214  {
215  face.reset(new Quad8);
216 
217  face->set_node(0) = this->get_node(2);
218  face->set_node(1) = this->get_node(0);
219  face->set_node(2) = this->get_node(3);
220  face->set_node(3) = this->get_node(5);
221  face->set_node(4) = this->get_node(8);
222  face->set_node(5) = this->get_node(9);
223  face->set_node(6) = this->get_node(14);
224  face->set_node(7) = this->get_node(11);
225 
226  break;
227  }
228  case 4: // the triangular face at z=1
229  {
230  face.reset(new Tri6);
231 
232  face->set_node(0) = this->get_node(3);
233  face->set_node(1) = this->get_node(4);
234  face->set_node(2) = this->get_node(5);
235  face->set_node(3) = this->get_node(12);
236  face->set_node(4) = this->get_node(13);
237  face->set_node(5) = this->get_node(14);
238 
239  break;
240  }
241  default:
242  {
243  libmesh_error();
244  }
245  }
246 
247  face->subdomain_id() = this->subdomain_id();
248  return face;
249  }
250 
251  // We'll never get here.
252  libmesh_error();
253  AutoPtr<Elem> ap(NULL); return ap;
254 }
255 
256 
257 AutoPtr<Elem> Prism15::build_edge (const unsigned int i) const
258 {
259  libmesh_assert_less (i, this->n_edges());
260 
261  return AutoPtr<Elem>(new SideEdge<Edge3,Prism15>(this,i));
262 }
263 
264 
265 void Prism15::connectivity(const unsigned int libmesh_dbg_var(sc),
266  const IOPackage iop,
267  std::vector<dof_id_type>& conn) const
268 {
270  libmesh_assert_less (sc, this->n_sub_elem());
271  libmesh_assert_not_equal_to (iop, INVALID_IO_PACKAGE);
272 
273  switch (iop)
274  {
275  case TECPLOT:
276  {
277  conn.resize(8);
278  conn[0] = this->node(0)+1;
279  conn[1] = this->node(1)+1;
280  conn[2] = this->node(2)+1;
281  conn[3] = this->node(2)+1;
282  conn[4] = this->node(3)+1;
283  conn[5] = this->node(4)+1;
284  conn[6] = this->node(5)+1;
285  conn[7] = this->node(5)+1;
286  return;
287  }
288 
289  case VTK:
290  {
291  /*
292  conn.resize(6);
293  conn[0] = this->node(0);
294  conn[1] = this->node(2);
295  conn[2] = this->node(1);
296  conn[3] = this->node(3);
297  conn[4] = this->node(5);
298  conn[5] = this->node(4);
299  */
300 
301  // VTK's VTK_QUADRATIC_WEDGE first 9 nodes match, then their
302  // middle and top layers of mid-edge nodes are reversed from
303  // LibMesh's.
304  conn.resize(15);
305  for (unsigned i=0; i<9; ++i)
306  conn[i] = this->node(i);
307 
308  // top "ring" of mid-edge nodes
309  conn[9] = this->node(12);
310  conn[10] = this->node(13);
311  conn[11] = this->node(14);
312 
313  // middle "ring" of mid-edge nodes
314  conn[12] = this->node(9);
315  conn[13] = this->node(10);
316  conn[14] = this->node(11);
317 
318 
319  return;
320  }
321 
322  default:
323  libmesh_error();
324  }
325 
326  libmesh_error();
327 
328 }
329 
330 
331 
332 
333 unsigned short int Prism15::second_order_adjacent_vertex (const unsigned int n,
334  const unsigned int v) const
335 {
336  libmesh_assert_greater_equal (n, this->n_vertices());
337  libmesh_assert_less (n, this->n_nodes());
338  libmesh_assert_less (v, 2);
339  return _second_order_adjacent_vertices[n-this->n_vertices()][v];
340 }
341 
342 
343 
344 std::pair<unsigned short int, unsigned short int>
345 Prism15::second_order_child_vertex (const unsigned int n) const
346 {
347  libmesh_assert_greater_equal (n, this->n_vertices());
348  libmesh_assert_less (n, this->n_nodes());
349 
350  return std::pair<unsigned short int, unsigned short int>
353 }
354 
355 } // namespace libMesh

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

Hosted By:
SourceForge.net Logo