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