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