libMesh
cell_prism15.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_prism15.h"
24 #include "libmesh/edge_edge3.h"
25 #include "libmesh/face_quad8.h"
26 #include "libmesh/face_tri6.h"
27 
28 namespace libMesh
29 {
30 
31 
32 
33 // ------------------------------------------------------------
34 // Prism15 class static member initializations
35 const unsigned int Prism15::side_nodes_map[5][8] =
36  {
37  {0, 2, 1, 8, 7, 6, 99, 99}, // Side 0
38  {0, 1, 4, 3, 6, 10, 12, 9}, // Side 1
39  {1, 2, 5, 4, 7, 11, 13, 10}, // Side 2
40  {2, 0, 3, 5, 8, 9, 14, 11}, // Side 3
41  {3, 4, 5, 12, 13, 14, 99, 99} // Side 4
42  };
43 
44 const unsigned int Prism15::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 // Prism15 class member functions
60 
61 bool Prism15::is_vertex(const unsigned int i) const
62 {
63  if (i < 6)
64  return true;
65  return false;
66 }
67 
68 bool Prism15::is_edge(const unsigned int i) const
69 {
70  if (i < 6)
71  return false;
72  return true;
73 }
74 
75 bool Prism15::is_face(const unsigned int) const
76 {
77  return false;
78 }
79 
80 bool Prism15::is_node_on_side(const unsigned int n,
81  const unsigned int s) const
82 {
83  libmesh_assert_less (s, n_sides());
84  for (unsigned int i = 0; i != 8; ++i)
85  if (side_nodes_map[s][i] == n)
86  return true;
87  return false;
88 }
89 
90 bool Prism15::is_node_on_edge(const unsigned int n,
91  const unsigned int e) const
92 {
93  libmesh_assert_less (e, n_edges());
94  for (unsigned int i = 0; i != 3; ++i)
95  if (edge_nodes_map[e][i] == n)
96  return true;
97  return false;
98 }
99 
100 
101 
102 bool Prism15::has_affine_map() const
103 {
104  // Make sure z edges are affine
105  Point v = this->point(3) - this->point(0);
106  if (!v.relative_fuzzy_equals(this->point(4) - this->point(1)) ||
107  !v.relative_fuzzy_equals(this->point(5) - this->point(2)))
108  return false;
109  // Make sure edges are straight
110  v /= 2;
111  if (!v.relative_fuzzy_equals(this->point(9) - this->point(0)) ||
112  !v.relative_fuzzy_equals(this->point(10) - this->point(1)) ||
113  !v.relative_fuzzy_equals(this->point(11) - this->point(2)))
114  return false;
115  v = (this->point(1) - this->point(0))/2;
116  if (!v.relative_fuzzy_equals(this->point(6) - this->point(0)) ||
117  !v.relative_fuzzy_equals(this->point(12) - this->point(3)))
118  return false;
119  v = (this->point(2) - this->point(0))/2;
120  if (!v.relative_fuzzy_equals(this->point(8) - this->point(0)) ||
121  !v.relative_fuzzy_equals(this->point(14) - this->point(3)))
122  return false;
123  v = (this->point(2) - this->point(1))/2;
124  if (!v.relative_fuzzy_equals(this->point(7) - this->point(1)) ||
125  !v.relative_fuzzy_equals(this->point(13) - this->point(4)))
126  return false;
127  return true;
128 }
129 
130 
131 
132 unsigned int Prism15::which_node_am_i(unsigned int side,
133  unsigned int side_node) const
134 {
135  libmesh_assert_less (side, this->n_sides());
136 
137  // Never more than 8 nodes per side.
138  libmesh_assert_less(side_node, 8);
139 
140  // Some sides have 6 nodes.
141  libmesh_assert(!(side==0 || side==4) || side_node < 6);
142 
143  return Prism15::side_nodes_map[side][side_node];
144 }
145 
146 
147 
148 UniquePtr<Elem> Prism15::build_side_ptr (const unsigned int i,
149  bool proxy)
150 {
151  libmesh_assert_less (i, this->n_sides());
152 
153  if (proxy)
154  {
155  switch (i)
156  {
157  case 0: // the triangular face at z=-1
158  case 4:
159  return UniquePtr<Elem>(new Side<Tri6,Prism15>(this,i));
160 
161  case 1:
162  case 2:
163  case 3:
164  return UniquePtr<Elem>(new Side<Quad8,Prism15>(this,i));
165 
166  default:
167  libmesh_error_msg("Invalid side i = " << i);
168  }
169  }
170 
171  else
172  {
173  // Create NULL pointer to be initialized, returned later.
174  Elem * face = libmesh_nullptr;
175 
176  switch (i)
177  {
178  case 0: // the triangular face at z=-1
179  case 4: // the triangular face at z=1
180  {
181  face = new Tri6;
182  break;
183  }
184  case 1: // the quad face at y=0
185  case 2: // the other quad face
186  case 3: // the quad face at x=0
187  {
188  face = new Quad8;
189  break;
190  }
191  default:
192  libmesh_error_msg("Invalid side i = " << i);
193  }
194 
195  face->subdomain_id() = this->subdomain_id();
196 
197  // Set the nodes
198  for (unsigned n=0; n<face->n_nodes(); ++n)
199  face->set_node(n) = this->node_ptr(Prism15::side_nodes_map[i][n]);
200 
201  return UniquePtr<Elem>(face);
202  }
203 
204  libmesh_error_msg("We'll never get here!");
205  return UniquePtr<Elem>();
206 }
207 
208 
209 UniquePtr<Elem> Prism15::build_edge_ptr (const unsigned int i)
210 {
211  libmesh_assert_less (i, this->n_edges());
212 
213  return UniquePtr<Elem>(new SideEdge<Edge3,Prism15>(this,i));
214 }
215 
216 
217 void Prism15::connectivity(const unsigned int libmesh_dbg_var(sc),
218  const IOPackage iop,
219  std::vector<dof_id_type> & conn) const
220 {
221  libmesh_assert(_nodes);
222  libmesh_assert_less (sc, this->n_sub_elem());
223  libmesh_assert_not_equal_to (iop, INVALID_IO_PACKAGE);
224 
225  switch (iop)
226  {
227  case TECPLOT:
228  {
229  conn.resize(8);
230  conn[0] = this->node_id(0)+1;
231  conn[1] = this->node_id(1)+1;
232  conn[2] = this->node_id(2)+1;
233  conn[3] = this->node_id(2)+1;
234  conn[4] = this->node_id(3)+1;
235  conn[5] = this->node_id(4)+1;
236  conn[6] = this->node_id(5)+1;
237  conn[7] = this->node_id(5)+1;
238  return;
239  }
240 
241  case VTK:
242  {
243  /*
244  conn.resize(6);
245  conn[0] = this->node_id(0);
246  conn[1] = this->node_id(2);
247  conn[2] = this->node_id(1);
248  conn[3] = this->node_id(3);
249  conn[4] = this->node_id(5);
250  conn[5] = this->node_id(4);
251  */
252 
253  // VTK's VTK_QUADRATIC_WEDGE first 9 nodes match, then their
254  // middle and top layers of mid-edge nodes are reversed from
255  // LibMesh's.
256  conn.resize(15);
257  for (unsigned i=0; i<9; ++i)
258  conn[i] = this->node_id(i);
259 
260  // top "ring" of mid-edge nodes
261  conn[9] = this->node_id(12);
262  conn[10] = this->node_id(13);
263  conn[11] = this->node_id(14);
264 
265  // middle "ring" of mid-edge nodes
266  conn[12] = this->node_id(9);
267  conn[13] = this->node_id(10);
268  conn[14] = this->node_id(11);
269 
270 
271  return;
272  }
273 
274  default:
275  libmesh_error_msg("Unsupported IO package " << iop);
276  }
277 }
278 
279 
280 
281 
282 unsigned short int Prism15::second_order_adjacent_vertex (const unsigned int n,
283  const unsigned int v) const
284 {
285  libmesh_assert_greater_equal (n, this->n_vertices());
286  libmesh_assert_less (n, this->n_nodes());
287  libmesh_assert_less (v, 2);
288  return _second_order_adjacent_vertices[n-this->n_vertices()][v];
289 }
290 
291 
292 
293 std::pair<unsigned short int, unsigned short int>
294 Prism15::second_order_child_vertex (const unsigned int n) const
295 {
296  libmesh_assert_greater_equal (n, this->n_vertices());
297  libmesh_assert_less (n, this->n_nodes());
298 
299  return std::pair<unsigned short int, unsigned short int>
300  (_second_order_vertex_child_number[n],
301  _second_order_vertex_child_index[n]);
302 }
303 
304 
305 
306 Real Prism15::volume () const
307 {
308  // Make copies of our points. It makes the subsequent calculations a bit
309  // shorter and avoids dereferencing the same pointer multiple times.
310  Point
311  x0 = point(0), x1 = point(1), x2 = point(2), x3 = point(3), x4 = point(4),
312  x5 = point(5), x6 = point(6), x7 = point(7), x8 = point(8), x9 = point(9),
313  x10 = point(10), x11 = point(11), x12 = point(12), x13 = point(13), x14 = point(14);
314 
315  // Terms are copied directly from a Python script.
316  Point dx_dxi[10] =
317  {
318  -x0 - x1 + x10 + 2*x12 - x3 - x4 + 2*x6 - x9,
319  3*x0/2 + x1/2 + 2*x12 - 3*x3/2 - x4/2 - 2*x6,
320  -x0/2 + x1/2 - x10 - x3/2 + x4/2 + x9,
321  2*x0 - 2*x12 + 2*x13 - 2*x14 + 2*x3 - 2*x6 + 2*x7 - 2*x8,
322  -2*x0 - 2*x12 + 2*x13 - 2*x14 + 2*x3 + 2*x6 - 2*x7 + 2*x8,
323  Point(0,0,0),
324  2*x0 + 2*x1 - 4*x12 + 2*x3 + 2*x4 - 4*x6,
325  -2*x0 - 2*x1 - 4*x12 + 2*x3 + 2*x4 + 4*x6,
326  Point(0,0,0),
327  Point(0,0,0)
328  };
329 
330  Point dx_deta[10] =
331  {
332  -x0 + x11 + 2*x14 - x2 - x3 - x5 + 2*x8 - x9,
333  3*x0/2 + 2*x14 + x2/2 - 3*x3/2 - x5/2 - 2*x8,
334  -x0/2 - x11 + x2/2 - x3/2 + x5/2 + x9,
335  2*x0 - 4*x14 + 2*x2 + 2*x3 + 2*x5 - 4*x8,
336  -2*x0 - 4*x14 - 2*x2 + 2*x3 + 2*x5 + 4*x8,
337  Point(0,0,0),
338  2*x0 - 2*x12 + 2*x13 - 2*x14 + 2*x3 - 2*x6 + 2*x7 - 2*x8,
339  -2*x0 - 2*x12 + 2*x13 - 2*x14 + 2*x3 + 2*x6 - 2*x7 + 2*x8,
340  Point(0,0,0),
341  Point(0,0,0)
342  };
343 
344  Point dx_dzeta[10] =
345  {
346  -x0/2 + x3/2,
347  x0 + x3 - 2*x9,
348  Point(0,0,0),
349  3*x0/2 + 2*x14 + x2/2 - 3*x3/2 - x5/2 - 2*x8,
350  -x0 - 2*x11 + x2 - x3 + x5 + 2*x9,
351  -x0 - 2*x14 - x2 + x3 + x5 + 2*x8,
352  3*x0/2 + x1/2 + 2*x12 - 3*x3/2 - x4/2 - 2*x6,
353  -x0 + x1 - 2*x10 - x3 + x4 + 2*x9,
354  -2*x0 - 2*x12 + 2*x13 - 2*x14 + 2*x3 + 2*x6 - 2*x7 + 2*x8,
355  -x0 - x1 - 2*x12 + x3 + x4 + 2*x6
356  };
357 
358  // The quadrature rule for the Prism15 is a tensor product between a
359  // FOURTH-order TRI3 rule (in xi, eta) and a FIFTH-order EDGE2 rule
360  // in zeta.
361 
362  // Number of points in the 2D quadrature rule.
363  const int N2D = 6;
364 
365  // Parameters of the 2D rule
366  static const Real
367  w1 = 1.1169079483900573284750350421656140e-01L,
368  w2 = 5.4975871827660933819163162450105264e-02L,
369  a1 = 4.4594849091596488631832925388305199e-01L,
370  a2 = 9.1576213509770743459571463402201508e-02L;
371 
372  // Points and weights of the 2D rule
373  static const Real w2D[N2D] = {w1, w1, w1, w2, w2, w2};
374 
375  // Quadrature point locations raised to powers. xi[0][2] is
376  // quadrature point 0, squared, xi[1][1] is quadrature point 1 to the
377  // first power, etc. This lets us avoid calling std::pow inside the
378  // loops below.
379  static const Real xi[N2D][3] =
380  {
381  // ^0 ^1 ^2
382  { 1., a1, a1*a1},
383  { 1., 1-2*a1, (1-2*a1)*(1-2*a1)},
384  { 1., a1, a1*a1},
385  { 1., a2, a2*a2},
386  { 1., 1-2*a2, (1-2*a2)*(1-2*a2)},
387  { 1., a2, a2*a2}
388  };
389 
390  static const Real eta[N2D][3] =
391  {
392  // ^0 ^1 ^2
393  { 1., a1, a1*a1},
394  { 1., a1, a1*a1},
395  { 1., 1-2*a1, (1-2*a1)*(1-2*a1)},
396  { 1., a2, a2*a2},
397  { 1., a2, a2*a2},
398  { 1., 1-2*a2, (1-2*a2)*(1-2*a2)}
399  };
400 
401  // Number of points in the 1D quadrature rule.
402  const int N1D = 3;
403 
404  // Points and weights of the 1D quadrature rule.
405  static const Real w1D[N1D] = {5./9, 8./9, 5./9};
406 
407  const Real zeta[N1D][3] =
408  {
409  //^0 ^1 ^2
410  { 1., -std::sqrt(15)/5., 15./25},
411  { 1., 0., 0.},
412  { 1., std::sqrt(15)/5., 15./25}
413  };
414 
415  // The integer exponents for each term.
416  static const int exponents[10][3] =
417  {
418  {0, 0, 0},
419  {0, 0, 1},
420  {0, 0, 2},
421  {0, 1, 0},
422  {0, 1, 1},
423  {0, 2, 0},
424  {1, 0, 0},
425  {1, 0, 1},
426  {1, 1, 0},
427  {2, 0, 0}
428  };
429 
430  Real vol = 0.;
431  for (int i=0; i<N2D; ++i)
432  for (int j=0; j<N1D; ++j)
433  {
434  // Compute dx_dxi, dx_deta, dx_dzeta at the current quadrature point.
435  Point dx_dxi_q, dx_deta_q, dx_dzeta_q;
436  for (int c=0; c<10; ++c)
437  {
438  Real coeff =
439  xi[i][exponents[c][0]]*
440  eta[i][exponents[c][1]]*
441  zeta[j][exponents[c][2]];
442 
443  dx_dxi_q += coeff * dx_dxi[c];
444  dx_deta_q += coeff * dx_deta[c];
445  dx_dzeta_q += coeff * dx_dzeta[c];
446  }
447 
448  // Compute scalar triple product, multiply by weight, and accumulate volume.
449  vol += w2D[i] * w1D[j] * triple_product(dx_dxi_q, dx_deta_q, dx_dzeta_q);
450  }
451 
452  return vol;
453 }
454 
455 
456 
457 #ifdef LIBMESH_ENABLE_AMR
458 
459 const float Prism15::_embedding_matrix[8][15][15] =
460  {
461  // Embedding matrix for child 0
462  {
463  // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
464  { 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, // 0
465  { 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0 }, // 1
466  { 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0 }, // 2
467  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0 }, // 3
468  { -0.25, -0.25, 0, -0.25, -0.25, 0, 0.5, 0, 0, 0.5, 0.5, 0, 0.5, 0, 0 }, // 4
469  { -0.25, 0, -0.25, -0.25, 0, -0.25, 0, 0, 0.5, 0.5, 0, 0.5, 0, 0, 0.5 }, // 5
470  { 0.375, -0.125, 0, 0, 0, 0, 0.75, 0, 0, 0, 0, 0, 0, 0, 0 }, // 6
471  { 0, -0.125, -0.125, 0, 0, 0, 0.5, 0.25, 0.5, 0, 0, 0, 0, 0, 0 }, // 7
472  { 0.375, 0, -0.125, 0, 0, 0, 0, 0, 0.75, 0, 0, 0, 0, 0, 0 }, // 8
473  { 0.375, 0, 0, -0.125, 0, 0, 0, 0, 0, 0.75, 0, 0, 0, 0, 0 }, // 9
474  { -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.75, 0, 0, 0.375, 0.375, 0, 0.25, 0, 0 }, // 10
475  { -0.1875, 0, -0.1875, -0.1875, 0, -0.1875, 0, 0, 0.75, 0.375, 0, 0.375, 0, 0, 0.25 }, // 11
476  { -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.375, 0, 0, 0.75, 0.25, 0, 0.375, 0, 0 }, // 12
477  { -0.25, -0.1875, -0.1875, -0.25, -0.1875, -0.1875, 0.25, 0.125, 0.25, 0.5, 0.25, 0.25, 0.25, 0.125, 0.25 }, // 13
478  { -0.1875, 0, -0.1875, -0.1875, 0, -0.1875, 0, 0, 0.375, 0.75, 0, 0.25, 0, 0, 0.375 } // 14
479  },
480 
481  // Embedding matrix for child 1
482  {
483  // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
484  { 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0 }, // 0
485  { 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, // 1
486  { 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0 }, // 2
487  { -0.25, -0.25, 0, -0.25, -0.25, 0, 0.5, 0, 0, 0.5, 0.5, 0, 0.5, 0, 0 }, // 3
488  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0 }, // 4
489  { 0, -0.25, -0.25, 0, -0.25, -0.25, 0, 0.5, 0, 0, 0.5, 0.5, 0, 0.5, 0 }, // 5
490  { -0.125, 0.375, 0, 0, 0, 0, 0.75, 0, 0, 0, 0, 0, 0, 0, 0 }, // 6
491  { 0, 0.375, -0.125, 0, 0, 0, 0, 0.75, 0, 0, 0, 0, 0, 0, 0 }, // 7
492  { -0.125, 0, -0.125, 0, 0, 0, 0.5, 0.5, 0.25, 0, 0, 0, 0, 0, 0 }, // 8
493  { -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.75, 0, 0, 0.375, 0.375, 0, 0.25, 0, 0 }, // 9
494  { 0, 0.375, 0, 0, -0.125, 0, 0, 0, 0, 0, 0.75, 0, 0, 0, 0 }, // 10
495  { 0, -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.75, 0, 0, 0.375, 0.375, 0, 0.25, 0 }, // 11
496  { -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.375, 0, 0, 0.25, 0.75, 0, 0.375, 0, 0 }, // 12
497  { 0, -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.375, 0, 0, 0.75, 0.25, 0, 0.375, 0 }, // 13
498  { -0.1875, -0.25, -0.1875, -0.1875, -0.25, -0.1875, 0.25, 0.25, 0.125, 0.25, 0.5, 0.25, 0.25, 0.25, 0.125 } // 14
499  },
500 
501  // Embedding matrix for child 2
502  {
503  // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
504  { 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0 }, // 0
505  { 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0 }, // 1
506  { 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, // 2
507  { -0.25, 0, -0.25, -0.25, 0, -0.25, 0, 0, 0.5, 0.5, 0, 0.5, 0, 0, 0.5 }, // 3
508  { 0, -0.25, -0.25, 0, -0.25, -0.25, 0, 0.5, 0, 0, 0.5, 0.5, 0, 0.5, 0 }, // 4
509  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0 }, // 5
510  { -0.125, -0.125, 0, 0, 0, 0, 0.25, 0.5, 0.5, 0, 0, 0, 0, 0, 0 }, // 6
511  { 0, -0.125, 0.375, 0, 0, 0, 0, 0.75, 0, 0, 0, 0, 0, 0, 0 }, // 7
512  { -0.125, 0, 0.375, 0, 0, 0, 0, 0, 0.75, 0, 0, 0, 0, 0, 0 }, // 8
513  { -0.1875, 0, -0.1875, -0.1875, 0, -0.1875, 0, 0, 0.75, 0.375, 0, 0.375, 0, 0, 0.25 }, // 9
514  { 0, -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.75, 0, 0, 0.375, 0.375, 0, 0.25, 0 }, // 10
515  { 0, 0, 0.375, 0, 0, -0.125, 0, 0, 0, 0, 0, 0.75, 0, 0, 0 }, // 11
516  { -0.1875, -0.1875, -0.25, -0.1875, -0.1875, -0.25, 0.125, 0.25, 0.25, 0.25, 0.25, 0.5, 0.125, 0.25, 0.25 }, // 12
517  { 0, -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.375, 0, 0, 0.25, 0.75, 0, 0.375, 0 }, // 13
518  { -0.1875, 0, -0.1875, -0.1875, 0, -0.1875, 0, 0, 0.375, 0.25, 0, 0.75, 0, 0, 0.375 } // 14
519  },
520 
521  // Embedding matrix for child 3
522  {
523  // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
524  { 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0 }, // 0
525  { 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0 }, // 1
526  { 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0 }, // 2
527  { -0.25, -0.25, 0, -0.25, -0.25, 0, 0.5, 0, 0, 0.5, 0.5, 0, 0.5, 0, 0 }, // 3
528  { 0, -0.25, -0.25, 0, -0.25, -0.25, 0, 0.5, 0, 0, 0.5, 0.5, 0, 0.5, 0 }, // 4
529  { -0.25, 0, -0.25, -0.25, 0, -0.25, 0, 0, 0.5, 0.5, 0, 0.5, 0, 0, 0.5 }, // 5
530  { -0.125, 0, -0.125, 0, 0, 0, 0.5, 0.5, 0.25, 0, 0, 0, 0, 0, 0 }, // 6
531  { -0.125, -0.125, 0, 0, 0, 0, 0.25, 0.5, 0.5, 0, 0, 0, 0, 0, 0 }, // 7
532  { 0, -0.125, -0.125, 0, 0, 0, 0.5, 0.25, 0.5, 0, 0, 0, 0, 0, 0 }, // 8
533  { -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.75, 0, 0, 0.375, 0.375, 0, 0.25, 0, 0 }, // 9
534  { 0, -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.75, 0, 0, 0.375, 0.375, 0, 0.25, 0 }, // 10
535  { -0.1875, 0, -0.1875, -0.1875, 0, -0.1875, 0, 0, 0.75, 0.375, 0, 0.375, 0, 0, 0.25 }, // 11
536  { -0.1875, -0.25, -0.1875, -0.1875, -0.25, -0.1875, 0.25, 0.25, 0.125, 0.25, 0.5, 0.25, 0.25, 0.25, 0.125 }, // 12
537  { -0.1875, -0.1875, -0.25, -0.1875, -0.1875, -0.25, 0.125, 0.25, 0.25, 0.25, 0.25, 0.5, 0.125, 0.25, 0.25 }, // 13
538  { -0.25, -0.1875, -0.1875, -0.25, -0.1875, -0.1875, 0.25, 0.125, 0.25, 0.5, 0.25, 0.25, 0.25, 0.125, 0.25 } // 14
539  },
540 
541  // Embedding matrix for child 4
542  {
543  // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
544  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0 }, // 0
545  { -0.25, -0.25, 0, -0.25, -0.25, 0, 0.5, 0, 0, 0.5, 0.5, 0, 0.5, 0, 0 }, // 1
546  { -0.25, 0, -0.25, -0.25, 0, -0.25, 0, 0, 0.5, 0.5, 0, 0.5, 0, 0, 0.5 }, // 2
547  { 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, // 3
548  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0 }, // 4
549  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1 }, // 5
550  { -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.375, 0, 0, 0.75, 0.25, 0, 0.375, 0, 0 }, // 6
551  { -0.25, -0.1875, -0.1875, -0.25, -0.1875, -0.1875, 0.25, 0.125, 0.25, 0.5, 0.25, 0.25, 0.25, 0.125, 0.25 }, // 7
552  { -0.1875, 0, -0.1875, -0.1875, 0, -0.1875, 0, 0, 0.375, 0.75, 0, 0.25, 0, 0, 0.375 }, // 8
553  { -0.125, 0, 0, 0.375, 0, 0, 0, 0, 0, 0.75, 0, 0, 0, 0, 0 }, // 9
554  { -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.25, 0, 0, 0.375, 0.375, 0, 0.75, 0, 0 }, // 10
555  { -0.1875, 0, -0.1875, -0.1875, 0, -0.1875, 0, 0, 0.25, 0.375, 0, 0.375, 0, 0, 0.75 }, // 11
556  { 0, 0, 0, 0.375, -0.125, 0, 0, 0, 0, 0, 0, 0, 0.75, 0, 0 }, // 12
557  { 0, 0, 0, 0, -0.125, -0.125, 0, 0, 0, 0, 0, 0, 0.5, 0.25, 0.5 }, // 13
558  { 0, 0, 0, 0.375, 0, -0.125, 0, 0, 0, 0, 0, 0, 0, 0, 0.75 } // 14
559  },
560 
561  // Embedding matrix for child 5
562  {
563  // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
564  { -0.25, -0.25, 0, -0.25, -0.25, 0, 0.5, 0, 0, 0.5, 0.5, 0, 0.5, 0, 0 }, // 0
565  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0 }, // 1
566  { 0, -0.25, -0.25, 0, -0.25, -0.25, 0, 0.5, 0, 0, 0.5, 0.5, 0, 0.5, 0 }, // 2
567  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0 }, // 3
568  { 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, // 4
569  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0 }, // 5
570  { -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.375, 0, 0, 0.25, 0.75, 0, 0.375, 0, 0 }, // 6
571  { 0, -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.375, 0, 0, 0.75, 0.25, 0, 0.375, 0 }, // 7
572  { -0.1875, -0.25, -0.1875, -0.1875, -0.25, -0.1875, 0.25, 0.25, 0.125, 0.25, 0.5, 0.25, 0.25, 0.25, 0.125 }, // 8
573  { -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.25, 0, 0, 0.375, 0.375, 0, 0.75, 0, 0 }, // 9
574  { 0, -0.125, 0, 0, 0.375, 0, 0, 0, 0, 0, 0.75, 0, 0, 0, 0 }, // 10
575  { 0, -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.25, 0, 0, 0.375, 0.375, 0, 0.75, 0 }, // 11
576  { 0, 0, 0, -0.125, 0.375, 0, 0, 0, 0, 0, 0, 0, 0.75, 0, 0 }, // 12
577  { 0, 0, 0, 0, 0.375, -0.125, 0, 0, 0, 0, 0, 0, 0, 0.75, 0 }, // 13
578  { 0, 0, 0, -0.125, 0, -0.125, 0, 0, 0, 0, 0, 0, 0.5, 0.5, 0.25 } // 14
579  },
580 
581  // Embedding matrix for child 6
582  {
583  // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
584  { -0.25, 0, -0.25, -0.25, 0, -0.25, 0, 0, 0.5, 0.5, 0, 0.5, 0, 0, 0.5 }, // 0
585  { 0, -0.25, -0.25, 0, -0.25, -0.25, 0, 0.5, 0, 0, 0.5, 0.5, 0, 0.5, 0 }, // 1
586  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0 }, // 2
587  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1 }, // 3
588  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0 }, // 4
589  { 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, // 5
590  { -0.1875, -0.1875, -0.25, -0.1875, -0.1875, -0.25, 0.125, 0.25, 0.25, 0.25, 0.25, 0.5, 0.125, 0.25, 0.25 }, // 6
591  { 0, -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.375, 0, 0, 0.25, 0.75, 0, 0.375, 0 }, // 7
592  { -0.1875, 0, -0.1875, -0.1875, 0, -0.1875, 0, 0, 0.375, 0.25, 0, 0.75, 0, 0, 0.375 }, // 8
593  { -0.1875, 0, -0.1875, -0.1875, 0, -0.1875, 0, 0, 0.25, 0.375, 0, 0.375, 0, 0, 0.75 }, // 9
594  { 0, -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.25, 0, 0, 0.375, 0.375, 0, 0.75, 0 }, // 10
595  { 0, 0, -0.125, 0, 0, 0.375, 0, 0, 0, 0, 0, 0.75, 0, 0, 0 }, // 11
596  { 0, 0, 0, -0.125, -0.125, 0, 0, 0, 0, 0, 0, 0, 0.25, 0.5, 0.5 }, // 12
597  { 0, 0, 0, 0, -0.125, 0.375, 0, 0, 0, 0, 0, 0, 0, 0.75, 0 }, // 13
598  { 0, 0, 0, -0.125, 0, 0.375, 0, 0, 0, 0, 0, 0, 0, 0, 0.75 } // 14
599  },
600 
601  // Embedding matrix for child 7
602  {
603  // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
604  { -0.25, -0.25, 0, -0.25, -0.25, 0, 0.5, 0, 0, 0.5, 0.5, 0, 0.5, 0, 0 }, // 0
605  { 0, -0.25, -0.25, 0, -0.25, -0.25, 0, 0.5, 0, 0, 0.5, 0.5, 0, 0.5, 0 }, // 1
606  { -0.25, 0, -0.25, -0.25, 0, -0.25, 0, 0, 0.5, 0.5, 0, 0.5, 0, 0, 0.5 }, // 2
607  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0 }, // 3
608  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0 }, // 4
609  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1 }, // 5
610  { -0.1875, -0.25, -0.1875, -0.1875, -0.25, -0.1875, 0.25, 0.25, 0.125, 0.25, 0.5, 0.25, 0.25, 0.25, 0.125 }, // 6
611  { -0.1875, -0.1875, -0.25, -0.1875, -0.1875, -0.25, 0.125, 0.25, 0.25, 0.25, 0.25, 0.5, 0.125, 0.25, 0.25 }, // 7
612  { -0.25, -0.1875, -0.1875, -0.25, -0.1875, -0.1875, 0.25, 0.125, 0.25, 0.5, 0.25, 0.25, 0.25, 0.125, 0.25 }, // 8
613  { -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.25, 0, 0, 0.375, 0.375, 0, 0.75, 0, 0 }, // 9
614  { 0, -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.25, 0, 0, 0.375, 0.375, 0, 0.75, 0 }, // 10
615  { -0.1875, 0, -0.1875, -0.1875, 0, -0.1875, 0, 0, 0.25, 0.375, 0, 0.375, 0, 0, 0.75 }, // 11
616  { 0, 0, 0, -0.125, 0, -0.125, 0, 0, 0, 0, 0, 0, 0.5, 0.5, 0.25 }, // 12
617  { 0, 0, 0, -0.125, -0.125, 0, 0, 0, 0, 0, 0, 0, 0.25, 0.5, 0.5 }, // 13
618  { 0, 0, 0, 0, -0.125, -0.125, 0, 0, 0, 0, 0, 0, 0.5, 0.25, 0.5 } // 14
619  }
620  };
621 
622 #endif
623 
624 } // namespace libMesh
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
PetscErrorCode Vec Mat libmesh_dbg_var(j)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
IOPackage
libMesh interfaces with several different software packages for the purposes of creating, reading, and writing mesh files.