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
54// ------------------------------------------------------------
55// Prism21 class member functions
56
57bool Prism21::is_vertex(const unsigned int i) const
58{
59 if (i < 6)
60 return true;
61 return false;
62}
63
64bool Prism21::is_edge(const unsigned int i) const
65{
66 if (i < 6)
67 return false;
68 if (i > 14)
69 return false;
70 return true;
71}
72
73bool Prism21::is_face(const unsigned int i) const
74{
75 if (i > 19)
76 return false;
77 if (i > 14)
78 return true;
79 return false;
80}
81
82bool Prism21::is_node_on_side(const unsigned int n,
83 const unsigned int s) const
84{
85 libmesh_assert_less (s, n_sides());
86 return std::find(std::begin(side_nodes_map[s]),
87 std::end(side_nodes_map[s]),
88 n) != std::end(side_nodes_map[s]);
89}
90
91std::vector<unsigned>
92Prism21::nodes_on_side(const unsigned int s) const
93{
94 libmesh_assert_less(s, n_sides());
95 auto trim = (s > 0 && s < 4) ? 0 : 2;
96 return {std::begin(side_nodes_map[s]), std::end(side_nodes_map[s]) - trim};
97}
98
99std::vector<unsigned>
100Prism21::nodes_on_edge(const unsigned int e) const
101{
102 libmesh_assert_less(e, n_edges());
103 return {std::begin(edge_nodes_map[e]), std::end(edge_nodes_map[e])};
104}
105
106bool Prism21::is_node_on_edge(const unsigned int n,
107 const unsigned int e) const
108{
109 libmesh_assert_less (e, n_edges());
110 return std::find(std::begin(edge_nodes_map[e]),
111 std::end(edge_nodes_map[e]),
112 n) != std::end(edge_nodes_map[e]);
113}
114
115
116
117bool Prism21::has_affine_map() const
118{
119 // Make sure z edges are affine
120 Point v = this->point(3) - this->point(0);
121 if (!v.relative_fuzzy_equals(this->point(4) - this->point(1), affine_tol) ||
122 !v.relative_fuzzy_equals(this->point(5) - this->point(2), affine_tol))
123 return false;
124
125 // Make sure edges are straight
126 v /= 2;
127 if (!v.relative_fuzzy_equals(this->point(9) - this->point(0), affine_tol) ||
128 !v.relative_fuzzy_equals(this->point(10) - this->point(1), affine_tol) ||
129 !v.relative_fuzzy_equals(this->point(11) - this->point(2), affine_tol) ||
130 !v.relative_fuzzy_equals(this->point(15) - this->point(6), affine_tol) ||
131 !v.relative_fuzzy_equals(this->point(16) - this->point(7), affine_tol) ||
132 !v.relative_fuzzy_equals(this->point(17) - this->point(8), affine_tol))
133 return false;
134 v = (this->point(1) - this->point(0))/2;
135 if (!v.relative_fuzzy_equals(this->point(6) - this->point(0), affine_tol) ||
136 !v.relative_fuzzy_equals(this->point(12) - this->point(3), affine_tol))
137 return false;
138 v = (this->point(2) - this->point(0))/2;
139 if (!v.relative_fuzzy_equals(this->point(8) - this->point(0), affine_tol) ||
140 !v.relative_fuzzy_equals(this->point(14) - this->point(3), affine_tol))
141 return false;
142 v = (this->point(2) - this->point(1))/2;
143 if (!v.relative_fuzzy_equals(this->point(7) - this->point(1), affine_tol) ||
144 !v.relative_fuzzy_equals(this->point(13) - this->point(4), affine_tol))
145 return false;
146
147 // Make sure triangle face midpoints are centered
148 v = this->point(2) + this->point(1) - 2*this->point(0);
149 if (!v.relative_fuzzy_equals((this->point(18)-this->point(0))*3))
150 return false;
151 v = this->point(5) + this->point(4) - 2*this->point(3);
152 if (!v.relative_fuzzy_equals((this->point(19)-this->point(3))*3))
153 return false;
154 v = this->point(11) + this->point(10) - 2*this->point(9);
155 if (!v.relative_fuzzy_equals((this->point(20)-this->point(9))*3))
156 return false;
157
158 return true;
159}
160
161
162
163Order Prism21::default_order() const
164{
165 return THIRD;
166}
167
168dof_id_type Prism21::key (const unsigned int s) const
169{
170 libmesh_assert_less (s, this->n_sides());
171
172 switch (s)
173 {
174 case 0: // the triangular face at z=0
175 {
176 return Elem::compute_key (this->node_id(18));
177 }
178 case 1: // the quad face at y=0
179 {
180 return Elem::compute_key (this->node_id(15));
181 }
182 case 2: // the other quad face
183 {
184 return Elem::compute_key (this->node_id(16));
185 }
186 case 3: // the quad face at x=0
187 {
188 return Elem::compute_key (this->node_id(17));
189 }
190 case 4: // the triangular face at z=1
191 {
192 return Elem::compute_key (this->node_id(19));
193 }
194 default:
195 libmesh_error_msg("Invalid side " << s);
196 }
197}
198
199
200
201unsigned int Prism21::local_side_node(unsigned int side,
202 unsigned int side_node) const
203{
204 libmesh_assert_less (side, this->n_sides());
205
206 // Never more than 9 nodes per side.
207 libmesh_assert_less(side_node, Prism21::nodes_per_side);
208
209 // Some sides have 7 nodes.
210 libmesh_assert(!(side==0 || side==4) || side_node < 7);
211
212 return Prism21::side_nodes_map[side][side_node];
213}
214
215
216
217unsigned int Prism21::local_edge_node(unsigned int edge,
218 unsigned int edge_node) const
219{
220 libmesh_assert_less(edge, this->n_edges());
221 libmesh_assert_less(edge_node, Prism21::nodes_per_edge);
222
223 return Prism21::edge_nodes_map[edge][edge_node];
224}
225
226
227
228std::unique_ptr<Elem> Prism21::build_side_ptr (const unsigned int i)
229{
230 libmesh_assert_less (i, this->n_sides());
231
232 std::unique_ptr<Elem> face;
233
234 switch (i)
235 {
236 case 0: // the triangular face at z=-1
237 case 4: // the triangular face at z=1
238 {
239 face = std::make_unique<Tri7>();
240 break;
241 }
242 case 1: // the quad face at y=0
243 case 2: // the other quad face
244 case 3: // the quad face at x=0
245 {
246 face = std::make_unique<Quad9>();
247 break;
248 }
249 default:
250 libmesh_error_msg("Invalid side i = " << i);
251 }
252
253 // Set the nodes
254 for (auto n : face->node_index_range())
255 face->set_node(n, this->node_ptr(Prism21::side_nodes_map[i][n]));
256
257 face->set_interior_parent(this);
258 face->inherit_data_from(*this);
259
260 return face;
261}
262
263
264
265void Prism21::build_side_ptr (std::unique_ptr<Elem> & side,
266 const unsigned int i)
267{
268 libmesh_assert_less (i, this->n_sides());
269
270 switch (i)
271 {
272 case 0: // the triangular face at z=-1
273 case 4: // the triangular face at z=1
274 {
275 if (!side.get() || side->type() != TRI7)
276 {
277 side = this->build_side_ptr(i);
278 return;
279 }
280 break;
281 }
282
283 case 1: // the quad face at y=0
284 case 2: // the other quad face
285 case 3: // the quad face at x=0
286 {
287 if (!side.get() || side->type() != QUAD9)
288 {
289 side = this->build_side_ptr(i);
290 return;
291 }
292 break;
293 }
294
295 default:
296 libmesh_error_msg("Invalid side i = " << i);
297 }
298
299 side->inherit_data_from(*this);
300
301 // Set the nodes
302 for (auto n : side->node_index_range())
303 side->set_node(n, this->node_ptr(Prism21::side_nodes_map[i][n]));
304}
305
306
307
308std::unique_ptr<Elem> Prism21::build_edge_ptr (const unsigned int i)
309{
310 return this->simple_build_edge_ptr<Edge3,Prism21>(i);
311}
312
313
314
315void Prism21::build_edge_ptr (std::unique_ptr<Elem> & edge, const unsigned int i)
316{
317 this->simple_build_edge_ptr<Prism21>(edge, i, EDGE3);
318}
319
320
321
322void Prism21::connectivity(const unsigned int /*sc*/,
323 const IOPackage /*iop*/,
324 std::vector<dof_id_type> & /*conn*/) const
325{
326 libmesh_not_implemented(); // FIXME RHS
327
328/*
329 libmesh_assert(_nodes);
330 libmesh_assert_less (sc, this->n_sub_elem());
331 libmesh_assert_not_equal_to (iop, INVALID_IO_PACKAGE);
332
333 switch (iop)
334 {
335 case TECPLOT:
336 {
337 conn.resize(8);
338 switch (sc)
339 {
340
341 case 0:
342 {
343 conn[0] = this->node_id(0)+1;
344 conn[1] = this->node_id(6)+1;
345 conn[2] = this->node_id(8)+1;
346 conn[3] = this->node_id(8)+1;
347 conn[4] = this->node_id(9)+1;
348 conn[5] = this->node_id(15)+1;
349 conn[6] = this->node_id(17)+1;
350 conn[7] = this->node_id(17)+1;
351
352 return;
353 }
354
355 case 1:
356 {
357 conn[0] = this->node_id(6)+1;
358 conn[1] = this->node_id(1)+1;
359 conn[2] = this->node_id(7)+1;
360 conn[3] = this->node_id(7)+1;
361 conn[4] = this->node_id(15)+1;
362 conn[5] = this->node_id(10)+1;
363 conn[6] = this->node_id(16)+1;
364 conn[7] = this->node_id(16)+1;
365
366 return;
367 }
368
369 case 2:
370 {
371 conn[0] = this->node_id(8)+1;
372 conn[1] = this->node_id(7)+1;
373 conn[2] = this->node_id(2)+1;
374 conn[3] = this->node_id(2)+1;
375 conn[4] = this->node_id(17)+1;
376 conn[5] = this->node_id(16)+1;
377 conn[6] = this->node_id(11)+1;
378 conn[7] = this->node_id(11)+1;
379
380 return;
381 }
382
383 case 3:
384 {
385 conn[0] = this->node_id(6)+1;
386 conn[1] = this->node_id(7)+1;
387 conn[2] = this->node_id(8)+1;
388 conn[3] = this->node_id(8)+1;
389 conn[4] = this->node_id(15)+1;
390 conn[5] = this->node_id(16)+1;
391 conn[6] = this->node_id(17)+1;
392 conn[7] = this->node_id(17)+1;
393
394 return;
395 }
396
397 case 4:
398 {
399 conn[0] = this->node_id(9)+1;
400 conn[1] = this->node_id(15)+1;
401 conn[2] = this->node_id(17)+1;
402 conn[3] = this->node_id(17)+1;
403 conn[4] = this->node_id(3)+1;
404 conn[5] = this->node_id(12)+1;
405 conn[6] = this->node_id(14)+1;
406 conn[7] = this->node_id(14)+1;
407
408 return;
409 }
410
411 case 5:
412 {
413 conn[0] = this->node_id(15)+1;
414 conn[1] = this->node_id(10)+1;
415 conn[2] = this->node_id(16)+1;
416 conn[3] = this->node_id(16)+1;
417 conn[4] = this->node_id(12)+1;
418 conn[5] = this->node_id(4)+1;
419 conn[6] = this->node_id(13)+1;
420 conn[7] = this->node_id(13)+1;
421
422 return;
423 }
424
425 case 6:
426 {
427 conn[0] = this->node_id(17)+1;
428 conn[1] = this->node_id(16)+1;
429 conn[2] = this->node_id(11)+1;
430 conn[3] = this->node_id(11)+1;
431 conn[4] = this->node_id(14)+1;
432 conn[5] = this->node_id(13)+1;
433 conn[6] = this->node_id(5)+1;
434 conn[7] = this->node_id(5)+1;
435
436 return;
437 }
438
439 case 7:
440 {
441 conn[0] = this->node_id(15)+1;
442 conn[1] = this->node_id(16)+1;
443 conn[2] = this->node_id(17)+1;
444 conn[3] = this->node_id(17)+1;
445 conn[4] = this->node_id(12)+1;
446 conn[5] = this->node_id(13)+1;
447 conn[6] = this->node_id(14)+1;
448 conn[7] = this->node_id(14)+1;
449
450 return;
451 }
452
453 default:
454 libmesh_error_msg("Invalid sc = " << sc);
455 }
456
457 }
458
459 case VTK:
460 {
461 // VTK now supports VTK_BIQUADRATIC_QUADRATIC_WEDGE directly
462 const unsigned int conn_size = 18;
463 conn.resize(conn_size);
464
465 // VTK's VTK_BIQUADRATIC_QUADRATIC_WEDGE first 9 (vertex) and
466 // last 3 (mid-face) nodes match. The middle and top layers
467 // of mid-edge nodes are reversed from LibMesh's.
468 for (auto i : index_range(conn))
469 conn[i] = this->node_id(i);
470
471 // top "ring" of mid-edge nodes
472 conn[9] = this->node_id(12);
473 conn[10] = this->node_id(13);
474 conn[11] = this->node_id(14);
475
476 // middle "ring" of mid-edge nodes
477 conn[12] = this->node_id(9);
478 conn[13] = this->node_id(10);
479 conn[14] = this->node_id(11);
480
481 return;
482 }
483
484 default:
485 libmesh_error_msg("Unsupported IO package " << iop);
486 }
487*/
488}
489
490
491
492
493unsigned int Prism21::n_second_order_adjacent_vertices (const unsigned int n) const
494{
495 switch (n)
496 {
497 case 6:
498 case 7:
499 case 8:
500 case 9:
501 case 10:
502 case 11:
503 case 12:
504 case 13:
505 case 14:
506 return 2;
507
508 case 15:
509 case 16:
510 case 17:
511 return 4;
512
513 case 18:
514 case 19:
515 return 3;
516
517 case 20:
518 return 6;
519
520 default:
521 libmesh_error_msg("Invalid node n = " << n);
522 }
523}
524
525
526
527
528
529unsigned short int Prism21::second_order_adjacent_vertex (const unsigned int n,
530 const unsigned int v) const
531{
532 libmesh_assert_greater_equal (n, this->n_vertices());
533 libmesh_assert_less (n, this->n_nodes());
534
535 switch (n)
536 {
537 /*
538 * These nodes are unique to \p Prism20,
539 * let our _remaining_... matrix handle
540 * this.
541 */
542 case 15:
543 case 16:
544 case 17:
545 case 18:
546 case 19:
547 {
548 libmesh_assert_less (v, 4);
549 return _remaining_second_order_adjacent_vertices[n-15][v];
550 }
551
552 /*
553 * for the bubble node the return value is simply v.
554 * Why? -- the user asks for the v-th adjacent vertex,
555 * from \p n_second_order_adjacent_vertices() there
556 * are 6 adjacent vertices, and these happen to be
557 * 0..5
558 */
559 case 20:
560 {
561 libmesh_assert_less (v, 6);
562 return static_cast<unsigned short int>(v);
563 }
564
565 /*
566 * All other second-order nodes (6,...,14) are
567 * identical with Prism15 and are therefore
568 * delegated to the _second_order matrix of
569 * \p Prism
570 */
571 default:
572 {
573 libmesh_assert_less (v, 2);
574 return _second_order_adjacent_vertices[n-this->n_vertices()][v];
575 }
576
577 }
578
579 return static_cast<unsigned short int>(-1);
580}
581
582
583
584const unsigned short int Prism21::_remaining_second_order_adjacent_vertices[5][4] =
585 {
586 { 0, 1, 3, 4}, // vertices adjacent to node 15
587 { 1, 2, 4, 5}, // vertices adjacent to node 16
588 { 0, 2, 3, 5}, // vertices adjacent to node 17
589 { 0, 1, 2, 99}, // vertices adjacent to node 18
590 { 3, 4, 5, 99} // vertices adjacent to node 19
591 };
592
593
594
595std::pair<unsigned short int, unsigned short int>
596Prism21::second_order_child_vertex (const unsigned int n) const
597{
598 libmesh_assert_greater_equal (n, this->n_vertices());
599 libmesh_assert_less (n, this->n_nodes());
600
601 return std::pair<unsigned short int, unsigned short int>
602 (_second_order_vertex_child_number[n],
603 _second_order_vertex_child_index[n]);
604}
605
606
607
608#ifdef LIBMESH_ENABLE_AMR
609
610// FIXME RHS
611
612const Real Prism21::_embedding_matrix[Prism21::num_children][Prism21::num_nodes][Prism21::num_nodes] =
613 {
614 // embedding matrix for child 0
615 {
616 // 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20
617 { 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // 0
618 { 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // 1
619 { 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // 2
620 { 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // 3
621 { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0}, // 4
622 { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0}, // 5
623 { 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
624 { 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
625 { 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
626 { 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
627 { 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
628 { 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
629 { 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
630 { 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
631 { 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
632 { 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
633 { 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
634 { 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
635 { 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
636 { 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
637 { 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
638 },
639 // embedding matrix for child 1
640 {
641 // 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20
642 { 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // 0
643 { 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // 1
644 { 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // 2
645 { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0}, // 3
646 { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // 4
647 { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0}, // 5
648 { -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
649 { 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
650 { -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
651 { 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
652 { 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
653 { 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
654 { 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
655 { 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
656 { 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
657 { -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
658 { 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
659 { -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
660 { -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
661 { 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
662 { -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
663 },
664 // embedding matrix for child 2
665 {
666 // 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20
667 { 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // 0
668 { 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // 1
669 { 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // 2
670 { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0}, // 3
671 { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0}, // 4
672 { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // 5
673 { -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
674 { 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
675 { -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
676 { 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
677 { 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
678 { 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
679 { 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
680 { 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
681 { 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
682 { -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
683 { 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
684 { -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
685 { -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
686 { 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
687 { -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
688 },
689 // embedding matrix for child 3
690 {
691 // 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20
692 { 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // 0
693 { 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // 1
694 { 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // 2
695 { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0}, // 3
696 { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0}, // 4
697 { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0}, // 5
698 { -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
699 { -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
700 { 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
701 { 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
702 { 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
703 { 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
704 { 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
705 { 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
706 { 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
707 { -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
708 { -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
709 { 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
710 { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0}, // 18
711 { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1}, // 19
712 { 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
713 },
714 // embedding matrix for child 4
715 {
716 // 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20
717 { 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // 0
718 { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0}, // 1
719 { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0}, // 2
720 { 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // 3
721 { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0}, // 4
722 { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0}, // 5
723 { 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
724 { 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
725 { 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
726 { -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
727 { 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
728 { 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
729 { 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
730 { 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
731 { 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
732 { -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
733 { -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
734 { -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
735 { 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
736 { 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
737 { -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
738 },
739 // embedding matrix for child 5
740 {
741 // 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20
742 { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0}, // 0
743 { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // 1
744 { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0}, // 2
745 { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0}, // 3
746 { 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // 4
747 { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0}, // 5
748 { 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
749 { 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
750 { 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
751 { 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
752 { 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
753 { 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
754 { 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
755 { 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
756 { 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
757 { 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
758 { 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
759 { 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
760 { 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
761 { 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
762 { 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
763 },
764 // embedding matrix for child 6
765 {
766 // 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20
767 { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0}, // 0
768 { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0}, // 1
769 { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // 2
770 { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0}, // 3
771 { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0}, // 4
772 { 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // 5
773 { 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
774 { 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
775 { 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
776 { 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
777 { 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
778 { 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
779 { 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
780 { 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
781 { 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
782 { 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
783 { 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
784 { 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
785 { 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
786 { 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
787 { 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
788 },
789 // embedding matrix for child 7
790 {
791 // 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20
792 { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0}, // 0
793 { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0}, // 1
794 { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0}, // 2
795 { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0}, // 3
796 { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0}, // 4
797 { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0}, // 5
798 { 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
799 { 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
800 { 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
801 { 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
802 { 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
803 { 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
804 { 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
805 { 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
806 { 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
807 { 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
808 { 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
809 { -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
810 { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1}, // 18
811 { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0}, // 19
812 { 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
813 }
814 };
815
816const std::vector<std::pair<unsigned char, unsigned char>>
817 Prism21::_parent_bracketing_nodes[Prism21::num_children][Prism21::num_nodes] =
818 {
819 // Child 0
820 { {},{{0,1}},{{0,2}},{{0,3}},{{0,4},{1,3},{6,12},{9,10}},{{0,5},{2,3},{8,14},{9,11}},
821 {{0,6}},{{6,8}},{{0,8}},{{0,9}},{{6,15}},{{8,17}},{{9,15}},{{15,17}},{{9,17}},
822 {{0,15},{6,9}},{{6,17},{8,15}},{{0,17},{8,9}},{{0,18}},{{9,20}},{{0,20},{9,18}} },
823 // Child 1
824 { {{0,1}}, {},{{1,2}},{{0,4},{1,3},{6,12},{9,10}},{{1,4}},{{1,5},{2,4},{7,13},{10,11}},
825 {{1,6}},{{1,7}},{{6,7}},{{6,15}},{{1,10}},{{7,16}},{{10,15}},{{10,16}},{{15,16}},
826 {{1,15},{6,10}},{{1,16},{7,10}},{{6,16},{7,15}},{{1,18}},{{10,20}},{{1,20},{10,18}} },
827 // Child 2
828 { {{0,2}},{{1,2}}, {},{{0,5},{2,3},{8,14},{9,11}},{{1,5},{2,4},{7,13},{10,11}},{{2,5}},
829 {{7,8}},{{2,7}},{{2,8}},{{8,17}},{{7,16}},{{2,11}},{{16,17}},{{11,16}},{{11,17}},
830 {{7,17},{8,16}},{{2,16},{7,11}},{{2,17},{8,11}},{{2,18}},{{11,20}},{{2,20},{11,18}} },
831 // Child 3
832 { {{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}},
833 {{6,7}},{{7,8}},{{6,8}},{{6,15}},{{7,16}},{{8,17}},{{15,16}},{{16,17}},{{15,17}},
834 {{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}} },
835 // Child 4
836 { {{0,3}},{{0,4},{1,3},{6,11},{9,10}},{{0,5},{2,3},{8,14},{9,11}}, {},{{3,4}},{{3,5}},
837 {{9,15}},{{15,17}},{{9,17}},{{3,9}},{{12,15}},{{14,17}},{{3,12}},{{12,14}},{{3,14}},
838 {{3,15},{9,12}},{{12,17},{14,15}},{{3,17},{9,14}},{{9,20}},{{3,19}},{{9,19},{3,20}} },
839 // Child 5
840 { {{0,4},{1,3},{6,12},{9,10}},{{1,4}},{{1,5},{2,4},{7,13},{10,11}},{{3,4}}, {},{{4,5}},
841 {{10,15}},{{10,16}},{{15,16}},{{12,15}},{{4,10}},{{13,16}},{{4,12}},{{4,13}},{{12,13}},
842 {{4,15},{10,12}},{{4,16},{10,13}},{{12,16},{13,15}},{{10,20}},{{4,19}},{{10,19},{4,20}} },
843 // Child 6
844 { {{0,5},{2,3},{8,14},{9,11}},{{1,5},{2,4},{7,13},{10,11}},{{2,5}},{{3,5}},{{4,5}}, {},
845 {{16,17}},{{11,16}},{{11,17}},{{14,17}},{{13,16}},{{5,11}},{{13,14}},{{5,13}},{{5,14}},
846 {{13,17},{14,16}},{{5,16},{11,13}},{{5,17},{11,14}},{{11,20}},{{5,19}},{{11,19},{5,20}} },
847 // Child 7
848 { {{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}},
849 {{15,16}},{{16,17}},{{15,17}},{{12,15}},{{13,16}},{{14,17}},{{12,13}},{{13,14}},{{12,14}},
850 {{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}} }
851 };
852
853#endif
854
855
856void
857Prism21::permute(unsigned int perm_num)
858{
859 libmesh_assert_less (perm_num, 6);
860 const unsigned int side = perm_num % 2;
861 const unsigned int rotate = perm_num / 2;
862
863 for (unsigned int i = 0; i != rotate; ++i)
864 {
865 swap3nodes(0,1,2);
866 swap3nodes(3,4,5);
867 swap3nodes(6,7,8);
868 swap3nodes(9,10,11);
869 swap3nodes(12,13,14);
870 swap3nodes(15,16,17);
871 swap3neighbors(1,2,3);
872 }
873
874 switch (side) {
875 case 0:
876 break;
877 case 1:
878 swap2nodes(1,3);
879 swap2nodes(0,4);
880 swap2nodes(2,5);
881 swap2nodes(6,12);
882 swap2nodes(9,10);
883 swap2nodes(7,14);
884 swap2nodes(8,13);
885 swap2nodes(16,17);
886 swap2nodes(18,19);
887 swap2neighbors(0,4);
888 swap2neighbors(2,3);
889 break;
890 default:
891 libmesh_error();
892 }
893}
894
895
896void
897Prism21::flip(BoundaryInfo * boundary_info)
898{
899 libmesh_assert(boundary_info);
900
901 swap2nodes(0,1);
902 swap2nodes(3,4);
903 swap2nodes(7,8);
904 swap2nodes(9,10);
905 swap2nodes(13,14);
906 swap2nodes(16,17);
907 swap2neighbors(2,3);
908 swap2boundarysides(2,3,boundary_info);
909 swap2boundaryedges(0,1,boundary_info);
910 swap2boundaryedges(3,4,boundary_info);
911 swap2boundaryedges(7,8,boundary_info);
912}
913
914
915unsigned int Prism21::center_node_on_side(const unsigned short side) const
916{
917 libmesh_assert_less (side, Prism21::num_sides);
918 if (side >= 1 && side <= 3)
919 return side + 14;
920 if (side == 4)
921 return 19;
922 return 18;
923}
924
925
926ElemType
927Prism21::side_type (const unsigned int s) const
928{
929 libmesh_assert_less (s, 5);
930 if (s == 0 || s == 4)
931 return TRI7;
932 return QUAD9;
933}
934
935
936} // 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:2191
libMesh::Elem::point
const Point & point(const unsigned int i) const
Definition elem.h:2583
libMesh::Elem::swap2nodes
void swap2nodes(unsigned int n1, unsigned int n2)
Swaps two node_ptrs.
Definition elem.h:2222
libMesh::Elem::swap3neighbors
void swap3neighbors(unsigned int n1, unsigned int n2, unsigned int n3)
Swaps three neighbor_ptrs, "rotating" them.
Definition elem.h:2263
libMesh::Elem::swap2neighbors
void swap2neighbors(unsigned int n1, unsigned int n2)
Swaps two neighbor_ptrs.
Definition elem.h:2232
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:3982
libMesh::Elem::compute_key
static dof_id_type compute_key(dof_id_type n0)
Definition elem.h:3432
libMesh::Elem::swap3nodes
void swap3nodes(unsigned int n1, unsigned int n2, unsigned int n3)
Swaps three node_ptrs, "rotating" them.
Definition elem.h:2254
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:3998
libMesh::Elem::node_id
dof_id_type node_id(const unsigned int i) const
Definition elem.h:2605
libMesh::Point
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition point.h:41
libMesh::Prism21::edge_nodes_map
static constexpr const unsigned int(& edge_nodes_map)[num_edges][nodes_per_edge]
Definition cell_prism21.h:256
libMesh::Prism21::nodes_per_edge
static const int nodes_per_edge
Definition cell_prism21.h:247
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:309
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:897
libMesh::Prism21::nodes_per_side
static const int nodes_per_side
Definition cell_prism21.h:246
libMesh::Prism21::n_second_order_adjacent_vertices
virtual unsigned int n_second_order_adjacent_vertices(const unsigned int) const override
Definition cell_prism21.C:493
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:82
libMesh::Prism21::nodes_on_side
virtual std::vector< unsigned int > nodes_on_side(const unsigned int s) const override
Definition cell_prism21.C:92
libMesh::Prism21::is_edge
virtual bool is_edge(const unsigned int i) const override
Definition cell_prism21.C:64
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:322
libMesh::Prism21::default_order
virtual Order default_order() const override
Definition cell_prism21.C:163
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:529
libMesh::Prism21::local_side_node
virtual unsigned int local_side_node(unsigned int side, unsigned int side_node) const override
Definition cell_prism21.C:201
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:308
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:596
libMesh::Prism21::is_vertex
virtual bool is_vertex(const unsigned int i) const override
Definition cell_prism21.C:57
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:857
libMesh::Prism21::center_node_on_side
unsigned int center_node_on_side(const unsigned short side) const override final
Definition cell_prism21.C:915
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:324
libMesh::Prism21::key
virtual dof_id_type key() const
Don't hide Elem::key() defined in the base class.
Definition elem.C:678
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:302
libMesh::Prism21::nodes_on_edge
virtual std::vector< unsigned int > nodes_on_edge(const unsigned int e) const override
Definition cell_prism21.C:100
libMesh::Prism21::n_nodes
virtual unsigned int n_nodes() const override
Definition cell_prism21.h:110
libMesh::Prism21::is_face
virtual bool is_face(const unsigned int i) const override
Definition cell_prism21.C:73
libMesh::Prism21::has_affine_map
virtual bool has_affine_map() const override
Definition cell_prism21.C:117
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:106
libMesh::Prism21::side_type
ElemType side_type(const unsigned int s) const override final
Definition cell_prism21.C:927
libMesh::Prism21::local_edge_node
virtual unsigned int local_edge_node(unsigned int edge, unsigned int edge_node) const override
Definition cell_prism21.C:217
libMesh::Prism21::side_nodes_map
static constexpr const unsigned int(& side_nodes_map)[num_sides][nodes_per_side]
Definition cell_prism21.h:261
libMesh::Prism21::num_nodes
static const int num_nodes
Geometric constants for Prism21.
Definition cell_prism21.h:245
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:228
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:213
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:208
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:218
libMesh::TypeVector::relative_fuzzy_equals
bool relative_fuzzy_equals(const TypeVector< T > &rhs, Real tol=TOLERANCE) const
Definition type_vector.h:1052
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:148
Generated on Wed Nov 11 2026 20:15:53 for libMesh by doxygen 1.9.8