libMesh
cell_prism6.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 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}, // Edge 0
56  {1, 2}, // Edge 1
57  {0, 2}, // Edge 2
58  {0, 3}, // Edge 3
59  {1, 4}, // Edge 4
60  {2, 5}, // Edge 5
61  {3, 4}, // Edge 6
62  {4, 5}, // Edge 7
63  {3, 5} // Edge 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 
107 bool Prism6::has_affine_map() const
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 UniquePtr<Elem> Prism6::build_side_ptr (const unsigned int i,
120  bool proxy)
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  return UniquePtr<Elem>(new Side<Tri3,Prism6>(this,i));
131 
132  case 1:
133  case 2:
134  case 3:
135  return UniquePtr<Elem>(new Side<Quad4,Prism6>(this,i));
136 
137  default:
138  libmesh_error_msg("Invalid side i = " << i);
139  }
140  }
141 
142  else
143  {
144  // Create NULL pointer to be initialized, returned later.
145  Elem * face = libmesh_nullptr;
146 
147  switch (i)
148  {
149  case 0: // the triangular face at z=-1
150  case 4: // the triangular face at z=1
151  {
152  face = new Tri3;
153  break;
154  }
155  case 1: // the quad face at y=0
156  case 2: // the other quad face
157  case 3: // the quad face at x=0
158  {
159  face = new Quad4;
160  break;
161  }
162  default:
163  libmesh_error_msg("Invalid side i = " << i);
164  }
165 
166  face->subdomain_id() = this->subdomain_id();
167 
168  // Set the nodes
169  for (unsigned n=0; n<face->n_nodes(); ++n)
170  face->set_node(n) = this->node_ptr(Prism6::side_nodes_map[i][n]);
171 
172  return UniquePtr<Elem>(face);
173  }
174 
175  libmesh_error_msg("We'll never get here!");
176  return UniquePtr<Elem>();
177 }
178 
179 
180 
181 UniquePtr<Elem> Prism6::build_edge_ptr (const unsigned int i)
182 {
183  libmesh_assert_less (i, this->n_edges());
184 
185  return UniquePtr<Elem>(new SideEdge<Edge2,Prism6>(this,i));
186 }
187 
188 
189 
190 void Prism6::connectivity(const unsigned int libmesh_dbg_var(sc),
191  const IOPackage iop,
192  std::vector<dof_id_type> & conn) const
193 {
194  libmesh_assert(_nodes);
195  libmesh_assert_less (sc, this->n_sub_elem());
196  libmesh_assert_not_equal_to (iop, INVALID_IO_PACKAGE);
197 
198  switch (iop)
199  {
200  case TECPLOT:
201  {
202  conn.resize(8);
203  conn[0] = this->node_id(0)+1;
204  conn[1] = this->node_id(1)+1;
205  conn[2] = this->node_id(2)+1;
206  conn[3] = this->node_id(2)+1;
207  conn[4] = this->node_id(3)+1;
208  conn[5] = this->node_id(4)+1;
209  conn[6] = this->node_id(5)+1;
210  conn[7] = this->node_id(5)+1;
211  return;
212  }
213 
214  case VTK:
215  {
216  conn.resize(6);
217  conn[0] = this->node_id(0);
218  conn[1] = this->node_id(2);
219  conn[2] = this->node_id(1);
220  conn[3] = this->node_id(3);
221  conn[4] = this->node_id(5);
222  conn[5] = this->node_id(4);
223  return;
224  }
225 
226  default:
227  libmesh_error_msg("Unsupported IO package " << iop);
228  }
229 }
230 
231 
232 
233 #ifdef LIBMESH_ENABLE_AMR
234 
235 const float Prism6::_embedding_matrix[8][6][6] =
236  {
237  // embedding matrix for child 0
238  {
239  // 0 1 2 3 4 5
240  { 1.0, 0.0, 0.0, 0.0, 0.0, 0.0}, // 0
241  { 0.5, 0.5, 0.0, 0.0, 0.0, 0.0}, // 1
242  { 0.5, 0.0, 0.5, 0.0, 0.0, 0.0}, // 2
243  { 0.5, 0.0, 0.0, 0.5, 0.0, 0.0}, // 3
244  { .25, .25, 0.0, .25, .25, 0.0}, // 4
245  { .25, 0.0, .25, .25, 0.0, .25} // 5
246  },
247 
248  // embedding matrix for child 1
249  {
250  // 0 1 2 3 4 5
251  { 0.5, 0.5, 0.0, 0.0, 0.0, 0.0}, // 0
252  { 0.0, 1.0, 0.0, 0.0, 0.0, 0.0}, // 1
253  { 0.0, 0.5, 0.5, 0.0, 0.0, 0.0}, // 2
254  { .25, .25, 0.0, .25, .25, 0.0}, // 3
255  { 0.0, 0.5, 0.0, 0.0, 0.5, 0.0}, // 4
256  { 0.0, .25, .25, 0.0, .25, .25} // 5
257  },
258 
259  // embedding matrix for child 2
260  {
261  // 0 1 2 3 4 5
262  { 0.5, 0.0, 0.5, 0.0, 0.0, 0.0}, // 0
263  { 0.0, 0.5, 0.5, 0.0, 0.0, 0.0}, // 1
264  { 0.0, 0.0, 1.0, 0.0, 0.0, 0.0}, // 2
265  { .25, 0.0, .25, .25, 0.0, .25}, // 3
266  { 0.0, .25, .25, 0.0, .25, .25}, // 4
267  { 0.0, 0.0, 0.5, 0.0, 0.0, 0.5} // 5
268  },
269 
270  // embedding matrix for child 3
271  {
272  // 0 1 2 3 4 5
273  { 0.5, 0.5, 0.0, 0.0, 0.0, 0.0}, // 0
274  { 0.0, 0.5, 0.5, 0.0, 0.0, 0.0}, // 1
275  { 0.5, 0.0, 0.5, 0.0, 0.0, 0.0}, // 2
276  { .25, .25, 0.0, .25, .25, 0.0}, // 3
277  { 0.0, .25, .25, 0.0, .25, .25}, // 4
278  { .25, 0.0, .25, .25, 0.0, .25} // 5
279  },
280 
281  // embedding matrix for child 4
282  {
283  // 0 1 2 3 4 5
284  { 0.5, 0.0, 0.0, 0.5, 0.0, 0.0}, // 0
285  { .25, .25, 0.0, .25, .25, 0.0}, // 1
286  { .25, 0.0, .25, .25, 0.0, .25}, // 2
287  { 0.0, 0.0, 0.0, 1.0, 0.0, 0.0}, // 3
288  { 0.0, 0.0, 0.0, 0.5, 0.5, 0.0}, // 4
289  { 0.0, 0.0, 0.0, 0.5, 0.0, 0.5} // 5
290  },
291 
292  // embedding matrix for child 5
293  {
294  // 0 1 2 3 4 5
295  { .25, .25, 0.0, .25, .25, 0.0}, // 0
296  { 0.0, 0.5, 0.0, 0.0, 0.5, 0.0}, // 1
297  { 0.0, .25, .25, 0.0, .25, .25}, // 2
298  { 0.0, 0.0, 0.0, 0.5, 0.5, 0.0}, // 3
299  { 0.0, 0.0, 0.0, 0.0, 1.0, 0.0}, // 4
300  { 0.0, 0.0, 0.0, 0.0, 0.5, 0.5} // 5
301  },
302 
303  // embedding matrix for child 6
304  {
305  // 0 1 2 3 4 5
306  { .25, 0.0, .25, .25, 0.0, .25}, // 0
307  { 0.0, .25, .25, 0.0, .25, .25}, // 1
308  { 0.0, 0.0, 0.5, 0.0, 0.0, 0.5}, // 2
309  { 0.0, 0.0, 0.0, 0.5, 0.0, 0.5}, // 3
310  { 0.0, 0.0, 0.0, 0.0, 0.5, 0.5}, // 4
311  { 0.0, 0.0, 0.0, 0.0, 0.0, 1.0} // 5
312  },
313 
314  // embedding matrix for child 7
315  {
316  // 0 1 2 3 4 5
317  { .25, .25, 0.0, .25, .25, 0.0}, // 0
318  { 0.0, .25, .25, 0.0, .25, .25}, // 1
319  { .25, 0.0, .25, .25, 0.0, .25}, // 2
320  { 0.0, 0.0, 0.0, 0.5, 0.5, 0.0}, // 3
321  { 0.0, 0.0, 0.0, 0.0, 0.5, 0.5}, // 4
322  { 0.0, 0.0, 0.0, 0.5, 0.0, 0.5} // 5
323  }
324  };
325 
326 #endif
327 
328 
329 
330 Real Prism6::volume () const
331 {
332  // Make copies of our points. It makes the subsequent calculations a bit
333  // shorter and avoids dereferencing the same pointer multiple times.
334  Point
335  x0 = point(0), x1 = point(1), x2 = point(2),
336  x3 = point(3), x4 = point(4), x5 = point(5);
337 
338  // constant and zeta terms only. These are copied directly from a
339  // Python script.
340  Point dx_dxi[2] =
341  {
342  -x0/2 + x1/2 - x3/2 + x4/2, // constant
343  x0/2 - x1/2 - x3/2 + x4/2, // zeta
344  };
345 
346  // constant and zeta terms only. These are copied directly from a
347  // Python script.
348  Point dx_deta[2] =
349  {
350  -x0/2 + x2/2 - x3/2 + x5/2, // constant
351  x0/2 - x2/2 - x3/2 + x5/2, // zeta
352  };
353 
354  // Constant, xi, and eta terms
355  Point dx_dzeta[3] =
356  {
357  -x0/2 + x3/2, // constant
358  x0/2 - x2/2 - x3/2 + x5/2, // eta
359  x0/2 - x1/2 - x3/2 + x4/2 // xi
360  };
361 
362  // The quadrature rule the Prism6 is a tensor product between a
363  // four-point TRI3 rule (in xi, eta) and a two-point EDGE2 rule (in
364  // zeta) which is capable of integrating cubics exactly.
365 
366  // Number of points in the 2D quadrature rule.
367  const int N2D = 4;
368 
369  static const Real w2D[N2D] =
370  {
371  1.5902069087198858469718450103758e-01L,
372  9.0979309128011415302815498962418e-02L,
373  1.5902069087198858469718450103758e-01L,
374  9.0979309128011415302815498962418e-02L
375  };
376 
377  static const Real xi[N2D] =
378  {
379  1.5505102572168219018027159252941e-01L,
380  6.4494897427831780981972840747059e-01L,
381  1.5505102572168219018027159252941e-01L,
382  6.4494897427831780981972840747059e-01L
383  };
384 
385  static const Real eta[N2D] =
386  {
387  1.7855872826361642311703513337422e-01L,
388  7.5031110222608118177475598324603e-02L,
389  6.6639024601470138670269327409637e-01L,
390  2.8001991549907407200279599420481e-01L
391  };
392 
393  // Number of points in the 1D quadrature rule. The weights of the
394  // 1D quadrature rule are equal to 1.
395  const int N1D = 2;
396 
397  // Points of the 1D quadrature rule
398  static const Real zeta[N1D] =
399  {
400  -std::sqrt(3.)/3,
401  std::sqrt(3.)/3.
402  };
403 
404  Real vol = 0.;
405  for (int i=0; i<N2D; ++i)
406  {
407  // dx_dzeta depends only on the 2D quadrature rule points.
408  Point dx_dzeta_q = dx_dzeta[0] + eta[i]*dx_dzeta[1] + xi[i]*dx_dzeta[2];
409 
410  for (int j=0; j<N1D; ++j)
411  {
412  // dx_dxi and dx_deta only depend on the 1D quadrature rule points.
413  Point
414  dx_dxi_q = dx_dxi[0] + zeta[j]*dx_dxi[1],
415  dx_deta_q = dx_deta[0] + zeta[j]*dx_deta[1];
416 
417  // Compute scalar triple product, multiply by weight, and accumulate volume.
418  vol += w2D[i] * triple_product(dx_dxi_q, dx_deta_q, dx_dzeta_q);
419  }
420  }
421 
422  return vol;
423 }
424 
425 } // namespace libMesh
const class libmesh_nullptr_t libmesh_nullptr
The libMesh namespace provides an interface to certain functionality in the library.
libmesh_assert(j)
T triple_product(const TypeVector< T > &a, const TypeVector< T > &b, const TypeVector< T > &c)
Definition: type_vector.h:1051
PetscErrorCode Vec Mat libmesh_dbg_var(j)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
IOPackage
libMesh interfaces with several different software packages for the purposes of creating, reading, and writing mesh files.