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