libMesh
cell_prism15.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2019 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 // Local includes
20 #include "libmesh/side.h"
21 #include "libmesh/cell_prism15.h"
22 #include "libmesh/edge_edge3.h"
23 #include "libmesh/face_quad8.h"
24 #include "libmesh/face_tri6.h"
25 #include "libmesh/enum_io_package.h"
26 #include "libmesh/enum_order.h"
27 
28 namespace libMesh
29 {
30 
31 
32 
33 // ------------------------------------------------------------
34 // Prism15 class static member initializations
35 const int Prism15::num_nodes;
36 const int Prism15::num_sides;
37 const int Prism15::num_edges;
38 const int Prism15::num_children;
39 const int Prism15::nodes_per_side;
40 const int Prism15::nodes_per_edge;
41 
42 const unsigned int Prism15::side_nodes_map[Prism15::num_sides][Prism15::nodes_per_side] =
43  {
44  {0, 2, 1, 8, 7, 6, 99, 99}, // Side 0
45  {0, 1, 4, 3, 6, 10, 12, 9}, // Side 1
46  {1, 2, 5, 4, 7, 11, 13, 10}, // Side 2
47  {2, 0, 3, 5, 8, 9, 14, 11}, // Side 3
48  {3, 4, 5, 12, 13, 14, 99, 99} // Side 4
49  };
50 
51 const unsigned int Prism15::edge_nodes_map[Prism15::num_edges][Prism15::nodes_per_edge] =
52  {
53  {0, 1, 6}, // Edge 0
54  {1, 2, 7}, // Edge 1
55  {0, 2, 8}, // Edge 2
56  {0, 3, 9}, // Edge 3
57  {1, 4, 10}, // Edge 4
58  {2, 5, 11}, // Edge 5
59  {3, 4, 12}, // Edge 6
60  {4, 5, 13}, // Edge 7
61  {3, 5, 14} // Edge 8
62  };
63 
64 
65 // ------------------------------------------------------------
66 // Prism15 class member functions
67 
68 bool Prism15::is_vertex(const unsigned int i) const
69 {
70  if (i < 6)
71  return true;
72  return false;
73 }
74 
75 bool Prism15::is_edge(const unsigned int i) const
76 {
77  if (i < 6)
78  return false;
79  return true;
80 }
81 
82 bool Prism15::is_face(const unsigned int) const
83 {
84  return false;
85 }
86 
87 bool Prism15::is_node_on_side(const unsigned int n,
88  const unsigned int s) const
89 {
90  libmesh_assert_less (s, n_sides());
91  return std::find(std::begin(side_nodes_map[s]),
92  std::end(side_nodes_map[s]),
93  n) != std::end(side_nodes_map[s]);
94 }
95 
96 std::vector<unsigned int>
97 Prism15::nodes_on_side(const unsigned int s) const
98 {
99  libmesh_assert_less(s, n_sides());
100  auto trim = (s > 0 && s < 4) ? 0 : 2;
101  return {std::begin(side_nodes_map[s]), std::end(side_nodes_map[s]) - trim};
102 }
103 
104 bool Prism15::is_node_on_edge(const unsigned int n,
105  const unsigned int e) const
106 {
107  libmesh_assert_less (e, n_edges());
108  return std::find(std::begin(edge_nodes_map[e]),
109  std::end(edge_nodes_map[e]),
110  n) != std::end(edge_nodes_map[e]);
111 }
112 
113 
114 
115 bool Prism15::has_affine_map() const
116 {
117  // Make sure z edges are affine
118  Point v = this->point(3) - this->point(0);
119  if (!v.relative_fuzzy_equals(this->point(4) - this->point(1)) ||
120  !v.relative_fuzzy_equals(this->point(5) - this->point(2)))
121  return false;
122  // Make sure edges are straight
123  v /= 2;
124  if (!v.relative_fuzzy_equals(this->point(9) - this->point(0)) ||
125  !v.relative_fuzzy_equals(this->point(10) - this->point(1)) ||
126  !v.relative_fuzzy_equals(this->point(11) - this->point(2)))
127  return false;
128  v = (this->point(1) - this->point(0))/2;
129  if (!v.relative_fuzzy_equals(this->point(6) - this->point(0)) ||
130  !v.relative_fuzzy_equals(this->point(12) - this->point(3)))
131  return false;
132  v = (this->point(2) - this->point(0))/2;
133  if (!v.relative_fuzzy_equals(this->point(8) - this->point(0)) ||
134  !v.relative_fuzzy_equals(this->point(14) - this->point(3)))
135  return false;
136  v = (this->point(2) - this->point(1))/2;
137  if (!v.relative_fuzzy_equals(this->point(7) - this->point(1)) ||
138  !v.relative_fuzzy_equals(this->point(13) - this->point(4)))
139  return false;
140  return true;
141 }
142 
143 
144 
145 Order Prism15::default_order() const
146 {
147  return SECOND;
148 }
149 
150 
151 
152 unsigned int Prism15::which_node_am_i(unsigned int side,
153  unsigned int side_node) const
154 {
155  libmesh_assert_less (side, this->n_sides());
156 
157  // Never more than 8 nodes per side.
158  libmesh_assert_less(side_node, Prism15::nodes_per_side);
159 
160  // Some sides have 6 nodes.
161  libmesh_assert(!(side==0 || side==4) || side_node < 6);
162 
163  return Prism15::side_nodes_map[side][side_node];
164 }
165 
166 
167 
168 std::unique_ptr<Elem> Prism15::build_side_ptr (const unsigned int i,
169  bool proxy)
170 {
171  libmesh_assert_less (i, this->n_sides());
172 
173  if (proxy)
174  {
175  switch (i)
176  {
177  case 0: // the triangular face at z=-1
178  case 4:
179  return libmesh_make_unique<Side<Tri6,Prism15>>(this,i);
180 
181  case 1:
182  case 2:
183  case 3:
184  return libmesh_make_unique<Side<Quad8,Prism15>>(this,i);
185 
186  default:
187  libmesh_error_msg("Invalid side i = " << i);
188  }
189  }
190 
191  else
192  {
193  // Return value
194  std::unique_ptr<Elem> face;
195 
196  switch (i)
197  {
198  case 0: // the triangular face at z=-1
199  case 4: // the triangular face at z=1
200  {
201  face = libmesh_make_unique<Tri6>();
202  break;
203  }
204  case 1: // the quad face at y=0
205  case 2: // the other quad face
206  case 3: // the quad face at x=0
207  {
208  face = libmesh_make_unique<Quad8>();
209  break;
210  }
211  default:
212  libmesh_error_msg("Invalid side i = " << i);
213  }
214 
215  face->subdomain_id() = this->subdomain_id();
216 
217  // Set the nodes
218  for (auto n : face->node_index_range())
219  face->set_node(n) = this->node_ptr(Prism15::side_nodes_map[i][n]);
220 
221  return face;
222  }
223 }
224 
225 
226 void Prism15::build_side_ptr (std::unique_ptr<Elem> & side,
227  const unsigned int i)
228 {
229  libmesh_assert_less (i, this->n_sides());
230 
231  switch (i)
232  {
233  case 0: // the triangular face at z=-1
234  case 4: // the triangular face at z=1
235  {
236  if (!side.get() || side->type() != TRI6)
237  {
238  side = this->build_side_ptr(i, false);
239  return;
240  }
241  break;
242  }
243 
244  case 1: // the quad face at y=0
245  case 2: // the other quad face
246  case 3: // the quad face at x=0
247  {
248  if (!side.get() || side->type() != QUAD8)
249  {
250  side = this->build_side_ptr(i, false);
251  return;
252  }
253  break;
254  }
255 
256  default:
257  libmesh_error_msg("Invalid side i = " << i);
258  }
259 
260  side->subdomain_id() = this->subdomain_id();
261 
262  // Set the nodes
263  for (auto n : side->node_index_range())
264  side->set_node(n) = this->node_ptr(Prism15::side_nodes_map[i][n]);
265 }
266 
267 
268 
269 std::unique_ptr<Elem> Prism15::build_edge_ptr (const unsigned int i)
270 {
271  libmesh_assert_less (i, this->n_edges());
272 
273  return libmesh_make_unique<SideEdge<Edge3,Prism15>>(this,i);
274 }
275 
276 
277 void Prism15::connectivity(const unsigned int libmesh_dbg_var(sc),
278  const IOPackage iop,
279  std::vector<dof_id_type> & conn) const
280 {
281  libmesh_assert(_nodes);
282  libmesh_assert_less (sc, this->n_sub_elem());
283  libmesh_assert_not_equal_to (iop, INVALID_IO_PACKAGE);
284 
285  switch (iop)
286  {
287  case TECPLOT:
288  {
289  conn.resize(8);
290  conn[0] = this->node_id(0)+1;
291  conn[1] = this->node_id(1)+1;
292  conn[2] = this->node_id(2)+1;
293  conn[3] = this->node_id(2)+1;
294  conn[4] = this->node_id(3)+1;
295  conn[5] = this->node_id(4)+1;
296  conn[6] = this->node_id(5)+1;
297  conn[7] = this->node_id(5)+1;
298  return;
299  }
300 
301  case VTK:
302  {
303  /*
304  conn.resize(6);
305  conn[0] = this->node_id(0);
306  conn[1] = this->node_id(2);
307  conn[2] = this->node_id(1);
308  conn[3] = this->node_id(3);
309  conn[4] = this->node_id(5);
310  conn[5] = this->node_id(4);
311  */
312 
313  // VTK's VTK_QUADRATIC_WEDGE first 9 nodes match, then their
314  // middle and top layers of mid-edge nodes are reversed from
315  // LibMesh's.
316  conn.resize(15);
317  for (unsigned i=0; i<9; ++i)
318  conn[i] = this->node_id(i);
319 
320  // top "ring" of mid-edge nodes
321  conn[9] = this->node_id(12);
322  conn[10] = this->node_id(13);
323  conn[11] = this->node_id(14);
324 
325  // middle "ring" of mid-edge nodes
326  conn[12] = this->node_id(9);
327  conn[13] = this->node_id(10);
328  conn[14] = this->node_id(11);
329 
330 
331  return;
332  }
333 
334  default:
335  libmesh_error_msg("Unsupported IO package " << iop);
336  }
337 }
338 
339 
340 
341 
342 unsigned short int Prism15::second_order_adjacent_vertex (const unsigned int n,
343  const unsigned int v) const
344 {
345  libmesh_assert_greater_equal (n, this->n_vertices());
346  libmesh_assert_less (n, this->n_nodes());
347  libmesh_assert_less (v, 2);
348  return _second_order_adjacent_vertices[n-this->n_vertices()][v];
349 }
350 
351 
352 
353 std::pair<unsigned short int, unsigned short int>
354 Prism15::second_order_child_vertex (const unsigned int n) const
355 {
356  libmesh_assert_greater_equal (n, this->n_vertices());
357  libmesh_assert_less (n, this->n_nodes());
358 
359  return std::pair<unsigned short int, unsigned short int>
360  (_second_order_vertex_child_number[n],
361  _second_order_vertex_child_index[n]);
362 }
363 
364 
365 
366 Real Prism15::volume () const
367 {
368  // This specialization is good for Lagrange mappings only
369  if (this->mapping_type() != LAGRANGE_MAP)
370  return this->Elem::volume();
371 
372  // Make copies of our points. It makes the subsequent calculations a bit
373  // shorter and avoids dereferencing the same pointer multiple times.
374  Point
375  x0 = point(0), x1 = point(1), x2 = point(2), x3 = point(3), x4 = point(4),
376  x5 = point(5), x6 = point(6), x7 = point(7), x8 = point(8), x9 = point(9),
377  x10 = point(10), x11 = point(11), x12 = point(12), x13 = point(13), x14 = point(14);
378 
379  // Terms are copied directly from a Python script.
380  Point dx_dxi[10] =
381  {
382  -x0 - x1 + x10 + 2*x12 - x3 - x4 + 2*x6 - x9,
383  3*x0/2 + x1/2 + 2*x12 - 3*x3/2 - x4/2 - 2*x6,
384  -x0/2 + x1/2 - x10 - x3/2 + x4/2 + x9,
385  2*x0 - 2*x12 + 2*x13 - 2*x14 + 2*x3 - 2*x6 + 2*x7 - 2*x8,
386  -2*x0 - 2*x12 + 2*x13 - 2*x14 + 2*x3 + 2*x6 - 2*x7 + 2*x8,
387  Point(0,0,0),
388  2*x0 + 2*x1 - 4*x12 + 2*x3 + 2*x4 - 4*x6,
389  -2*x0 - 2*x1 - 4*x12 + 2*x3 + 2*x4 + 4*x6,
390  Point(0,0,0),
391  Point(0,0,0)
392  };
393 
394  Point dx_deta[10] =
395  {
396  -x0 + x11 + 2*x14 - x2 - x3 - x5 + 2*x8 - x9,
397  3*x0/2 + 2*x14 + x2/2 - 3*x3/2 - x5/2 - 2*x8,
398  -x0/2 - x11 + x2/2 - x3/2 + x5/2 + x9,
399  2*x0 - 4*x14 + 2*x2 + 2*x3 + 2*x5 - 4*x8,
400  -2*x0 - 4*x14 - 2*x2 + 2*x3 + 2*x5 + 4*x8,
401  Point(0,0,0),
402  2*x0 - 2*x12 + 2*x13 - 2*x14 + 2*x3 - 2*x6 + 2*x7 - 2*x8,
403  -2*x0 - 2*x12 + 2*x13 - 2*x14 + 2*x3 + 2*x6 - 2*x7 + 2*x8,
404  Point(0,0,0),
405  Point(0,0,0)
406  };
407 
408  Point dx_dzeta[10] =
409  {
410  -x0/2 + x3/2,
411  x0 + x3 - 2*x9,
412  Point(0,0,0),
413  3*x0/2 + 2*x14 + x2/2 - 3*x3/2 - x5/2 - 2*x8,
414  -x0 - 2*x11 + x2 - x3 + x5 + 2*x9,
415  -x0 - 2*x14 - x2 + x3 + x5 + 2*x8,
416  3*x0/2 + x1/2 + 2*x12 - 3*x3/2 - x4/2 - 2*x6,
417  -x0 + x1 - 2*x10 - x3 + x4 + 2*x9,
418  -2*x0 - 2*x12 + 2*x13 - 2*x14 + 2*x3 + 2*x6 - 2*x7 + 2*x8,
419  -x0 - x1 - 2*x12 + x3 + x4 + 2*x6
420  };
421 
422  // The quadrature rule for the Prism15 is a tensor product between a
423  // FOURTH-order TRI3 rule (in xi, eta) and a FIFTH-order EDGE2 rule
424  // in zeta.
425 
426  // Number of points in the 2D quadrature rule.
427  const int N2D = 6;
428 
429  // Parameters of the 2D rule
430  static const Real
431  w1 = Real(1.1169079483900573284750350421656140e-01L),
432  w2 = Real(5.4975871827660933819163162450105264e-02L),
433  a1 = Real(4.4594849091596488631832925388305199e-01L),
434  a2 = Real(9.1576213509770743459571463402201508e-02L);
435 
436  // Points and weights of the 2D rule
437  static const Real w2D[N2D] = {w1, w1, w1, w2, w2, w2};
438 
439  // Quadrature point locations raised to powers. xi[0][2] is
440  // quadrature point 0, squared, xi[1][1] is quadrature point 1 to the
441  // first power, etc. This lets us avoid calling std::pow inside the
442  // loops below.
443  static const Real xi[N2D][3] =
444  {
445  // ^0 ^1 ^2
446  { 1., a1, a1*a1},
447  { 1., 1-2*a1, (1-2*a1)*(1-2*a1)},
448  { 1., a1, a1*a1},
449  { 1., a2, a2*a2},
450  { 1., 1-2*a2, (1-2*a2)*(1-2*a2)},
451  { 1., a2, a2*a2}
452  };
453 
454  static const Real eta[N2D][3] =
455  {
456  // ^0 ^1 ^2
457  { 1., a1, a1*a1},
458  { 1., a1, a1*a1},
459  { 1., 1-2*a1, (1-2*a1)*(1-2*a1)},
460  { 1., a2, a2*a2},
461  { 1., a2, a2*a2},
462  { 1., 1-2*a2, (1-2*a2)*(1-2*a2)}
463  };
464 
465  // Number of points in the 1D quadrature rule.
466  const int N1D = 3;
467 
468  // Points and weights of the 1D quadrature rule.
469  static const Real w1D[N1D] = {5./9, 8./9, 5./9};
470 
471  const Real zeta[N1D][3] =
472  {
473  //^0 ^1 ^2
474  { 1., -std::sqrt(15)/5., 15./25},
475  { 1., 0., 0.},
476  { 1., std::sqrt(15)/5., 15./25}
477  };
478 
479  // The integer exponents for each term.
480  static const int exponents[10][3] =
481  {
482  {0, 0, 0},
483  {0, 0, 1},
484  {0, 0, 2},
485  {0, 1, 0},
486  {0, 1, 1},
487  {0, 2, 0},
488  {1, 0, 0},
489  {1, 0, 1},
490  {1, 1, 0},
491  {2, 0, 0}
492  };
493 
494  Real vol = 0.;
495  for (int i=0; i<N2D; ++i)
496  for (int j=0; j<N1D; ++j)
497  {
498  // Compute dx_dxi, dx_deta, dx_dzeta at the current quadrature point.
499  Point dx_dxi_q, dx_deta_q, dx_dzeta_q;
500  for (int c=0; c<10; ++c)
501  {
502  Real coeff =
503  xi[i][exponents[c][0]]*
504  eta[i][exponents[c][1]]*
505  zeta[j][exponents[c][2]];
506 
507  dx_dxi_q += coeff * dx_dxi[c];
508  dx_deta_q += coeff * dx_deta[c];
509  dx_dzeta_q += coeff * dx_dzeta[c];
510  }
511 
512  // Compute scalar triple product, multiply by weight, and accumulate volume.
513  vol += w2D[i] * w1D[j] * triple_product(dx_dxi_q, dx_deta_q, dx_dzeta_q);
514  }
515 
516  return vol;
517 }
518 
519 
520 
521 #ifdef LIBMESH_ENABLE_AMR
522 
523 const float Prism15::_embedding_matrix[Prism15::num_children][Prism15::num_nodes][Prism15::num_nodes] =
524  {
525  // Embedding matrix for child 0
526  {
527  // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
528  { 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, // 0
529  { 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0 }, // 1
530  { 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0 }, // 2
531  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0 }, // 3
532  { -0.25, -0.25, 0, -0.25, -0.25, 0, 0.5, 0, 0, 0.5, 0.5, 0, 0.5, 0, 0 }, // 4
533  { -0.25, 0, -0.25, -0.25, 0, -0.25, 0, 0, 0.5, 0.5, 0, 0.5, 0, 0, 0.5 }, // 5
534  { 0.375, -0.125, 0, 0, 0, 0, 0.75, 0, 0, 0, 0, 0, 0, 0, 0 }, // 6
535  { 0, -0.125, -0.125, 0, 0, 0, 0.5, 0.25, 0.5, 0, 0, 0, 0, 0, 0 }, // 7
536  { 0.375, 0, -0.125, 0, 0, 0, 0, 0, 0.75, 0, 0, 0, 0, 0, 0 }, // 8
537  { 0.375, 0, 0, -0.125, 0, 0, 0, 0, 0, 0.75, 0, 0, 0, 0, 0 }, // 9
538  { -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.75, 0, 0, 0.375, 0.375, 0, 0.25, 0, 0 }, // 10
539  { -0.1875, 0, -0.1875, -0.1875, 0, -0.1875, 0, 0, 0.75, 0.375, 0, 0.375, 0, 0, 0.25 }, // 11
540  { -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.375, 0, 0, 0.75, 0.25, 0, 0.375, 0, 0 }, // 12
541  { -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
542  { -0.1875, 0, -0.1875, -0.1875, 0, -0.1875, 0, 0, 0.375, 0.75, 0, 0.25, 0, 0, 0.375 } // 14
543  },
544 
545  // Embedding matrix for child 1
546  {
547  // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
548  { 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0 }, // 0
549  { 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, // 1
550  { 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0 }, // 2
551  { -0.25, -0.25, 0, -0.25, -0.25, 0, 0.5, 0, 0, 0.5, 0.5, 0, 0.5, 0, 0 }, // 3
552  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0 }, // 4
553  { 0, -0.25, -0.25, 0, -0.25, -0.25, 0, 0.5, 0, 0, 0.5, 0.5, 0, 0.5, 0 }, // 5
554  { -0.125, 0.375, 0, 0, 0, 0, 0.75, 0, 0, 0, 0, 0, 0, 0, 0 }, // 6
555  { 0, 0.375, -0.125, 0, 0, 0, 0, 0.75, 0, 0, 0, 0, 0, 0, 0 }, // 7
556  { -0.125, 0, -0.125, 0, 0, 0, 0.5, 0.5, 0.25, 0, 0, 0, 0, 0, 0 }, // 8
557  { -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.75, 0, 0, 0.375, 0.375, 0, 0.25, 0, 0 }, // 9
558  { 0, 0.375, 0, 0, -0.125, 0, 0, 0, 0, 0, 0.75, 0, 0, 0, 0 }, // 10
559  { 0, -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.75, 0, 0, 0.375, 0.375, 0, 0.25, 0 }, // 11
560  { -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.375, 0, 0, 0.25, 0.75, 0, 0.375, 0, 0 }, // 12
561  { 0, -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.375, 0, 0, 0.75, 0.25, 0, 0.375, 0 }, // 13
562  { -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
563  },
564 
565  // Embedding matrix for child 2
566  {
567  // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
568  { 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0 }, // 0
569  { 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0 }, // 1
570  { 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, // 2
571  { -0.25, 0, -0.25, -0.25, 0, -0.25, 0, 0, 0.5, 0.5, 0, 0.5, 0, 0, 0.5 }, // 3
572  { 0, -0.25, -0.25, 0, -0.25, -0.25, 0, 0.5, 0, 0, 0.5, 0.5, 0, 0.5, 0 }, // 4
573  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0 }, // 5
574  { -0.125, -0.125, 0, 0, 0, 0, 0.25, 0.5, 0.5, 0, 0, 0, 0, 0, 0 }, // 6
575  { 0, -0.125, 0.375, 0, 0, 0, 0, 0.75, 0, 0, 0, 0, 0, 0, 0 }, // 7
576  { -0.125, 0, 0.375, 0, 0, 0, 0, 0, 0.75, 0, 0, 0, 0, 0, 0 }, // 8
577  { -0.1875, 0, -0.1875, -0.1875, 0, -0.1875, 0, 0, 0.75, 0.375, 0, 0.375, 0, 0, 0.25 }, // 9
578  { 0, -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.75, 0, 0, 0.375, 0.375, 0, 0.25, 0 }, // 10
579  { 0, 0, 0.375, 0, 0, -0.125, 0, 0, 0, 0, 0, 0.75, 0, 0, 0 }, // 11
580  { -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
581  { 0, -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.375, 0, 0, 0.25, 0.75, 0, 0.375, 0 }, // 13
582  { -0.1875, 0, -0.1875, -0.1875, 0, -0.1875, 0, 0, 0.375, 0.25, 0, 0.75, 0, 0, 0.375 } // 14
583  },
584 
585  // Embedding matrix for child 3
586  {
587  // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
588  { 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0 }, // 0
589  { 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0 }, // 1
590  { 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0 }, // 2
591  { -0.25, -0.25, 0, -0.25, -0.25, 0, 0.5, 0, 0, 0.5, 0.5, 0, 0.5, 0, 0 }, // 3
592  { 0, -0.25, -0.25, 0, -0.25, -0.25, 0, 0.5, 0, 0, 0.5, 0.5, 0, 0.5, 0 }, // 4
593  { -0.25, 0, -0.25, -0.25, 0, -0.25, 0, 0, 0.5, 0.5, 0, 0.5, 0, 0, 0.5 }, // 5
594  { -0.125, 0, -0.125, 0, 0, 0, 0.5, 0.5, 0.25, 0, 0, 0, 0, 0, 0 }, // 6
595  { -0.125, -0.125, 0, 0, 0, 0, 0.25, 0.5, 0.5, 0, 0, 0, 0, 0, 0 }, // 7
596  { 0, -0.125, -0.125, 0, 0, 0, 0.5, 0.25, 0.5, 0, 0, 0, 0, 0, 0 }, // 8
597  { -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.75, 0, 0, 0.375, 0.375, 0, 0.25, 0, 0 }, // 9
598  { 0, -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.75, 0, 0, 0.375, 0.375, 0, 0.25, 0 }, // 10
599  { -0.1875, 0, -0.1875, -0.1875, 0, -0.1875, 0, 0, 0.75, 0.375, 0, 0.375, 0, 0, 0.25 }, // 11
600  { -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
601  { -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
602  { -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
603  },
604 
605  // Embedding matrix for child 4
606  {
607  // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
608  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0 }, // 0
609  { -0.25, -0.25, 0, -0.25, -0.25, 0, 0.5, 0, 0, 0.5, 0.5, 0, 0.5, 0, 0 }, // 1
610  { -0.25, 0, -0.25, -0.25, 0, -0.25, 0, 0, 0.5, 0.5, 0, 0.5, 0, 0, 0.5 }, // 2
611  { 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, // 3
612  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0 }, // 4
613  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1 }, // 5
614  { -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.375, 0, 0, 0.75, 0.25, 0, 0.375, 0, 0 }, // 6
615  { -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
616  { -0.1875, 0, -0.1875, -0.1875, 0, -0.1875, 0, 0, 0.375, 0.75, 0, 0.25, 0, 0, 0.375 }, // 8
617  { -0.125, 0, 0, 0.375, 0, 0, 0, 0, 0, 0.75, 0, 0, 0, 0, 0 }, // 9
618  { -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.25, 0, 0, 0.375, 0.375, 0, 0.75, 0, 0 }, // 10
619  { -0.1875, 0, -0.1875, -0.1875, 0, -0.1875, 0, 0, 0.25, 0.375, 0, 0.375, 0, 0, 0.75 }, // 11
620  { 0, 0, 0, 0.375, -0.125, 0, 0, 0, 0, 0, 0, 0, 0.75, 0, 0 }, // 12
621  { 0, 0, 0, 0, -0.125, -0.125, 0, 0, 0, 0, 0, 0, 0.5, 0.25, 0.5 }, // 13
622  { 0, 0, 0, 0.375, 0, -0.125, 0, 0, 0, 0, 0, 0, 0, 0, 0.75 } // 14
623  },
624 
625  // Embedding matrix for child 5
626  {
627  // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
628  { -0.25, -0.25, 0, -0.25, -0.25, 0, 0.5, 0, 0, 0.5, 0.5, 0, 0.5, 0, 0 }, // 0
629  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0 }, // 1
630  { 0, -0.25, -0.25, 0, -0.25, -0.25, 0, 0.5, 0, 0, 0.5, 0.5, 0, 0.5, 0 }, // 2
631  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0 }, // 3
632  { 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, // 4
633  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0 }, // 5
634  { -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.375, 0, 0, 0.25, 0.75, 0, 0.375, 0, 0 }, // 6
635  { 0, -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.375, 0, 0, 0.75, 0.25, 0, 0.375, 0 }, // 7
636  { -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
637  { -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.25, 0, 0, 0.375, 0.375, 0, 0.75, 0, 0 }, // 9
638  { 0, -0.125, 0, 0, 0.375, 0, 0, 0, 0, 0, 0.75, 0, 0, 0, 0 }, // 10
639  { 0, -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.25, 0, 0, 0.375, 0.375, 0, 0.75, 0 }, // 11
640  { 0, 0, 0, -0.125, 0.375, 0, 0, 0, 0, 0, 0, 0, 0.75, 0, 0 }, // 12
641  { 0, 0, 0, 0, 0.375, -0.125, 0, 0, 0, 0, 0, 0, 0, 0.75, 0 }, // 13
642  { 0, 0, 0, -0.125, 0, -0.125, 0, 0, 0, 0, 0, 0, 0.5, 0.5, 0.25 } // 14
643  },
644 
645  // Embedding matrix for child 6
646  {
647  // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
648  { -0.25, 0, -0.25, -0.25, 0, -0.25, 0, 0, 0.5, 0.5, 0, 0.5, 0, 0, 0.5 }, // 0
649  { 0, -0.25, -0.25, 0, -0.25, -0.25, 0, 0.5, 0, 0, 0.5, 0.5, 0, 0.5, 0 }, // 1
650  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0 }, // 2
651  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1 }, // 3
652  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0 }, // 4
653  { 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, // 5
654  { -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
655  { 0, -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.375, 0, 0, 0.25, 0.75, 0, 0.375, 0 }, // 7
656  { -0.1875, 0, -0.1875, -0.1875, 0, -0.1875, 0, 0, 0.375, 0.25, 0, 0.75, 0, 0, 0.375 }, // 8
657  { -0.1875, 0, -0.1875, -0.1875, 0, -0.1875, 0, 0, 0.25, 0.375, 0, 0.375, 0, 0, 0.75 }, // 9
658  { 0, -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.25, 0, 0, 0.375, 0.375, 0, 0.75, 0 }, // 10
659  { 0, 0, -0.125, 0, 0, 0.375, 0, 0, 0, 0, 0, 0.75, 0, 0, 0 }, // 11
660  { 0, 0, 0, -0.125, -0.125, 0, 0, 0, 0, 0, 0, 0, 0.25, 0.5, 0.5 }, // 12
661  { 0, 0, 0, 0, -0.125, 0.375, 0, 0, 0, 0, 0, 0, 0, 0.75, 0 }, // 13
662  { 0, 0, 0, -0.125, 0, 0.375, 0, 0, 0, 0, 0, 0, 0, 0, 0.75 } // 14
663  },
664 
665  // Embedding matrix for child 7
666  {
667  // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
668  { -0.25, -0.25, 0, -0.25, -0.25, 0, 0.5, 0, 0, 0.5, 0.5, 0, 0.5, 0, 0 }, // 0
669  { 0, -0.25, -0.25, 0, -0.25, -0.25, 0, 0.5, 0, 0, 0.5, 0.5, 0, 0.5, 0 }, // 1
670  { -0.25, 0, -0.25, -0.25, 0, -0.25, 0, 0, 0.5, 0.5, 0, 0.5, 0, 0, 0.5 }, // 2
671  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0 }, // 3
672  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0 }, // 4
673  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1 }, // 5
674  { -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
675  { -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
676  { -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
677  { -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.25, 0, 0, 0.375, 0.375, 0, 0.75, 0, 0 }, // 9
678  { 0, -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.25, 0, 0, 0.375, 0.375, 0, 0.75, 0 }, // 10
679  { -0.1875, 0, -0.1875, -0.1875, 0, -0.1875, 0, 0, 0.25, 0.375, 0, 0.375, 0, 0, 0.75 }, // 11
680  { 0, 0, 0, -0.125, 0, -0.125, 0, 0, 0, 0, 0, 0, 0.5, 0.5, 0.25 }, // 12
681  { 0, 0, 0, -0.125, -0.125, 0, 0, 0, 0, 0, 0, 0, 0.25, 0.5, 0.5 }, // 13
682  { 0, 0, 0, 0, -0.125, -0.125, 0, 0, 0, 0, 0, 0, 0.5, 0.25, 0.5 } // 14
683  }
684  };
685 
686 #endif
687 
688 } // namespace libMesh
libMesh::Prism::n_vertices
virtual unsigned int n_vertices() const override final
Definition: cell_prism.h:84
libMesh::Prism15::second_order_adjacent_vertex
virtual unsigned short int second_order_adjacent_vertex(const unsigned int n, const unsigned int v) const override
Definition: cell_prism15.C:342
libMesh::IOPackage
IOPackage
libMesh interfaces with several different software packages for the purposes of creating,...
Definition: enum_io_package.h:37
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:164
libMesh::Prism15::which_node_am_i
virtual unsigned int which_node_am_i(unsigned int side, unsigned int side_node) const override
Definition: cell_prism15.C:152
libMesh::Prism15::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_prism15.h:215
libMesh
The libMesh namespace provides an interface to certain functionality in the library.
Definition: factoryfunction.C:55
libMesh::Order
Order
Definition: enum_order.h:40
libMesh::INVALID_IO_PACKAGE
Definition: enum_io_package.h:48
end
IterBase * end
Also have a polymorphic pointer to the end object, this prevents iterating past the end.
Definition: variant_filter_iterator.h:343
std::sqrt
MetaPhysicL::DualNumber< T, D > sqrt(const MetaPhysicL::DualNumber< T, D > &in)
libMesh::Prism::n_edges
virtual unsigned int n_edges() const override final
Definition: cell_prism.h:89
libMesh::Prism15::num_nodes
static const int num_nodes
Geometric constants for Prism15.
Definition: cell_prism15.h:198
libMesh::triple_product
T triple_product(const TypeVector< T > &a, const TypeVector< T > &b, const TypeVector< T > &c)
Definition: type_vector.h:1106
libMesh::SECOND
Definition: enum_order.h:43
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:159
libMesh::Prism15::is_edge
virtual bool is_edge(const unsigned int i) const override
Definition: cell_prism15.C:75
libMesh::Elem::_nodes
Node ** _nodes
Pointers to the nodes we are connected to.
Definition: elem.h:1743
libMesh::Elem::point
const Point & point(const unsigned int i) const
Definition: elem.h:1955
libMesh::Prism15::volume
virtual Real volume() const override
A specialization for computing the volume of a Prism15.
Definition: cell_prism15.C:366
libMesh::Prism15::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_prism15.h:209
libMesh::libmesh_assert
libmesh_assert(ctx)
libMesh::Prism15::n_sub_elem
virtual unsigned int n_sub_elem() const override
Definition: cell_prism15.h:100
libMesh::Elem::volume
virtual Real volume() const
Definition: elem.C:2617
libMesh::Prism15::is_node_on_side
virtual bool is_node_on_side(const unsigned int n, const unsigned int s) const override
Definition: cell_prism15.C:87
libMesh::Prism15::default_order
virtual Order default_order() const override
Definition: cell_prism15.C:145
libMesh::Point
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
libMesh::Prism15::_embedding_matrix
static const float _embedding_matrix[num_children][num_nodes][num_nodes]
Matrix that computes new nodal locations/solution values from current nodes/solution.
Definition: cell_prism15.h:245
libMesh::Elem::mapping_type
ElemMappingType mapping_type() const
Definition: elem.h:2524
libMesh::Prism15::nodes_per_edge
static const int nodes_per_edge
Definition: cell_prism15.h:203
libMesh::Prism15::n_nodes
virtual unsigned int n_nodes() const override
Definition: cell_prism15.h:95
libMesh::Prism15::num_edges
static const int num_edges
Definition: cell_prism15.h:200
libMesh::Prism15::is_vertex
virtual bool is_vertex(const unsigned int i) const override
Definition: cell_prism15.C:68
libMesh::TRI6
Definition: enum_elem_type.h:40
libMesh::LAGRANGE_MAP
Definition: enum_elem_type.h:83
libMesh::Prism::n_sides
virtual unsigned int n_sides() const override final
Definition: cell_prism.h:79
libMesh::Prism15::nodes_on_side
virtual std::vector< unsigned int > nodes_on_side(const unsigned int s) const override
Definition: cell_prism15.C:97
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:154
libMesh::Elem::subdomain_id
subdomain_id_type subdomain_id() const
Definition: elem.h:2069
libMesh::Prism15::connectivity
virtual void connectivity(const unsigned int sc, const IOPackage iop, std::vector< dof_id_type > &conn) const override
Definition: cell_prism15.C:277
libMesh::Prism15::is_node_on_edge
virtual bool is_node_on_edge(const unsigned int n, const unsigned int e) const override
Definition: cell_prism15.C:104
libMesh::Prism15::num_sides
static const int num_sides
Definition: cell_prism15.h:199
libMesh::VTK
Definition: enum_io_package.h:42
libMesh::Prism15::num_children
static const int num_children
Definition: cell_prism15.h:201
libMesh::Prism15::has_affine_map
virtual bool has_affine_map() const override
Definition: cell_prism15.C:115
libMesh::Prism15::build_side_ptr
virtual std::unique_ptr< Elem > build_side_ptr(const unsigned int i, bool proxy=true) override
Builds a QUAD8 or TRI6 built coincident with face i.
Definition: cell_prism15.C:168
libMesh::Elem::node_id
dof_id_type node_id(const unsigned int i) const
Definition: elem.h:1977
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:121
libMesh::Prism15::is_face
virtual bool is_face(const unsigned int i) const override
Definition: cell_prism15.C:82
libMesh::Prism15::nodes_per_side
static const int nodes_per_side
Definition: cell_prism15.h:202
libMesh::Prism15::build_edge_ptr
virtual std::unique_ptr< Elem > build_edge_ptr(const unsigned int i) override
Builds a EDGE3 or INFEDGE2 coincident with edge i.
Definition: cell_prism15.C:269
libMesh::Elem::node_ptr
const Node * node_ptr(const unsigned int i) const
Definition: elem.h:2009
libMesh::Prism15::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_prism15.C:354
libMesh::TypeVector::relative_fuzzy_equals
bool relative_fuzzy_equals(const TypeVector< T > &rhs, Real tol=TOLERANCE) const
Definition: type_vector.h:1042
libMesh::TECPLOT
Definition: enum_io_package.h:39
libMesh::QUAD8
Definition: enum_elem_type.h:42
Generated on Sat Jan 25 2020 12:06:59 for libMesh by   doxygen 1.8.16