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