libMesh
Loading...
Searching...
No Matches
cell_prism21.C
Go to the documentation of this file.
1// The libMesh Finite Element Library.
2// Copyright (C) 2002-2026 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
29namespace {
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
43namespace libMesh
44{
45
46
47
48// ------------------------------------------------------------
49// Prism21 class static member initializations
50const int Prism21::num_nodes;
51const int Prism21::nodes_per_side;
52const int Prism21::nodes_per_edge;
53
54const 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
63const 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
79bool Prism21::is_vertex(const unsigned int i) const
80{
81 if (i < 6)
82 return true;
83 return false;
84}
85
86bool 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
95bool 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
104bool 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
113std::vector<unsigned>
114Prism21::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
121std::vector<unsigned>
122Prism21::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
128bool 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
139bool 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
185Order Prism21::default_order() const
186{
187 return THIRD;
188}
189
190dof_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
223unsigned 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
239unsigned 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
250std::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
287void 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
330std::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
337void 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
344void 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
515unsigned 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
551unsigned 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
606const 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
617std::pair<unsigned short int, unsigned short int>
618Prism21::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
634const 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
838const 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
878void
879Prism21::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
918void
919Prism21::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
937unsigned 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
948ElemType
949Prism21::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::BoundaryInfo
The BoundaryInfo class contains information relevant to boundary conditions including storing faces,...
Definition boundary_info.h:58
libMesh::Elem::affine_tol
static constexpr Real affine_tol
Default tolerance to use in has_affine_map().
Definition elem.h:2070
libMesh::Elem::point
const Point & point(const unsigned int i) const
Definition elem.h:2462
libMesh::Elem::swap2nodes
void swap2nodes(unsigned int n1, unsigned int n2)
Swaps two node_ptrs.
Definition elem.h:2101
libMesh::Elem::swap3neighbors
void swap3neighbors(unsigned int n1, unsigned int n2, unsigned int n3)
Swaps three neighbor_ptrs, "rotating" them.
Definition elem.h:2142
libMesh::Elem::swap2neighbors
void swap2neighbors(unsigned int n1, unsigned int n2)
Swaps two neighbor_ptrs.
Definition elem.h:2111
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:3579
libMesh::Elem::compute_key
static dof_id_type compute_key(dof_id_type n0)
Definition elem.h:3311
libMesh::Elem::swap3nodes
void swap3nodes(unsigned int n1, unsigned int n2, unsigned int n3)
Swaps three node_ptrs, "rotating" them.
Definition elem.h:2133
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:3595
libMesh::Elem::node_id
dof_id_type node_id(const unsigned int i) const
Definition elem.h:2484
libMesh::Point
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition point.h:40
libMesh::Prism21::nodes_per_edge
static const int nodes_per_edge
Definition cell_prism21.h:245
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::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::Prism21::nodes_per_side
static const int nodes_per_side
Definition cell_prism21.h:244
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::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::nodes_on_side
virtual std::vector< unsigned int > nodes_on_side(const unsigned int s) const override
Definition cell_prism21.C:114
libMesh::Prism21::is_edge
virtual bool is_edge(const unsigned int i) const override
Definition cell_prism21.C:86
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::default_order
virtual Order default_order() const override
Definition cell_prism21.C:185
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_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::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::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::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::Prism21::is_vertex
virtual bool is_vertex(const unsigned int i) const override
Definition cell_prism21.C:79
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::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::Prism21::key
virtual dof_id_type key() const
Don't hide Elem::key() defined in the base class.
Definition elem.C:738
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::Prism21::nodes_on_edge
virtual std::vector< unsigned int > nodes_on_edge(const unsigned int e) const override
Definition cell_prism21.C:122
libMesh::Prism21::n_nodes
virtual unsigned int n_nodes() const override
Definition cell_prism21.h:108
libMesh::Prism21::is_face
virtual bool is_face(const unsigned int i) const override
Definition cell_prism21.C:95
libMesh::Prism21::has_affine_map
virtual bool has_affine_map() const override
Definition cell_prism21.C:139
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::side_type
ElemType side_type(const unsigned int s) const override final
Definition cell_prism21.C:949
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::Prism21::num_nodes
static const int num_nodes
Geometric constants for Prism21.
Definition cell_prism21.h:243
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::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::Prism::n_vertices
virtual unsigned int n_vertices() const override final
Definition cell_prism.h:91
libMesh::Prism::n_edges
virtual unsigned int n_edges() const override final
Definition cell_prism.h:96
libMesh::Prism::num_children
static const int num_children
Definition cell_prism.h:75
libMesh::Prism::n_sides
virtual unsigned int n_sides() const override final
Definition cell_prism.h:86
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::Prism::num_edges
static const int num_edges
Definition cell_prism.h:74
libMesh::Prism::num_sides
static const int num_sides
Geometric constants for all Prisms.
Definition cell_prism.h:73
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::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::TypeVector::relative_fuzzy_equals
bool relative_fuzzy_equals(const TypeVector< T > &rhs, Real tol=TOLERANCE) const
Definition type_vector.h:984
libMesh
The libMesh namespace provides an interface to certain functionality in the library.
libMesh::IOPackage
IOPackage
libMesh interfaces with several different software packages for the purposes of creating,...
Definition enum_io_package.h:38
libMesh::ElemType
ElemType
Defines an enum for geometric element types.
Definition enum_elem_type.h:33
libMesh::libmesh_assert
libmesh_assert(ctx)
libMesh::dof_id_type
uint8_t dof_id_type
Definition id_types.h:67
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition libmesh_common.h:144
Generated on Fri Aug 21 2026 13:32:29 for libMesh by doxygen 1.9.8