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