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