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_prism18.h"
21 : #include "libmesh/edge_edge3.h"
22 : #include "libmesh/face_quad9.h"
23 : #include "libmesh/face_tri6.h"
24 : #include "libmesh/enum_io_package.h"
25 : #include "libmesh/enum_order.h"
26 : #include "libmesh/int_range.h"
27 :
28 : namespace libMesh
29 : {
30 :
31 :
32 :
33 : // ------------------------------------------------------------
34 : // Prism18 class static member initializations
35 : const int Prism18::num_nodes;
36 : const int Prism18::nodes_per_side;
37 : const int Prism18::nodes_per_edge;
38 :
39 : const unsigned int Prism18::side_nodes_map[Prism18::num_sides][Prism18::nodes_per_side] =
40 : {
41 : {0, 2, 1, 8, 7, 6, 99, 99, 99}, // Side 0
42 : {0, 1, 4, 3, 6, 10, 12, 9, 15}, // Side 1
43 : {1, 2, 5, 4, 7, 11, 13, 10, 16}, // Side 2
44 : {2, 0, 3, 5, 8, 9, 14, 11, 17}, // Side 3
45 : {3, 4, 5, 12, 13, 14, 99, 99, 99} // Side 4
46 : };
47 :
48 : const unsigned int Prism18::edge_nodes_map[Prism18::num_edges][Prism18::nodes_per_edge] =
49 : {
50 : {0, 1, 6}, // Edge 0
51 : {1, 2, 7}, // Edge 1
52 : {0, 2, 8}, // Edge 2
53 : {0, 3, 9}, // Edge 3
54 : {1, 4, 10}, // Edge 4
55 : {2, 5, 11}, // Edge 5
56 : {3, 4, 12}, // Edge 6
57 : {4, 5, 13}, // Edge 7
58 : {3, 5, 14} // Edge 8
59 : };
60 :
61 : // ------------------------------------------------------------
62 : // Prism18 class member functions
63 :
64 2041158 : bool Prism18::is_vertex(const unsigned int i) const
65 : {
66 2041158 : if (i < 6)
67 643555 : return true;
68 140097 : return false;
69 : }
70 :
71 13824 : bool Prism18::is_edge(const unsigned int i) const
72 : {
73 13824 : if (i < 6)
74 0 : return false;
75 13824 : if (i > 14)
76 3456 : return false;
77 864 : return true;
78 : }
79 :
80 909293 : bool Prism18::is_face(const unsigned int i) const
81 : {
82 909293 : if (i > 14)
83 117218 : return true;
84 117315 : return false;
85 : }
86 :
87 166610 : bool Prism18::is_node_on_side(const unsigned int n,
88 : const unsigned int s) const
89 : {
90 20365 : libmesh_assert_less (s, n_sides());
91 20365 : return std::find(std::begin(side_nodes_map[s]),
92 166610 : std::end(side_nodes_map[s]),
93 166610 : n) != std::end(side_nodes_map[s]);
94 : }
95 :
96 : std::vector<unsigned>
97 3592909 : Prism18::nodes_on_side(const unsigned int s) const
98 : {
99 320006 : libmesh_assert_less(s, n_sides());
100 3592909 : auto trim = (s > 0 && s < 4) ? 0 : 3;
101 3592909 : return {std::begin(side_nodes_map[s]), std::end(side_nodes_map[s]) - trim};
102 : }
103 :
104 : std::vector<unsigned>
105 11007 : Prism18::nodes_on_edge(const unsigned int e) const
106 : {
107 882 : libmesh_assert_less(e, n_edges());
108 11889 : return {std::begin(edge_nodes_map[e]), std::end(edge_nodes_map[e])};
109 : }
110 :
111 50382 : bool Prism18::is_node_on_edge(const unsigned int n,
112 : const unsigned int e) const
113 : {
114 3672 : libmesh_assert_less (e, n_edges());
115 3672 : return std::find(std::begin(edge_nodes_map[e]),
116 50382 : std::end(edge_nodes_map[e]),
117 50382 : n) != std::end(edge_nodes_map[e]);
118 : }
119 :
120 :
121 :
122 581117 : bool Prism18::has_affine_map() const
123 : {
124 : // Make sure z edges are affine
125 115758 : Point v = this->point(3) - this->point(0);
126 696875 : if (!v.relative_fuzzy_equals(this->point(4) - this->point(1), affine_tol) ||
127 696875 : !v.relative_fuzzy_equals(this->point(5) - this->point(2), affine_tol))
128 0 : return false;
129 : // Make sure edges are straight
130 57879 : v /= 2;
131 638996 : if (!v.relative_fuzzy_equals(this->point(9) - this->point(0), affine_tol) ||
132 638996 : !v.relative_fuzzy_equals(this->point(10) - this->point(1), affine_tol) ||
133 638996 : !v.relative_fuzzy_equals(this->point(11) - this->point(2), affine_tol) ||
134 638996 : !v.relative_fuzzy_equals(this->point(15) - this->point(6), affine_tol) ||
135 1277992 : !v.relative_fuzzy_equals(this->point(16) - this->point(7), affine_tol) ||
136 638996 : !v.relative_fuzzy_equals(this->point(17) - this->point(8), affine_tol))
137 0 : return false;
138 638996 : v = (this->point(1) - this->point(0))/2;
139 696875 : if (!v.relative_fuzzy_equals(this->point(6) - this->point(0), affine_tol) ||
140 638996 : !v.relative_fuzzy_equals(this->point(12) - this->point(3), affine_tol))
141 0 : return false;
142 638996 : v = (this->point(2) - this->point(0))/2;
143 696875 : if (!v.relative_fuzzy_equals(this->point(8) - this->point(0), affine_tol) ||
144 638996 : !v.relative_fuzzy_equals(this->point(14) - this->point(3), affine_tol))
145 0 : return false;
146 638996 : v = (this->point(2) - this->point(1))/2;
147 696875 : if (!v.relative_fuzzy_equals(this->point(7) - this->point(1), affine_tol) ||
148 696875 : !v.relative_fuzzy_equals(this->point(13) - this->point(4), affine_tol))
149 0 : return false;
150 57879 : return true;
151 : }
152 :
153 :
154 :
155 18736407 : Order Prism18::default_order() const
156 : {
157 18736407 : return SECOND;
158 : }
159 :
160 0 : dof_id_type Prism18::key (const unsigned int s) const
161 : {
162 0 : libmesh_assert_less (s, this->n_sides());
163 :
164 0 : switch (s)
165 : {
166 0 : case 0: // the triangular face at z=0
167 : {
168 0 : return Prism::key(0);
169 : }
170 0 : case 1: // the quad face at y=0
171 : {
172 0 : return Elem::compute_key (this->node_id(15));
173 : }
174 0 : case 2: // the other quad face
175 : {
176 0 : return Elem::compute_key (this->node_id(16));
177 : }
178 0 : case 3: // the quad face at x=0
179 : {
180 0 : return Elem::compute_key (this->node_id(17));
181 : }
182 0 : case 4: // the triangular face at z=1
183 : {
184 0 : return Prism::key(4);
185 : }
186 0 : default:
187 0 : libmesh_error_msg("Invalid side " << s);
188 : }
189 : }
190 :
191 :
192 :
193 322856 : unsigned int Prism18::local_side_node(unsigned int side,
194 : unsigned int side_node) const
195 : {
196 63350 : libmesh_assert_less (side, this->n_sides());
197 :
198 : // Never more than 9 nodes per side.
199 63350 : libmesh_assert_less(side_node, Prism18::nodes_per_side);
200 :
201 : // Some sides have 6 nodes.
202 63350 : libmesh_assert(!(side==0 || side==4) || side_node < 6);
203 :
204 322856 : return Prism18::side_nodes_map[side][side_node];
205 : }
206 :
207 :
208 :
209 7216929 : unsigned int Prism18::local_edge_node(unsigned int edge,
210 : unsigned int edge_node) const
211 : {
212 642564 : libmesh_assert_less(edge, this->n_edges());
213 642564 : libmesh_assert_less(edge_node, Prism18::nodes_per_edge);
214 :
215 7216929 : return Prism18::edge_nodes_map[edge][edge_node];
216 : }
217 :
218 :
219 :
220 524343 : std::unique_ptr<Elem> Prism18::build_side_ptr (const unsigned int i)
221 : {
222 40384 : libmesh_assert_less (i, this->n_sides());
223 :
224 524343 : std::unique_ptr<Elem> face;
225 :
226 524343 : switch (i)
227 : {
228 294365 : case 0: // the triangular face at z=-1
229 : case 4: // the triangular face at z=1
230 : {
231 294365 : face = std::make_unique<Tri6>();
232 294365 : break;
233 : }
234 229978 : case 1: // the quad face at y=0
235 : case 2: // the other quad face
236 : case 3: // the quad face at x=0
237 : {
238 229978 : face = std::make_unique<Quad9>();
239 229978 : break;
240 : }
241 0 : default:
242 0 : libmesh_error_msg("Invalid side i = " << i);
243 : }
244 :
245 : // Set the nodes
246 4360335 : for (auto n : face->node_index_range())
247 4135914 : face->set_node(n, this->node_ptr(Prism18::side_nodes_map[i][n]));
248 :
249 524343 : face->set_interior_parent(this);
250 483959 : face->inherit_data_from(*this);
251 :
252 524343 : return face;
253 0 : }
254 :
255 :
256 :
257 436595 : void Prism18::build_side_ptr (std::unique_ptr<Elem> & side,
258 : const unsigned int i)
259 : {
260 25258 : libmesh_assert_less (i, this->n_sides());
261 :
262 436595 : switch (i)
263 : {
264 7292 : case 0: // the triangular face at z=-1
265 : case 4: // the triangular face at z=1
266 : {
267 111900 : if (!side.get() || side->type() != TRI6)
268 : {
269 207808 : side = this->build_side_ptr(i);
270 111132 : return;
271 : }
272 64 : break;
273 : }
274 :
275 17966 : case 1: // the quad face at y=0
276 : case 2: // the other quad face
277 : case 3: // the quad face at x=0
278 : {
279 324695 : if (!side.get() || side->type() != QUAD9)
280 : {
281 207588 : side = this->build_side_ptr(i);
282 110370 : return;
283 : }
284 11390 : break;
285 : }
286 :
287 0 : default:
288 0 : libmesh_error_msg("Invalid side i = " << i);
289 : }
290 :
291 203639 : side->inherit_data_from(*this);
292 :
293 : // Set the nodes
294 2148626 : for (auto n : side->node_index_range())
295 2036427 : side->set_node(n, this->node_ptr(Prism18::side_nodes_map[i][n]));
296 : }
297 :
298 :
299 :
300 859747 : std::unique_ptr<Elem> Prism18::build_edge_ptr (const unsigned int i)
301 : {
302 859747 : return this->simple_build_edge_ptr<Edge3,Prism18>(i);
303 : }
304 :
305 :
306 :
307 0 : void Prism18::build_edge_ptr (std::unique_ptr<Elem> & edge, const unsigned int i)
308 : {
309 0 : this->simple_build_edge_ptr<Prism18>(edge, i, EDGE3);
310 0 : }
311 :
312 :
313 :
314 0 : void Prism18::connectivity(const unsigned int sc,
315 : const IOPackage iop,
316 : std::vector<dof_id_type> & conn) const
317 : {
318 0 : libmesh_assert(_nodes);
319 0 : libmesh_assert_less (sc, this->n_sub_elem());
320 0 : libmesh_assert_not_equal_to (iop, INVALID_IO_PACKAGE);
321 :
322 0 : switch (iop)
323 : {
324 0 : case TECPLOT:
325 : {
326 0 : conn.resize(8);
327 0 : switch (sc)
328 : {
329 :
330 0 : case 0:
331 : {
332 0 : conn[0] = this->node_id(0)+1;
333 0 : conn[1] = this->node_id(6)+1;
334 0 : conn[2] = this->node_id(8)+1;
335 0 : conn[3] = this->node_id(8)+1;
336 0 : conn[4] = this->node_id(9)+1;
337 0 : conn[5] = this->node_id(15)+1;
338 0 : conn[6] = this->node_id(17)+1;
339 0 : conn[7] = this->node_id(17)+1;
340 :
341 0 : return;
342 : }
343 :
344 0 : case 1:
345 : {
346 0 : conn[0] = this->node_id(6)+1;
347 0 : conn[1] = this->node_id(1)+1;
348 0 : conn[2] = this->node_id(7)+1;
349 0 : conn[3] = this->node_id(7)+1;
350 0 : conn[4] = this->node_id(15)+1;
351 0 : conn[5] = this->node_id(10)+1;
352 0 : conn[6] = this->node_id(16)+1;
353 0 : conn[7] = this->node_id(16)+1;
354 :
355 0 : return;
356 : }
357 :
358 0 : case 2:
359 : {
360 0 : conn[0] = this->node_id(8)+1;
361 0 : conn[1] = this->node_id(7)+1;
362 0 : conn[2] = this->node_id(2)+1;
363 0 : conn[3] = this->node_id(2)+1;
364 0 : conn[4] = this->node_id(17)+1;
365 0 : conn[5] = this->node_id(16)+1;
366 0 : conn[6] = this->node_id(11)+1;
367 0 : conn[7] = this->node_id(11)+1;
368 :
369 0 : return;
370 : }
371 :
372 0 : case 3:
373 : {
374 0 : conn[0] = this->node_id(6)+1;
375 0 : conn[1] = this->node_id(7)+1;
376 0 : conn[2] = this->node_id(8)+1;
377 0 : conn[3] = this->node_id(8)+1;
378 0 : conn[4] = this->node_id(15)+1;
379 0 : conn[5] = this->node_id(16)+1;
380 0 : conn[6] = this->node_id(17)+1;
381 0 : conn[7] = this->node_id(17)+1;
382 :
383 0 : return;
384 : }
385 :
386 0 : case 4:
387 : {
388 0 : conn[0] = this->node_id(9)+1;
389 0 : conn[1] = this->node_id(15)+1;
390 0 : conn[2] = this->node_id(17)+1;
391 0 : conn[3] = this->node_id(17)+1;
392 0 : conn[4] = this->node_id(3)+1;
393 0 : conn[5] = this->node_id(12)+1;
394 0 : conn[6] = this->node_id(14)+1;
395 0 : conn[7] = this->node_id(14)+1;
396 :
397 0 : return;
398 : }
399 :
400 0 : case 5:
401 : {
402 0 : conn[0] = this->node_id(15)+1;
403 0 : conn[1] = this->node_id(10)+1;
404 0 : conn[2] = this->node_id(16)+1;
405 0 : conn[3] = this->node_id(16)+1;
406 0 : conn[4] = this->node_id(12)+1;
407 0 : conn[5] = this->node_id(4)+1;
408 0 : conn[6] = this->node_id(13)+1;
409 0 : conn[7] = this->node_id(13)+1;
410 :
411 0 : return;
412 : }
413 :
414 0 : case 6:
415 : {
416 0 : conn[0] = this->node_id(17)+1;
417 0 : conn[1] = this->node_id(16)+1;
418 0 : conn[2] = this->node_id(11)+1;
419 0 : conn[3] = this->node_id(11)+1;
420 0 : conn[4] = this->node_id(14)+1;
421 0 : conn[5] = this->node_id(13)+1;
422 0 : conn[6] = this->node_id(5)+1;
423 0 : conn[7] = this->node_id(5)+1;
424 :
425 0 : return;
426 : }
427 :
428 0 : case 7:
429 : {
430 0 : conn[0] = this->node_id(15)+1;
431 0 : conn[1] = this->node_id(16)+1;
432 0 : conn[2] = this->node_id(17)+1;
433 0 : conn[3] = this->node_id(17)+1;
434 0 : conn[4] = this->node_id(12)+1;
435 0 : conn[5] = this->node_id(13)+1;
436 0 : conn[6] = this->node_id(14)+1;
437 0 : conn[7] = this->node_id(14)+1;
438 :
439 0 : return;
440 : }
441 :
442 0 : default:
443 0 : libmesh_error_msg("Invalid sc = " << sc);
444 : }
445 :
446 : }
447 :
448 0 : case VTK:
449 : {
450 : // VTK now supports VTK_BIQUADRATIC_QUADRATIC_WEDGE directly
451 0 : const unsigned int conn_size = 18;
452 0 : conn.resize(conn_size);
453 :
454 : // VTK's VTK_BIQUADRATIC_QUADRATIC_WEDGE first 9 (vertex) and
455 : // last 3 (mid-face) nodes match. The middle and top layers
456 : // of mid-edge nodes are reversed from LibMesh's.
457 0 : for (auto i : index_range(conn))
458 0 : conn[i] = this->node_id(i);
459 :
460 : // top "ring" of mid-edge nodes
461 0 : conn[9] = this->node_id(12);
462 0 : conn[10] = this->node_id(13);
463 0 : conn[11] = this->node_id(14);
464 :
465 : // middle "ring" of mid-edge nodes
466 0 : conn[12] = this->node_id(9);
467 0 : conn[13] = this->node_id(10);
468 0 : conn[14] = this->node_id(11);
469 :
470 0 : return;
471 : }
472 :
473 0 : default:
474 0 : libmesh_error_msg("Unsupported IO package " << iop);
475 : }
476 : }
477 :
478 :
479 :
480 :
481 2592 : unsigned int Prism18::n_second_order_adjacent_vertices (const unsigned int n) const
482 : {
483 2592 : switch (n)
484 : {
485 108 : case 6:
486 : case 7:
487 : case 8:
488 : case 9:
489 : case 10:
490 : case 11:
491 : case 12:
492 : case 13:
493 : case 14:
494 108 : return 2;
495 :
496 648 : case 15:
497 : case 16:
498 : case 17:
499 648 : return 4;
500 :
501 0 : default:
502 0 : libmesh_error_msg("Invalid node n = " << n);
503 : }
504 : }
505 :
506 :
507 :
508 :
509 :
510 6480 : unsigned short int Prism18::second_order_adjacent_vertex (const unsigned int n,
511 : const unsigned int v) const
512 : {
513 360 : libmesh_assert_greater_equal (n, this->n_vertices());
514 360 : libmesh_assert_less (n, this->n_nodes());
515 :
516 6480 : switch (n)
517 : {
518 : /*
519 : * These nodes are unique to \p Prism18,
520 : * let our _remaining_... matrix handle
521 : * this.
522 : */
523 2592 : case 15:
524 : case 16:
525 : case 17:
526 : {
527 144 : libmesh_assert_less (v, 4);
528 2592 : return _remaining_second_order_adjacent_vertices[n-15][v];
529 : }
530 :
531 : /*
532 : * All other second-order nodes (6,...,14) are
533 : * identical with Prism15 and are therefore
534 : * delegated to the _second_order matrix of
535 : * \p Prism
536 : */
537 3888 : default:
538 : {
539 216 : libmesh_assert_less (v, 2);
540 3888 : return _second_order_adjacent_vertices[n-this->n_vertices()][v];
541 : }
542 :
543 : }
544 :
545 : // libmesh_error_msg("We'll never get here!"); // static checkers agree
546 : return static_cast<unsigned short int>(-1);
547 : }
548 :
549 :
550 :
551 : const unsigned short int Prism18::_remaining_second_order_adjacent_vertices[3][4] =
552 : {
553 : { 0, 1, 3, 4}, // vertices adjacent to node 15
554 : { 1, 2, 4, 5}, // vertices adjacent to node 16
555 : { 0, 2, 3, 5} // vertices adjacent to node 17
556 : };
557 :
558 :
559 :
560 : std::pair<unsigned short int, unsigned short int>
561 0 : Prism18::second_order_child_vertex (const unsigned int n) const
562 : {
563 0 : libmesh_assert_greater_equal (n, this->n_vertices());
564 0 : libmesh_assert_less (n, this->n_nodes());
565 :
566 0 : return std::pair<unsigned short int, unsigned short int>
567 0 : (_second_order_vertex_child_number[n],
568 0 : _second_order_vertex_child_index[n]);
569 : }
570 :
571 :
572 :
573 153 : Real Prism18::volume () const
574 : {
575 : // This specialization is good for Lagrange mappings only
576 153 : if (this->mapping_type() != LAGRANGE_MAP)
577 0 : return this->Elem::volume();
578 :
579 : // Make copies of our points. It makes the subsequent calculations a bit
580 : // shorter and avoids dereferencing the same pointer multiple times.
581 : Point
582 459 : x0 = point(0), x1 = point(1), x2 = point(2), x3 = point(3), x4 = point(4), x5 = point(5),
583 408 : x6 = point(6), x7 = point(7), x8 = point(8), x9 = point(9), x10 = point(10), x11 = point(11),
584 408 : x12 = point(12), x13 = point(13), x14 = point(14), x15 = point(15), x16 = point(16), x17 = point(17);
585 :
586 : // The number of components in the dx_dxi, dx_deta, and dx_dzeta arrays.
587 51 : const int n_components = 16;
588 :
589 : // Terms are copied directly from a Python script.
590 : Point dx_dxi[n_components] =
591 : {
592 51 : -x10 + 4*x15 - 3*x9,
593 51 : 3*x0/2 + x1/2 + 2*x12 - 3*x3/2 - x4/2 - 2*x6,
594 51 : -3*x0/2 - x1/2 + x10 + 2*x12 - 4*x15 - 3*x3/2 - x4/2 + 2*x6 + 3*x9,
595 51 : -4*x15 + 4*x16 - 4*x17 + 4*x9,
596 51 : -2*x0 - 2*x12 + 2*x13 - 2*x14 + 2*x3 + 2*x6 - 2*x7 + 2*x8,
597 51 : 2*x0 - 2*x12 + 2*x13 - 2*x14 + 4*x15 - 4*x16 + 4*x17 + 2*x3 - 2*x6 + 2*x7 - 2*x8 - 4*x9,
598 : Point(0,0,0),
599 : Point(0,0,0),
600 51 : 4*x10 - 8*x15 + 4*x9,
601 51 : -2*x0 - 2*x1 - 4*x12 + 2*x3 + 2*x4 + 4*x6,
602 51 : 2*x0 + 2*x1 - 4*x10 - 4*x12 + 8*x15 + 2*x3 + 2*x4 - 4*x6 - 4*x9,
603 : Point(0,0,0),
604 : Point(0,0,0),
605 : Point(0,0,0),
606 : Point(0,0,0),
607 : Point(0,0,0)
608 510 : };
609 :
610 : Point dx_deta[n_components] =
611 : {
612 51 : -x11 + 4*x17 - 3*x9,
613 51 : 3*x0/2 + 2*x14 + x2/2 - 3*x3/2 - x5/2 - 2*x8,
614 51 : -3*x0/2 + x11 + 2*x14 - 4*x17 - x2/2 - 3*x3/2 - x5/2 + 2*x8 + 3*x9,
615 51 : 4*x11 - 8*x17 + 4*x9,
616 51 : -2*x0 - 4*x14 - 2*x2 + 2*x3 + 2*x5 + 4*x8,
617 51 : 2*x0 - 4*x11 - 4*x14 + 8*x17 + 2*x2 + 2*x3 + 2*x5 - 4*x8 - 4*x9,
618 : Point(0,0,0),
619 : Point(0,0,0),
620 51 : -4*x15 + 4*x16 - 4*x17 + 4*x9,
621 51 : -2*x0 - 2*x12 + 2*x13 - 2*x14 + 2*x3 + 2*x6 - 2*x7 + 2*x8,
622 51 : 2*x0 - 2*x12 + 2*x13 - 2*x14 + 4*x15 - 4*x16 + 4*x17 + 2*x3 - 2*x6 + 2*x7 - 2*x8 - 4*x9,
623 : Point(0,0,0),
624 : Point(0,0,0),
625 : Point(0,0,0),
626 : Point(0,0,0),
627 : Point(0,0,0)
628 510 : };
629 :
630 : Point dx_dzeta[n_components] =
631 : {
632 51 : -x0/2 + x3/2,
633 51 : x0 + x3 - 2*x9,
634 : Point(0,0,0),
635 51 : 3*x0/2 + 2*x14 + x2/2 - 3*x3/2 - x5/2 - 2*x8,
636 51 : -3*x0 + 2*x11 + 4*x14 - 8*x17 - x2 - 3*x3 - x5 + 4*x8 + 6*x9,
637 : Point(0,0,0),
638 51 : -x0 - 2*x14 - x2 + x3 + x5 + 2*x8,
639 51 : 2*x0 - 4*x11 - 4*x14 + 8*x17 + 2*x2 + 2*x3 + 2*x5 - 4*x8 - 4*x9,
640 51 : 3*x0/2 + x1/2 + 2*x12 - 3*x3/2 - x4/2 - 2*x6,
641 51 : -3*x0 - x1 + 2*x10 + 4*x12 - 8*x15 - 3*x3 - x4 + 4*x6 + 6*x9,
642 : Point(0,0,0),
643 51 : -2*x0 - 2*x12 + 2*x13 - 2*x14 + 2*x3 + 2*x6 - 2*x7 + 2*x8,
644 51 : 4*x0 - 4*x12 + 4*x13 - 4*x14 + 8*x15 - 8*x16 + 8*x17 + 4*x3 - 4*x6 + 4*x7 - 4*x8 - 8*x9,
645 : Point(0,0,0),
646 51 : -x0 - x1 - 2*x12 + x3 + x4 + 2*x6,
647 51 : 2*x0 + 2*x1 - 4*x10 - 4*x12 + 8*x15 + 2*x3 + 2*x4 - 4*x6 - 4*x9
648 663 : };
649 :
650 : // The quadrature rule for the Prism18 is a tensor product between a
651 : // FOURTH-order TRI rule (in xi, eta) and a FIFTH-order EDGE rule
652 : // in zeta.
653 :
654 : // Number of points in the 2D quadrature rule.
655 51 : const int N2D = 6;
656 :
657 : // Parameters of the 2D rule
658 : static const Real
659 : w1 = 1.1169079483900573284750350421656140e-01_R,
660 : w2 = 5.4975871827660933819163162450105264e-02_R,
661 : a1 = 4.4594849091596488631832925388305199e-01_R,
662 : a2 = 9.1576213509770743459571463402201508e-02_R;
663 :
664 : // Points and weights of the 2D rule
665 : static const Real w2D[N2D] = {w1, w1, w1, w2, w2, w2};
666 :
667 : // Quadrature point locations raised to powers. xi[0][2] is
668 : // quadrature point 0, squared, xi[1][1] is quadrature point 1 to the
669 : // first power, etc. This lets us avoid calling std::pow inside the
670 : // loops below.
671 : static const Real xi[N2D][3] =
672 : {
673 : // ^0 ^1 ^2
674 : { 1., a1, a1*a1},
675 : { 1., 1-2*a1, (1-2*a1)*(1-2*a1)},
676 : { 1., a1, a1*a1},
677 : { 1., a2, a2*a2},
678 : { 1., 1-2*a2, (1-2*a2)*(1-2*a2)},
679 : { 1., a2, a2*a2}
680 : };
681 :
682 : static const Real eta[N2D][3] =
683 : {
684 : // ^0 ^1 ^2
685 : { 1., a1, a1*a1},
686 : { 1., a1, a1*a1},
687 : { 1., 1-2*a1, (1-2*a1)*(1-2*a1)},
688 : { 1., a2, a2*a2},
689 : { 1., a2, a2*a2},
690 : { 1., 1-2*a2, (1-2*a2)*(1-2*a2)}
691 : };
692 :
693 : // Number of points in the 1D quadrature rule.
694 51 : const int N1D = 3;
695 :
696 : // Points and weights of the 1D quadrature rule.
697 : static const Real w1D[N1D] = {5./9, 8./9, 5./9};
698 :
699 153 : const Real zeta[N1D][3] =
700 : {
701 : //^0 ^1 ^2
702 : { 1., -std::sqrt(15)/5., 15./25},
703 : { 1., 0., 0.},
704 : { 1., std::sqrt(15)/5., 15./25}
705 : };
706 :
707 : // The integer exponents for each term.
708 : static const int exponents[n_components][3] =
709 : {
710 : {0, 0, 0},
711 : {0, 0, 1},
712 : {0, 0, 2},
713 : {0, 1, 0},
714 : {0, 1, 1},
715 : {0, 1, 2},
716 : {0, 2, 0},
717 : {0, 2, 1},
718 : {1, 0, 0},
719 : {1, 0, 1},
720 : {1, 0, 2},
721 : {1, 1, 0},
722 : {1, 1, 1},
723 : {1, 2, 0},
724 : {2, 0, 0},
725 : {2, 0, 1}
726 : };
727 :
728 51 : Real vol = 0.;
729 1071 : for (int i=0; i<N2D; ++i)
730 3672 : for (int j=0; j<N1D; ++j)
731 : {
732 : // Compute dx_dxi, dx_deta, dx_dzeta at the current quadrature point.
733 918 : Point dx_dxi_q, dx_deta_q, dx_dzeta_q;
734 46818 : for (int c=0; c<n_components; ++c)
735 : {
736 14688 : Real coeff =
737 73440 : xi[i][exponents[c][0]]*
738 73440 : eta[i][exponents[c][1]]*
739 44064 : zeta[j][exponents[c][2]];
740 :
741 14688 : dx_dxi_q += coeff * dx_dxi[c];
742 14688 : dx_deta_q += coeff * dx_deta[c];
743 14688 : dx_dzeta_q += coeff * dx_dzeta[c];
744 : }
745 :
746 : // Compute scalar triple product, multiply by weight, and accumulate volume.
747 3672 : vol += w2D[i] * w1D[j] * triple_product(dx_dxi_q, dx_deta_q, dx_dzeta_q);
748 : }
749 :
750 51 : return vol;
751 : }
752 :
753 :
754 :
755 :
756 : #ifdef LIBMESH_ENABLE_AMR
757 :
758 : const Real Prism18::_embedding_matrix[Prism18::num_children][Prism18::num_nodes][Prism18::num_nodes] =
759 : {
760 : // embedding matrix for child 0
761 : {
762 : // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
763 : { 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}, // 0
764 : { 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}, // 1
765 : { 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0.}, // 2
766 : { 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0.}, // 3
767 : { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0.}, // 4
768 : { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1.}, // 5
769 : { 0.375, -0.125, 0., 0., 0., 0., 0.75, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}, // 6
770 : { 0., -0.125, -0.125, 0., 0., 0., 0.5, 0.25, 0.5, 0., 0., 0., 0., 0., 0., 0., 0., 0.}, // 7
771 : { 0.375, 0., -0.125, 0., 0., 0., 0., 0., 0.75, 0., 0., 0., 0., 0., 0., 0., 0., 0.}, // 8
772 : { 0.375, 0., 0., -0.125, 0., 0., 0., 0., 0., 0.75, 0., 0., 0., 0., 0., 0., 0., 0.}, // 9
773 : { 0., 0., 0., 0., 0., 0., 0.375, 0., 0., 0., 0., 0., -0.125, 0., 0., 0.75, 0., 0.}, // 10
774 : { 0., 0., 0., 0., 0., 0., 0., 0., 0.375, 0., 0., 0., 0., 0., -0.125, 0., 0., 0.75}, // 11
775 : { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.375, -0.125, 0., 0., 0., 0., 0.75, 0., 0.}, // 12
776 : { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., -0.125, -0.125, 0., 0., 0., 0.5, 0.25, 0.5}, // 13
777 : { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.375, 0., -0.125, 0., 0., 0., 0., 0., 0.75}, // 14
778 : { 0.140625,-0.046875, 0.,-0.046875, 0.015625, 0., 0.28125, 0., 0., 0.28125, -0.09375, 0., -0.09375, 0., 0., 0.5625, 0., 0.}, // 15
779 : { 0.,-0.046875,-0.046875, 0., 0.015625, 0.015625, 0.1875, 0.09375, 0.1875, 0., -0.09375, -0.09375, -0.0625, -0.03125, -0.0625, 0.375, 0.1875, 0.375}, // 16
780 : { 0.140625, 0.,-0.046875,-0.046875, 0., 0.015625, 0., 0., 0.28125, 0.28125, 0., -0.09375, 0., 0., -0.09375, 0., 0., 0.5625} // 17
781 : },
782 :
783 : // embedding matrix for child 1
784 : {
785 : // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
786 : { 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}, // 0
787 : { 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}, // 1
788 : { 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}, // 2
789 : { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0.}, // 3
790 : { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0.}, // 4
791 : { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0.}, // 5
792 : { -0.125, 0.375, 0., 0., 0., 0., 0.75, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}, // 6
793 : { 0., 0.375, -0.125, 0., 0., 0., 0., 0.75, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}, // 7
794 : { -0.125, 0., -0.125, 0., 0., 0., 0.5, 0.5, 0.25, 0., 0., 0., 0., 0., 0., 0., 0., 0.}, // 8
795 : { 0., 0., 0., 0., 0., 0., 0.375, 0., 0., 0., 0., 0., -0.125, 0., 0., 0.75, 0., 0.}, // 9
796 : { 0., 0.375, 0., 0., -0.125, 0., 0., 0., 0., 0., 0.75, 0., 0., 0., 0., 0., 0., 0.}, // 10
797 : { 0., 0., 0., 0., 0., 0., 0., 0.375, 0., 0., 0., 0., 0., -0.125, 0., 0., 0.75, 0.}, // 11
798 : { 0., 0., 0., 0., 0., 0., 0., 0., 0., -0.125, 0.375, 0., 0., 0., 0., 0.75, 0., 0.}, // 12
799 : { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.375, -0.125, 0., 0., 0., 0., 0.75, 0.}, // 13
800 : { 0., 0., 0., 0., 0., 0., 0., 0., 0., -0.125, 0., -0.125, 0., 0., 0., 0.5, 0.5, 0.25}, // 14
801 : {-0.046875, 0.140625, 0., 0.015625,-0.046875, 0., 0.28125, 0., 0., -0.09375, 0.28125, 0., -0.09375, 0., 0., 0.5625, 0., 0.}, // 15
802 : { 0., 0.140625,-0.046875, 0.,-0.046875, 0.015625, 0., 0.28125, 0., 0., 0.28125, -0.09375, 0., -0.09375, 0., 0., 0.5625, 0.}, // 16
803 : {-0.046875, 0.,-0.046875, 0.015625, 0., 0.015625, 0.1875, 0.1875, 0.09375, -0.09375, 0., -0.09375, -0.0625, -0.0625, -0.03125, 0.375, 0.375, 0.1875} // 17
804 : },
805 :
806 : // embedding matrix for child 2
807 : {
808 : // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
809 : { 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0.}, // 0
810 : { 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}, // 1
811 : { 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}, // 2
812 : { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1.}, // 3
813 : { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0.}, // 4
814 : { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0.}, // 5
815 : { -0.125, -0.125, 0., 0., 0., 0., 0.25, 0.5, 0.5, 0., 0., 0., 0., 0., 0., 0., 0., 0.}, // 6
816 : { 0., -0.125, 0.375, 0., 0., 0., 0., 0.75, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}, // 7
817 : { -0.125, 0., 0.375, 0., 0., 0., 0., 0., 0.75, 0., 0., 0., 0., 0., 0., 0., 0., 0.}, // 8
818 : { 0., 0., 0., 0., 0., 0., 0., 0., 0.375, 0., 0., 0., 0., 0., -0.125, 0., 0., 0.75}, // 9
819 : { 0., 0., 0., 0., 0., 0., 0., 0.375, 0., 0., 0., 0., 0., -0.125, 0., 0., 0.75, 0.}, // 10
820 : { 0., 0., 0.375, 0., 0., -0.125, 0., 0., 0., 0., 0., 0.75, 0., 0., 0., 0., 0., 0.}, // 11
821 : { 0., 0., 0., 0., 0., 0., 0., 0., 0., -0.125, -0.125, 0., 0., 0., 0., 0.25, 0.5, 0.5}, // 12
822 : { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., -0.125, 0.375, 0., 0., 0., 0., 0.75, 0.}, // 13
823 : { 0., 0., 0., 0., 0., 0., 0., 0., 0., -0.125, 0., 0.375, 0., 0., 0., 0., 0., 0.75}, // 14
824 : {-0.046875,-0.046875, 0., 0.015625, 0.015625, 0., 0.09375, 0.1875, 0.1875, -0.09375, -0.09375, 0., -0.03125, -0.0625, -0.0625, 0.1875, 0.375, 0.375}, // 15
825 : { 0.,-0.046875, 0.140625, 0., 0.015625,-0.046875, 0., 0.28125, 0., 0., -0.09375, 0.28125, 0., -0.09375, 0., 0., 0.5625, 0.}, // 16
826 : {-0.046875, 0., 0.140625, 0.015625, 0.,-0.046875, 0., 0., 0.28125, -0.09375, 0., 0.28125, 0., 0., -0.09375, 0., 0., 0.5625} // 17
827 : },
828 :
829 : // embedding matrix for child 3
830 : {
831 : // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
832 : { 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}, // 0
833 : { 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}, // 1
834 : { 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0.}, // 2
835 : { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0.}, // 3
836 : { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0.}, // 4
837 : { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1.}, // 5
838 : { -0.125, 0., -0.125, 0., 0., 0., 0.5, 0.5, 0.25, 0., 0., 0., 0., 0., 0., 0., 0., 0.}, // 6
839 : { -0.125, -0.125, 0., 0., 0., 0., 0.25, 0.5, 0.5, 0., 0., 0., 0., 0., 0., 0., 0., 0.}, // 7
840 : { 0., -0.125, -0.125, 0., 0., 0., 0.5, 0.25, 0.5, 0., 0., 0., 0., 0., 0., 0., 0., 0.}, // 8
841 : { 0., 0., 0., 0., 0., 0., 0.375, 0., 0., 0., 0., 0., -0.125, 0., 0., 0.75, 0., 0.}, // 9
842 : { 0., 0., 0., 0., 0., 0., 0., 0.375, 0., 0., 0., 0., 0., -0.125, 0., 0., 0.75, 0.}, // 10
843 : { 0., 0., 0., 0., 0., 0., 0., 0., 0.375, 0., 0., 0., 0., 0., -0.125, 0., 0., 0.75}, // 11
844 : { 0., 0., 0., 0., 0., 0., 0., 0., 0., -0.125, 0., -0.125, 0., 0., 0., 0.5, 0.5, 0.25}, // 12
845 : { 0., 0., 0., 0., 0., 0., 0., 0., 0., -0.125, -0.125, 0., 0., 0., 0., 0.25, 0.5, 0.5}, // 13
846 : { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., -0.125, -0.125, 0., 0., 0., 0.5, 0.25, 0.5}, // 14
847 : {-0.046875, 0.,-0.046875, 0.015625, 0., 0.015625, 0.1875, 0.1875, 0.09375, -0.09375, 0., -0.09375, -0.0625, -0.0625, -0.03125, 0.375, 0.375, 0.1875}, // 15
848 : {-0.046875,-0.046875, 0., 0.015625, 0.015625, 0., 0.09375, 0.1875, 0.1875, -0.09375, -0.09375, 0., -0.03125, -0.0625, -0.0625, 0.1875, 0.375, 0.375}, // 16
849 : { 0.,-0.046875,-0.046875, 0., 0.015625, 0.015625, 0.1875, 0.09375, 0.1875, 0., -0.09375, -0.09375, -0.0625, -0.03125, -0.0625, 0.375, 0.1875, 0.375} // 17
850 : },
851 :
852 : // embedding matrix for child 4
853 : {
854 : // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
855 : { 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0.}, // 0
856 : { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0.}, // 1
857 : { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1.}, // 2
858 : { 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}, // 3
859 : { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0.}, // 4
860 : { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0.}, // 5
861 : { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.375, -0.125, 0., 0., 0., 0., 0.75, 0., 0.}, // 6
862 : { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., -0.125, -0.125, 0., 0., 0., 0.5, 0.25, 0.5}, // 7
863 : { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.375, 0., -0.125, 0., 0., 0., 0., 0., 0.75}, // 8
864 : { -0.125, 0., 0., 0.375, 0., 0., 0., 0., 0., 0.75, 0., 0., 0., 0., 0., 0., 0., 0.}, // 9
865 : { 0., 0., 0., 0., 0., 0., -0.125, 0., 0., 0., 0., 0., 0.375, 0., 0., 0.75, 0., 0.}, // 10
866 : { 0., 0., 0., 0., 0., 0., 0., 0., -0.125, 0., 0., 0., 0., 0., 0.375, 0., 0., 0.75}, // 11
867 : { 0., 0., 0., 0.375, -0.125, 0., 0., 0., 0., 0., 0., 0., 0.75, 0., 0., 0., 0., 0.}, // 12
868 : { 0., 0., 0., 0., -0.125, -0.125, 0., 0., 0., 0., 0., 0., 0.5, 0.25, 0.5, 0., 0., 0.}, // 13
869 : { 0., 0., 0., 0.375, 0., -0.125, 0., 0., 0., 0., 0., 0., 0., 0., 0.75, 0., 0., 0.}, // 14
870 : {-0.046875, 0.015625, 0., 0.140625,-0.046875, 0., -0.09375, 0., 0., 0.28125, -0.09375, 0., 0.28125, 0., 0., 0.5625, 0., 0.}, // 15
871 : { 0., 0.015625, 0.015625, 0.,-0.046875,-0.046875, -0.0625, -0.03125, -0.0625, 0., -0.09375, -0.09375, 0.1875, 0.09375, 0.1875, 0.375, 0.1875, 0.375}, // 16
872 : {-0.046875, 0., 0.015625, 0.140625, 0.,-0.046875, 0., 0., -0.09375, 0.28125, 0., -0.09375, 0., 0., 0.28125, 0., 0., 0.5625} // 17
873 : },
874 :
875 : // embedding matrix for child 5
876 : {
877 : // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
878 : { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0.}, // 0
879 : { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0.}, // 1
880 : { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0.}, // 2
881 : { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0.}, // 3
882 : { 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}, // 4
883 : { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0.}, // 5
884 : { 0., 0., 0., 0., 0., 0., 0., 0., 0., -0.125, 0.375, 0., 0., 0., 0., 0.75, 0., 0.}, // 6
885 : { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.375, -0.125, 0., 0., 0., 0., 0.75, 0.}, // 7
886 : { 0., 0., 0., 0., 0., 0., 0., 0., 0., -0.125, 0., -0.125, 0., 0., 0., 0.5, 0.5, 0.25}, // 8
887 : { 0., 0., 0., 0., 0., 0., -0.125, 0., 0., 0., 0., 0., 0.375, 0., 0., 0.75, 0., 0.}, // 9
888 : { 0., -0.125, 0., 0., 0.375, 0., 0., 0., 0., 0., 0.75, 0., 0., 0., 0., 0., 0., 0.}, // 10
889 : { 0., 0., 0., 0., 0., 0., 0., -0.125, 0., 0., 0., 0., 0., 0.375, 0., 0., 0.75, 0.}, // 11
890 : { 0., 0., 0., -0.125, 0.375, 0., 0., 0., 0., 0., 0., 0., 0.75, 0., 0., 0., 0., 0.}, // 12
891 : { 0., 0., 0., 0., 0.375, -0.125, 0., 0., 0., 0., 0., 0., 0., 0.75, 0., 0., 0., 0.}, // 13
892 : { 0., 0., 0., -0.125, 0., -0.125, 0., 0., 0., 0., 0., 0., 0.5, 0.5, 0.25, 0., 0., 0.}, // 14
893 : { 0.015625,-0.046875, 0.,-0.046875, 0.140625, 0., -0.09375, 0., 0., -0.09375, 0.28125, 0., 0.28125, 0., 0., 0.5625, 0., 0.}, // 15
894 : { 0.,-0.046875, 0.015625, 0., 0.140625,-0.046875, 0., -0.09375, 0., 0., 0.28125, -0.09375, 0., 0.28125, 0., 0., 0.5625, 0.}, // 16
895 : { 0.015625, 0., 0.015625,-0.046875, 0.,-0.046875, -0.0625, -0.0625, -0.03125, -0.09375, 0., -0.09375, 0.1875, 0.1875, 0.09375, 0.375, 0.375, 0.1875} // 17
896 : },
897 :
898 : // embedding matrix for child 6
899 : {
900 : // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
901 : { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1.}, // 0
902 : { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0.}, // 1
903 : { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0.}, // 2
904 : { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0.}, // 3
905 : { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0.}, // 4
906 : { 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}, // 5
907 : { 0., 0., 0., 0., 0., 0., 0., 0., 0., -0.125, -0.125, 0., 0., 0., 0., 0.25, 0.5, 0.5}, // 6
908 : { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., -0.125, 0.375, 0., 0., 0., 0., 0.75, 0.}, // 7
909 : { 0., 0., 0., 0., 0., 0., 0., 0., 0., -0.125, 0., 0.375, 0., 0., 0., 0., 0., 0.75}, // 8
910 : { 0., 0., 0., 0., 0., 0., 0., 0., -0.125, 0., 0., 0., 0., 0., 0.375, 0., 0., 0.75}, // 9
911 : { 0., 0., 0., 0., 0., 0., 0., -0.125, 0., 0., 0., 0., 0., 0.375, 0., 0., 0.75, 0.}, // 10
912 : { 0., 0., -0.125, 0., 0., 0.375, 0., 0., 0., 0., 0., 0.75, 0., 0., 0., 0., 0., 0.}, // 11
913 : { 0., 0., 0., -0.125, -0.125, 0., 0., 0., 0., 0., 0., 0., 0.25, 0.5, 0.5, 0., 0., 0.}, // 12
914 : { 0., 0., 0., 0., -0.125, 0.375, 0., 0., 0., 0., 0., 0., 0., 0.75, 0., 0., 0., 0.}, // 13
915 : { 0., 0., 0., -0.125, 0., 0.375, 0., 0., 0., 0., 0., 0., 0., 0., 0.75, 0., 0., 0.}, // 14
916 : { 0.015625, 0.015625, 0.,-0.046875,-0.046875, 0., -0.03125, -0.0625, -0.0625, -0.09375, -0.09375, 0., 0.09375, 0.1875, 0.1875, 0.1875, 0.375, 0.375}, // 15
917 : { 0., 0.015625,-0.046875, 0.,-0.046875, 0.140625, 0., -0.09375, 0., 0., -0.09375, 0.28125, 0., 0.28125, 0., 0., 0.5625, 0.}, // 16
918 : { 0.015625, 0.,-0.046875,-0.046875, 0., 0.140625, 0., 0., -0.09375, -0.09375, 0., 0.28125, 0., 0., 0.28125, 0., 0., 0.5625} // 17
919 : },
920 :
921 : // embedding matrix for child 7
922 : {
923 : // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
924 : { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0.}, // 0
925 : { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0.}, // 1
926 : { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1.}, // 2
927 : { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0.}, // 3
928 : { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0.}, // 4
929 : { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0.}, // 5
930 : { 0., 0., 0., 0., 0., 0., 0., 0., 0., -0.125, 0., -0.125, 0., 0., 0., 0.5, 0.5, 0.25}, // 6
931 : { 0., 0., 0., 0., 0., 0., 0., 0., 0., -0.125, -0.125, 0., 0., 0., 0., 0.25, 0.5, 0.5}, // 7
932 : { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., -0.125, -0.125, 0., 0., 0., 0.5, 0.25, 0.5}, // 8
933 : { 0., 0., 0., 0., 0., 0., -0.125, 0., 0., 0., 0., 0., 0.375, 0., 0., 0.75, 0., 0.}, // 9
934 : { 0., 0., 0., 0., 0., 0., 0., -0.125, 0., 0., 0., 0., 0., 0.375, 0., 0., 0.75, 0.}, // 10
935 : { 0., 0., 0., 0., 0., 0., 0., 0., -0.125, 0., 0., 0., 0., 0., 0.375, 0., 0., 0.75}, // 11
936 : { 0., 0., 0., -0.125, 0., -0.125, 0., 0., 0., 0., 0., 0., 0.5, 0.5, 0.25, 0., 0., 0.}, // 12
937 : { 0., 0., 0., -0.125, -0.125, 0., 0., 0., 0., 0., 0., 0., 0.25, 0.5, 0.5, 0., 0., 0.}, // 13
938 : { 0., 0., 0., 0., -0.125, -0.125, 0., 0., 0., 0., 0., 0., 0.5, 0.25, 0.5, 0., 0., 0.}, // 14
939 : { 0.015625, 0., 0.015625,-0.046875, 0.,-0.046875, -0.0625, -0.0625, -0.03125, -0.09375, 0., -0.09375, 0.1875, 0.1875, 0.09375, 0.375, 0.375, 0.1875}, // 15
940 : { 0.015625, 0.015625, 0.,-0.046875,-0.046875, 0., -0.03125, -0.0625, -0.0625, -0.09375, -0.09375, 0., 0.09375, 0.1875, 0.1875, 0.1875, 0.375, 0.375}, // 16
941 : { 0., 0.015625, 0.015625, 0.,-0.046875,-0.046875, -0.0625, -0.03125, -0.0625, 0., -0.09375, -0.09375, 0.1875, 0.09375, 0.1875, 0.375, 0.1875, 0.375} // 17
942 : }
943 : };
944 :
945 : #endif
946 :
947 :
948 : void
949 5508 : Prism18::permute(unsigned int perm_num)
950 : {
951 492 : libmesh_assert_less (perm_num, 6);
952 5508 : const unsigned int side = perm_num % 2;
953 5508 : const unsigned int rotate = perm_num / 2;
954 :
955 13194 : for (unsigned int i = 0; i != rotate; ++i)
956 : {
957 7686 : swap3nodes(0,1,2);
958 6996 : swap3nodes(3,4,5);
959 6996 : swap3nodes(6,7,8);
960 6996 : swap3nodes(9,10,11);
961 6996 : swap3nodes(12,13,14);
962 6996 : swap3nodes(15,16,17);
963 6996 : swap3neighbors(1,2,3);
964 : }
965 :
966 5508 : switch (side) {
967 48 : case 0:
968 48 : break;
969 4932 : case 1:
970 4932 : swap2nodes(1,3);
971 4932 : swap2nodes(0,4);
972 4932 : swap2nodes(2,5);
973 4932 : swap2nodes(6,12);
974 4932 : swap2nodes(9,10);
975 4932 : swap2nodes(7,14);
976 4932 : swap2nodes(8,13);
977 4932 : swap2nodes(16,17);
978 444 : swap2neighbors(0,4);
979 444 : swap2neighbors(2,3);
980 444 : break;
981 0 : default:
982 0 : libmesh_error();
983 : }
984 5508 : }
985 :
986 :
987 : void
988 576 : Prism18::flip(BoundaryInfo * boundary_info)
989 : {
990 48 : libmesh_assert(boundary_info);
991 :
992 576 : swap2nodes(0,1);
993 576 : swap2nodes(3,4);
994 576 : swap2nodes(7,8);
995 576 : swap2nodes(9,10);
996 576 : swap2nodes(13,14);
997 576 : swap2nodes(16,17);
998 48 : swap2neighbors(2,3);
999 576 : swap2boundarysides(2,3,boundary_info);
1000 576 : swap2boundaryedges(0,1,boundary_info);
1001 576 : swap2boundaryedges(3,4,boundary_info);
1002 576 : swap2boundaryedges(7,8,boundary_info);
1003 576 : }
1004 :
1005 :
1006 960 : unsigned int Prism18::center_node_on_side(const unsigned short side) const
1007 : {
1008 80 : libmesh_assert_less (side, Prism18::num_sides);
1009 960 : return (side >= 1 && side <= 3) ? side + 14 : invalid_uint;
1010 : }
1011 :
1012 :
1013 : ElemType
1014 2880 : Prism18::side_type (const unsigned int s) const
1015 : {
1016 240 : libmesh_assert_less (s, 5);
1017 2880 : if (s == 0 || s == 4)
1018 1152 : return TRI6;
1019 144 : return QUAD9;
1020 : }
1021 :
1022 :
1023 : } // namespace libMesh
|
