cell_pyramid5.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_pyramid5.h"
24 #include "libmesh/edge_edge2.h"
25 #include "libmesh/face_tri3.h"
26 #include "libmesh/face_quad4.h"
27 
28 namespace libMesh
29 {
30 
31 
32 
33 
34 // ------------------------------------------------------------
35 // Pyramid5 class static member initializations
36 const unsigned int Pyramid5::side_nodes_map[5][4] =
37 {
38  {0, 1, 4, 99}, // Side 0
39  {1, 2, 4, 99}, // Side 1
40  {2, 3, 4, 99}, // Side 2
41  {3, 0, 4, 99}, // Side 3
42  {0, 3, 2, 1} // Side 4
43 };
44 
45 const unsigned int Pyramid5::edge_nodes_map[8][2] =
46 {
47  {0, 1}, // Side 0
48  {1, 2}, // Side 1
49  {2, 3}, // Side 2
50  {0, 3}, // Side 3
51  {0, 4}, // Side 4
52  {1, 4}, // Side 5
53  {2, 4}, // Side 6
54  {3, 4} // Side 7
55 };
56 
57 
58 
59 // ------------------------------------------------------------
60 // Pyramid5 class member functions
61 
62 bool Pyramid5::is_vertex(const unsigned int) const
63 {
64  return true;
65 }
66 
67 bool Pyramid5::is_edge(const unsigned int) const
68 {
69  return false;
70 }
71 
72 bool Pyramid5::is_face(const unsigned int) const
73 {
74  return false;
75 }
76 
77 bool Pyramid5::is_node_on_side(const unsigned int n,
78  const unsigned int s) const
79 {
80  libmesh_assert_less (s, n_sides());
81  for (unsigned int i = 0; i != 4; ++i)
82  if (side_nodes_map[s][i] == n)
83  return true;
84  return false;
85 }
86 
87 bool Pyramid5::is_node_on_edge(const unsigned int n,
88  const unsigned int e) const
89 {
90  libmesh_assert_less (e, n_edges());
91  for (unsigned int i = 0; i != 2; ++i)
92  if (edge_nodes_map[e][i] == n)
93  return true;
94  return false;
95 }
96 
97 
98 
100 {
101 // Point v = this->point(3) - this->point(0);
102 // return (v.relative_fuzzy_equals(this->point(2) - this->point(1)));
103  return false;
104 }
105 
106 
107 
108 AutoPtr<Elem> Pyramid5::build_side (const unsigned int i,
109  bool proxy) const
110 {
111  libmesh_assert_less (i, this->n_sides());
112 
113  if (proxy)
114  {
115  switch (i)
116  {
117  case 0:
118  case 1:
119  case 2:
120  case 3:
121  {
122  AutoPtr<Elem> face(new Side<Tri3,Pyramid5>(this,i));
123  return face;
124  }
125 
126  case 4:
127  {
128  AutoPtr<Elem> face(new Side<Quad4,Pyramid5>(this,i));
129  return face;
130  }
131 
132  default:
133  {
134  libmesh_error();
135  }
136  }
137  }
138 
139  else
140  {
141  // Create NULL pointer to be initialized, returned later.
142  AutoPtr<Elem> face(NULL);
143 
144  switch (i)
145  {
146  case 0: // triangular face 1
147  {
148  face.reset(new Tri3);
149 
150  face->set_node(0) = this->get_node(0);
151  face->set_node(1) = this->get_node(1);
152  face->set_node(2) = this->get_node(4);
153 
154  break;
155  }
156  case 1: // triangular face 2
157  {
158  face.reset(new Tri3);
159 
160  face->set_node(0) = this->get_node(1);
161  face->set_node(1) = this->get_node(2);
162  face->set_node(2) = this->get_node(4);
163 
164  break;
165  }
166  case 2: // triangular face 3
167  {
168  face.reset(new Tri3);
169 
170  face->set_node(0) = this->get_node(2);
171  face->set_node(1) = this->get_node(3);
172  face->set_node(2) = this->get_node(4);
173 
174  break;
175  }
176  case 3: // triangular face 4
177  {
178  face.reset(new Tri3);
179 
180  face->set_node(0) = this->get_node(3);
181  face->set_node(1) = this->get_node(0);
182  face->set_node(2) = this->get_node(4);
183 
184  break;
185  }
186  case 4: // the quad face at z=0
187  {
188  face.reset(new Quad4);
189 
190  face->set_node(0) = this->get_node(0);
191  face->set_node(1) = this->get_node(3);
192  face->set_node(2) = this->get_node(2);
193  face->set_node(3) = this->get_node(1);
194 
195  break;
196  }
197  default:
198  {
199  libmesh_error();
200  }
201  }
202 
203  face->subdomain_id() = this->subdomain_id();
204  return face;
205  }
206 
207 
208  // We'll never get here.
209  libmesh_error();
210  AutoPtr<Elem> ap(NULL); return ap;
211 }
212 
213 
214 
215 AutoPtr<Elem> Pyramid5::build_edge (const unsigned int i) const
216 {
217  libmesh_assert_less (i, this->n_edges());
218 
219  return AutoPtr<Elem>(new SideEdge<Edge2,Pyramid5>(this,i));
220 }
221 
222 
223 
224 void Pyramid5::connectivity(const unsigned int libmesh_dbg_var(sc),
225  const IOPackage iop,
226  std::vector<dof_id_type>& conn) const
227 {
229  libmesh_assert_less (sc, this->n_sub_elem());
230  libmesh_assert_not_equal_to (iop, INVALID_IO_PACKAGE);
231 
232  switch (iop)
233  {
234  case TECPLOT:
235  {
236  conn.resize(8);
237  conn[0] = this->node(0)+1;
238  conn[1] = this->node(1)+1;
239  conn[2] = this->node(2)+1;
240  conn[3] = this->node(3)+1;
241  conn[4] = this->node(4)+1;
242  conn[5] = this->node(4)+1;
243  conn[6] = this->node(4)+1;
244  conn[7] = this->node(4)+1;
245  return;
246  }
247 
248  case VTK:
249  {
250  conn.resize(5);
251  conn[0] = this->node(3);
252  conn[1] = this->node(2);
253  conn[2] = this->node(1);
254  conn[3] = this->node(0);
255  conn[4] = this->node(4);
256  return;
257  }
258 
259  default:
260  libmesh_error();
261  }
262 
263  libmesh_error();
264 }
265 
266 
268 {
269  // The pyramid with a bilinear base has volume given by the
270  // formula in: "Calculation of the Volume of a General Hexahedron
271  // for Flow Predictions", AIAA Journal v.23, no.6, 1984, p.954-
272  Node* node0 = this->get_node(0);
273  Node* node1 = this->get_node(1);
274  Node* node2 = this->get_node(2);
275  Node* node3 = this->get_node(3);
276  Node* node4 = this->get_node(4);
277 
278  // Construct Various edge and diagonal vectors
279  Point v40 ( *node0 - *node4 );
280  Point v13 ( *node3 - *node1 );
281  Point v02 ( *node2 - *node0 );
282  Point v03 ( *node3 - *node0 );
283  Point v01 ( *node1 - *node0 );
284 
285  // Finally, ready to return the volume!
286  return (1./6.)*(v40*(v13.cross(v02))) + (1./12.)*(v02*(v01.cross(v03)));
287 }
288 
289 } // namespace libMesh

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

Hosted By:
SourceForge.net Logo