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/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 : // ------------------------------------------------------------
55 : // Prism21 class member functions
56 :
57 4958012 : bool Prism21::is_vertex(const unsigned int i) const
58 : {
59 4958012 : if (i < 6)
60 1430258 : return true;
61 242730 : return false;
62 : }
63 :
64 135664 : bool Prism21::is_edge(const unsigned int i) const
65 : {
66 135664 : if (i < 6)
67 0 : return false;
68 135664 : if (i > 14)
69 54448 : return false;
70 6768 : return true;
71 : }
72 :
73 2361331 : bool Prism21::is_face(const unsigned int i) const
74 : {
75 2361331 : if (i > 19)
76 1248 : return false;
77 2353043 : if (i > 14)
78 400583 : return true;
79 158520 : return false;
80 : }
81 :
82 1032459 : bool Prism21::is_node_on_side(const unsigned int n,
83 : const unsigned int s) const
84 : {
85 81298 : libmesh_assert_less (s, n_sides());
86 81298 : return std::find(std::begin(side_nodes_map[s]),
87 1032459 : std::end(side_nodes_map[s]),
88 1032459 : n) != std::end(side_nodes_map[s]);
89 : }
90 :
91 : std::vector<unsigned>
92 264159259 : Prism21::nodes_on_side(const unsigned int s) const
93 : {
94 22851167 : libmesh_assert_less(s, n_sides());
95 264159259 : auto trim = (s > 0 && s < 4) ? 0 : 2;
96 264159259 : return {std::begin(side_nodes_map[s]), std::end(side_nodes_map[s]) - trim};
97 : }
98 :
99 : std::vector<unsigned>
100 38655 : Prism21::nodes_on_edge(const unsigned int e) const
101 : {
102 3186 : libmesh_assert_less(e, n_edges());
103 41841 : return {std::begin(edge_nodes_map[e]), std::end(edge_nodes_map[e])};
104 : }
105 :
106 561915 : bool Prism21::is_node_on_edge(const unsigned int n,
107 : const unsigned int e) const
108 : {
109 40668 : libmesh_assert_less (e, n_edges());
110 40668 : return std::find(std::begin(edge_nodes_map[e]),
111 561915 : std::end(edge_nodes_map[e]),
112 561915 : n) != std::end(edge_nodes_map[e]);
113 : }
114 :
115 :
116 :
117 799403 : bool Prism21::has_affine_map() const
118 : {
119 : // Make sure z edges are affine
120 137462 : Point v = this->point(3) - this->point(0);
121 936861 : if (!v.relative_fuzzy_equals(this->point(4) - this->point(1), affine_tol) ||
122 936792 : !v.relative_fuzzy_equals(this->point(5) - this->point(2), affine_tol))
123 71 : return false;
124 :
125 : // Make sure edges are straight
126 68731 : v /= 2;
127 867267 : if (!v.relative_fuzzy_equals(this->point(9) - this->point(0), affine_tol) ||
128 867123 : !v.relative_fuzzy_equals(this->point(10) - this->point(1), affine_tol) ||
129 867123 : !v.relative_fuzzy_equals(this->point(11) - this->point(2), affine_tol) ||
130 867123 : !v.relative_fuzzy_equals(this->point(15) - this->point(6), affine_tol) ||
131 1735114 : !v.relative_fuzzy_equals(this->point(16) - this->point(7), affine_tol) ||
132 867199 : !v.relative_fuzzy_equals(this->point(17) - this->point(8), affine_tol))
133 864 : return false;
134 867123 : v = (this->point(1) - this->point(0))/2;
135 935782 : if (!v.relative_fuzzy_equals(this->point(6) - this->point(0), affine_tol) ||
136 867127 : !v.relative_fuzzy_equals(this->point(12) - this->point(3), affine_tol))
137 0 : return false;
138 867123 : v = (this->point(2) - this->point(0))/2;
139 935782 : if (!v.relative_fuzzy_equals(this->point(8) - this->point(0), affine_tol) ||
140 867127 : !v.relative_fuzzy_equals(this->point(14) - this->point(3), affine_tol))
141 0 : return false;
142 867123 : v = (this->point(2) - this->point(1))/2;
143 935782 : if (!v.relative_fuzzy_equals(this->point(7) - this->point(1), affine_tol) ||
144 867127 : !v.relative_fuzzy_equals(this->point(13) - this->point(4), affine_tol))
145 0 : return false;
146 :
147 : // Make sure triangle face midpoints are centered
148 867123 : v = this->point(2) + this->point(1) - 2*this->point(0);
149 798468 : if (!v.relative_fuzzy_equals((this->point(18)-this->point(0))*3))
150 0 : return false;
151 867123 : v = this->point(5) + this->point(4) - 2*this->point(3);
152 798468 : if (!v.relative_fuzzy_equals((this->point(19)-this->point(3))*3))
153 0 : return false;
154 867123 : v = this->point(11) + this->point(10) - 2*this->point(9);
155 798468 : if (!v.relative_fuzzy_equals((this->point(20)-this->point(9))*3))
156 0 : return false;
157 :
158 68659 : return true;
159 : }
160 :
161 :
162 :
163 26941070 : Order Prism21::default_order() const
164 : {
165 26941070 : return THIRD;
166 : }
167 :
168 0 : dof_id_type Prism21::key (const unsigned int s) const
169 : {
170 0 : libmesh_assert_less (s, this->n_sides());
171 :
172 0 : switch (s)
173 : {
174 0 : case 0: // the triangular face at z=0
175 : {
176 0 : return Elem::compute_key (this->node_id(18));
177 : }
178 0 : case 1: // the quad face at y=0
179 : {
180 0 : return Elem::compute_key (this->node_id(15));
181 : }
182 0 : case 2: // the other quad face
183 : {
184 0 : return Elem::compute_key (this->node_id(16));
185 : }
186 0 : case 3: // the quad face at x=0
187 : {
188 0 : return Elem::compute_key (this->node_id(17));
189 : }
190 0 : case 4: // the triangular face at z=1
191 : {
192 0 : return Elem::compute_key (this->node_id(19));
193 : }
194 0 : default:
195 0 : libmesh_error_msg("Invalid side " << s);
196 : }
197 : }
198 :
199 :
200 :
201 1248699168 : unsigned int Prism21::local_side_node(unsigned int side,
202 : unsigned int side_node) const
203 : {
204 108711714 : libmesh_assert_less (side, this->n_sides());
205 :
206 : // Never more than 9 nodes per side.
207 108711714 : libmesh_assert_less(side_node, Prism21::nodes_per_side);
208 :
209 : // Some sides have 7 nodes.
210 108711714 : libmesh_assert(!(side==0 || side==4) || side_node < 7);
211 :
212 1248699168 : return Prism21::side_nodes_map[side][side_node];
213 : }
214 :
215 :
216 :
217 99010737 : unsigned int Prism21::local_edge_node(unsigned int edge,
218 : unsigned int edge_node) const
219 : {
220 8830542 : libmesh_assert_less(edge, this->n_edges());
221 8830542 : libmesh_assert_less(edge_node, Prism21::nodes_per_edge);
222 :
223 99010737 : return Prism21::edge_nodes_map[edge][edge_node];
224 : }
225 :
226 :
227 :
228 869460 : std::unique_ptr<Elem> Prism21::build_side_ptr (const unsigned int i)
229 : {
230 34319 : libmesh_assert_less (i, this->n_sides());
231 :
232 869460 : std::unique_ptr<Elem> face;
233 :
234 869460 : switch (i)
235 : {
236 505915 : case 0: // the triangular face at z=-1
237 : case 4: // the triangular face at z=1
238 : {
239 505915 : face = std::make_unique<Tri7>();
240 505915 : break;
241 : }
242 363545 : 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 363545 : face = std::make_unique<Quad9>();
247 363545 : break;
248 : }
249 0 : default:
250 0 : libmesh_error_msg("Invalid side i = " << i);
251 : }
252 :
253 : // Set the nodes
254 7682770 : for (auto n : face->node_index_range())
255 7088389 : face->set_node(n, this->node_ptr(Prism21::side_nodes_map[i][n]));
256 :
257 869460 : face->set_interior_parent(this);
258 835141 : face->inherit_data_from(*this);
259 :
260 869460 : return face;
261 0 : }
262 :
263 :
264 :
265 823607 : void Prism21::build_side_ptr (std::unique_ptr<Elem> & side,
266 : const unsigned int i)
267 : {
268 19230 : libmesh_assert_less (i, this->n_sides());
269 :
270 823607 : switch (i)
271 : {
272 4745 : case 0: // the triangular face at z=-1
273 : case 4: // the triangular face at z=1
274 : {
275 205801 : if (!side.get() || side->type() != TRI7)
276 : {
277 400702 : side = this->build_side_ptr(i);
278 205032 : return;
279 : }
280 64 : break;
281 : }
282 :
283 14485 : 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 617806 : if (!side.get() || side->type() != QUAD9)
288 : {
289 406410 : side = this->build_side_ptr(i);
290 208394 : return;
291 : }
292 9296 : break;
293 : }
294 :
295 0 : default:
296 0 : libmesh_error_msg("Invalid side i = " << i);
297 : }
298 :
299 400821 : side->inherit_data_from(*this);
300 :
301 : // Set the nodes
302 4100272 : for (auto n : side->node_index_range())
303 3774203 : side->set_node(n, this->node_ptr(Prism21::side_nodes_map[i][n]));
304 : }
305 :
306 :
307 :
308 869633 : std::unique_ptr<Elem> Prism21::build_edge_ptr (const unsigned int i)
309 : {
310 869633 : return this->simple_build_edge_ptr<Edge3,Prism21>(i);
311 : }
312 :
313 :
314 :
315 0 : void Prism21::build_edge_ptr (std::unique_ptr<Elem> & edge, const unsigned int i)
316 : {
317 0 : this->simple_build_edge_ptr<Prism21>(edge, i, EDGE3);
318 0 : }
319 :
320 :
321 :
322 0 : void Prism21::connectivity(const unsigned int /*sc*/,
323 : const IOPackage /*iop*/,
324 : std::vector<dof_id_type> & /*conn*/) const
325 : {
326 0 : 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 :
493 239355 : unsigned int Prism21::n_second_order_adjacent_vertices (const unsigned int n) const
494 : {
495 11490 : switch (n)
496 : {
497 6894 : 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 6894 : return 2;
507 :
508 2298 : case 15:
509 : case 16:
510 : case 17:
511 2298 : return 4;
512 :
513 1532 : case 18:
514 : case 19:
515 1532 : return 3;
516 :
517 766 : case 20:
518 766 : return 6;
519 :
520 0 : default:
521 0 : libmesh_error_msg("Invalid node n = " << n);
522 : }
523 : }
524 :
525 :
526 :
527 :
528 :
529 670194 : unsigned short int Prism21::second_order_adjacent_vertex (const unsigned int n,
530 : const unsigned int v) const
531 : {
532 32172 : libmesh_assert_greater_equal (n, this->n_vertices());
533 32172 : libmesh_assert_less (n, this->n_nodes());
534 :
535 670194 : switch (n)
536 : {
537 : /*
538 : * These nodes are unique to \p Prism20,
539 : * let our _remaining_... matrix handle
540 : * this.
541 : */
542 287226 : case 15:
543 : case 16:
544 : case 17:
545 : case 18:
546 : case 19:
547 : {
548 13788 : libmesh_assert_less (v, 4);
549 287226 : 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 95742 : case 20:
560 : {
561 4596 : libmesh_assert_less (v, 6);
562 95742 : 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 287226 : default:
572 : {
573 13788 : libmesh_assert_less (v, 2);
574 287226 : 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 :
584 : const 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 :
595 : std::pair<unsigned short int, unsigned short int>
596 0 : Prism21::second_order_child_vertex (const unsigned int n) const
597 : {
598 0 : libmesh_assert_greater_equal (n, this->n_vertices());
599 0 : libmesh_assert_less (n, this->n_nodes());
600 :
601 0 : return std::pair<unsigned short int, unsigned short int>
602 0 : (_second_order_vertex_child_number[n],
603 0 : _second_order_vertex_child_index[n]);
604 : }
605 :
606 :
607 :
608 : #ifdef LIBMESH_ENABLE_AMR
609 :
610 : // FIXME RHS
611 :
612 : const 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 :
816 : const 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 :
856 : void
857 49810 : Prism21::permute(unsigned int perm_num)
858 : {
859 2504 : libmesh_assert_less (perm_num, 6);
860 49810 : const unsigned int side = perm_num % 2;
861 49810 : const unsigned int rotate = perm_num / 2;
862 :
863 107802 : for (unsigned int i = 0; i != rotate; ++i)
864 : {
865 57992 : swap3nodes(0,1,2);
866 55026 : swap3nodes(3,4,5);
867 55026 : swap3nodes(6,7,8);
868 55026 : swap3nodes(9,10,11);
869 55026 : swap3nodes(12,13,14);
870 55026 : swap3nodes(15,16,17);
871 55026 : swap3neighbors(1,2,3);
872 : }
873 :
874 49810 : switch (side) {
875 890 : case 0:
876 890 : break;
877 29837 : case 1:
878 29837 : swap2nodes(1,3);
879 29837 : swap2nodes(0,4);
880 29837 : swap2nodes(2,5);
881 29837 : swap2nodes(6,12);
882 29837 : swap2nodes(9,10);
883 29837 : swap2nodes(7,14);
884 29837 : swap2nodes(8,13);
885 29837 : swap2nodes(16,17);
886 29837 : swap2nodes(18,19);
887 1614 : swap2neighbors(0,4);
888 1614 : swap2neighbors(2,3);
889 1614 : break;
890 0 : default:
891 0 : libmesh_error();
892 : }
893 49810 : }
894 :
895 :
896 : void
897 576 : Prism21::flip(BoundaryInfo * boundary_info)
898 : {
899 48 : libmesh_assert(boundary_info);
900 :
901 576 : swap2nodes(0,1);
902 576 : swap2nodes(3,4);
903 576 : swap2nodes(7,8);
904 576 : swap2nodes(9,10);
905 576 : swap2nodes(13,14);
906 576 : swap2nodes(16,17);
907 48 : swap2neighbors(2,3);
908 576 : swap2boundarysides(2,3,boundary_info);
909 576 : swap2boundaryedges(0,1,boundary_info);
910 576 : swap2boundaryedges(3,4,boundary_info);
911 576 : swap2boundaryedges(7,8,boundary_info);
912 576 : }
913 :
914 :
915 960 : unsigned int Prism21::center_node_on_side(const unsigned short side) const
916 : {
917 80 : libmesh_assert_less (side, Prism21::num_sides);
918 960 : if (side >= 1 && side <= 3)
919 576 : return side + 14;
920 384 : if (side == 4)
921 192 : return 19;
922 16 : return 18;
923 : }
924 :
925 :
926 : ElemType
927 528119886 : Prism21::side_type (const unsigned int s) const
928 : {
929 45685802 : libmesh_assert_less (s, 5);
930 528119886 : if (s == 0 || s == 4)
931 104644692 : return TRI7;
932 36690542 : return QUAD9;
933 : }
934 :
935 :
936 : } // namespace libMesh
|
