Loading [MathJax]/extensions/tex2jax.js
libMesh
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
cell_prism15.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2025 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 // Local includes
20 #include "libmesh/cell_prism15.h"
21 #include "libmesh/edge_edge3.h"
22 #include "libmesh/face_quad8.h"
23 #include "libmesh/face_tri6.h"
24 #include "libmesh/enum_io_package.h"
25 #include "libmesh/enum_order.h"
26 
27 namespace libMesh
28 {
29 
30 
31 
32 // ------------------------------------------------------------
33 // Prism15 class static member initializations
34 const int Prism15::num_nodes;
35 const int Prism15::nodes_per_side;
36 const int Prism15::nodes_per_edge;
37 
39  {
40  {0, 2, 1, 8, 7, 6, 99, 99}, // Side 0
41  {0, 1, 4, 3, 6, 10, 12, 9}, // Side 1
42  {1, 2, 5, 4, 7, 11, 13, 10}, // Side 2
43  {2, 0, 3, 5, 8, 9, 14, 11}, // Side 3
44  {3, 4, 5, 12, 13, 14, 99, 99} // Side 4
45  };
46 
48  {
49  {0, 1, 6}, // Edge 0
50  {1, 2, 7}, // Edge 1
51  {0, 2, 8}, // Edge 2
52  {0, 3, 9}, // Edge 3
53  {1, 4, 10}, // Edge 4
54  {2, 5, 11}, // Edge 5
55  {3, 4, 12}, // Edge 6
56  {4, 5, 13}, // Edge 7
57  {3, 5, 14} // Edge 8
58  };
59 
60 // ------------------------------------------------------------
61 // Prism15 class member functions
62 
63 bool Prism15::is_vertex(const unsigned int i) const
64 {
65  if (i < 6)
66  return true;
67  return false;
68 }
69 
70 bool Prism15::is_edge(const unsigned int i) const
71 {
72  if (i < 6)
73  return false;
74  return true;
75 }
76 
77 bool Prism15::is_face(const unsigned int) const
78 {
79  return false;
80 }
81 
82 bool Prism15::is_node_on_side(const unsigned int n,
83  const unsigned int s) const
84 {
85  libmesh_assert_less (s, n_sides());
86  return std::find(std::begin(side_nodes_map[s]),
87  std::end(side_nodes_map[s]),
88  n) != std::end(side_nodes_map[s]);
89 }
90 
91 std::vector<unsigned int>
92 Prism15::nodes_on_side(const unsigned int s) const
93 {
94  libmesh_assert_less(s, n_sides());
95  auto trim = (s > 0 && s < 4) ? 0 : 2;
96  return {std::begin(side_nodes_map[s]), std::end(side_nodes_map[s]) - trim};
97 }
98 
99 std::vector<unsigned>
100 Prism15::nodes_on_edge(const unsigned int e) const
101 {
102  libmesh_assert_less(e, n_edges());
103  return {std::begin(edge_nodes_map[e]), std::end(edge_nodes_map[e])};
104 }
105 
106 bool Prism15::is_node_on_edge(const unsigned int n,
107  const unsigned int e) const
108 {
109  libmesh_assert_less (e, n_edges());
110  return std::find(std::begin(edge_nodes_map[e]),
111  std::end(edge_nodes_map[e]),
112  n) != std::end(edge_nodes_map[e]);
113 }
114 
115 
116 
118 {
119  // Make sure z edges are affine
120  Point v = this->point(3) - this->point(0);
121  if (!v.relative_fuzzy_equals(this->point(4) - this->point(1), affine_tol) ||
122  !v.relative_fuzzy_equals(this->point(5) - this->point(2), affine_tol))
123  return false;
124  // Make sure edges are straight
125  v /= 2;
126  if (!v.relative_fuzzy_equals(this->point(9) - this->point(0), affine_tol) ||
127  !v.relative_fuzzy_equals(this->point(10) - this->point(1), affine_tol) ||
128  !v.relative_fuzzy_equals(this->point(11) - this->point(2), affine_tol))
129  return false;
130  v = (this->point(1) - this->point(0))/2;
131  if (!v.relative_fuzzy_equals(this->point(6) - this->point(0), affine_tol) ||
132  !v.relative_fuzzy_equals(this->point(12) - this->point(3), affine_tol))
133  return false;
134  v = (this->point(2) - this->point(0))/2;
135  if (!v.relative_fuzzy_equals(this->point(8) - this->point(0), affine_tol) ||
136  !v.relative_fuzzy_equals(this->point(14) - this->point(3), affine_tol))
137  return false;
138  v = (this->point(2) - this->point(1))/2;
139  if (!v.relative_fuzzy_equals(this->point(7) - this->point(1), affine_tol) ||
140  !v.relative_fuzzy_equals(this->point(13) - this->point(4), affine_tol))
141  return false;
142  return true;
143 }
144 
145 
146 
148 {
149  return SECOND;
150 }
151 
152 
153 
154 unsigned int Prism15::local_side_node(unsigned int side,
155  unsigned int side_node) const
156 {
157  libmesh_assert_less (side, this->n_sides());
158 
159  // Never more than 8 nodes per side.
160  libmesh_assert_less(side_node, Prism15::nodes_per_side);
161 
162  // Some sides have 6 nodes.
163  libmesh_assert(!(side==0 || side==4) || side_node < 6);
164 
165  return Prism15::side_nodes_map[side][side_node];
166 }
167 
168 
169 
170 unsigned int Prism15::local_edge_node(unsigned int edge,
171  unsigned int edge_node) const
172 {
173  libmesh_assert_less(edge, this->n_edges());
174  libmesh_assert_less(edge_node, Prism15::nodes_per_edge);
175 
176  return Prism15::edge_nodes_map[edge][edge_node];
177 }
178 
179 
180 
181 std::unique_ptr<Elem> Prism15::build_side_ptr (const unsigned int i)
182 {
183  libmesh_assert_less (i, this->n_sides());
184 
185  std::unique_ptr<Elem> face;
186 
187  switch (i)
188  {
189  case 0: // the triangular face at z=-1
190  case 4: // the triangular face at z=1
191  {
192  face = std::make_unique<Tri6>();
193  break;
194  }
195  case 1: // the quad face at y=0
196  case 2: // the other quad face
197  case 3: // the quad face at x=0
198  {
199  face = std::make_unique<Quad8>();
200  break;
201  }
202  default:
203  libmesh_error_msg("Invalid side i = " << i);
204  }
205 
206  // Set the nodes
207  for (auto n : face->node_index_range())
208  face->set_node(n, this->node_ptr(Prism15::side_nodes_map[i][n]));
209 
210  face->set_interior_parent(this);
211  face->inherit_data_from(*this);
212 
213  return face;
214 }
215 
216 
217 void Prism15::build_side_ptr (std::unique_ptr<Elem> & side,
218  const unsigned int i)
219 {
220  libmesh_assert_less (i, this->n_sides());
221 
222  switch (i)
223  {
224  case 0: // the triangular face at z=-1
225  case 4: // the triangular face at z=1
226  {
227  if (!side.get() || side->type() != TRI6)
228  {
229  side = this->build_side_ptr(i);
230  return;
231  }
232  break;
233  }
234 
235  case 1: // the quad face at y=0
236  case 2: // the other quad face
237  case 3: // the quad face at x=0
238  {
239  if (!side.get() || side->type() != QUAD8)
240  {
241  side = this->build_side_ptr(i);
242  return;
243  }
244  break;
245  }
246 
247  default:
248  libmesh_error_msg("Invalid side i = " << i);
249  }
250 
251  side->inherit_data_from(*this);
252 
253  // Set the nodes
254  for (auto n : side->node_index_range())
255  side->set_node(n, this->node_ptr(Prism15::side_nodes_map[i][n]));
256 }
257 
258 
259 
260 std::unique_ptr<Elem> Prism15::build_edge_ptr (const unsigned int i)
261 {
262  return this->simple_build_edge_ptr<Edge3,Prism15>(i);
263 }
264 
265 
266 
267 void Prism15::build_edge_ptr (std::unique_ptr<Elem> & edge, const unsigned int i)
268 {
269  this->simple_build_edge_ptr<Prism15>(edge, i, EDGE3);
270 }
271 
272 
273 
274 void Prism15::connectivity(const unsigned int libmesh_dbg_var(sc),
275  const IOPackage iop,
276  std::vector<dof_id_type> & conn) const
277 {
279  libmesh_assert_less (sc, this->n_sub_elem());
280  libmesh_assert_not_equal_to (iop, INVALID_IO_PACKAGE);
281 
282  switch (iop)
283  {
284  case TECPLOT:
285  {
286  conn.resize(8);
287  conn[0] = this->node_id(0)+1;
288  conn[1] = this->node_id(1)+1;
289  conn[2] = this->node_id(2)+1;
290  conn[3] = this->node_id(2)+1;
291  conn[4] = this->node_id(3)+1;
292  conn[5] = this->node_id(4)+1;
293  conn[6] = this->node_id(5)+1;
294  conn[7] = this->node_id(5)+1;
295  return;
296  }
297 
298  case VTK:
299  {
300  // VTK's VTK_QUADRATIC_WEDGE first 9 nodes match, then their
301  // middle and top layers of mid-edge nodes are reversed from
302  // LibMesh's.
303  conn.resize(15);
304  for (unsigned i=0; i<9; ++i)
305  conn[i] = this->node_id(i);
306 
307  // top "ring" of mid-edge nodes
308  conn[9] = this->node_id(12);
309  conn[10] = this->node_id(13);
310  conn[11] = this->node_id(14);
311 
312  // middle "ring" of mid-edge nodes
313  conn[12] = this->node_id(9);
314  conn[13] = this->node_id(10);
315  conn[14] = this->node_id(11);
316 
317  return;
318  }
319 
320  default:
321  libmesh_error_msg("Unsupported IO package " << iop);
322  }
323 }
324 
325 
326 
327 
328 unsigned short int Prism15::second_order_adjacent_vertex (const unsigned int n,
329  const unsigned int v) const
330 {
331  libmesh_assert_greater_equal (n, this->n_vertices());
332  libmesh_assert_less (n, this->n_nodes());
333  libmesh_assert_less (v, 2);
334  return _second_order_adjacent_vertices[n-this->n_vertices()][v];
335 }
336 
337 
338 
339 std::pair<unsigned short int, unsigned short int>
340 Prism15::second_order_child_vertex (const unsigned int n) const
341 {
342  libmesh_assert_greater_equal (n, this->n_vertices());
343  libmesh_assert_less (n, this->n_nodes());
344 
345  return std::pair<unsigned short int, unsigned short int>
348 }
349 
350 
351 
353 {
354  // This specialization is good for Lagrange mappings only
355  if (this->mapping_type() != LAGRANGE_MAP)
356  return this->Elem::volume();
357 
358  // Make copies of our points. It makes the subsequent calculations a bit
359  // shorter and avoids dereferencing the same pointer multiple times.
360  Point
361  x0 = point(0), x1 = point(1), x2 = point(2), x3 = point(3), x4 = point(4),
362  x5 = point(5), x6 = point(6), x7 = point(7), x8 = point(8), x9 = point(9),
363  x10 = point(10), x11 = point(11), x12 = point(12), x13 = point(13), x14 = point(14);
364 
365  // Terms are copied directly from a Python script.
366  Point dx_dxi[10] =
367  {
368  -x0 - x1 + x10 + 2*x12 - x3 - x4 + 2*x6 - x9,
369  3*x0/2 + x1/2 + 2*x12 - 3*x3/2 - x4/2 - 2*x6,
370  -x0/2 + x1/2 - x10 - x3/2 + x4/2 + x9,
371  2*x0 - 2*x12 + 2*x13 - 2*x14 + 2*x3 - 2*x6 + 2*x7 - 2*x8,
372  -2*x0 - 2*x12 + 2*x13 - 2*x14 + 2*x3 + 2*x6 - 2*x7 + 2*x8,
373  Point(0,0,0),
374  2*x0 + 2*x1 - 4*x12 + 2*x3 + 2*x4 - 4*x6,
375  -2*x0 - 2*x1 - 4*x12 + 2*x3 + 2*x4 + 4*x6,
376  Point(0,0,0),
377  Point(0,0,0)
378  };
379 
380  Point dx_deta[10] =
381  {
382  -x0 + x11 + 2*x14 - x2 - x3 - x5 + 2*x8 - x9,
383  3*x0/2 + 2*x14 + x2/2 - 3*x3/2 - x5/2 - 2*x8,
384  -x0/2 - x11 + x2/2 - x3/2 + x5/2 + x9,
385  2*x0 - 4*x14 + 2*x2 + 2*x3 + 2*x5 - 4*x8,
386  -2*x0 - 4*x14 - 2*x2 + 2*x3 + 2*x5 + 4*x8,
387  Point(0,0,0),
388  2*x0 - 2*x12 + 2*x13 - 2*x14 + 2*x3 - 2*x6 + 2*x7 - 2*x8,
389  -2*x0 - 2*x12 + 2*x13 - 2*x14 + 2*x3 + 2*x6 - 2*x7 + 2*x8,
390  Point(0,0,0),
391  Point(0,0,0)
392  };
393 
394  Point dx_dzeta[10] =
395  {
396  -x0/2 + x3/2,
397  x0 + x3 - 2*x9,
398  Point(0,0,0),
399  3*x0/2 + 2*x14 + x2/2 - 3*x3/2 - x5/2 - 2*x8,
400  -x0 - 2*x11 + x2 - x3 + x5 + 2*x9,
401  -x0 - 2*x14 - x2 + x3 + x5 + 2*x8,
402  3*x0/2 + x1/2 + 2*x12 - 3*x3/2 - x4/2 - 2*x6,
403  -x0 + x1 - 2*x10 - x3 + x4 + 2*x9,
404  -2*x0 - 2*x12 + 2*x13 - 2*x14 + 2*x3 + 2*x6 - 2*x7 + 2*x8,
405  -x0 - x1 - 2*x12 + x3 + x4 + 2*x6
406  };
407 
408  // The quadrature rule for the Prism15 is a tensor product between a
409  // FOURTH-order TRI3 rule (in xi, eta) and a FIFTH-order EDGE2 rule
410  // in zeta.
411 
412  // Number of points in the 2D quadrature rule.
413  const int N2D = 6;
414 
415  // Parameters of the 2D rule
416  static const Real
417  w1 = 1.1169079483900573284750350421656140e-01_R,
418  w2 = 5.4975871827660933819163162450105264e-02_R,
419  a1 = 4.4594849091596488631832925388305199e-01_R,
420  a2 = 9.1576213509770743459571463402201508e-02_R;
421 
422  // Points and weights of the 2D rule
423  static const Real w2D[N2D] = {w1, w1, w1, w2, w2, w2};
424 
425  // Quadrature point locations raised to powers. xi[0][2] is
426  // quadrature point 0, squared, xi[1][1] is quadrature point 1 to the
427  // first power, etc. This lets us avoid calling std::pow inside the
428  // loops below.
429  static const Real xi[N2D][3] =
430  {
431  // ^0 ^1 ^2
432  { 1., a1, a1*a1},
433  { 1., 1-2*a1, (1-2*a1)*(1-2*a1)},
434  { 1., a1, a1*a1},
435  { 1., a2, a2*a2},
436  { 1., 1-2*a2, (1-2*a2)*(1-2*a2)},
437  { 1., a2, a2*a2}
438  };
439 
440  static const Real eta[N2D][3] =
441  {
442  // ^0 ^1 ^2
443  { 1., a1, a1*a1},
444  { 1., a1, a1*a1},
445  { 1., 1-2*a1, (1-2*a1)*(1-2*a1)},
446  { 1., a2, a2*a2},
447  { 1., a2, a2*a2},
448  { 1., 1-2*a2, (1-2*a2)*(1-2*a2)}
449  };
450 
451  // Number of points in the 1D quadrature rule.
452  const int N1D = 3;
453 
454  // Points and weights of the 1D quadrature rule.
455  static const Real w1D[N1D] = {5./9, 8./9, 5./9};
456 
457  const Real zeta[N1D][3] =
458  {
459  //^0 ^1 ^2
460  { 1., -std::sqrt(15)/5., 15./25},
461  { 1., 0., 0.},
462  { 1., std::sqrt(15)/5., 15./25}
463  };
464 
465  // The integer exponents for each term.
466  static const int exponents[10][3] =
467  {
468  {0, 0, 0},
469  {0, 0, 1},
470  {0, 0, 2},
471  {0, 1, 0},
472  {0, 1, 1},
473  {0, 2, 0},
474  {1, 0, 0},
475  {1, 0, 1},
476  {1, 1, 0},
477  {2, 0, 0}
478  };
479 
480  Real vol = 0.;
481  for (int i=0; i<N2D; ++i)
482  for (int j=0; j<N1D; ++j)
483  {
484  // Compute dx_dxi, dx_deta, dx_dzeta at the current quadrature point.
485  Point dx_dxi_q, dx_deta_q, dx_dzeta_q;
486  for (int c=0; c<10; ++c)
487  {
488  Real coeff =
489  xi[i][exponents[c][0]]*
490  eta[i][exponents[c][1]]*
491  zeta[j][exponents[c][2]];
492 
493  dx_dxi_q += coeff * dx_dxi[c];
494  dx_deta_q += coeff * dx_deta[c];
495  dx_dzeta_q += coeff * dx_dzeta[c];
496  }
497 
498  // Compute scalar triple product, multiply by weight, and accumulate volume.
499  vol += w2D[i] * w1D[j] * triple_product(dx_dxi_q, dx_deta_q, dx_dzeta_q);
500  }
501 
502  return vol;
503 }
504 
505 
506 
507 #ifdef LIBMESH_ENABLE_AMR
508 
510  {
511  // Embedding matrix for child 0
512  {
513  // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
514  { 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, // 0
515  { 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0 }, // 1
516  { 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0 }, // 2
517  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0 }, // 3
518  { -0.25, -0.25, 0, -0.25, -0.25, 0, 0.5, 0, 0, 0.5, 0.5, 0, 0.5, 0, 0 }, // 4
519  { -0.25, 0, -0.25, -0.25, 0, -0.25, 0, 0, 0.5, 0.5, 0, 0.5, 0, 0, 0.5 }, // 5
520  { 0.375, -0.125, 0, 0, 0, 0, 0.75, 0, 0, 0, 0, 0, 0, 0, 0 }, // 6
521  { 0, -0.125, -0.125, 0, 0, 0, 0.5, 0.25, 0.5, 0, 0, 0, 0, 0, 0 }, // 7
522  { 0.375, 0, -0.125, 0, 0, 0, 0, 0, 0.75, 0, 0, 0, 0, 0, 0 }, // 8
523  { 0.375, 0, 0, -0.125, 0, 0, 0, 0, 0, 0.75, 0, 0, 0, 0, 0 }, // 9
524  { -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.75, 0, 0, 0.375, 0.375, 0, 0.25, 0, 0 }, // 10
525  { -0.1875, 0, -0.1875, -0.1875, 0, -0.1875, 0, 0, 0.75, 0.375, 0, 0.375, 0, 0, 0.25 }, // 11
526  { -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.375, 0, 0, 0.75, 0.25, 0, 0.375, 0, 0 }, // 12
527  { -0.25, -0.1875, -0.1875, -0.25, -0.1875, -0.1875, 0.25, 0.125, 0.25, 0.5, 0.25, 0.25, 0.25, 0.125, 0.25 }, // 13
528  { -0.1875, 0, -0.1875, -0.1875, 0, -0.1875, 0, 0, 0.375, 0.75, 0, 0.25, 0, 0, 0.375 } // 14
529  },
530 
531  // Embedding matrix for child 1
532  {
533  // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
534  { 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0 }, // 0
535  { 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, // 1
536  { 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0 }, // 2
537  { -0.25, -0.25, 0, -0.25, -0.25, 0, 0.5, 0, 0, 0.5, 0.5, 0, 0.5, 0, 0 }, // 3
538  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0 }, // 4
539  { 0, -0.25, -0.25, 0, -0.25, -0.25, 0, 0.5, 0, 0, 0.5, 0.5, 0, 0.5, 0 }, // 5
540  { -0.125, 0.375, 0, 0, 0, 0, 0.75, 0, 0, 0, 0, 0, 0, 0, 0 }, // 6
541  { 0, 0.375, -0.125, 0, 0, 0, 0, 0.75, 0, 0, 0, 0, 0, 0, 0 }, // 7
542  { -0.125, 0, -0.125, 0, 0, 0, 0.5, 0.5, 0.25, 0, 0, 0, 0, 0, 0 }, // 8
543  { -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.75, 0, 0, 0.375, 0.375, 0, 0.25, 0, 0 }, // 9
544  { 0, 0.375, 0, 0, -0.125, 0, 0, 0, 0, 0, 0.75, 0, 0, 0, 0 }, // 10
545  { 0, -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.75, 0, 0, 0.375, 0.375, 0, 0.25, 0 }, // 11
546  { -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.375, 0, 0, 0.25, 0.75, 0, 0.375, 0, 0 }, // 12
547  { 0, -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.375, 0, 0, 0.75, 0.25, 0, 0.375, 0 }, // 13
548  { -0.1875, -0.25, -0.1875, -0.1875, -0.25, -0.1875, 0.25, 0.25, 0.125, 0.25, 0.5, 0.25, 0.25, 0.25, 0.125 } // 14
549  },
550 
551  // Embedding matrix for child 2
552  {
553  // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
554  { 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0 }, // 0
555  { 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0 }, // 1
556  { 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, // 2
557  { -0.25, 0, -0.25, -0.25, 0, -0.25, 0, 0, 0.5, 0.5, 0, 0.5, 0, 0, 0.5 }, // 3
558  { 0, -0.25, -0.25, 0, -0.25, -0.25, 0, 0.5, 0, 0, 0.5, 0.5, 0, 0.5, 0 }, // 4
559  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0 }, // 5
560  { -0.125, -0.125, 0, 0, 0, 0, 0.25, 0.5, 0.5, 0, 0, 0, 0, 0, 0 }, // 6
561  { 0, -0.125, 0.375, 0, 0, 0, 0, 0.75, 0, 0, 0, 0, 0, 0, 0 }, // 7
562  { -0.125, 0, 0.375, 0, 0, 0, 0, 0, 0.75, 0, 0, 0, 0, 0, 0 }, // 8
563  { -0.1875, 0, -0.1875, -0.1875, 0, -0.1875, 0, 0, 0.75, 0.375, 0, 0.375, 0, 0, 0.25 }, // 9
564  { 0, -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.75, 0, 0, 0.375, 0.375, 0, 0.25, 0 }, // 10
565  { 0, 0, 0.375, 0, 0, -0.125, 0, 0, 0, 0, 0, 0.75, 0, 0, 0 }, // 11
566  { -0.1875, -0.1875, -0.25, -0.1875, -0.1875, -0.25, 0.125, 0.25, 0.25, 0.25, 0.25, 0.5, 0.125, 0.25, 0.25 }, // 12
567  { 0, -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.375, 0, 0, 0.25, 0.75, 0, 0.375, 0 }, // 13
568  { -0.1875, 0, -0.1875, -0.1875, 0, -0.1875, 0, 0, 0.375, 0.25, 0, 0.75, 0, 0, 0.375 } // 14
569  },
570 
571  // Embedding matrix for child 3
572  {
573  // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
574  { 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0 }, // 0
575  { 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0 }, // 1
576  { 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0 }, // 2
577  { -0.25, -0.25, 0, -0.25, -0.25, 0, 0.5, 0, 0, 0.5, 0.5, 0, 0.5, 0, 0 }, // 3
578  { 0, -0.25, -0.25, 0, -0.25, -0.25, 0, 0.5, 0, 0, 0.5, 0.5, 0, 0.5, 0 }, // 4
579  { -0.25, 0, -0.25, -0.25, 0, -0.25, 0, 0, 0.5, 0.5, 0, 0.5, 0, 0, 0.5 }, // 5
580  { -0.125, 0, -0.125, 0, 0, 0, 0.5, 0.5, 0.25, 0, 0, 0, 0, 0, 0 }, // 6
581  { -0.125, -0.125, 0, 0, 0, 0, 0.25, 0.5, 0.5, 0, 0, 0, 0, 0, 0 }, // 7
582  { 0, -0.125, -0.125, 0, 0, 0, 0.5, 0.25, 0.5, 0, 0, 0, 0, 0, 0 }, // 8
583  { -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.75, 0, 0, 0.375, 0.375, 0, 0.25, 0, 0 }, // 9
584  { 0, -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.75, 0, 0, 0.375, 0.375, 0, 0.25, 0 }, // 10
585  { -0.1875, 0, -0.1875, -0.1875, 0, -0.1875, 0, 0, 0.75, 0.375, 0, 0.375, 0, 0, 0.25 }, // 11
586  { -0.1875, -0.25, -0.1875, -0.1875, -0.25, -0.1875, 0.25, 0.25, 0.125, 0.25, 0.5, 0.25, 0.25, 0.25, 0.125 }, // 12
587  { -0.1875, -0.1875, -0.25, -0.1875, -0.1875, -0.25, 0.125, 0.25, 0.25, 0.25, 0.25, 0.5, 0.125, 0.25, 0.25 }, // 13
588  { -0.25, -0.1875, -0.1875, -0.25, -0.1875, -0.1875, 0.25, 0.125, 0.25, 0.5, 0.25, 0.25, 0.25, 0.125, 0.25 } // 14
589  },
590 
591  // Embedding matrix for child 4
592  {
593  // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
594  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0 }, // 0
595  { -0.25, -0.25, 0, -0.25, -0.25, 0, 0.5, 0, 0, 0.5, 0.5, 0, 0.5, 0, 0 }, // 1
596  { -0.25, 0, -0.25, -0.25, 0, -0.25, 0, 0, 0.5, 0.5, 0, 0.5, 0, 0, 0.5 }, // 2
597  { 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, // 3
598  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0 }, // 4
599  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1 }, // 5
600  { -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.375, 0, 0, 0.75, 0.25, 0, 0.375, 0, 0 }, // 6
601  { -0.25, -0.1875, -0.1875, -0.25, -0.1875, -0.1875, 0.25, 0.125, 0.25, 0.5, 0.25, 0.25, 0.25, 0.125, 0.25 }, // 7
602  { -0.1875, 0, -0.1875, -0.1875, 0, -0.1875, 0, 0, 0.375, 0.75, 0, 0.25, 0, 0, 0.375 }, // 8
603  { -0.125, 0, 0, 0.375, 0, 0, 0, 0, 0, 0.75, 0, 0, 0, 0, 0 }, // 9
604  { -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.25, 0, 0, 0.375, 0.375, 0, 0.75, 0, 0 }, // 10
605  { -0.1875, 0, -0.1875, -0.1875, 0, -0.1875, 0, 0, 0.25, 0.375, 0, 0.375, 0, 0, 0.75 }, // 11
606  { 0, 0, 0, 0.375, -0.125, 0, 0, 0, 0, 0, 0, 0, 0.75, 0, 0 }, // 12
607  { 0, 0, 0, 0, -0.125, -0.125, 0, 0, 0, 0, 0, 0, 0.5, 0.25, 0.5 }, // 13
608  { 0, 0, 0, 0.375, 0, -0.125, 0, 0, 0, 0, 0, 0, 0, 0, 0.75 } // 14
609  },
610 
611  // Embedding matrix for child 5
612  {
613  // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
614  { -0.25, -0.25, 0, -0.25, -0.25, 0, 0.5, 0, 0, 0.5, 0.5, 0, 0.5, 0, 0 }, // 0
615  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0 }, // 1
616  { 0, -0.25, -0.25, 0, -0.25, -0.25, 0, 0.5, 0, 0, 0.5, 0.5, 0, 0.5, 0 }, // 2
617  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0 }, // 3
618  { 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, // 4
619  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0 }, // 5
620  { -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.375, 0, 0, 0.25, 0.75, 0, 0.375, 0, 0 }, // 6
621  { 0, -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.375, 0, 0, 0.75, 0.25, 0, 0.375, 0 }, // 7
622  { -0.1875, -0.25, -0.1875, -0.1875, -0.25, -0.1875, 0.25, 0.25, 0.125, 0.25, 0.5, 0.25, 0.25, 0.25, 0.125 }, // 8
623  { -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.25, 0, 0, 0.375, 0.375, 0, 0.75, 0, 0 }, // 9
624  { 0, -0.125, 0, 0, 0.375, 0, 0, 0, 0, 0, 0.75, 0, 0, 0, 0 }, // 10
625  { 0, -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.25, 0, 0, 0.375, 0.375, 0, 0.75, 0 }, // 11
626  { 0, 0, 0, -0.125, 0.375, 0, 0, 0, 0, 0, 0, 0, 0.75, 0, 0 }, // 12
627  { 0, 0, 0, 0, 0.375, -0.125, 0, 0, 0, 0, 0, 0, 0, 0.75, 0 }, // 13
628  { 0, 0, 0, -0.125, 0, -0.125, 0, 0, 0, 0, 0, 0, 0.5, 0.5, 0.25 } // 14
629  },
630 
631  // Embedding matrix for child 6
632  {
633  // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
634  { -0.25, 0, -0.25, -0.25, 0, -0.25, 0, 0, 0.5, 0.5, 0, 0.5, 0, 0, 0.5 }, // 0
635  { 0, -0.25, -0.25, 0, -0.25, -0.25, 0, 0.5, 0, 0, 0.5, 0.5, 0, 0.5, 0 }, // 1
636  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0 }, // 2
637  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1 }, // 3
638  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0 }, // 4
639  { 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, // 5
640  { -0.1875, -0.1875, -0.25, -0.1875, -0.1875, -0.25, 0.125, 0.25, 0.25, 0.25, 0.25, 0.5, 0.125, 0.25, 0.25 }, // 6
641  { 0, -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.375, 0, 0, 0.25, 0.75, 0, 0.375, 0 }, // 7
642  { -0.1875, 0, -0.1875, -0.1875, 0, -0.1875, 0, 0, 0.375, 0.25, 0, 0.75, 0, 0, 0.375 }, // 8
643  { -0.1875, 0, -0.1875, -0.1875, 0, -0.1875, 0, 0, 0.25, 0.375, 0, 0.375, 0, 0, 0.75 }, // 9
644  { 0, -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.25, 0, 0, 0.375, 0.375, 0, 0.75, 0 }, // 10
645  { 0, 0, -0.125, 0, 0, 0.375, 0, 0, 0, 0, 0, 0.75, 0, 0, 0 }, // 11
646  { 0, 0, 0, -0.125, -0.125, 0, 0, 0, 0, 0, 0, 0, 0.25, 0.5, 0.5 }, // 12
647  { 0, 0, 0, 0, -0.125, 0.375, 0, 0, 0, 0, 0, 0, 0, 0.75, 0 }, // 13
648  { 0, 0, 0, -0.125, 0, 0.375, 0, 0, 0, 0, 0, 0, 0, 0, 0.75 } // 14
649  },
650 
651  // Embedding matrix for child 7
652  {
653  // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
654  { -0.25, -0.25, 0, -0.25, -0.25, 0, 0.5, 0, 0, 0.5, 0.5, 0, 0.5, 0, 0 }, // 0
655  { 0, -0.25, -0.25, 0, -0.25, -0.25, 0, 0.5, 0, 0, 0.5, 0.5, 0, 0.5, 0 }, // 1
656  { -0.25, 0, -0.25, -0.25, 0, -0.25, 0, 0, 0.5, 0.5, 0, 0.5, 0, 0, 0.5 }, // 2
657  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0 }, // 3
658  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0 }, // 4
659  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1 }, // 5
660  { -0.1875, -0.25, -0.1875, -0.1875, -0.25, -0.1875, 0.25, 0.25, 0.125, 0.25, 0.5, 0.25, 0.25, 0.25, 0.125 }, // 6
661  { -0.1875, -0.1875, -0.25, -0.1875, -0.1875, -0.25, 0.125, 0.25, 0.25, 0.25, 0.25, 0.5, 0.125, 0.25, 0.25 }, // 7
662  { -0.25, -0.1875, -0.1875, -0.25, -0.1875, -0.1875, 0.25, 0.125, 0.25, 0.5, 0.25, 0.25, 0.25, 0.125, 0.25 }, // 8
663  { -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.25, 0, 0, 0.375, 0.375, 0, 0.75, 0, 0 }, // 9
664  { 0, -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.25, 0, 0, 0.375, 0.375, 0, 0.75, 0 }, // 10
665  { -0.1875, 0, -0.1875, -0.1875, 0, -0.1875, 0, 0, 0.25, 0.375, 0, 0.375, 0, 0, 0.75 }, // 11
666  { 0, 0, 0, -0.125, 0, -0.125, 0, 0, 0, 0, 0, 0, 0.5, 0.5, 0.25 }, // 12
667  { 0, 0, 0, -0.125, -0.125, 0, 0, 0, 0, 0, 0, 0, 0.25, 0.5, 0.5 }, // 13
668  { 0, 0, 0, 0, -0.125, -0.125, 0, 0, 0, 0, 0, 0, 0.5, 0.25, 0.5 } // 14
669  }
670  };
671 
672 #endif
673 
674 
675 void
676 Prism15::permute(unsigned int perm_num)
677 {
678  libmesh_assert_less (perm_num, 6);
679  const unsigned int side = perm_num % 2;
680  const unsigned int rotate = perm_num / 2;
681 
682  for (unsigned int i = 0; i != rotate; ++i)
683  {
684  swap3nodes(0,1,2);
685  swap3nodes(3,4,5);
686  swap3nodes(6,7,8);
687  swap3nodes(9,10,11);
688  swap3nodes(12,13,14);
689  swap3neighbors(1,2,3);
690  }
691 
692  switch (side) {
693  case 0:
694  break;
695  case 1:
696  swap2nodes(1,3);
697  swap2nodes(0,4);
698  swap2nodes(2,5);
699  swap2nodes(6,12);
700  swap2nodes(9,10);
701  swap2nodes(7,14);
702  swap2nodes(8,13);
703  swap2neighbors(0,4);
704  swap2neighbors(2,3);
705  break;
706  default:
707  libmesh_error();
708  }
709 
710 }
711 
712 
713 void
714 Prism15::flip(BoundaryInfo * boundary_info)
715 {
716  libmesh_assert(boundary_info);
717 
718  swap2nodes(0,1);
719  swap2nodes(3,4);
720  swap2nodes(7,8);
721  swap2nodes(9,10);
722  swap2nodes(13,14);
723  swap2neighbors(2,3);
724  swap2boundarysides(2,3,boundary_info);
725  swap2boundaryedges(0,1,boundary_info);
726  swap2boundaryedges(3,4,boundary_info);
727  swap2boundaryedges(7,8,boundary_info);
728 }
729 
730 
731 ElemType
732 Prism15::side_type (const unsigned int s) const
733 {
734  libmesh_assert_less (s, 5);
735  if (s == 0 || s == 4)
736  return TRI6;
737  return QUAD8;
738 }
739 
740 } // namespace libMesh
ElemType
Defines an enum for geometric element types.
void swap2boundaryedges(unsigned short e1, unsigned short e2, BoundaryInfo *boundary_info) const
Swaps two edges in boundary_info, if it is non-null.
Definition: elem.C:3550
ElemType side_type(const unsigned int s) const override final
Definition: cell_prism15.C:732
virtual bool has_affine_map() const override
Definition: cell_prism15.C:117
Order
defines an enum for polynomial orders.
Definition: enum_order.h:40
Node ** _nodes
Pointers to the nodes we are connected to.
Definition: elem.h:2245
static const int num_edges
Definition: cell_prism.h:74
static const unsigned int edge_nodes_map[num_edges][nodes_per_edge]
This maps the node of the edge to element node numbers.
Definition: cell_prism15.h:234
static const unsigned int side_nodes_map[num_sides][nodes_per_side]
This maps the node of the side to element node numbers.
Definition: cell_prism15.h:228
virtual unsigned int n_nodes() const override
Definition: cell_prism15.h:101
virtual unsigned int n_sides() const override final
Definition: cell_prism.h:86
IOPackage
libMesh interfaces with several different software packages for the purposes of creating, reading, and writing mesh files.
virtual Order default_order() const override
Definition: cell_prism15.C:147
void swap2boundarysides(unsigned short s1, unsigned short s2, BoundaryInfo *boundary_info) const
Swaps two sides in boundary_info, if it is non-null.
Definition: elem.C:3534
virtual std::pair< unsigned short int, unsigned short int > second_order_child_vertex(const unsigned int n) const override
Definition: cell_prism15.C:340
virtual void connectivity(const unsigned int sc, const IOPackage iop, std::vector< dof_id_type > &conn) const override
Definition: cell_prism15.C:274
virtual bool is_node_on_edge(const unsigned int n, const unsigned int e) const override
Definition: cell_prism15.C:106
The libMesh namespace provides an interface to certain functionality in the library.
static const int num_children
Definition: cell_prism.h:75
virtual bool is_vertex(const unsigned int i) const override
Definition: cell_prism15.C:63
virtual bool is_edge(const unsigned int i) const override
Definition: cell_prism15.C:70
virtual void permute(unsigned int perm_num) override final
Permutes the element (by swapping node and neighbor pointers) according to the specified index...
Definition: cell_prism15.C:676
static const unsigned short int _second_order_adjacent_vertices[9][2]
Matrix that tells which vertices define the location of mid-side (or second-order) nodes...
Definition: cell_prism.h:199
void swap3neighbors(unsigned int n1, unsigned int n2, unsigned int n3)
Swaps three neighbor_ptrs, "rotating" them.
Definition: elem.h:2133
virtual bool is_face(const unsigned int i) const override
Definition: cell_prism15.C:77
virtual unsigned short int second_order_adjacent_vertex(const unsigned int n, const unsigned int v) const override
Definition: cell_prism15.C:328
void swap3nodes(unsigned int n1, unsigned int n2, unsigned int n3)
Swaps three node_ptrs, "rotating" them.
Definition: elem.h:2124
T triple_product(const TypeVector< T > &a, const TypeVector< T > &b, const TypeVector< T > &c)
Definition: type_vector.h:1029
virtual void flip(BoundaryInfo *) override final
Flips the element (by swapping node and neighbor pointers) to have a mapping Jacobian of opposite sig...
Definition: cell_prism15.C:714
virtual unsigned int local_side_node(unsigned int side, unsigned int side_node) const override
Definition: cell_prism15.C:154
virtual std::vector< unsigned int > nodes_on_edge(const unsigned int e) const override
Definition: cell_prism15.C:100
ElemMappingType mapping_type() const
Definition: elem.h:3120
void swap2nodes(unsigned int n1, unsigned int n2)
Swaps two node_ptrs.
Definition: elem.h:2092
virtual unsigned int n_sub_elem() const override
Definition: cell_prism15.h:106
static const unsigned short int _second_order_vertex_child_index[18]
Vector that names the child vertex index for each second order node.
Definition: cell_prism.h:209
virtual std::vector< unsigned int > nodes_on_side(const unsigned int s) const override
Definition: cell_prism15.C:92
The BoundaryInfo class contains information relevant to boundary conditions including storing faces...
Definition: boundary_info.h:57
libmesh_assert(ctx)
static constexpr Real affine_tol
Default tolerance to use in has_affine_map().
Definition: elem.h:2061
virtual unsigned int n_vertices() const override final
Definition: cell_prism.h:91
RealTensorValue rotate(MeshBase &mesh, const Real phi, const Real theta=0., const Real psi=0.)
Rotates the mesh in the xy plane.
static const unsigned short int _second_order_vertex_child_number[18]
Vector that names a child sharing each second order node.
Definition: cell_prism.h:204
void swap2neighbors(unsigned int n1, unsigned int n2)
Swaps two neighbor_ptrs.
Definition: elem.h:2102
virtual bool is_node_on_side(const unsigned int n, const unsigned int s) const override
Definition: cell_prism15.C:82
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual unsigned int n_edges() const override final
Definition: cell_prism.h:96
static const int num_sides
Geometric constants for all Prisms.
Definition: cell_prism.h:73
static const int num_nodes
Geometric constants for Prism15.
Definition: cell_prism15.h:220
virtual Real volume() const
Definition: elem.C:3429
static const int nodes_per_edge
Definition: cell_prism15.h:222
virtual std::unique_ptr< Elem > build_side_ptr(const unsigned int i) override
Builds a QUAD8 or TRI6 built coincident with face i.
Definition: cell_prism15.C:181
static const Real _embedding_matrix[num_children][num_nodes][num_nodes]
Matrix that computes new nodal locations/solution values from current nodes/solution.
Definition: cell_prism15.h:270
static const int nodes_per_side
Definition: cell_prism15.h:221
virtual std::unique_ptr< Elem > build_edge_ptr(const unsigned int i) override
Builds a EDGE3 or INFEDGE2 coincident with edge i.
Definition: cell_prism15.C:260
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
dof_id_type node_id(const unsigned int i) const
Definition: elem.h:2475
const Point & point(const unsigned int i) const
Definition: elem.h:2453
bool relative_fuzzy_equals(const TypeVector< T > &rhs, Real tol=TOLERANCE) const
Definition: type_vector.h:981
virtual Real volume() const override
A specialization for computing the volume of a Prism15.
Definition: cell_prism15.C:352
virtual unsigned int local_edge_node(unsigned int edge, unsigned int edge_node) const override
Definition: cell_prism15.C:170