libMesh
cell_prism18.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_prism18.h"
21 #include "libmesh/edge_edge3.h"
22 #include "libmesh/face_quad9.h"
23 #include "libmesh/face_tri6.h"
24 #include "libmesh/enum_io_package.h"
25 #include "libmesh/enum_order.h"
26 #include "libmesh/int_range.h"
27 
28 namespace libMesh
29 {
30 
31 
32 
33 // ------------------------------------------------------------
34 // Prism18 class static member initializations
35 const int Prism18::num_nodes;
36 const int Prism18::nodes_per_side;
37 const int Prism18::nodes_per_edge;
38 
39 const unsigned int Prism18::side_nodes_map[Prism18::num_sides][Prism18::nodes_per_side] =
40  {
41  {0, 2, 1, 8, 7, 6, 99, 99, 99}, // Side 0
42  {0, 1, 4, 3, 6, 10, 12, 9, 15}, // Side 1
43  {1, 2, 5, 4, 7, 11, 13, 10, 16}, // Side 2
44  {2, 0, 3, 5, 8, 9, 14, 11, 17}, // Side 3
45  {3, 4, 5, 12, 13, 14, 99, 99, 99} // Side 4
46  };
47 
48 const unsigned int Prism18::edge_nodes_map[Prism18::num_edges][Prism18::nodes_per_edge] =
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 // Prism18 class member functions
63 
64 bool Prism18::is_vertex(const unsigned int i) const
65 {
66  if (i < 6)
67  return true;
68  return false;
69 }
70 
71 bool Prism18::is_edge(const unsigned int i) const
72 {
73  if (i < 6)
74  return false;
75  if (i > 14)
76  return false;
77  return true;
78 }
79 
80 bool Prism18::is_face(const unsigned int i) const
81 {
82  if (i > 14)
83  return true;
84  return false;
85 }
86 
87 bool Prism18::is_node_on_side(const unsigned int n,
88  const unsigned int s) const
89 {
90  libmesh_assert_less (s, n_sides());
91  return std::find(std::begin(side_nodes_map[s]),
92  std::end(side_nodes_map[s]),
93  n) != std::end(side_nodes_map[s]);
94 }
95 
96 std::vector<unsigned>
97 Prism18::nodes_on_side(const unsigned int s) const
98 {
99  libmesh_assert_less(s, n_sides());
100  auto trim = (s > 0 && s < 4) ? 0 : 3;
101  return {std::begin(side_nodes_map[s]), std::end(side_nodes_map[s]) - trim};
102 }
103 
104 std::vector<unsigned>
105 Prism18::nodes_on_edge(const unsigned int e) const
106 {
107  libmesh_assert_less(e, n_edges());
108  return {std::begin(edge_nodes_map[e]), std::end(edge_nodes_map[e])};
109 }
110 
111 bool Prism18::is_node_on_edge(const unsigned int n,
112  const unsigned int e) const
113 {
114  libmesh_assert_less (e, n_edges());
115  return std::find(std::begin(edge_nodes_map[e]),
116  std::end(edge_nodes_map[e]),
117  n) != std::end(edge_nodes_map[e]);
118 }
119 
120 
121 
122 bool Prism18::has_affine_map() const
123 {
124  // Make sure z edges are affine
125  Point v = this->point(3) - this->point(0);
126  if (!v.relative_fuzzy_equals(this->point(4) - this->point(1), affine_tol) ||
127  !v.relative_fuzzy_equals(this->point(5) - this->point(2), affine_tol))
128  return false;
129  // Make sure edges are straight
130  v /= 2;
131  if (!v.relative_fuzzy_equals(this->point(9) - this->point(0), affine_tol) ||
132  !v.relative_fuzzy_equals(this->point(10) - this->point(1), affine_tol) ||
133  !v.relative_fuzzy_equals(this->point(11) - this->point(2), affine_tol) ||
134  !v.relative_fuzzy_equals(this->point(15) - this->point(6), affine_tol) ||
135  !v.relative_fuzzy_equals(this->point(16) - this->point(7), affine_tol) ||
136  !v.relative_fuzzy_equals(this->point(17) - this->point(8), affine_tol))
137  return false;
138  v = (this->point(1) - this->point(0))/2;
139  if (!v.relative_fuzzy_equals(this->point(6) - this->point(0), affine_tol) ||
140  !v.relative_fuzzy_equals(this->point(12) - this->point(3), affine_tol))
141  return false;
142  v = (this->point(2) - this->point(0))/2;
143  if (!v.relative_fuzzy_equals(this->point(8) - this->point(0), affine_tol) ||
144  !v.relative_fuzzy_equals(this->point(14) - this->point(3), affine_tol))
145  return false;
146  v = (this->point(2) - this->point(1))/2;
147  if (!v.relative_fuzzy_equals(this->point(7) - this->point(1), affine_tol) ||
148  !v.relative_fuzzy_equals(this->point(13) - this->point(4), affine_tol))
149  return false;
150  return true;
151 }
152 
153 
154 
155 Order Prism18::default_order() const
156 {
157  return SECOND;
158 }
159 
160 dof_id_type Prism18::key (const unsigned int s) const
161 {
162  libmesh_assert_less (s, this->n_sides());
163 
164  switch (s)
165  {
166  case 0: // the triangular face at z=0
167  {
168  return Prism::key(0);
169  }
170  case 1: // the quad face at y=0
171  {
172  return Elem::compute_key (this->node_id(15));
173  }
174  case 2: // the other quad face
175  {
176  return Elem::compute_key (this->node_id(16));
177  }
178  case 3: // the quad face at x=0
179  {
180  return Elem::compute_key (this->node_id(17));
181  }
182  case 4: // the triangular face at z=1
183  {
184  return Prism::key(4);
185  }
186  default:
187  libmesh_error_msg("Invalid side " << s);
188  }
189 }
190 
191 
192 
193 unsigned int Prism18::local_side_node(unsigned int side,
194  unsigned int side_node) const
195 {
196  libmesh_assert_less (side, this->n_sides());
197 
198  // Never more than 9 nodes per side.
199  libmesh_assert_less(side_node, Prism18::nodes_per_side);
200 
201  // Some sides have 6 nodes.
202  libmesh_assert(!(side==0 || side==4) || side_node < 6);
203 
204  return Prism18::side_nodes_map[side][side_node];
205 }
206 
207 
208 
209 unsigned int Prism18::local_edge_node(unsigned int edge,
210  unsigned int edge_node) const
211 {
212  libmesh_assert_less(edge, this->n_edges());
213  libmesh_assert_less(edge_node, Prism18::nodes_per_edge);
214 
215  return Prism18::edge_nodes_map[edge][edge_node];
216 }
217 
218 
219 
220 std::unique_ptr<Elem> Prism18::build_side_ptr (const unsigned int i)
221 {
222  libmesh_assert_less (i, this->n_sides());
223 
224  std::unique_ptr<Elem> face;
225 
226  switch (i)
227  {
228  case 0: // the triangular face at z=-1
229  case 4: // the triangular face at z=1
230  {
231  face = std::make_unique<Tri6>();
232  break;
233  }
234  case 1: // the quad face at y=0
235  case 2: // the other quad face
236  case 3: // the quad face at x=0
237  {
238  face = std::make_unique<Quad9>();
239  break;
240  }
241  default:
242  libmesh_error_msg("Invalid side i = " << i);
243  }
244 
245  // Set the nodes
246  for (auto n : face->node_index_range())
247  face->set_node(n, this->node_ptr(Prism18::side_nodes_map[i][n]));
248 
249  face->set_interior_parent(this);
250  face->inherit_data_from(*this);
251 
252  return face;
253 }
254 
255 
256 
257 void Prism18::build_side_ptr (std::unique_ptr<Elem> & side,
258  const unsigned int i)
259 {
260  libmesh_assert_less (i, this->n_sides());
261 
262  switch (i)
263  {
264  case 0: // the triangular face at z=-1
265  case 4: // the triangular face at z=1
266  {
267  if (!side.get() || side->type() != TRI6)
268  {
269  side = this->build_side_ptr(i);
270  return;
271  }
272  break;
273  }
274 
275  case 1: // the quad face at y=0
276  case 2: // the other quad face
277  case 3: // the quad face at x=0
278  {
279  if (!side.get() || side->type() != QUAD9)
280  {
281  side = this->build_side_ptr(i);
282  return;
283  }
284  break;
285  }
286 
287  default:
288  libmesh_error_msg("Invalid side i = " << i);
289  }
290 
291  side->inherit_data_from(*this);
292 
293  // Set the nodes
294  for (auto n : side->node_index_range())
295  side->set_node(n, this->node_ptr(Prism18::side_nodes_map[i][n]));
296 }
297 
298 
299 
300 std::unique_ptr<Elem> Prism18::build_edge_ptr (const unsigned int i)
301 {
302  return this->simple_build_edge_ptr<Edge3,Prism18>(i);
303 }
304 
305 
306 
307 void Prism18::build_edge_ptr (std::unique_ptr<Elem> & edge, const unsigned int i)
308 {
309  this->simple_build_edge_ptr<Prism18>(edge, i, EDGE3);
310 }
311 
312 
313 
314 void Prism18::connectivity(const unsigned int sc,
315  const IOPackage iop,
316  std::vector<dof_id_type> & conn) const
317 {
318  libmesh_assert(_nodes);
319  libmesh_assert_less (sc, this->n_sub_elem());
320  libmesh_assert_not_equal_to (iop, INVALID_IO_PACKAGE);
321 
322  switch (iop)
323  {
324  case TECPLOT:
325  {
326  conn.resize(8);
327  switch (sc)
328  {
329 
330  case 0:
331  {
332  conn[0] = this->node_id(0)+1;
333  conn[1] = this->node_id(6)+1;
334  conn[2] = this->node_id(8)+1;
335  conn[3] = this->node_id(8)+1;
336  conn[4] = this->node_id(9)+1;
337  conn[5] = this->node_id(15)+1;
338  conn[6] = this->node_id(17)+1;
339  conn[7] = this->node_id(17)+1;
340 
341  return;
342  }
343 
344  case 1:
345  {
346  conn[0] = this->node_id(6)+1;
347  conn[1] = this->node_id(1)+1;
348  conn[2] = this->node_id(7)+1;
349  conn[3] = this->node_id(7)+1;
350  conn[4] = this->node_id(15)+1;
351  conn[5] = this->node_id(10)+1;
352  conn[6] = this->node_id(16)+1;
353  conn[7] = this->node_id(16)+1;
354 
355  return;
356  }
357 
358  case 2:
359  {
360  conn[0] = this->node_id(8)+1;
361  conn[1] = this->node_id(7)+1;
362  conn[2] = this->node_id(2)+1;
363  conn[3] = this->node_id(2)+1;
364  conn[4] = this->node_id(17)+1;
365  conn[5] = this->node_id(16)+1;
366  conn[6] = this->node_id(11)+1;
367  conn[7] = this->node_id(11)+1;
368 
369  return;
370  }
371 
372  case 3:
373  {
374  conn[0] = this->node_id(6)+1;
375  conn[1] = this->node_id(7)+1;
376  conn[2] = this->node_id(8)+1;
377  conn[3] = this->node_id(8)+1;
378  conn[4] = this->node_id(15)+1;
379  conn[5] = this->node_id(16)+1;
380  conn[6] = this->node_id(17)+1;
381  conn[7] = this->node_id(17)+1;
382 
383  return;
384  }
385 
386  case 4:
387  {
388  conn[0] = this->node_id(9)+1;
389  conn[1] = this->node_id(15)+1;
390  conn[2] = this->node_id(17)+1;
391  conn[3] = this->node_id(17)+1;
392  conn[4] = this->node_id(3)+1;
393  conn[5] = this->node_id(12)+1;
394  conn[6] = this->node_id(14)+1;
395  conn[7] = this->node_id(14)+1;
396 
397  return;
398  }
399 
400  case 5:
401  {
402  conn[0] = this->node_id(15)+1;
403  conn[1] = this->node_id(10)+1;
404  conn[2] = this->node_id(16)+1;
405  conn[3] = this->node_id(16)+1;
406  conn[4] = this->node_id(12)+1;
407  conn[5] = this->node_id(4)+1;
408  conn[6] = this->node_id(13)+1;
409  conn[7] = this->node_id(13)+1;
410 
411  return;
412  }
413 
414  case 6:
415  {
416  conn[0] = this->node_id(17)+1;
417  conn[1] = this->node_id(16)+1;
418  conn[2] = this->node_id(11)+1;
419  conn[3] = this->node_id(11)+1;
420  conn[4] = this->node_id(14)+1;
421  conn[5] = this->node_id(13)+1;
422  conn[6] = this->node_id(5)+1;
423  conn[7] = this->node_id(5)+1;
424 
425  return;
426  }
427 
428  case 7:
429  {
430  conn[0] = this->node_id(15)+1;
431  conn[1] = this->node_id(16)+1;
432  conn[2] = this->node_id(17)+1;
433  conn[3] = this->node_id(17)+1;
434  conn[4] = this->node_id(12)+1;
435  conn[5] = this->node_id(13)+1;
436  conn[6] = this->node_id(14)+1;
437  conn[7] = this->node_id(14)+1;
438 
439  return;
440  }
441 
442  default:
443  libmesh_error_msg("Invalid sc = " << sc);
444  }
445 
446  }
447 
448  case VTK:
449  {
450  // VTK now supports VTK_BIQUADRATIC_QUADRATIC_WEDGE directly
451  const unsigned int conn_size = 18;
452  conn.resize(conn_size);
453 
454  // VTK's VTK_BIQUADRATIC_QUADRATIC_WEDGE first 9 (vertex) and
455  // last 3 (mid-face) nodes match. The middle and top layers
456  // of mid-edge nodes are reversed from LibMesh's.
457  for (auto i : index_range(conn))
458  conn[i] = this->node_id(i);
459 
460  // top "ring" of mid-edge nodes
461  conn[9] = this->node_id(12);
462  conn[10] = this->node_id(13);
463  conn[11] = this->node_id(14);
464 
465  // middle "ring" of mid-edge nodes
466  conn[12] = this->node_id(9);
467  conn[13] = this->node_id(10);
468  conn[14] = this->node_id(11);
469 
470  return;
471  }
472 
473  default:
474  libmesh_error_msg("Unsupported IO package " << iop);
475  }
476 }
477 
478 
479 
480 
481 unsigned int Prism18::n_second_order_adjacent_vertices (const unsigned int n) const
482 {
483  switch (n)
484  {
485  case 6:
486  case 7:
487  case 8:
488  case 9:
489  case 10:
490  case 11:
491  case 12:
492  case 13:
493  case 14:
494  return 2;
495 
496  case 15:
497  case 16:
498  case 17:
499  return 4;
500 
501  default:
502  libmesh_error_msg("Invalid node n = " << n);
503  }
504 }
505 
506 
507 
508 
509 
510 unsigned short int Prism18::second_order_adjacent_vertex (const unsigned int n,
511  const unsigned int v) const
512 {
513  libmesh_assert_greater_equal (n, this->n_vertices());
514  libmesh_assert_less (n, this->n_nodes());
515 
516  switch (n)
517  {
518  /*
519  * These nodes are unique to \p Prism18,
520  * let our _remaining_... matrix handle
521  * this.
522  */
523  case 15:
524  case 16:
525  case 17:
526  {
527  libmesh_assert_less (v, 4);
528  return _remaining_second_order_adjacent_vertices[n-15][v];
529  }
530 
531  /*
532  * All other second-order nodes (6,...,14) are
533  * identical with Prism15 and are therefore
534  * delegated to the _second_order matrix of
535  * \p Prism
536  */
537  default:
538  {
539  libmesh_assert_less (v, 2);
540  return _second_order_adjacent_vertices[n-this->n_vertices()][v];
541  }
542 
543  }
544 
545  // libmesh_error_msg("We'll never get here!"); // static checkers agree
546  return static_cast<unsigned short int>(-1);
547 }
548 
549 
550 
551 const unsigned short int Prism18::_remaining_second_order_adjacent_vertices[3][4] =
552  {
553  { 0, 1, 3, 4}, // vertices adjacent to node 15
554  { 1, 2, 4, 5}, // vertices adjacent to node 16
555  { 0, 2, 3, 5} // vertices adjacent to node 17
556  };
557 
558 
559 
560 std::pair<unsigned short int, unsigned short int>
561 Prism18::second_order_child_vertex (const unsigned int n) const
562 {
563  libmesh_assert_greater_equal (n, this->n_vertices());
564  libmesh_assert_less (n, this->n_nodes());
565 
566  return std::pair<unsigned short int, unsigned short int>
567  (_second_order_vertex_child_number[n],
568  _second_order_vertex_child_index[n]);
569 }
570 
571 
572 
573 Real Prism18::volume () const
574 {
575  // This specialization is good for Lagrange mappings only
576  if (this->mapping_type() != LAGRANGE_MAP)
577  return this->Elem::volume();
578 
579  // Make copies of our points. It makes the subsequent calculations a bit
580  // shorter and avoids dereferencing the same pointer multiple times.
581  Point
582  x0 = point(0), x1 = point(1), x2 = point(2), x3 = point(3), x4 = point(4), x5 = point(5),
583  x6 = point(6), x7 = point(7), x8 = point(8), x9 = point(9), x10 = point(10), x11 = point(11),
584  x12 = point(12), x13 = point(13), x14 = point(14), x15 = point(15), x16 = point(16), x17 = point(17);
585 
586  // The number of components in the dx_dxi, dx_deta, and dx_dzeta arrays.
587  const int n_components = 16;
588 
589  // Terms are copied directly from a Python script.
590  Point dx_dxi[n_components] =
591  {
592  -x10 + 4*x15 - 3*x9,
593  3*x0/2 + x1/2 + 2*x12 - 3*x3/2 - x4/2 - 2*x6,
594  -3*x0/2 - x1/2 + x10 + 2*x12 - 4*x15 - 3*x3/2 - x4/2 + 2*x6 + 3*x9,
595  -4*x15 + 4*x16 - 4*x17 + 4*x9,
596  -2*x0 - 2*x12 + 2*x13 - 2*x14 + 2*x3 + 2*x6 - 2*x7 + 2*x8,
597  2*x0 - 2*x12 + 2*x13 - 2*x14 + 4*x15 - 4*x16 + 4*x17 + 2*x3 - 2*x6 + 2*x7 - 2*x8 - 4*x9,
598  Point(0,0,0),
599  Point(0,0,0),
600  4*x10 - 8*x15 + 4*x9,
601  -2*x0 - 2*x1 - 4*x12 + 2*x3 + 2*x4 + 4*x6,
602  2*x0 + 2*x1 - 4*x10 - 4*x12 + 8*x15 + 2*x3 + 2*x4 - 4*x6 - 4*x9,
603  Point(0,0,0),
604  Point(0,0,0),
605  Point(0,0,0),
606  Point(0,0,0),
607  Point(0,0,0)
608  };
609 
610  Point dx_deta[n_components] =
611  {
612  -x11 + 4*x17 - 3*x9,
613  3*x0/2 + 2*x14 + x2/2 - 3*x3/2 - x5/2 - 2*x8,
614  -3*x0/2 + x11 + 2*x14 - 4*x17 - x2/2 - 3*x3/2 - x5/2 + 2*x8 + 3*x9,
615  4*x11 - 8*x17 + 4*x9,
616  -2*x0 - 4*x14 - 2*x2 + 2*x3 + 2*x5 + 4*x8,
617  2*x0 - 4*x11 - 4*x14 + 8*x17 + 2*x2 + 2*x3 + 2*x5 - 4*x8 - 4*x9,
618  Point(0,0,0),
619  Point(0,0,0),
620  -4*x15 + 4*x16 - 4*x17 + 4*x9,
621  -2*x0 - 2*x12 + 2*x13 - 2*x14 + 2*x3 + 2*x6 - 2*x7 + 2*x8,
622  2*x0 - 2*x12 + 2*x13 - 2*x14 + 4*x15 - 4*x16 + 4*x17 + 2*x3 - 2*x6 + 2*x7 - 2*x8 - 4*x9,
623  Point(0,0,0),
624  Point(0,0,0),
625  Point(0,0,0),
626  Point(0,0,0),
627  Point(0,0,0)
628  };
629 
630  Point dx_dzeta[n_components] =
631  {
632  -x0/2 + x3/2,
633  x0 + x3 - 2*x9,
634  Point(0,0,0),
635  3*x0/2 + 2*x14 + x2/2 - 3*x3/2 - x5/2 - 2*x8,
636  -3*x0 + 2*x11 + 4*x14 - 8*x17 - x2 - 3*x3 - x5 + 4*x8 + 6*x9,
637  Point(0,0,0),
638  -x0 - 2*x14 - x2 + x3 + x5 + 2*x8,
639  2*x0 - 4*x11 - 4*x14 + 8*x17 + 2*x2 + 2*x3 + 2*x5 - 4*x8 - 4*x9,
640  3*x0/2 + x1/2 + 2*x12 - 3*x3/2 - x4/2 - 2*x6,
641  -3*x0 - x1 + 2*x10 + 4*x12 - 8*x15 - 3*x3 - x4 + 4*x6 + 6*x9,
642  Point(0,0,0),
643  -2*x0 - 2*x12 + 2*x13 - 2*x14 + 2*x3 + 2*x6 - 2*x7 + 2*x8,
644  4*x0 - 4*x12 + 4*x13 - 4*x14 + 8*x15 - 8*x16 + 8*x17 + 4*x3 - 4*x6 + 4*x7 - 4*x8 - 8*x9,
645  Point(0,0,0),
646  -x0 - x1 - 2*x12 + x3 + x4 + 2*x6,
647  2*x0 + 2*x1 - 4*x10 - 4*x12 + 8*x15 + 2*x3 + 2*x4 - 4*x6 - 4*x9
648  };
649 
650  // The quadrature rule for the Prism18 is a tensor product between a
651  // FOURTH-order TRI rule (in xi, eta) and a FIFTH-order EDGE rule
652  // in zeta.
653 
654  // Number of points in the 2D quadrature rule.
655  const int N2D = 6;
656 
657  // Parameters of the 2D rule
658  static const Real
659  w1 = 1.1169079483900573284750350421656140e-01_R,
660  w2 = 5.4975871827660933819163162450105264e-02_R,
661  a1 = 4.4594849091596488631832925388305199e-01_R,
662  a2 = 9.1576213509770743459571463402201508e-02_R;
663 
664  // Points and weights of the 2D rule
665  static const Real w2D[N2D] = {w1, w1, w1, w2, w2, w2};
666 
667  // Quadrature point locations raised to powers. xi[0][2] is
668  // quadrature point 0, squared, xi[1][1] is quadrature point 1 to the
669  // first power, etc. This lets us avoid calling std::pow inside the
670  // loops below.
671  static const Real xi[N2D][3] =
672  {
673  // ^0 ^1 ^2
674  { 1., a1, a1*a1},
675  { 1., 1-2*a1, (1-2*a1)*(1-2*a1)},
676  { 1., a1, a1*a1},
677  { 1., a2, a2*a2},
678  { 1., 1-2*a2, (1-2*a2)*(1-2*a2)},
679  { 1., a2, a2*a2}
680  };
681 
682  static const Real eta[N2D][3] =
683  {
684  // ^0 ^1 ^2
685  { 1., a1, a1*a1},
686  { 1., a1, a1*a1},
687  { 1., 1-2*a1, (1-2*a1)*(1-2*a1)},
688  { 1., a2, a2*a2},
689  { 1., a2, a2*a2},
690  { 1., 1-2*a2, (1-2*a2)*(1-2*a2)}
691  };
692 
693  // Number of points in the 1D quadrature rule.
694  const int N1D = 3;
695 
696  // Points and weights of the 1D quadrature rule.
697  static const Real w1D[N1D] = {5./9, 8./9, 5./9};
698 
699  const Real zeta[N1D][3] =
700  {
701  //^0 ^1 ^2
702  { 1., -std::sqrt(15)/5., 15./25},
703  { 1., 0., 0.},
704  { 1., std::sqrt(15)/5., 15./25}
705  };
706 
707  // The integer exponents for each term.
708  static const int exponents[n_components][3] =
709  {
710  {0, 0, 0},
711  {0, 0, 1},
712  {0, 0, 2},
713  {0, 1, 0},
714  {0, 1, 1},
715  {0, 1, 2},
716  {0, 2, 0},
717  {0, 2, 1},
718  {1, 0, 0},
719  {1, 0, 1},
720  {1, 0, 2},
721  {1, 1, 0},
722  {1, 1, 1},
723  {1, 2, 0},
724  {2, 0, 0},
725  {2, 0, 1}
726  };
727 
728  Real vol = 0.;
729  for (int i=0; i<N2D; ++i)
730  for (int j=0; j<N1D; ++j)
731  {
732  // Compute dx_dxi, dx_deta, dx_dzeta at the current quadrature point.
733  Point dx_dxi_q, dx_deta_q, dx_dzeta_q;
734  for (int c=0; c<n_components; ++c)
735  {
736  Real coeff =
737  xi[i][exponents[c][0]]*
738  eta[i][exponents[c][1]]*
739  zeta[j][exponents[c][2]];
740 
741  dx_dxi_q += coeff * dx_dxi[c];
742  dx_deta_q += coeff * dx_deta[c];
743  dx_dzeta_q += coeff * dx_dzeta[c];
744  }
745 
746  // Compute scalar triple product, multiply by weight, and accumulate volume.
747  vol += w2D[i] * w1D[j] * triple_product(dx_dxi_q, dx_deta_q, dx_dzeta_q);
748  }
749 
750  return vol;
751 }
752 
753 
754 
755 
756 #ifdef LIBMESH_ENABLE_AMR
757 
758 const Real Prism18::_embedding_matrix[Prism18::num_children][Prism18::num_nodes][Prism18::num_nodes] =
759  {
760  // embedding matrix for child 0
761  {
762  // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
763  { 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}, // 0
764  { 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}, // 1
765  { 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0.}, // 2
766  { 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0.}, // 3
767  { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0.}, // 4
768  { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1.}, // 5
769  { 0.375, -0.125, 0., 0., 0., 0., 0.75, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}, // 6
770  { 0., -0.125, -0.125, 0., 0., 0., 0.5, 0.25, 0.5, 0., 0., 0., 0., 0., 0., 0., 0., 0.}, // 7
771  { 0.375, 0., -0.125, 0., 0., 0., 0., 0., 0.75, 0., 0., 0., 0., 0., 0., 0., 0., 0.}, // 8
772  { 0.375, 0., 0., -0.125, 0., 0., 0., 0., 0., 0.75, 0., 0., 0., 0., 0., 0., 0., 0.}, // 9
773  { 0., 0., 0., 0., 0., 0., 0.375, 0., 0., 0., 0., 0., -0.125, 0., 0., 0.75, 0., 0.}, // 10
774  { 0., 0., 0., 0., 0., 0., 0., 0., 0.375, 0., 0., 0., 0., 0., -0.125, 0., 0., 0.75}, // 11
775  { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.375, -0.125, 0., 0., 0., 0., 0.75, 0., 0.}, // 12
776  { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., -0.125, -0.125, 0., 0., 0., 0.5, 0.25, 0.5}, // 13
777  { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.375, 0., -0.125, 0., 0., 0., 0., 0., 0.75}, // 14
778  { 0.140625,-0.046875, 0.,-0.046875, 0.015625, 0., 0.28125, 0., 0., 0.28125, -0.09375, 0., -0.09375, 0., 0., 0.5625, 0., 0.}, // 15
779  { 0.,-0.046875,-0.046875, 0., 0.015625, 0.015625, 0.1875, 0.09375, 0.1875, 0., -0.09375, -0.09375, -0.0625, -0.03125, -0.0625, 0.375, 0.1875, 0.375}, // 16
780  { 0.140625, 0.,-0.046875,-0.046875, 0., 0.015625, 0., 0., 0.28125, 0.28125, 0., -0.09375, 0., 0., -0.09375, 0., 0., 0.5625} // 17
781  },
782 
783  // embedding matrix for child 1
784  {
785  // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
786  { 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}, // 0
787  { 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}, // 1
788  { 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}, // 2
789  { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0.}, // 3
790  { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0.}, // 4
791  { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0.}, // 5
792  { -0.125, 0.375, 0., 0., 0., 0., 0.75, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}, // 6
793  { 0., 0.375, -0.125, 0., 0., 0., 0., 0.75, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}, // 7
794  { -0.125, 0., -0.125, 0., 0., 0., 0.5, 0.5, 0.25, 0., 0., 0., 0., 0., 0., 0., 0., 0.}, // 8
795  { 0., 0., 0., 0., 0., 0., 0.375, 0., 0., 0., 0., 0., -0.125, 0., 0., 0.75, 0., 0.}, // 9
796  { 0., 0.375, 0., 0., -0.125, 0., 0., 0., 0., 0., 0.75, 0., 0., 0., 0., 0., 0., 0.}, // 10
797  { 0., 0., 0., 0., 0., 0., 0., 0.375, 0., 0., 0., 0., 0., -0.125, 0., 0., 0.75, 0.}, // 11
798  { 0., 0., 0., 0., 0., 0., 0., 0., 0., -0.125, 0.375, 0., 0., 0., 0., 0.75, 0., 0.}, // 12
799  { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.375, -0.125, 0., 0., 0., 0., 0.75, 0.}, // 13
800  { 0., 0., 0., 0., 0., 0., 0., 0., 0., -0.125, 0., -0.125, 0., 0., 0., 0.5, 0.5, 0.25}, // 14
801  {-0.046875, 0.140625, 0., 0.015625,-0.046875, 0., 0.28125, 0., 0., -0.09375, 0.28125, 0., -0.09375, 0., 0., 0.5625, 0., 0.}, // 15
802  { 0., 0.140625,-0.046875, 0.,-0.046875, 0.015625, 0., 0.28125, 0., 0., 0.28125, -0.09375, 0., -0.09375, 0., 0., 0.5625, 0.}, // 16
803  {-0.046875, 0.,-0.046875, 0.015625, 0., 0.015625, 0.1875, 0.1875, 0.09375, -0.09375, 0., -0.09375, -0.0625, -0.0625, -0.03125, 0.375, 0.375, 0.1875} // 17
804  },
805 
806  // embedding matrix for child 2
807  {
808  // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
809  { 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0.}, // 0
810  { 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}, // 1
811  { 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}, // 2
812  { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1.}, // 3
813  { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0.}, // 4
814  { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0.}, // 5
815  { -0.125, -0.125, 0., 0., 0., 0., 0.25, 0.5, 0.5, 0., 0., 0., 0., 0., 0., 0., 0., 0.}, // 6
816  { 0., -0.125, 0.375, 0., 0., 0., 0., 0.75, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}, // 7
817  { -0.125, 0., 0.375, 0., 0., 0., 0., 0., 0.75, 0., 0., 0., 0., 0., 0., 0., 0., 0.}, // 8
818  { 0., 0., 0., 0., 0., 0., 0., 0., 0.375, 0., 0., 0., 0., 0., -0.125, 0., 0., 0.75}, // 9
819  { 0., 0., 0., 0., 0., 0., 0., 0.375, 0., 0., 0., 0., 0., -0.125, 0., 0., 0.75, 0.}, // 10
820  { 0., 0., 0.375, 0., 0., -0.125, 0., 0., 0., 0., 0., 0.75, 0., 0., 0., 0., 0., 0.}, // 11
821  { 0., 0., 0., 0., 0., 0., 0., 0., 0., -0.125, -0.125, 0., 0., 0., 0., 0.25, 0.5, 0.5}, // 12
822  { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., -0.125, 0.375, 0., 0., 0., 0., 0.75, 0.}, // 13
823  { 0., 0., 0., 0., 0., 0., 0., 0., 0., -0.125, 0., 0.375, 0., 0., 0., 0., 0., 0.75}, // 14
824  {-0.046875,-0.046875, 0., 0.015625, 0.015625, 0., 0.09375, 0.1875, 0.1875, -0.09375, -0.09375, 0., -0.03125, -0.0625, -0.0625, 0.1875, 0.375, 0.375}, // 15
825  { 0.,-0.046875, 0.140625, 0., 0.015625,-0.046875, 0., 0.28125, 0., 0., -0.09375, 0.28125, 0., -0.09375, 0., 0., 0.5625, 0.}, // 16
826  {-0.046875, 0., 0.140625, 0.015625, 0.,-0.046875, 0., 0., 0.28125, -0.09375, 0., 0.28125, 0., 0., -0.09375, 0., 0., 0.5625} // 17
827  },
828 
829  // embedding matrix for child 3
830  {
831  // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
832  { 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}, // 0
833  { 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}, // 1
834  { 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0.}, // 2
835  { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0.}, // 3
836  { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0.}, // 4
837  { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1.}, // 5
838  { -0.125, 0., -0.125, 0., 0., 0., 0.5, 0.5, 0.25, 0., 0., 0., 0., 0., 0., 0., 0., 0.}, // 6
839  { -0.125, -0.125, 0., 0., 0., 0., 0.25, 0.5, 0.5, 0., 0., 0., 0., 0., 0., 0., 0., 0.}, // 7
840  { 0., -0.125, -0.125, 0., 0., 0., 0.5, 0.25, 0.5, 0., 0., 0., 0., 0., 0., 0., 0., 0.}, // 8
841  { 0., 0., 0., 0., 0., 0., 0.375, 0., 0., 0., 0., 0., -0.125, 0., 0., 0.75, 0., 0.}, // 9
842  { 0., 0., 0., 0., 0., 0., 0., 0.375, 0., 0., 0., 0., 0., -0.125, 0., 0., 0.75, 0.}, // 10
843  { 0., 0., 0., 0., 0., 0., 0., 0., 0.375, 0., 0., 0., 0., 0., -0.125, 0., 0., 0.75}, // 11
844  { 0., 0., 0., 0., 0., 0., 0., 0., 0., -0.125, 0., -0.125, 0., 0., 0., 0.5, 0.5, 0.25}, // 12
845  { 0., 0., 0., 0., 0., 0., 0., 0., 0., -0.125, -0.125, 0., 0., 0., 0., 0.25, 0.5, 0.5}, // 13
846  { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., -0.125, -0.125, 0., 0., 0., 0.5, 0.25, 0.5}, // 14
847  {-0.046875, 0.,-0.046875, 0.015625, 0., 0.015625, 0.1875, 0.1875, 0.09375, -0.09375, 0., -0.09375, -0.0625, -0.0625, -0.03125, 0.375, 0.375, 0.1875}, // 15
848  {-0.046875,-0.046875, 0., 0.015625, 0.015625, 0., 0.09375, 0.1875, 0.1875, -0.09375, -0.09375, 0., -0.03125, -0.0625, -0.0625, 0.1875, 0.375, 0.375}, // 16
849  { 0.,-0.046875,-0.046875, 0., 0.015625, 0.015625, 0.1875, 0.09375, 0.1875, 0., -0.09375, -0.09375, -0.0625, -0.03125, -0.0625, 0.375, 0.1875, 0.375} // 17
850  },
851 
852  // embedding matrix for child 4
853  {
854  // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
855  { 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0.}, // 0
856  { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0.}, // 1
857  { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1.}, // 2
858  { 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}, // 3
859  { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0.}, // 4
860  { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0.}, // 5
861  { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.375, -0.125, 0., 0., 0., 0., 0.75, 0., 0.}, // 6
862  { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., -0.125, -0.125, 0., 0., 0., 0.5, 0.25, 0.5}, // 7
863  { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.375, 0., -0.125, 0., 0., 0., 0., 0., 0.75}, // 8
864  { -0.125, 0., 0., 0.375, 0., 0., 0., 0., 0., 0.75, 0., 0., 0., 0., 0., 0., 0., 0.}, // 9
865  { 0., 0., 0., 0., 0., 0., -0.125, 0., 0., 0., 0., 0., 0.375, 0., 0., 0.75, 0., 0.}, // 10
866  { 0., 0., 0., 0., 0., 0., 0., 0., -0.125, 0., 0., 0., 0., 0., 0.375, 0., 0., 0.75}, // 11
867  { 0., 0., 0., 0.375, -0.125, 0., 0., 0., 0., 0., 0., 0., 0.75, 0., 0., 0., 0., 0.}, // 12
868  { 0., 0., 0., 0., -0.125, -0.125, 0., 0., 0., 0., 0., 0., 0.5, 0.25, 0.5, 0., 0., 0.}, // 13
869  { 0., 0., 0., 0.375, 0., -0.125, 0., 0., 0., 0., 0., 0., 0., 0., 0.75, 0., 0., 0.}, // 14
870  {-0.046875, 0.015625, 0., 0.140625,-0.046875, 0., -0.09375, 0., 0., 0.28125, -0.09375, 0., 0.28125, 0., 0., 0.5625, 0., 0.}, // 15
871  { 0., 0.015625, 0.015625, 0.,-0.046875,-0.046875, -0.0625, -0.03125, -0.0625, 0., -0.09375, -0.09375, 0.1875, 0.09375, 0.1875, 0.375, 0.1875, 0.375}, // 16
872  {-0.046875, 0., 0.015625, 0.140625, 0.,-0.046875, 0., 0., -0.09375, 0.28125, 0., -0.09375, 0., 0., 0.28125, 0., 0., 0.5625} // 17
873  },
874 
875  // embedding matrix for child 5
876  {
877  // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
878  { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0.}, // 0
879  { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0.}, // 1
880  { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0.}, // 2
881  { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0.}, // 3
882  { 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}, // 4
883  { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0.}, // 5
884  { 0., 0., 0., 0., 0., 0., 0., 0., 0., -0.125, 0.375, 0., 0., 0., 0., 0.75, 0., 0.}, // 6
885  { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.375, -0.125, 0., 0., 0., 0., 0.75, 0.}, // 7
886  { 0., 0., 0., 0., 0., 0., 0., 0., 0., -0.125, 0., -0.125, 0., 0., 0., 0.5, 0.5, 0.25}, // 8
887  { 0., 0., 0., 0., 0., 0., -0.125, 0., 0., 0., 0., 0., 0.375, 0., 0., 0.75, 0., 0.}, // 9
888  { 0., -0.125, 0., 0., 0.375, 0., 0., 0., 0., 0., 0.75, 0., 0., 0., 0., 0., 0., 0.}, // 10
889  { 0., 0., 0., 0., 0., 0., 0., -0.125, 0., 0., 0., 0., 0., 0.375, 0., 0., 0.75, 0.}, // 11
890  { 0., 0., 0., -0.125, 0.375, 0., 0., 0., 0., 0., 0., 0., 0.75, 0., 0., 0., 0., 0.}, // 12
891  { 0., 0., 0., 0., 0.375, -0.125, 0., 0., 0., 0., 0., 0., 0., 0.75, 0., 0., 0., 0.}, // 13
892  { 0., 0., 0., -0.125, 0., -0.125, 0., 0., 0., 0., 0., 0., 0.5, 0.5, 0.25, 0., 0., 0.}, // 14
893  { 0.015625,-0.046875, 0.,-0.046875, 0.140625, 0., -0.09375, 0., 0., -0.09375, 0.28125, 0., 0.28125, 0., 0., 0.5625, 0., 0.}, // 15
894  { 0.,-0.046875, 0.015625, 0., 0.140625,-0.046875, 0., -0.09375, 0., 0., 0.28125, -0.09375, 0., 0.28125, 0., 0., 0.5625, 0.}, // 16
895  { 0.015625, 0., 0.015625,-0.046875, 0.,-0.046875, -0.0625, -0.0625, -0.03125, -0.09375, 0., -0.09375, 0.1875, 0.1875, 0.09375, 0.375, 0.375, 0.1875} // 17
896  },
897 
898  // embedding matrix for child 6
899  {
900  // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
901  { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1.}, // 0
902  { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0.}, // 1
903  { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0.}, // 2
904  { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0.}, // 3
905  { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0.}, // 4
906  { 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}, // 5
907  { 0., 0., 0., 0., 0., 0., 0., 0., 0., -0.125, -0.125, 0., 0., 0., 0., 0.25, 0.5, 0.5}, // 6
908  { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., -0.125, 0.375, 0., 0., 0., 0., 0.75, 0.}, // 7
909  { 0., 0., 0., 0., 0., 0., 0., 0., 0., -0.125, 0., 0.375, 0., 0., 0., 0., 0., 0.75}, // 8
910  { 0., 0., 0., 0., 0., 0., 0., 0., -0.125, 0., 0., 0., 0., 0., 0.375, 0., 0., 0.75}, // 9
911  { 0., 0., 0., 0., 0., 0., 0., -0.125, 0., 0., 0., 0., 0., 0.375, 0., 0., 0.75, 0.}, // 10
912  { 0., 0., -0.125, 0., 0., 0.375, 0., 0., 0., 0., 0., 0.75, 0., 0., 0., 0., 0., 0.}, // 11
913  { 0., 0., 0., -0.125, -0.125, 0., 0., 0., 0., 0., 0., 0., 0.25, 0.5, 0.5, 0., 0., 0.}, // 12
914  { 0., 0., 0., 0., -0.125, 0.375, 0., 0., 0., 0., 0., 0., 0., 0.75, 0., 0., 0., 0.}, // 13
915  { 0., 0., 0., -0.125, 0., 0.375, 0., 0., 0., 0., 0., 0., 0., 0., 0.75, 0., 0., 0.}, // 14
916  { 0.015625, 0.015625, 0.,-0.046875,-0.046875, 0., -0.03125, -0.0625, -0.0625, -0.09375, -0.09375, 0., 0.09375, 0.1875, 0.1875, 0.1875, 0.375, 0.375}, // 15
917  { 0., 0.015625,-0.046875, 0.,-0.046875, 0.140625, 0., -0.09375, 0., 0., -0.09375, 0.28125, 0., 0.28125, 0., 0., 0.5625, 0.}, // 16
918  { 0.015625, 0.,-0.046875,-0.046875, 0., 0.140625, 0., 0., -0.09375, -0.09375, 0., 0.28125, 0., 0., 0.28125, 0., 0., 0.5625} // 17
919  },
920 
921  // embedding matrix for child 7
922  {
923  // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
924  { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0.}, // 0
925  { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0.}, // 1
926  { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1.}, // 2
927  { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0.}, // 3
928  { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0.}, // 4
929  { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0.}, // 5
930  { 0., 0., 0., 0., 0., 0., 0., 0., 0., -0.125, 0., -0.125, 0., 0., 0., 0.5, 0.5, 0.25}, // 6
931  { 0., 0., 0., 0., 0., 0., 0., 0., 0., -0.125, -0.125, 0., 0., 0., 0., 0.25, 0.5, 0.5}, // 7
932  { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., -0.125, -0.125, 0., 0., 0., 0.5, 0.25, 0.5}, // 8
933  { 0., 0., 0., 0., 0., 0., -0.125, 0., 0., 0., 0., 0., 0.375, 0., 0., 0.75, 0., 0.}, // 9
934  { 0., 0., 0., 0., 0., 0., 0., -0.125, 0., 0., 0., 0., 0., 0.375, 0., 0., 0.75, 0.}, // 10
935  { 0., 0., 0., 0., 0., 0., 0., 0., -0.125, 0., 0., 0., 0., 0., 0.375, 0., 0., 0.75}, // 11
936  { 0., 0., 0., -0.125, 0., -0.125, 0., 0., 0., 0., 0., 0., 0.5, 0.5, 0.25, 0., 0., 0.}, // 12
937  { 0., 0., 0., -0.125, -0.125, 0., 0., 0., 0., 0., 0., 0., 0.25, 0.5, 0.5, 0., 0., 0.}, // 13
938  { 0., 0., 0., 0., -0.125, -0.125, 0., 0., 0., 0., 0., 0., 0.5, 0.25, 0.5, 0., 0., 0.}, // 14
939  { 0.015625, 0., 0.015625,-0.046875, 0.,-0.046875, -0.0625, -0.0625, -0.03125, -0.09375, 0., -0.09375, 0.1875, 0.1875, 0.09375, 0.375, 0.375, 0.1875}, // 15
940  { 0.015625, 0.015625, 0.,-0.046875,-0.046875, 0., -0.03125, -0.0625, -0.0625, -0.09375, -0.09375, 0., 0.09375, 0.1875, 0.1875, 0.1875, 0.375, 0.375}, // 16
941  { 0., 0.015625, 0.015625, 0.,-0.046875,-0.046875, -0.0625, -0.03125, -0.0625, 0., -0.09375, -0.09375, 0.1875, 0.09375, 0.1875, 0.375, 0.1875, 0.375} // 17
942  }
943  };
944 
945 #endif
946 
947 
948 void
949 Prism18::permute(unsigned int perm_num)
950 {
951  libmesh_assert_less (perm_num, 6);
952  const unsigned int side = perm_num % 2;
953  const unsigned int rotate = perm_num / 2;
954 
955  for (unsigned int i = 0; i != rotate; ++i)
956  {
957  swap3nodes(0,1,2);
958  swap3nodes(3,4,5);
959  swap3nodes(6,7,8);
960  swap3nodes(9,10,11);
961  swap3nodes(12,13,14);
962  swap3nodes(15,16,17);
963  swap3neighbors(1,2,3);
964  }
965 
966  switch (side) {
967  case 0:
968  break;
969  case 1:
970  swap2nodes(1,3);
971  swap2nodes(0,4);
972  swap2nodes(2,5);
973  swap2nodes(6,12);
974  swap2nodes(9,10);
975  swap2nodes(7,14);
976  swap2nodes(8,13);
977  swap2nodes(16,17);
978  swap2neighbors(0,4);
979  swap2neighbors(2,3);
980  break;
981  default:
982  libmesh_error();
983  }
984 }
985 
986 
987 void
988 Prism18::flip(BoundaryInfo * boundary_info)
989 {
990  libmesh_assert(boundary_info);
991 
992  swap2nodes(0,1);
993  swap2nodes(3,4);
994  swap2nodes(7,8);
995  swap2nodes(9,10);
996  swap2nodes(13,14);
997  swap2nodes(16,17);
998  swap2neighbors(2,3);
999  swap2boundarysides(2,3,boundary_info);
1000  swap2boundaryedges(0,1,boundary_info);
1001  swap2boundaryedges(3,4,boundary_info);
1002  swap2boundaryedges(7,8,boundary_info);
1003 }
1004 
1005 
1006 unsigned int Prism18::center_node_on_side(const unsigned short side) const
1007 {
1008  libmesh_assert_less (side, Prism18::num_sides);
1009  return (side >= 1 && side <= 3) ? side + 14 : invalid_uint;
1010 }
1011 
1012 
1013 ElemType
1014 Prism18::side_type (const unsigned int s) const
1015 {
1016  libmesh_assert_less (s, 5);
1017  if (s == 0 || s == 4)
1018  return TRI6;
1019  return QUAD9;
1020 }
1021 
1022 
1023 } // namespace libMesh
libMesh::Prism18::side_type
ElemType side_type(const unsigned int s) const override final
Definition: cell_prism18.C:1014
libMesh::Prism18::is_node_on_side
virtual bool is_node_on_side(const unsigned int n, const unsigned int s) const override
Definition: cell_prism18.C:87
libMesh::ElemType
ElemType
Defines an enum for geometric element types.
Definition: enum_elem_type.h:33
libMesh::Elem::swap2boundaryedges
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:3565
libMesh::Prism18::_remaining_second_order_adjacent_vertices
static const unsigned short int _remaining_second_order_adjacent_vertices[3][4]
Matrix that tells which vertices define the location of mid-side (or second-order) nodes...
Definition: cell_prism18.h:301
libMesh::Order
Order
defines an enum for polynomial orders.
Definition: enum_order.h:40
libMesh::Elem::_nodes
Node ** _nodes
Pointers to the nodes we are connected to.
Definition: elem.h:2246
libMesh::Prism::num_edges
static const int num_edges
Definition: cell_prism.h:74
libMesh::invalid_uint
const unsigned int invalid_uint
A number which is used quite often to represent an invalid or uninitialized value for an unsigned int...
Definition: libmesh.h:310
libMesh::Prism18::center_node_on_side
unsigned int center_node_on_side(const unsigned short side) const override final
Definition: cell_prism18.C:1006
libMesh::Elem::key
virtual dof_id_type key() const
Definition: elem.C:751
libMesh::Prism::n_sides
virtual unsigned int n_sides() const override final
Definition: cell_prism.h:86
libMesh::IOPackage
IOPackage
libMesh interfaces with several different software packages for the purposes of creating, reading, and writing mesh files.
Definition: enum_io_package.h:37
libMesh::Elem::swap2boundarysides
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:3549
libMesh
The libMesh namespace provides an interface to certain functionality in the library.
libMesh::Prism18::default_order
virtual Order default_order() const override
Definition: cell_prism18.C:155
libMesh::Prism18::num_nodes
static const int num_nodes
Geometric constants for Prism18.
Definition: cell_prism18.h:235
libMesh::Prism18::local_edge_node
virtual unsigned int local_edge_node(unsigned int edge, unsigned int edge_node) const override
Definition: cell_prism18.C:209
libMesh::Prism::num_children
static const int num_children
Definition: cell_prism.h:75
libMesh::Prism18::connectivity
virtual void connectivity(const unsigned int sc, const IOPackage iop, std::vector< dof_id_type > &conn) const override
Definition: cell_prism18.C:314
libMesh::Prism::_second_order_adjacent_vertices
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
libMesh::Elem::swap3neighbors
void swap3neighbors(unsigned int n1, unsigned int n2, unsigned int n3)
Swaps three neighbor_ptrs, "rotating" them.
Definition: elem.h:2134
libMesh::Prism18::_embedding_matrix
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_prism18.h:287
libMesh::Elem::swap3nodes
void swap3nodes(unsigned int n1, unsigned int n2, unsigned int n3)
Swaps three node_ptrs, "rotating" them.
Definition: elem.h:2125
libMesh::triple_product
T triple_product(const TypeVector< T > &a, const TypeVector< T > &b, const TypeVector< T > &c)
Definition: type_vector.h:1029
libMesh::Elem::mapping_type
ElemMappingType mapping_type() const
Definition: elem.h:3121
libMesh::Elem::swap2nodes
void swap2nodes(unsigned int n1, unsigned int n2)
Swaps two node_ptrs.
Definition: elem.h:2093
libMesh::Prism18::is_node_on_edge
virtual bool is_node_on_edge(const unsigned int n, const unsigned int e) const override
Definition: cell_prism18.C:111
libMesh::Prism::_second_order_vertex_child_index
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
libMesh::Prism18::build_side_ptr
virtual std::unique_ptr< Elem > build_side_ptr(const unsigned int i) override
Builds a QUAD9 or TRI6 built coincident with face i.
Definition: cell_prism18.C:220
libMesh::Prism18::n_nodes
virtual unsigned int n_nodes() const override
Definition: cell_prism18.h:101
libMesh::BoundaryInfo
The BoundaryInfo class contains information relevant to boundary conditions including storing faces...
Definition: boundary_info.h:57
libMesh::libmesh_assert
libmesh_assert(ctx)
libMesh::Prism18::build_edge_ptr
virtual std::unique_ptr< Elem > build_edge_ptr(const unsigned int i) override
Builds a EDGE3 or INFEDGE2 built coincident with edge i.
Definition: cell_prism18.C:300
libMesh::Elem::affine_tol
static constexpr Real affine_tol
Default tolerance to use in has_affine_map().
Definition: elem.h:2062
libMesh::Prism18::has_affine_map
virtual bool has_affine_map() const override
Definition: cell_prism18.C:122
libMesh::Prism::n_vertices
virtual unsigned int n_vertices() const override final
Definition: cell_prism.h:91
libMesh::MeshTools::Modification::rotate
RealTensorValue rotate(MeshBase &mesh, const Real phi, const Real theta=0., const Real psi=0.)
Rotates the mesh in the xy plane.
Definition: mesh_modification.C:344
libMesh::Prism::_second_order_vertex_child_number
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
libMesh::Prism18::side_nodes_map
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_prism18.h:243
libMesh::Elem::swap2neighbors
void swap2neighbors(unsigned int n1, unsigned int n2)
Swaps two neighbor_ptrs.
Definition: elem.h:2103
libMesh::Prism18::nodes_on_edge
virtual std::vector< unsigned int > nodes_on_edge(const unsigned int e) const override
Definition: cell_prism18.C:105
libMesh::Prism18::second_order_adjacent_vertex
virtual unsigned short int second_order_adjacent_vertex(const unsigned int n, const unsigned int v) const override
Definition: cell_prism18.C:510
libMesh::Prism18::nodes_per_edge
static const int nodes_per_edge
Definition: cell_prism18.h:237
libMesh::Prism18::flip
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_prism18.C:988
libMesh::Prism18::edge_nodes_map
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_prism18.h:249
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:144
libMesh::Prism::n_edges
virtual unsigned int n_edges() const override final
Definition: cell_prism.h:96
libMesh::Prism18::is_vertex
virtual bool is_vertex(const unsigned int i) const override
Definition: cell_prism18.C:64
libMesh::Prism::num_sides
static const int num_sides
Geometric constants for all Prisms.
Definition: cell_prism.h:73
libMesh::Prism18::local_side_node
virtual unsigned int local_side_node(unsigned int side, unsigned int side_node) const override
Definition: cell_prism18.C:193
libMesh::Prism18::nodes_on_side
virtual std::vector< unsigned int > nodes_on_side(const unsigned int s) const override
Definition: cell_prism18.C:97
libMesh::Elem::volume
virtual Real volume() const
Definition: elem.C:3433
libMesh::Prism18::n_sub_elem
virtual unsigned int n_sub_elem() const override
Definition: cell_prism18.h:106
libMesh::Prism18::volume
virtual Real volume() const override
A specialization for computing the volume of a Prism18.
Definition: cell_prism18.C:573
libMesh::Elem::compute_key
static dof_id_type compute_key(dof_id_type n0)
Definition: elem.h:3295
libMesh::Prism18::nodes_per_side
static const int nodes_per_side
Definition: cell_prism18.h:236
libMesh::Prism18::n_second_order_adjacent_vertices
virtual unsigned int n_second_order_adjacent_vertices(const unsigned int) const override
Definition: cell_prism18.C:481
libMesh::Prism18::permute
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_prism18.C:949
libMesh::Point
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
libMesh::Elem::node_id
dof_id_type node_id(const unsigned int i) const
Definition: elem.h:2476
libMesh::Elem::point
const Point & point(const unsigned int i) const
Definition: elem.h:2454
libMesh::Prism18::second_order_child_vertex
virtual std::pair< unsigned short int, unsigned short int > second_order_child_vertex(const unsigned int n) const override
Definition: cell_prism18.C:561
libMesh::TypeVector::relative_fuzzy_equals
bool relative_fuzzy_equals(const TypeVector< T > &rhs, Real tol=TOLERANCE) const
Definition: type_vector.h:981
libMesh::index_range
auto index_range(const T &sizable)
Helper function that returns an IntRange<std::size_t> representing all the indices of the passed-in v...
Definition: int_range.h:117
libMesh::dof_id_type
uint8_t dof_id_type
Definition: id_types.h:67
libMesh::Prism18::is_edge
virtual bool is_edge(const unsigned int i) const override
Definition: cell_prism18.C:71
libMesh::Prism18::is_face
virtual bool is_face(const unsigned int i) const override
Definition: cell_prism18.C:80
Generated on Thu Oct 30 2025 06:31:51 for libMesh by   doxygen 1.8.14