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