libMesh
cell_prism21.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 // Local includes
20 #include "libmesh/cell_prism21.h"
21 #include "libmesh/edge_edge3.h"
22 #include "libmesh/face_quad9.h"
23 #include "libmesh/face_tri7.h"
24 #include "libmesh/enum_io_package.h"
25 #include "libmesh/enum_order.h"
26 #include "libmesh/int_range.h"
27 
28 #ifdef LIBMESH_ENABLE_AMR
29 namespace {
30  // Avoid downcasting when Real > double
31  constexpr libMesh::Real r6 = 6;
32  constexpr libMesh::Real r9 = 9;
33  constexpr libMesh::Real r12 = 12;
34  constexpr libMesh::Real r18 = 18;
35  constexpr libMesh::Real r24 = 24;
36  constexpr libMesh::Real r36 = 36;
37  constexpr libMesh::Real r48 = 48;
38  constexpr libMesh::Real r72 = 72;
39  constexpr libMesh::Real r144 = 144;
40 }
41 #endif
42 
43 namespace libMesh
44 {
45 
46 
47 
48 // ------------------------------------------------------------
49 // Prism21 class static member initializations
50 const int Prism21::num_nodes;
51 const int Prism21::nodes_per_side;
52 const int Prism21::nodes_per_edge;
53 
54 const unsigned int Prism21::side_nodes_map[Prism21::num_sides][Prism21::nodes_per_side] =
55  {
56  {0, 2, 1, 8, 7, 6, 18, 99, 99}, // Side 0
57  {0, 1, 4, 3, 6, 10, 12, 9, 15}, // Side 1
58  {1, 2, 5, 4, 7, 11, 13, 10, 16}, // Side 2
59  {2, 0, 3, 5, 8, 9, 14, 11, 17}, // Side 3
60  {3, 4, 5, 12, 13, 14, 19, 99, 99} // Side 4
61  };
62 
63 const unsigned int Prism21::edge_nodes_map[Prism21::num_edges][Prism21::nodes_per_edge] =
64  {
65  {0, 1, 6}, // Edge 0
66  {1, 2, 7}, // Edge 1
67  {0, 2, 8}, // Edge 2
68  {0, 3, 9}, // Edge 3
69  {1, 4, 10}, // Edge 4
70  {2, 5, 11}, // Edge 5
71  {3, 4, 12}, // Edge 6
72  {4, 5, 13}, // Edge 7
73  {3, 5, 14} // Edge 8
74  };
75 
76 // ------------------------------------------------------------
77 // Prism21 class member functions
78 
79 bool Prism21::is_vertex(const unsigned int i) const
80 {
81  if (i < 6)
82  return true;
83  return false;
84 }
85 
86 bool Prism21::is_edge(const unsigned int i) const
87 {
88  if (i < 6)
89  return false;
90  if (i > 14)
91  return false;
92  return true;
93 }
94 
95 bool Prism21::is_face(const unsigned int i) const
96 {
97  if (i > 19)
98  return false;
99  if (i > 14)
100  return true;
101  return false;
102 }
103 
104 bool Prism21::is_node_on_side(const unsigned int n,
105  const unsigned int s) const
106 {
107  libmesh_assert_less (s, n_sides());
108  return std::find(std::begin(side_nodes_map[s]),
109  std::end(side_nodes_map[s]),
110  n) != std::end(side_nodes_map[s]);
111 }
112 
113 std::vector<unsigned>
114 Prism21::nodes_on_side(const unsigned int s) const
115 {
116  libmesh_assert_less(s, n_sides());
117  auto trim = (s > 0 && s < 4) ? 0 : 2;
118  return {std::begin(side_nodes_map[s]), std::end(side_nodes_map[s]) - trim};
119 }
120 
121 std::vector<unsigned>
122 Prism21::nodes_on_edge(const unsigned int e) const
123 {
124  libmesh_assert_less(e, n_edges());
125  return {std::begin(edge_nodes_map[e]), std::end(edge_nodes_map[e])};
126 }
127 
128 bool Prism21::is_node_on_edge(const unsigned int n,
129  const unsigned int e) const
130 {
131  libmesh_assert_less (e, n_edges());
132  return std::find(std::begin(edge_nodes_map[e]),
133  std::end(edge_nodes_map[e]),
134  n) != std::end(edge_nodes_map[e]);
135 }
136 
137 
138 
139 bool Prism21::has_affine_map() const
140 {
141  // Make sure z edges are affine
142  Point v = this->point(3) - this->point(0);
143  if (!v.relative_fuzzy_equals(this->point(4) - this->point(1), affine_tol) ||
144  !v.relative_fuzzy_equals(this->point(5) - this->point(2), affine_tol))
145  return false;
146 
147  // Make sure edges are straight
148  v /= 2;
149  if (!v.relative_fuzzy_equals(this->point(9) - this->point(0), affine_tol) ||
150  !v.relative_fuzzy_equals(this->point(10) - this->point(1), affine_tol) ||
151  !v.relative_fuzzy_equals(this->point(11) - this->point(2), affine_tol) ||
152  !v.relative_fuzzy_equals(this->point(15) - this->point(6), affine_tol) ||
153  !v.relative_fuzzy_equals(this->point(16) - this->point(7), affine_tol) ||
154  !v.relative_fuzzy_equals(this->point(17) - this->point(8), affine_tol))
155  return false;
156  v = (this->point(1) - this->point(0))/2;
157  if (!v.relative_fuzzy_equals(this->point(6) - this->point(0), affine_tol) ||
158  !v.relative_fuzzy_equals(this->point(12) - this->point(3), affine_tol))
159  return false;
160  v = (this->point(2) - this->point(0))/2;
161  if (!v.relative_fuzzy_equals(this->point(8) - this->point(0), affine_tol) ||
162  !v.relative_fuzzy_equals(this->point(14) - this->point(3), affine_tol))
163  return false;
164  v = (this->point(2) - this->point(1))/2;
165  if (!v.relative_fuzzy_equals(this->point(7) - this->point(1), affine_tol) ||
166  !v.relative_fuzzy_equals(this->point(13) - this->point(4), affine_tol))
167  return false;
168 
169  // Make sure triangle face midpoints are centered
170  v = this->point(2) + this->point(1) - 2*this->point(0);
171  if (!v.relative_fuzzy_equals((this->point(18)-this->point(0))*3))
172  return false;
173  v = this->point(5) + this->point(4) - 2*this->point(3);
174  if (!v.relative_fuzzy_equals((this->point(19)-this->point(3))*3))
175  return false;
176  v = this->point(11) + this->point(10) - 2*this->point(9);
177  if (!v.relative_fuzzy_equals((this->point(20)-this->point(9))*3))
178  return false;
179 
180  return true;
181 }
182 
183 
184 
185 Order Prism21::default_order() const
186 {
187  return THIRD;
188 }
189 
190 dof_id_type Prism21::key (const unsigned int s) const
191 {
192  libmesh_assert_less (s, this->n_sides());
193 
194  switch (s)
195  {
196  case 0: // the triangular face at z=0
197  {
198  return Elem::compute_key (this->node_id(18));
199  }
200  case 1: // the quad face at y=0
201  {
202  return Elem::compute_key (this->node_id(15));
203  }
204  case 2: // the other quad face
205  {
206  return Elem::compute_key (this->node_id(16));
207  }
208  case 3: // the quad face at x=0
209  {
210  return Elem::compute_key (this->node_id(17));
211  }
212  case 4: // the triangular face at z=1
213  {
214  return Elem::compute_key (this->node_id(19));
215  }
216  default:
217  libmesh_error_msg("Invalid side " << s);
218  }
219 }
220 
221 
222 
223 unsigned int Prism21::local_side_node(unsigned int side,
224  unsigned int side_node) const
225 {
226  libmesh_assert_less (side, this->n_sides());
227 
228  // Never more than 9 nodes per side.
229  libmesh_assert_less(side_node, Prism21::nodes_per_side);
230 
231  // Some sides have 7 nodes.
232  libmesh_assert(!(side==0 || side==4) || side_node < 7);
233 
234  return Prism21::side_nodes_map[side][side_node];
235 }
236 
237 
238 
239 unsigned int Prism21::local_edge_node(unsigned int edge,
240  unsigned int edge_node) const
241 {
242  libmesh_assert_less(edge, this->n_edges());
243  libmesh_assert_less(edge_node, Prism21::nodes_per_edge);
244 
245  return Prism21::edge_nodes_map[edge][edge_node];
246 }
247 
248 
249 
250 std::unique_ptr<Elem> Prism21::build_side_ptr (const unsigned int i)
251 {
252  libmesh_assert_less (i, this->n_sides());
253 
254  std::unique_ptr<Elem> face;
255 
256  switch (i)
257  {
258  case 0: // the triangular face at z=-1
259  case 4: // the triangular face at z=1
260  {
261  face = std::make_unique<Tri7>();
262  break;
263  }
264  case 1: // the quad face at y=0
265  case 2: // the other quad face
266  case 3: // the quad face at x=0
267  {
268  face = std::make_unique<Quad9>();
269  break;
270  }
271  default:
272  libmesh_error_msg("Invalid side i = " << i);
273  }
274 
275  // Set the nodes
276  for (auto n : face->node_index_range())
277  face->set_node(n, this->node_ptr(Prism21::side_nodes_map[i][n]));
278 
279  face->set_interior_parent(this);
280  face->inherit_data_from(*this);
281 
282  return face;
283 }
284 
285 
286 
287 void Prism21::build_side_ptr (std::unique_ptr<Elem> & side,
288  const unsigned int i)
289 {
290  libmesh_assert_less (i, this->n_sides());
291 
292  switch (i)
293  {
294  case 0: // the triangular face at z=-1
295  case 4: // the triangular face at z=1
296  {
297  if (!side.get() || side->type() != TRI7)
298  {
299  side = this->build_side_ptr(i);
300  return;
301  }
302  break;
303  }
304 
305  case 1: // the quad face at y=0
306  case 2: // the other quad face
307  case 3: // the quad face at x=0
308  {
309  if (!side.get() || side->type() != QUAD9)
310  {
311  side = this->build_side_ptr(i);
312  return;
313  }
314  break;
315  }
316 
317  default:
318  libmesh_error_msg("Invalid side i = " << i);
319  }
320 
321  side->inherit_data_from(*this);
322 
323  // Set the nodes
324  for (auto n : side->node_index_range())
325  side->set_node(n, this->node_ptr(Prism21::side_nodes_map[i][n]));
326 }
327 
328 
329 
330 std::unique_ptr<Elem> Prism21::build_edge_ptr (const unsigned int i)
331 {
332  return this->simple_build_edge_ptr<Edge3,Prism21>(i);
333 }
334 
335 
336 
337 void Prism21::build_edge_ptr (std::unique_ptr<Elem> & edge, const unsigned int i)
338 {
339  this->simple_build_edge_ptr<Prism21>(edge, i, EDGE3);
340 }
341 
342 
343 
344 void Prism21::connectivity(const unsigned int /*sc*/,
345  const IOPackage /*iop*/,
346  std::vector<dof_id_type> & /*conn*/) const
347 {
348  libmesh_not_implemented(); // FIXME RHS
349 
350 /*
351  libmesh_assert(_nodes);
352  libmesh_assert_less (sc, this->n_sub_elem());
353  libmesh_assert_not_equal_to (iop, INVALID_IO_PACKAGE);
354 
355  switch (iop)
356  {
357  case TECPLOT:
358  {
359  conn.resize(8);
360  switch (sc)
361  {
362 
363  case 0:
364  {
365  conn[0] = this->node_id(0)+1;
366  conn[1] = this->node_id(6)+1;
367  conn[2] = this->node_id(8)+1;
368  conn[3] = this->node_id(8)+1;
369  conn[4] = this->node_id(9)+1;
370  conn[5] = this->node_id(15)+1;
371  conn[6] = this->node_id(17)+1;
372  conn[7] = this->node_id(17)+1;
373 
374  return;
375  }
376 
377  case 1:
378  {
379  conn[0] = this->node_id(6)+1;
380  conn[1] = this->node_id(1)+1;
381  conn[2] = this->node_id(7)+1;
382  conn[3] = this->node_id(7)+1;
383  conn[4] = this->node_id(15)+1;
384  conn[5] = this->node_id(10)+1;
385  conn[6] = this->node_id(16)+1;
386  conn[7] = this->node_id(16)+1;
387 
388  return;
389  }
390 
391  case 2:
392  {
393  conn[0] = this->node_id(8)+1;
394  conn[1] = this->node_id(7)+1;
395  conn[2] = this->node_id(2)+1;
396  conn[3] = this->node_id(2)+1;
397  conn[4] = this->node_id(17)+1;
398  conn[5] = this->node_id(16)+1;
399  conn[6] = this->node_id(11)+1;
400  conn[7] = this->node_id(11)+1;
401 
402  return;
403  }
404 
405  case 3:
406  {
407  conn[0] = this->node_id(6)+1;
408  conn[1] = this->node_id(7)+1;
409  conn[2] = this->node_id(8)+1;
410  conn[3] = this->node_id(8)+1;
411  conn[4] = this->node_id(15)+1;
412  conn[5] = this->node_id(16)+1;
413  conn[6] = this->node_id(17)+1;
414  conn[7] = this->node_id(17)+1;
415 
416  return;
417  }
418 
419  case 4:
420  {
421  conn[0] = this->node_id(9)+1;
422  conn[1] = this->node_id(15)+1;
423  conn[2] = this->node_id(17)+1;
424  conn[3] = this->node_id(17)+1;
425  conn[4] = this->node_id(3)+1;
426  conn[5] = this->node_id(12)+1;
427  conn[6] = this->node_id(14)+1;
428  conn[7] = this->node_id(14)+1;
429 
430  return;
431  }
432 
433  case 5:
434  {
435  conn[0] = this->node_id(15)+1;
436  conn[1] = this->node_id(10)+1;
437  conn[2] = this->node_id(16)+1;
438  conn[3] = this->node_id(16)+1;
439  conn[4] = this->node_id(12)+1;
440  conn[5] = this->node_id(4)+1;
441  conn[6] = this->node_id(13)+1;
442  conn[7] = this->node_id(13)+1;
443 
444  return;
445  }
446 
447  case 6:
448  {
449  conn[0] = this->node_id(17)+1;
450  conn[1] = this->node_id(16)+1;
451  conn[2] = this->node_id(11)+1;
452  conn[3] = this->node_id(11)+1;
453  conn[4] = this->node_id(14)+1;
454  conn[5] = this->node_id(13)+1;
455  conn[6] = this->node_id(5)+1;
456  conn[7] = this->node_id(5)+1;
457 
458  return;
459  }
460 
461  case 7:
462  {
463  conn[0] = this->node_id(15)+1;
464  conn[1] = this->node_id(16)+1;
465  conn[2] = this->node_id(17)+1;
466  conn[3] = this->node_id(17)+1;
467  conn[4] = this->node_id(12)+1;
468  conn[5] = this->node_id(13)+1;
469  conn[6] = this->node_id(14)+1;
470  conn[7] = this->node_id(14)+1;
471 
472  return;
473  }
474 
475  default:
476  libmesh_error_msg("Invalid sc = " << sc);
477  }
478 
479  }
480 
481  case VTK:
482  {
483  // VTK now supports VTK_BIQUADRATIC_QUADRATIC_WEDGE directly
484  const unsigned int conn_size = 18;
485  conn.resize(conn_size);
486 
487  // VTK's VTK_BIQUADRATIC_QUADRATIC_WEDGE first 9 (vertex) and
488  // last 3 (mid-face) nodes match. The middle and top layers
489  // of mid-edge nodes are reversed from LibMesh's.
490  for (auto i : index_range(conn))
491  conn[i] = this->node_id(i);
492 
493  // top "ring" of mid-edge nodes
494  conn[9] = this->node_id(12);
495  conn[10] = this->node_id(13);
496  conn[11] = this->node_id(14);
497 
498  // middle "ring" of mid-edge nodes
499  conn[12] = this->node_id(9);
500  conn[13] = this->node_id(10);
501  conn[14] = this->node_id(11);
502 
503  return;
504  }
505 
506  default:
507  libmesh_error_msg("Unsupported IO package " << iop);
508  }
509 */
510 }
511 
512 
513 
514 
515 unsigned int Prism21::n_second_order_adjacent_vertices (const unsigned int n) const
516 {
517  switch (n)
518  {
519  case 6:
520  case 7:
521  case 8:
522  case 9:
523  case 10:
524  case 11:
525  case 12:
526  case 13:
527  case 14:
528  return 2;
529 
530  case 15:
531  case 16:
532  case 17:
533  return 4;
534 
535  case 18:
536  case 19:
537  return 3;
538 
539  case 20:
540  return 6;
541 
542  default:
543  libmesh_error_msg("Invalid node n = " << n);
544  }
545 }
546 
547 
548 
549 
550 
551 unsigned short int Prism21::second_order_adjacent_vertex (const unsigned int n,
552  const unsigned int v) const
553 {
554  libmesh_assert_greater_equal (n, this->n_vertices());
555  libmesh_assert_less (n, this->n_nodes());
556 
557  switch (n)
558  {
559  /*
560  * These nodes are unique to \p Prism20,
561  * let our _remaining_... matrix handle
562  * this.
563  */
564  case 15:
565  case 16:
566  case 17:
567  case 18:
568  case 19:
569  {
570  libmesh_assert_less (v, 4);
571  return _remaining_second_order_adjacent_vertices[n-15][v];
572  }
573 
574  /*
575  * for the bubble node the return value is simply v.
576  * Why? -- the user asks for the v-th adjacent vertex,
577  * from \p n_second_order_adjacent_vertices() there
578  * are 6 adjacent vertices, and these happen to be
579  * 0..5
580  */
581  case 20:
582  {
583  libmesh_assert_less (v, 6);
584  return static_cast<unsigned short int>(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  return static_cast<unsigned short int>(-1);
602 }
603 
604 
605 
606 const unsigned short int Prism21::_remaining_second_order_adjacent_vertices[5][4] =
607  {
608  { 0, 1, 3, 4}, // vertices adjacent to node 15
609  { 1, 2, 4, 5}, // vertices adjacent to node 16
610  { 0, 2, 3, 5}, // vertices adjacent to node 17
611  { 0, 1, 2, 99}, // vertices adjacent to node 18
612  { 3, 4, 5, 99} // vertices adjacent to node 19
613  };
614 
615 
616 
617 std::pair<unsigned short int, unsigned short int>
618 Prism21::second_order_child_vertex (const unsigned int n) const
619 {
620  libmesh_assert_greater_equal (n, this->n_vertices());
621  libmesh_assert_less (n, this->n_nodes());
622 
623  return std::pair<unsigned short int, unsigned short int>
624  (_second_order_vertex_child_number[n],
625  _second_order_vertex_child_index[n]);
626 }
627 
628 
629 
630 #ifdef LIBMESH_ENABLE_AMR
631 
632 // FIXME RHS
633 
634 const Real Prism21::_embedding_matrix[Prism21::num_children][Prism21::num_nodes][Prism21::num_nodes] =
635  {
636  // embedding matrix for child 0
637  {
638  // 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20
639  { 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // 0
640  { 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // 1
641  { 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // 2
642  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // 3
643  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0}, // 4
644  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0}, // 5
645  { 3/8., -1/8., 0, 0, 0, 0, 3/4., 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // 6
646  { 3/32., -1/32., -1/32., 0, 0, 0, 1/8., -1/8., 1/8., 0, 0, 0, 0, 0, 0, 0, 0, 0, 27/32., 0, 0}, // 7
647  { 3/8., 0, -1/8., 0, 0, 0, 0, 0, 3/4., 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // 8
648  { 3/8., 0, 0, -1/8., 0, 0, 0, 0, 0, 3/4., 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // 9
649  { 0, 0, 0, 0, 0, 0, 3/8., 0, 0, 0, 0, 0, -1/8., 0, 0, 3/4., 0, 0, 0, 0, 0}, // 10
650  { 0, 0, 0, 0, 0, 0, 0, 0, 3/8., 0, 0, 0, 0, 0, -1/8., 0, 0, 3/4., 0, 0, 0}, // 11
651  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 3/8., -1/8., 0, 0, 0, 0, 3/4., 0, 0, 0, 0, 0}, // 12
652  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 3/32., -1/32., -1/32., 0, 0, 0, 1/8., -1/8., 1/8., 0, 0, 27/32.}, // 13
653  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 3/8., 0, -1/8., 0, 0, 0, 0, 0, 3/4., 0, 0, 0}, // 14
654  { 9/64., -3/64., 0, -3/64., 1/64., 0, 9/32., 0, 0, 9/32., -3/32., 0, -3/32., 0, 0, 9/16., 0, 0, 0, 0, 0}, // 15
655  { 9/256., -3/256., -3/256., -3/256., 1/256., 1/256., 3/64., -3/64., 3/64., 9/128., -3/128., -3/128., -1/64., 1/64., -1/64., 3/32., -3/32., 3/32., 81/256., -27/256., 81/128.}, // 16
656  { 9/64., 0, -3/64., -3/64., 0, 1/64., 0, 0, 9/32., 9/32., 0, -3/32., 0, 0, -3/32., 0, 0, 9/16., 0, 0, 0}, // 17
657  { 5/r18, -1/r18, -1/r18, 0, 0, 0, 2/r9, -1/r9, 2/r9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1/2., 0, 0}, // 18
658  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 5/r18, -1/r18, -1/r18, 0, 0, 0, 2/r9, -1/r9, 2/r9, 0, 0, 1/2.}, // 19
659  { 5/r48, -1/r48, -1/r48, -5/r144, 1/r144, 1/r144, 1/r12, -1/r24, 1/r12, 5/r24, -1/r24, -1/r24, -1/r36, 1/r72, -1/r36, 1/r6, -1/r12, 1/r6, 3/16., -1/16., 3/8.} // 20
660  },
661  // embedding matrix for child 1
662  {
663  // 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20
664  { 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // 0
665  { 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // 1
666  { 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // 2
667  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0}, // 3
668  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // 4
669  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0}, // 5
670  { -1/8., 3/8., 0, 0, 0, 0, 3/4., 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // 6
671  { 0, 3/8., -1/8., 0, 0, 0, 0, 3/4., 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // 7
672  { -1/32., 3/32., -1/32., 0, 0, 0, 1/8., 1/8., -1/8., 0, 0, 0, 0, 0, 0, 0, 0, 0, 27/32., 0, 0}, // 8
673  { 0, 0, 0, 0, 0, 0, 3/8., 0, 0, 0, 0, 0, -1/8., 0, 0, 3/4., 0, 0, 0, 0, 0}, // 9
674  { 0, 3/8., 0, 0, -1/8., 0, 0, 0, 0, 0, 3/4., 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // 10
675  { 0, 0, 0, 0, 0, 0, 0, 3/8., 0, 0, 0, 0, 0, -1/8., 0, 0, 3/4., 0, 0, 0, 0}, // 11
676  { 0, 0, 0, 0, 0, 0, 0, 0, 0, -1/8., 3/8., 0, 0, 0, 0, 3/4., 0, 0, 0, 0, 0}, // 12
677  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3/8., -1/8., 0, 0, 0, 0, 3/4., 0, 0, 0, 0}, // 13
678  { 0, 0, 0, 0, 0, 0, 0, 0, 0, -1/32., 3/32., -1/32., 0, 0, 0, 1/8., 1/8., -1/8., 0, 0, 27/32.}, // 14
679  { -3/64., 9/64., 0, 1/64., -3/64., 0, 9/32., 0, 0, -3/32., 9/32., 0, -3/32., 0, 0, 9/16., 0, 0, 0, 0, 0}, // 15
680  { 0, 9/64., -3/64., 0, -3/64., 1/64., 0, 9/32., 0, 0, 9/32., -3/32., 0, -3/32., 0, 0, 9/16., 0, 0, 0, 0}, // 16
681  { -3/256., 9/256., -3/256., 1/256., -3/256., 1/256., 3/64., 3/64., -3/64., -3/128., 9/128., -3/128., -1/64., -1/64., 1/64., 3/32., 3/32., -3/32., 81/256., -27/256., 81/128.}, // 17
682  { -1/r18, 5/r18, -1/r18, 0, 0, 0, 2/r9, 2/r9, -1/r9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1/2., 0, 0}, // 18
683  { 0, 0, 0, 0, 0, 0, 0, 0, 0, -1/r18, 5/r18, -1/r18, 0, 0, 0, 2/r9, 2/r9, -1/r9, 0, 0, 1/2.}, // 19
684  { -1/r48, 5/r48, -1/r48, 1/r144, -5/r144, 1/r144, 1/r12, 1/r12, -1/r24, -1/r24, 5/r24, -1/r24, -1/r36, -1/r36, 1/r72, 1/r6, 1/r6, -1/r12, 3/16., -1/16., 3/8.} // 20
685  },
686  // embedding matrix for child 2
687  {
688  // 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20
689  { 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // 0
690  { 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // 1
691  { 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // 2
692  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0}, // 3
693  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0}, // 4
694  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // 5
695  { -1/32., -1/32., 3/32., 0, 0, 0, -1/8., 1/8., 1/8., 0, 0, 0, 0, 0, 0, 0, 0, 0, 27/32., 0, 0}, // 6
696  { 0, -1/8., 3/8., 0, 0, 0, 0, 3/4., 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // 7
697  { -1/8., 0, 3/8., 0, 0, 0, 0, 0, 3/4., 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // 8
698  { 0, 0, 0, 0, 0, 0, 0, 0, 3/8., 0, 0, 0, 0, 0, -1/8., 0, 0, 3/4., 0, 0, 0}, // 9
699  { 0, 0, 0, 0, 0, 0, 0, 3/8., 0, 0, 0, 0, 0, -1/8., 0, 0, 3/4., 0, 0, 0, 0}, // 10
700  { 0, 0, 3/8., 0, 0, -1/8., 0, 0, 0, 0, 0, 3/4., 0, 0, 0, 0, 0, 0, 0, 0, 0}, // 11
701  { 0, 0, 0, 0, 0, 0, 0, 0, 0, -1/32., -1/32., 3/32., 0, 0, 0, -1/8., 1/8., 1/8., 0, 0, 27/32.}, // 12
702  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1/8., 3/8., 0, 0, 0, 0, 3/4., 0, 0, 0, 0}, // 13
703  { 0, 0, 0, 0, 0, 0, 0, 0, 0, -1/8., 0, 3/8., 0, 0, 0, 0, 0, 3/4., 0, 0, 0}, // 14
704  { -3/256., -3/256., 9/256., 1/256., 1/256., -3/256., -3/64., 3/64., 3/64., -3/128., -3/128., 9/128., 1/64., -1/64., -1/64., -3/32., 3/32., 3/32., 81/256., -27/256., 81/128.}, // 15
705  { 0, -3/64., 9/64., 0, 1/64., -3/64., 0, 9/32., 0, 0, -3/32., 9/32., 0, -3/32., 0, 0, 9/16., 0, 0, 0, 0}, // 16
706  { -3/64., 0, 9/64., 1/64., 0, -3/64., 0, 0, 9/32., -3/32., 0, 9/32., 0, 0, -3/32., 0, 0, 9/16., 0, 0, 0}, // 17
707  { -1/r18, -1/r18, 5/r18, 0, 0, 0, -1/r9, 2/r9, 2/r9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1/2., 0, 0}, // 18
708  { 0, 0, 0, 0, 0, 0, 0, 0, 0, -1/r18, -1/r18, 5/r18, 0, 0, 0, -1/r9, 2/r9, 2/r9, 0, 0, 1/2.}, // 19
709  { -1/r48, -1/r48, 5/r48, 1/r144, 1/r144, -5/r144, -1/r24, 1/r12, 1/r12, -1/r24, -1/r24, 5/r24, 1/r72, -1/r36, -1/r36, -1/r12, 1/r6, 1/r6, 3/16., -1/16., 3/8.} // 20
710  },
711  // embedding matrix for child 3
712  {
713  // 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20
714  { 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // 0
715  { 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // 1
716  { 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // 2
717  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0}, // 3
718  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0}, // 4
719  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0}, // 5
720  { -1/32., 3/32., -1/32., 0, 0, 0, 1/8., 1/8., -1/8., 0, 0, 0, 0, 0, 0, 0, 0, 0, 27/32., 0, 0}, // 6
721  { -1/32., -1/32., 3/32., 0, 0, 0, -1/8., 1/8., 1/8., 0, 0, 0, 0, 0, 0, 0, 0, 0, 27/32., 0, 0}, // 7
722  { 3/32., -1/32., -1/32., 0, 0, 0, 1/8., -1/8., 1/8., 0, 0, 0, 0, 0, 0, 0, 0, 0, 27/32., 0, 0}, // 8
723  { 0, 0, 0, 0, 0, 0, 3/8., 0, 0, 0, 0, 0, -1/8., 0, 0, 3/4., 0, 0, 0, 0, 0}, // 9
724  { 0, 0, 0, 0, 0, 0, 0, 3/8., 0, 0, 0, 0, 0, -1/8., 0, 0, 3/4., 0, 0, 0, 0}, // 10
725  { 0, 0, 0, 0, 0, 0, 0, 0, 3/8., 0, 0, 0, 0, 0, -1/8., 0, 0, 3/4., 0, 0, 0}, // 11
726  { 0, 0, 0, 0, 0, 0, 0, 0, 0, -1/32., 3/32., -1/32., 0, 0, 0, 1/8., 1/8., -1/8., 0, 0, 27/32.}, // 12
727  { 0, 0, 0, 0, 0, 0, 0, 0, 0, -1/32., -1/32., 3/32., 0, 0, 0, -1/8., 1/8., 1/8., 0, 0, 27/32.}, // 13
728  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 3/32., -1/32., -1/32., 0, 0, 0, 1/8., -1/8., 1/8., 0, 0, 27/32.}, // 14
729  { -3/256., 9/256., -3/256., 1/256., -3/256., 1/256., 3/64., 3/64., -3/64., -3/128., 9/128., -3/128., -1/64., -1/64., 1/64., 3/32., 3/32., -3/32., 81/256., -27/256., 81/128.}, // 15
730  { -3/256., -3/256., 9/256., 1/256., 1/256., -3/256., -3/64., 3/64., 3/64., -3/128., -3/128., 9/128., 1/64., -1/64., -1/64., -3/32., 3/32., 3/32., 81/256., -27/256., 81/128.}, // 16
731  { 9/256., -3/256., -3/256., -3/256., 1/256., 1/256., 3/64., -3/64., 3/64., 9/128., -3/128., -3/128., -1/64., 1/64., -1/64., 3/32., -3/32., 3/32., 81/256., -27/256., 81/128.}, // 17
732  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0}, // 18
733  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1}, // 19
734  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3/8., -1/8., 3/4.} // 20
735  },
736  // embedding matrix for child 4
737  {
738  // 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20
739  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // 0
740  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0}, // 1
741  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0}, // 2
742  { 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // 3
743  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0}, // 4
744  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0}, // 5
745  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 3/8., -1/8., 0, 0, 0, 0, 3/4., 0, 0, 0, 0, 0}, // 6
746  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 3/32., -1/32., -1/32., 0, 0, 0, 1/8., -1/8., 1/8., 0, 0, 27/32.}, // 7
747  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 3/8., 0, -1/8., 0, 0, 0, 0, 0, 3/4., 0, 0, 0}, // 8
748  { -1/8., 0, 0, 3/8., 0, 0, 0, 0, 0, 3/4., 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // 9
749  { 0, 0, 0, 0, 0, 0, -1/8., 0, 0, 0, 0, 0, 3/8., 0, 0, 3/4., 0, 0, 0, 0, 0}, // 10
750  { 0, 0, 0, 0, 0, 0, 0, 0, -1/8., 0, 0, 0, 0, 0, 3/8., 0, 0, 3/4., 0, 0, 0}, // 11
751  { 0, 0, 0, 3/8., -1/8., 0, 0, 0, 0, 0, 0, 0, 3/4., 0, 0, 0, 0, 0, 0, 0, 0}, // 12
752  { 0, 0, 0, 3/32., -1/32., -1/32., 0, 0, 0, 0, 0, 0, 1/8., -1/8., 1/8., 0, 0, 0, 0, 27/32., 0}, // 13
753  { 0, 0, 0, 3/8., 0, -1/8., 0, 0, 0, 0, 0, 0, 0, 0, 3/4., 0, 0, 0, 0, 0, 0}, // 14
754  { -3/64., 1/64., 0, 9/64., -3/64., 0, -3/32., 0, 0, 9/32., -3/32., 0, 9/32., 0, 0, 9/16., 0, 0, 0, 0, 0}, // 15
755  { -3/256., 1/256., 1/256., 9/256., -3/256., -3/256., -1/64., 1/64., -1/64., 9/128., -3/128., -3/128., 3/64., -3/64., 3/64., 3/32., -3/32., 3/32., -27/256., 81/256., 81/128.}, // 16
756  { -3/64., 0, 1/64., 9/64., 0, -3/64., 0, 0, -3/32., 9/32., 0, -3/32., 0, 0, 9/32., 0, 0, 9/16., 0, 0, 0}, // 17
757  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 5/r18, -1/r18, -1/r18, 0, 0, 0, 2/r9, -1/r9, 2/r9, 0, 0, 1/2.}, // 18
758  { 0, 0, 0, 5/r18, -1/r18, -1/r18, 0, 0, 0, 0, 0, 0, 2/r9, -1/r9, 2/r9, 0, 0, 0, 0, 1/2., 0}, // 19
759  { -5/r144, 1/r144, 1/r144, 5/r48, -1/r48, -1/r48, -1/r36, 1/r72, -1/r36, 5/r24, -1/r24, -1/r24, 1/r12, -1/r24, 1/r12, 1/r6, -1/r12, 1/r6, -1/16., 3/16., 3/8.} // 20
760  },
761  // embedding matrix for child 5
762  {
763  // 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20
764  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0}, // 0
765  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // 1
766  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0}, // 2
767  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0}, // 3
768  { 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // 4
769  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0}, // 5
770  { 0, 0, 0, 0, 0, 0, 0, 0, 0, -1/8., 3/8., 0, 0, 0, 0, 3/4., 0, 0, 0, 0, 0}, // 6
771  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3/8., -1/8., 0, 0, 0, 0, 3/4., 0, 0, 0, 0}, // 7
772  { 0, 0, 0, 0, 0, 0, 0, 0, 0, -1/32., 3/32., -1/32., 0, 0, 0, 1/8., 1/8., -1/8., 0, 0, 27/32.}, // 8
773  { 0, 0, 0, 0, 0, 0, -1/8., 0, 0, 0, 0, 0, 3/8., 0, 0, 3/4., 0, 0, 0, 0, 0}, // 9
774  { 0, -1/8., 0, 0, 3/8., 0, 0, 0, 0, 0, 3/4., 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // 10
775  { 0, 0, 0, 0, 0, 0, 0, -1/8., 0, 0, 0, 0, 0, 3/8., 0, 0, 3/4., 0, 0, 0, 0}, // 11
776  { 0, 0, 0, -1/8., 3/8., 0, 0, 0, 0, 0, 0, 0, 3/4., 0, 0, 0, 0, 0, 0, 0, 0}, // 12
777  { 0, 0, 0, 0, 3/8., -1/8., 0, 0, 0, 0, 0, 0, 0, 3/4., 0, 0, 0, 0, 0, 0, 0}, // 13
778  { 0, 0, 0, -1/32., 3/32., -1/32., 0, 0, 0, 0, 0, 0, 1/8., 1/8., -1/8., 0, 0, 0, 0, 27/32., 0}, // 14
779  { 1/64., -3/64., 0, -3/64., 9/64., 0, -3/32., 0, 0, -3/32., 9/32., 0, 9/32., 0, 0, 9/16., 0, 0, 0, 0, 0}, // 15
780  { 0, -3/64., 1/64., 0, 9/64., -3/64., 0, -3/32., 0, 0, 9/32., -3/32., 0, 9/32., 0, 0, 9/16., 0, 0, 0, 0}, // 16
781  { 1/256., -3/256., 1/256., -3/256., 9/256., -3/256., -1/64., -1/64., 1/64., -3/128., 9/128., -3/128., 3/64., 3/64., -3/64., 3/32., 3/32., -3/32., -27/256., 81/256., 81/128.}, // 17
782  { 0, 0, 0, 0, 0, 0, 0, 0, 0, -1/r18, 5/r18, -1/r18, 0, 0, 0, 2/r9, 2/r9, -1/r9, 0, 0, 1/2.}, // 18
783  { 0, 0, 0, -1/r18, 5/r18, -1/r18, 0, 0, 0, 0, 0, 0, 2/r9, 2/r9, -1/r9, 0, 0, 0, 0, 1/2., 0}, // 19
784  { 1/r144, -5/r144, 1/r144, -1/r48, 5/r48, -1/r48, -1/r36, -1/r36, 1/r72, -1/r24, 5/r24, -1/r24, 1/r12, 1/r12, -1/r24, 1/r6, 1/r6, -1/r12, -1/16., 3/16., 3/8.} // 20
785  },
786  // embedding matrix for child 6
787  {
788  // 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20
789  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0}, // 0
790  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0}, // 1
791  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // 2
792  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0}, // 3
793  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0}, // 4
794  { 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // 5
795  { 0, 0, 0, 0, 0, 0, 0, 0, 0, -1/32., -1/32., 3/32., 0, 0, 0, -1/8., 1/8., 1/8., 0, 0, 27/32.}, // 6
796  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1/8., 3/8., 0, 0, 0, 0, 3/4., 0, 0, 0, 0}, // 7
797  { 0, 0, 0, 0, 0, 0, 0, 0, 0, -1/8., 0, 3/8., 0, 0, 0, 0, 0, 3/4., 0, 0, 0}, // 8
798  { 0, 0, 0, 0, 0, 0, 0, 0, -1/8., 0, 0, 0, 0, 0, 3/8., 0, 0, 3/4., 0, 0, 0}, // 9
799  { 0, 0, 0, 0, 0, 0, 0, -1/8., 0, 0, 0, 0, 0, 3/8., 0, 0, 3/4., 0, 0, 0, 0}, // 10
800  { 0, 0, -1/8., 0, 0, 3/8., 0, 0, 0, 0, 0, 3/4., 0, 0, 0, 0, 0, 0, 0, 0, 0}, // 11
801  { 0, 0, 0, -1/32., -1/32., 3/32., 0, 0, 0, 0, 0, 0, -1/8., 1/8., 1/8., 0, 0, 0, 0, 27/32., 0}, // 12
802  { 0, 0, 0, 0, -1/8., 3/8., 0, 0, 0, 0, 0, 0, 0, 3/4., 0, 0, 0, 0, 0, 0, 0}, // 13
803  { 0, 0, 0, -1/8., 0, 3/8., 0, 0, 0, 0, 0, 0, 0, 0, 3/4., 0, 0, 0, 0, 0, 0}, // 14
804  { 1/256., 1/256., -3/256., -3/256., -3/256., 9/256., 1/64., -1/64., -1/64., -3/128., -3/128., 9/128., -3/64., 3/64., 3/64., -3/32., 3/32., 3/32., -27/256., 81/256., 81/128.}, // 15
805  { 0, 1/64., -3/64., 0, -3/64., 9/64., 0, -3/32., 0, 0, -3/32., 9/32., 0, 9/32., 0, 0, 9/16., 0, 0, 0, 0}, // 16
806  { 1/64., 0, -3/64., -3/64., 0, 9/64., 0, 0, -3/32., -3/32., 0, 9/32., 0, 0, 9/32., 0, 0, 9/16., 0, 0, 0}, // 17
807  { 0, 0, 0, 0, 0, 0, 0, 0, 0, -1/r18, -1/r18, 5/r18, 0, 0, 0, -1/r9, 2/r9, 2/r9, 0, 0, 1/2.}, // 18
808  { 0, 0, 0, -1/r18, -1/r18, 5/r18, 0, 0, 0, 0, 0, 0, -1/r9, 2/r9, 2/r9, 0, 0, 0, 0, 1/2., 0}, // 19
809  { 1/r144, 1/r144, -5/r144, -1/r48, -1/r48, 5/r48, 1/r72, -1/r36, -1/r36, -1/r24, -1/r24, 5/r24, -1/r24, 1/r12, 1/r12, -1/r12, 1/r6, 1/r6, -1/16., 3/16., 3/8.} // 20
810  },
811  // embedding matrix for child 7
812  {
813  // 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20
814  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0}, // 0
815  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0}, // 1
816  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0}, // 2
817  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0}, // 3
818  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0}, // 4
819  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0}, // 5
820  { 0, 0, 0, 0, 0, 0, 0, 0, 0, -1/32., 3/32., -1/32., 0, 0, 0, 1/8., 1/8., -1/8., 0, 0, 27/32.}, // 6
821  { 0, 0, 0, 0, 0, 0, 0, 0, 0, -1/32., -1/32., 3/32., 0, 0, 0, -1/8., 1/8., 1/8., 0, 0, 27/32.}, // 7
822  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 3/32., -1/32., -1/32., 0, 0, 0, 1/8., -1/8., 1/8., 0, 0, 27/32.}, // 8
823  { 0, 0, 0, 0, 0, 0, -1/8., 0, 0, 0, 0, 0, 3/8., 0, 0, 3/4., 0, 0, 0, 0, 0}, // 9
824  { 0, 0, 0, 0, 0, 0, 0, -1/8., 0, 0, 0, 0, 0, 3/8., 0, 0, 3/4., 0, 0, 0, 0}, // 10
825  { 0, 0, 0, 0, 0, 0, 0, 0, -1/8., 0, 0, 0, 0, 0, 3/8., 0, 0, 3/4., 0, 0, 0}, // 11
826  { 0, 0, 0, -1/32., 3/32., -1/32., 0, 0, 0, 0, 0, 0, 1/8., 1/8., -1/8., 0, 0, 0, 0, 27/32., 0}, // 12
827  { 0, 0, 0, -1/32., -1/32., 3/32., 0, 0, 0, 0, 0, 0, -1/8., 1/8., 1/8., 0, 0, 0, 0, 27/32., 0}, // 13
828  { 0, 0, 0, 3/32., -1/32., -1/32., 0, 0, 0, 0, 0, 0, 1/8., -1/8., 1/8., 0, 0, 0, 0, 27/32., 0}, // 14
829  { 1/256., -3/256., 1/256., -3/256., 9/256., -3/256., -1/64., -1/64., 1/64., -3/128., 9/128., -3/128., 3/64., 3/64., -3/64., 3/32., 3/32., -3/32., -27/256., 81/256., 81/128.}, // 15
830  { 1/256., 1/256., -3/256., -3/256., -3/256., 9/256., 1/64., -1/64., -1/64., -3/128., -3/128., 9/128., -3/64., 3/64., 3/64., -3/32., 3/32., 3/32., -27/256., 81/256., 81/128.}, // 16
831  { -3/256., 1/256., 1/256., 9/256., -3/256., -3/256., -1/64., 1/64., -1/64., 9/128., -3/128., -3/128., 3/64., -3/64., 3/64., 3/32., -3/32., 3/32., -27/256., 81/256., 81/128.}, // 17
832  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1}, // 18
833  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0}, // 19
834  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1/8., 3/8., 3/4.} // 20
835  }
836  };
837 
838 const std::vector<std::pair<unsigned char, unsigned char>>
839  Prism21::_parent_bracketing_nodes[Prism21::num_children][Prism21::num_nodes] =
840  {
841  // Child 0
842  { {},{{0,1}},{{0,2}},{{0,3}},{{0,4},{1,3},{6,12},{9,10}},{{0,5},{2,3},{8,14},{9,11}},
843  {{0,6}},{{6,8}},{{0,8}},{{0,9}},{{6,15}},{{8,17}},{{9,15}},{{15,17}},{{9,17}},
844  {{0,15},{6,9}},{{6,17},{8,15}},{{0,17},{8,9}},{{0,18}},{{9,20}},{{0,20},{9,18}} },
845  // Child 1
846  { {{0,1}}, {},{{1,2}},{{0,4},{1,3},{6,12},{9,10}},{{1,4}},{{1,5},{2,4},{7,13},{10,11}},
847  {{1,6}},{{1,7}},{{6,7}},{{6,15}},{{1,10}},{{7,16}},{{10,15}},{{10,16}},{{15,16}},
848  {{1,15},{6,10}},{{1,16},{7,10}},{{6,16},{7,15}},{{1,18}},{{10,20}},{{1,20},{10,18}} },
849  // Child 2
850  { {{0,2}},{{1,2}}, {},{{0,5},{2,3},{8,14},{9,11}},{{1,5},{2,4},{7,13},{10,11}},{{2,5}},
851  {{7,8}},{{2,7}},{{2,8}},{{8,17}},{{7,16}},{{2,11}},{{16,17}},{{11,16}},{{11,17}},
852  {{7,17},{8,16}},{{2,16},{7,11}},{{2,17},{8,11}},{{2,18}},{{11,20}},{{2,20},{11,18}} },
853  // Child 3
854  { {{0,1}},{{1,2}},{{0,2}},{{0,4},{1,3},{6,12},{9,10}},{{1,5},{2,4},{7,13},{10,11}},{{0,5},{2,3},{8,14},{9,11}},
855  {{6,7}},{{7,8}},{{6,8}},{{6,15}},{{7,16}},{{8,17}},{{15,16}},{{16,17}},{{15,17}},
856  {{6,16},{7,15}},{{7,17},{8,16}},{{6,17},{8,15}},{{0,7},{1,8},{2,6}},{{9,16},{10,17},{11,15}},{{0,16},{1,17},{2,15}} },
857  // Child 4
858  { {{0,3}},{{0,4},{1,3},{6,11},{9,10}},{{0,5},{2,3},{8,14},{9,11}}, {},{{3,4}},{{3,5}},
859  {{9,15}},{{15,17}},{{9,17}},{{3,9}},{{12,15}},{{14,17}},{{3,12}},{{12,14}},{{3,14}},
860  {{3,15},{9,12}},{{12,17},{14,15}},{{3,17},{9,14}},{{9,20}},{{3,19}},{{9,19},{3,20}} },
861  // Child 5
862  { {{0,4},{1,3},{6,12},{9,10}},{{1,4}},{{1,5},{2,4},{7,13},{10,11}},{{3,4}}, {},{{4,5}},
863  {{10,15}},{{10,16}},{{15,16}},{{12,15}},{{4,10}},{{13,16}},{{4,12}},{{4,13}},{{12,13}},
864  {{4,15},{10,12}},{{4,16},{10,13}},{{12,16},{13,15}},{{10,20}},{{4,19}},{{10,19},{4,20}} },
865  // Child 6
866  { {{0,5},{2,3},{8,14},{9,11}},{{1,5},{2,4},{7,13},{10,11}},{{2,5}},{{3,5}},{{4,5}}, {},
867  {{16,17}},{{11,16}},{{11,17}},{{14,17}},{{13,16}},{{5,11}},{{13,14}},{{5,13}},{{5,14}},
868  {{13,17},{14,16}},{{5,16},{11,13}},{{5,17},{11,14}},{{11,20}},{{5,19}},{{11,19},{5,20}} },
869  // Child 7
870  { {{0,4},{1,3},{6,12},{9,10}},{{1,5},{2,4},{7,13},{10,11}},{{0,5},{2,3},{8,14},{9,11}},{{3,4}},{{4,5}},{{3,5}},
871  {{15,16}},{{16,17}},{{15,17}},{{12,15}},{{13,16}},{{14,17}},{{12,13}},{{13,14}},{{12,14}},
872  {{12,16},{13,15}},{{13,17},{14,16}},{{12,17},{14,15}},{{9,16},{10,17},{11,15}},{{3,13},{4,14},{5,12}},{{3,16},{4,17},{5,15}} }
873  };
874 
875 #endif
876 
877 
878 void
879 Prism21::permute(unsigned int perm_num)
880 {
881  libmesh_assert_less (perm_num, 6);
882  const unsigned int side = perm_num % 2;
883  const unsigned int rotate = perm_num / 2;
884 
885  for (unsigned int i = 0; i != rotate; ++i)
886  {
887  swap3nodes(0,1,2);
888  swap3nodes(3,4,5);
889  swap3nodes(6,7,8);
890  swap3nodes(9,10,11);
891  swap3nodes(12,13,14);
892  swap3nodes(15,16,17);
893  swap3neighbors(1,2,3);
894  }
895 
896  switch (side) {
897  case 0:
898  break;
899  case 1:
900  swap2nodes(1,3);
901  swap2nodes(0,4);
902  swap2nodes(2,5);
903  swap2nodes(6,12);
904  swap2nodes(9,10);
905  swap2nodes(7,14);
906  swap2nodes(8,13);
907  swap2nodes(16,17);
908  swap2nodes(18,19);
909  swap2neighbors(0,4);
910  swap2neighbors(2,3);
911  break;
912  default:
913  libmesh_error();
914  }
915 }
916 
917 
918 void
919 Prism21::flip(BoundaryInfo * boundary_info)
920 {
921  libmesh_assert(boundary_info);
922 
923  swap2nodes(0,1);
924  swap2nodes(3,4);
925  swap2nodes(7,8);
926  swap2nodes(9,10);
927  swap2nodes(13,14);
928  swap2nodes(16,17);
929  swap2neighbors(2,3);
930  swap2boundarysides(2,3,boundary_info);
931  swap2boundaryedges(0,1,boundary_info);
932  swap2boundaryedges(3,4,boundary_info);
933  swap2boundaryedges(7,8,boundary_info);
934 }
935 
936 
937 unsigned int Prism21::center_node_on_side(const unsigned short side) const
938 {
939  libmesh_assert_less (side, Prism21::num_sides);
940  if (side >= 1 && side <= 3)
941  return side + 14;
942  if (side == 4)
943  return 19;
944  return 18;
945 }
946 
947 
948 ElemType
949 Prism21::side_type (const unsigned int s) const
950 {
951  libmesh_assert_less (s, 5);
952  if (s == 0 || s == 4)
953  return TRI7;
954  return QUAD9;
955 }
956 
957 
958 } // namespace libMesh
libMesh::Prism21::build_side_ptr
virtual std::unique_ptr< Elem > build_side_ptr(const unsigned int i) override
Builds a QUAD9 or TRI6 built coincident with face i.
Definition: cell_prism21.C:250
libMesh::ElemType
ElemType
Defines an enum for geometric element types.
Definition: enum_elem_type.h:33
libMesh::Elem::swap2boundaryedges
void swap2boundaryedges(unsigned short e1, unsigned short e2, BoundaryInfo *boundary_info) const
Swaps two edges in boundary_info, if it is non-null.
Definition: elem.C:3550
libMesh::Order
Order
defines an enum for polynomial orders.
Definition: enum_order.h:40
libMesh::Prism::num_edges
static const int num_edges
Definition: cell_prism.h:74
libMesh::Elem::key
virtual dof_id_type key() const
Definition: elem.C:753
libMesh::Prism21::local_edge_node
virtual unsigned int local_edge_node(unsigned int edge, unsigned int edge_node) const override
Definition: cell_prism21.C:239
libMesh::Prism::n_sides
virtual unsigned int n_sides() const override final
Definition: cell_prism.h:86
libMesh::Prism21::is_face
virtual bool is_face(const unsigned int i) const override
Definition: cell_prism21.C:95
libMesh::Prism21::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_prism21.h:251
libMesh::IOPackage
IOPackage
libMesh interfaces with several different software packages for the purposes of creating, reading, and writing mesh files.
Definition: enum_io_package.h:37
libMesh::Elem::swap2boundarysides
void swap2boundarysides(unsigned short s1, unsigned short s2, BoundaryInfo *boundary_info) const
Swaps two sides in boundary_info, if it is non-null.
Definition: elem.C:3534
libMesh::Prism21::is_edge
virtual bool is_edge(const unsigned int i) const override
Definition: cell_prism21.C:86
libMesh
The libMesh namespace provides an interface to certain functionality in the library.
libMesh::Prism21::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_prism21.h:257
libMesh::Prism21::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_prism21.C:618
libMesh::Prism::num_children
static const int num_children
Definition: cell_prism.h:75
libMesh::Prism21::nodes_per_edge
static const int nodes_per_edge
Definition: cell_prism21.h:245
libMesh::Prism21::nodes_per_side
static const int nodes_per_side
Definition: cell_prism21.h:244
libMesh::Prism::_second_order_adjacent_vertices
static const unsigned short int _second_order_adjacent_vertices[9][2]
Matrix that tells which vertices define the location of mid-side (or second-order) nodes...
Definition: cell_prism.h:199
libMesh::Elem::swap3neighbors
void swap3neighbors(unsigned int n1, unsigned int n2, unsigned int n3)
Swaps three neighbor_ptrs, "rotating" them.
Definition: elem.h:2133
libMesh::Prism21::n_nodes
virtual unsigned int n_nodes() const override
Definition: cell_prism21.h:108
libMesh::Elem::swap3nodes
void swap3nodes(unsigned int n1, unsigned int n2, unsigned int n3)
Swaps three node_ptrs, "rotating" them.
Definition: elem.h:2124
libMesh::Prism21::is_node_on_edge
virtual bool is_node_on_edge(const unsigned int n, const unsigned int e) const override
Definition: cell_prism21.C:128
libMesh::Prism21::flip
virtual void flip(BoundaryInfo *) override final
Flips the element (by swapping node and neighbor pointers) to have a mapping Jacobian of opposite sig...
Definition: cell_prism21.C:919
libMesh::Elem::swap2nodes
void swap2nodes(unsigned int n1, unsigned int n2)
Swaps two node_ptrs.
Definition: elem.h:2092
libMesh::Prism::_second_order_vertex_child_index
static const unsigned short int _second_order_vertex_child_index[18]
Vector that names the child vertex index for each second order node.
Definition: cell_prism.h:209
libMesh::Prism21::num_nodes
static const int num_nodes
Geometric constants for Prism21.
Definition: cell_prism21.h:243
libMesh::Prism21::is_node_on_side
virtual bool is_node_on_side(const unsigned int n, const unsigned int s) const override
Definition: cell_prism21.C:104
libMesh::Prism21::connectivity
virtual void connectivity(const unsigned int sc, const IOPackage iop, std::vector< dof_id_type > &conn) const override
Definition: cell_prism21.C:344
libMesh::Prism21::second_order_adjacent_vertex
virtual unsigned short int second_order_adjacent_vertex(const unsigned int n, const unsigned int v) const override
Definition: cell_prism21.C:551
libMesh::Prism21::has_affine_map
virtual bool has_affine_map() const override
Definition: cell_prism21.C:139
libMesh::Prism21::build_edge_ptr
virtual std::unique_ptr< Elem > build_edge_ptr(const unsigned int i) override
Builds a EDGE3 or INFEDGE2 built coincident with edge i.
Definition: cell_prism21.C:330
libMesh::Prism21::nodes_on_side
virtual std::vector< unsigned int > nodes_on_side(const unsigned int s) const override
Definition: cell_prism21.C:114
libMesh::Prism21::side_type
ElemType side_type(const unsigned int s) const override final
Definition: cell_prism21.C:949
libMesh::BoundaryInfo
The BoundaryInfo class contains information relevant to boundary conditions including storing faces...
Definition: boundary_info.h:57
libMesh::libmesh_assert
libmesh_assert(ctx)
libMesh::Prism21::n_second_order_adjacent_vertices
virtual unsigned int n_second_order_adjacent_vertices(const unsigned int) const override
Definition: cell_prism21.C:515
libMesh::Elem::affine_tol
static constexpr Real affine_tol
Default tolerance to use in has_affine_map().
Definition: elem.h:2061
libMesh::Prism21::default_order
virtual Order default_order() const override
Definition: cell_prism21.C:185
libMesh::Prism21::center_node_on_side
unsigned int center_node_on_side(const unsigned short side) const override final
Definition: cell_prism21.C:937
libMesh::Prism21::_remaining_second_order_adjacent_vertices
static const unsigned short int _remaining_second_order_adjacent_vertices[5][4]
Matrix that tells which vertices define the location of mid-side (or second-order) nodes...
Definition: cell_prism21.h:320
libMesh::Prism::n_vertices
virtual unsigned int n_vertices() const override final
Definition: cell_prism.h:91
libMesh::MeshTools::Modification::rotate
RealTensorValue rotate(MeshBase &mesh, const Real phi, const Real theta=0., const Real psi=0.)
Rotates the mesh in the xy plane.
Definition: mesh_modification.C:344
libMesh::Prism::_second_order_vertex_child_number
static const unsigned short int _second_order_vertex_child_number[18]
Vector that names a child sharing each second order node.
Definition: cell_prism.h:204
libMesh::Prism21::_parent_bracketing_nodes
static const std::vector< std::pair< unsigned char, unsigned char > > _parent_bracketing_nodes[num_children][num_nodes]
Pairs of nodes that bracket child nodes when doing mesh refinement.
Definition: cell_prism21.h:305
libMesh::Elem::swap2neighbors
void swap2neighbors(unsigned int n1, unsigned int n2)
Swaps two neighbor_ptrs.
Definition: elem.h:2102
libMesh::Prism21::_embedding_matrix
static const Real _embedding_matrix[num_children][num_nodes][num_nodes]
Matrix that computes new nodal locations/solution values from current nodes/solution.
Definition: cell_prism21.h:298
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:144
libMesh::Prism::n_edges
virtual unsigned int n_edges() const override final
Definition: cell_prism.h:96
libMesh::Prism::num_sides
static const int num_sides
Geometric constants for all Prisms.
Definition: cell_prism.h:73
libMesh::Prism21::permute
virtual void permute(unsigned int perm_num) override final
Permutes the element (by swapping node and neighbor pointers) according to the specified index...
Definition: cell_prism21.C:879
libMesh::Prism21::local_side_node
virtual unsigned int local_side_node(unsigned int side, unsigned int side_node) const override
Definition: cell_prism21.C:223
libMesh::Prism21::is_vertex
virtual bool is_vertex(const unsigned int i) const override
Definition: cell_prism21.C:79
libMesh::Elem::compute_key
static dof_id_type compute_key(dof_id_type n0)
Definition: elem.h:3294
libMesh::Prism21::nodes_on_edge
virtual std::vector< unsigned int > nodes_on_edge(const unsigned int e) const override
Definition: cell_prism21.C:122
libMesh::Point
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
libMesh::Elem::node_id
dof_id_type node_id(const unsigned int i) const
Definition: elem.h:2475
libMesh::Elem::point
const Point & point(const unsigned int i) const
Definition: elem.h:2453
libMesh::TypeVector::relative_fuzzy_equals
bool relative_fuzzy_equals(const TypeVector< T > &rhs, Real tol=TOLERANCE) const
Definition: type_vector.h:981
libMesh::dof_id_type
uint8_t dof_id_type
Definition: id_types.h:67
Generated on Thu Jul 17 2025 01:29:08 for libMesh by   doxygen 1.8.14