libMesh
Loading...
Searching...
No Matches
cell_prism18.C
Go to the documentation of this file.
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
28namespace libMesh
29{
30
31
32
33// ------------------------------------------------------------
34// Prism18 class static member initializations
35const int Prism18::num_nodes;
36const int Prism18::nodes_per_side;
37const int Prism18::nodes_per_edge;
38
39// ------------------------------------------------------------
40// Prism18 class member functions
41
42bool Prism18::is_vertex(const unsigned int i) const
43{
44 if (i < 6)
45 return true;
46 return false;
47}
48
49bool Prism18::is_edge(const unsigned int i) const
50{
51 if (i < 6)
52 return false;
53 if (i > 14)
54 return false;
55 return true;
56}
57
58bool Prism18::is_face(const unsigned int i) const
59{
60 if (i > 14)
61 return true;
62 return false;
63}
64
65bool Prism18::is_node_on_side(const unsigned int n,
66 const unsigned int s) const
67{
68 libmesh_assert_less (s, n_sides());
69 return std::find(std::begin(side_nodes_map[s]),
70 std::end(side_nodes_map[s]),
71 n) != std::end(side_nodes_map[s]);
72}
73
74std::vector<unsigned>
75Prism18::nodes_on_side(const unsigned int s) const
76{
77 libmesh_assert_less(s, n_sides());
78 auto trim = (s > 0 && s < 4) ? 0 : 3;
79 return {std::begin(side_nodes_map[s]), std::end(side_nodes_map[s]) - trim};
80}
81
82std::vector<unsigned>
83Prism18::nodes_on_edge(const unsigned int e) const
84{
85 libmesh_assert_less(e, n_edges());
86 return {std::begin(edge_nodes_map[e]), std::end(edge_nodes_map[e])};
87}
88
89bool Prism18::is_node_on_edge(const unsigned int n,
90 const unsigned int e) const
91{
92 libmesh_assert_less (e, n_edges());
93 return std::find(std::begin(edge_nodes_map[e]),
94 std::end(edge_nodes_map[e]),
95 n) != std::end(edge_nodes_map[e]);
96}
97
98
99
100bool Prism18::has_affine_map() const
101{
102 // Make sure z edges are affine
103 Point v = this->point(3) - this->point(0);
104 if (!v.relative_fuzzy_equals(this->point(4) - this->point(1), affine_tol) ||
105 !v.relative_fuzzy_equals(this->point(5) - this->point(2), affine_tol))
106 return false;
107 // Make sure edges are straight
108 v /= 2;
109 if (!v.relative_fuzzy_equals(this->point(9) - this->point(0), affine_tol) ||
110 !v.relative_fuzzy_equals(this->point(10) - this->point(1), affine_tol) ||
111 !v.relative_fuzzy_equals(this->point(11) - this->point(2), affine_tol) ||
112 !v.relative_fuzzy_equals(this->point(15) - this->point(6), affine_tol) ||
113 !v.relative_fuzzy_equals(this->point(16) - this->point(7), affine_tol) ||
114 !v.relative_fuzzy_equals(this->point(17) - this->point(8), affine_tol))
115 return false;
116 v = (this->point(1) - this->point(0))/2;
117 if (!v.relative_fuzzy_equals(this->point(6) - this->point(0), affine_tol) ||
118 !v.relative_fuzzy_equals(this->point(12) - this->point(3), affine_tol))
119 return false;
120 v = (this->point(2) - this->point(0))/2;
121 if (!v.relative_fuzzy_equals(this->point(8) - this->point(0), affine_tol) ||
122 !v.relative_fuzzy_equals(this->point(14) - this->point(3), affine_tol))
123 return false;
124 v = (this->point(2) - this->point(1))/2;
125 if (!v.relative_fuzzy_equals(this->point(7) - this->point(1), affine_tol) ||
126 !v.relative_fuzzy_equals(this->point(13) - this->point(4), affine_tol))
127 return false;
128 return true;
129}
130
131
132
133Order Prism18::default_order() const
134{
135 return SECOND;
136}
137
138dof_id_type Prism18::key (const unsigned int s) const
139{
140 libmesh_assert_less (s, this->n_sides());
141
142 switch (s)
143 {
144 case 0: // the triangular face at z=0
145 {
146 return Prism::key(0);
147 }
148 case 1: // the quad face at y=0
149 {
150 return Elem::compute_key (this->node_id(15));
151 }
152 case 2: // the other quad face
153 {
154 return Elem::compute_key (this->node_id(16));
155 }
156 case 3: // the quad face at x=0
157 {
158 return Elem::compute_key (this->node_id(17));
159 }
160 case 4: // the triangular face at z=1
161 {
162 return Prism::key(4);
163 }
164 default:
165 libmesh_error_msg("Invalid side " << s);
166 }
167}
168
169
170
171unsigned int Prism18::local_side_node(unsigned int side,
172 unsigned int side_node) const
173{
174 libmesh_assert_less (side, this->n_sides());
175
176 // Never more than 9 nodes per side.
177 libmesh_assert_less(side_node, Prism18::nodes_per_side);
178
179 // Some sides have 6 nodes.
180 libmesh_assert(!(side==0 || side==4) || side_node < 6);
181
182 return Prism18::side_nodes_map[side][side_node];
183}
184
185
186
187unsigned int Prism18::local_edge_node(unsigned int edge,
188 unsigned int edge_node) const
189{
190 libmesh_assert_less(edge, this->n_edges());
191 libmesh_assert_less(edge_node, Prism18::nodes_per_edge);
192
193 return Prism18::edge_nodes_map[edge][edge_node];
194}
195
196
197
198std::unique_ptr<Elem> Prism18::build_side_ptr (const unsigned int i)
199{
200 libmesh_assert_less (i, this->n_sides());
201
202 std::unique_ptr<Elem> face;
203
204 switch (i)
205 {
206 case 0: // the triangular face at z=-1
207 case 4: // the triangular face at z=1
208 {
209 face = std::make_unique<Tri6>();
210 break;
211 }
212 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 face = std::make_unique<Quad9>();
217 break;
218 }
219 default:
220 libmesh_error_msg("Invalid side i = " << i);
221 }
222
223 // Set the nodes
224 for (auto n : face->node_index_range())
225 face->set_node(n, this->node_ptr(Prism18::side_nodes_map[i][n]));
226
227 face->set_interior_parent(this);
228 face->inherit_data_from(*this);
229
230 return face;
231}
232
233
234
235void Prism18::build_side_ptr (std::unique_ptr<Elem> & side,
236 const unsigned int i)
237{
238 libmesh_assert_less (i, this->n_sides());
239
240 switch (i)
241 {
242 case 0: // the triangular face at z=-1
243 case 4: // the triangular face at z=1
244 {
245 if (!side.get() || side->type() != TRI6)
246 {
247 side = this->build_side_ptr(i);
248 return;
249 }
250 break;
251 }
252
253 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 if (!side.get() || side->type() != QUAD9)
258 {
259 side = this->build_side_ptr(i);
260 return;
261 }
262 break;
263 }
264
265 default:
266 libmesh_error_msg("Invalid side i = " << i);
267 }
268
269 side->inherit_data_from(*this);
270
271 // Set the nodes
272 for (auto n : side->node_index_range())
273 side->set_node(n, this->node_ptr(Prism18::side_nodes_map[i][n]));
274}
275
276
277
278std::unique_ptr<Elem> Prism18::build_edge_ptr (const unsigned int i)
279{
280 return this->simple_build_edge_ptr<Edge3,Prism18>(i);
281}
282
283
284
285void Prism18::build_edge_ptr (std::unique_ptr<Elem> & edge, const unsigned int i)
286{
287 this->simple_build_edge_ptr<Prism18>(edge, i, EDGE3);
288}
289
290
291
292void Prism18::connectivity(const unsigned int sc,
293 const IOPackage iop,
294 std::vector<dof_id_type> & conn) const
295{
296 libmesh_assert(_nodes);
297 libmesh_assert_less (sc, this->n_sub_elem());
298 libmesh_assert_not_equal_to (iop, INVALID_IO_PACKAGE);
299
300 switch (iop)
301 {
302 case TECPLOT:
303 {
304 conn.resize(8);
305 switch (sc)
306 {
307
308 case 0:
309 {
310 conn[0] = this->node_id(0)+1;
311 conn[1] = this->node_id(6)+1;
312 conn[2] = this->node_id(8)+1;
313 conn[3] = this->node_id(8)+1;
314 conn[4] = this->node_id(9)+1;
315 conn[5] = this->node_id(15)+1;
316 conn[6] = this->node_id(17)+1;
317 conn[7] = this->node_id(17)+1;
318
319 return;
320 }
321
322 case 1:
323 {
324 conn[0] = this->node_id(6)+1;
325 conn[1] = this->node_id(1)+1;
326 conn[2] = this->node_id(7)+1;
327 conn[3] = this->node_id(7)+1;
328 conn[4] = this->node_id(15)+1;
329 conn[5] = this->node_id(10)+1;
330 conn[6] = this->node_id(16)+1;
331 conn[7] = this->node_id(16)+1;
332
333 return;
334 }
335
336 case 2:
337 {
338 conn[0] = this->node_id(8)+1;
339 conn[1] = this->node_id(7)+1;
340 conn[2] = this->node_id(2)+1;
341 conn[3] = this->node_id(2)+1;
342 conn[4] = this->node_id(17)+1;
343 conn[5] = this->node_id(16)+1;
344 conn[6] = this->node_id(11)+1;
345 conn[7] = this->node_id(11)+1;
346
347 return;
348 }
349
350 case 3:
351 {
352 conn[0] = this->node_id(6)+1;
353 conn[1] = this->node_id(7)+1;
354 conn[2] = this->node_id(8)+1;
355 conn[3] = this->node_id(8)+1;
356 conn[4] = this->node_id(15)+1;
357 conn[5] = this->node_id(16)+1;
358 conn[6] = this->node_id(17)+1;
359 conn[7] = this->node_id(17)+1;
360
361 return;
362 }
363
364 case 4:
365 {
366 conn[0] = this->node_id(9)+1;
367 conn[1] = this->node_id(15)+1;
368 conn[2] = this->node_id(17)+1;
369 conn[3] = this->node_id(17)+1;
370 conn[4] = this->node_id(3)+1;
371 conn[5] = this->node_id(12)+1;
372 conn[6] = this->node_id(14)+1;
373 conn[7] = this->node_id(14)+1;
374
375 return;
376 }
377
378 case 5:
379 {
380 conn[0] = this->node_id(15)+1;
381 conn[1] = this->node_id(10)+1;
382 conn[2] = this->node_id(16)+1;
383 conn[3] = this->node_id(16)+1;
384 conn[4] = this->node_id(12)+1;
385 conn[5] = this->node_id(4)+1;
386 conn[6] = this->node_id(13)+1;
387 conn[7] = this->node_id(13)+1;
388
389 return;
390 }
391
392 case 6:
393 {
394 conn[0] = this->node_id(17)+1;
395 conn[1] = this->node_id(16)+1;
396 conn[2] = this->node_id(11)+1;
397 conn[3] = this->node_id(11)+1;
398 conn[4] = this->node_id(14)+1;
399 conn[5] = this->node_id(13)+1;
400 conn[6] = this->node_id(5)+1;
401 conn[7] = this->node_id(5)+1;
402
403 return;
404 }
405
406 case 7:
407 {
408 conn[0] = this->node_id(15)+1;
409 conn[1] = this->node_id(16)+1;
410 conn[2] = this->node_id(17)+1;
411 conn[3] = this->node_id(17)+1;
412 conn[4] = this->node_id(12)+1;
413 conn[5] = this->node_id(13)+1;
414 conn[6] = this->node_id(14)+1;
415 conn[7] = this->node_id(14)+1;
416
417 return;
418 }
419
420 default:
421 libmesh_error_msg("Invalid sc = " << sc);
422 }
423
424 }
425
426 case VTK:
427 {
428 // VTK now supports VTK_BIQUADRATIC_QUADRATIC_WEDGE directly
429 const unsigned int conn_size = 18;
430 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 for (auto i : index_range(conn))
436 conn[i] = this->node_id(i);
437
438 // top "ring" of mid-edge nodes
439 conn[9] = this->node_id(12);
440 conn[10] = this->node_id(13);
441 conn[11] = this->node_id(14);
442
443 // middle "ring" of mid-edge nodes
444 conn[12] = this->node_id(9);
445 conn[13] = this->node_id(10);
446 conn[14] = this->node_id(11);
447
448 return;
449 }
450
451 default:
452 libmesh_error_msg("Unsupported IO package " << iop);
453 }
454}
455
456
457
458
459unsigned int Prism18::n_second_order_adjacent_vertices (const unsigned int n) const
460{
461 switch (n)
462 {
463 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 return 2;
473
474 case 15:
475 case 16:
476 case 17:
477 return 4;
478
479 default:
480 libmesh_error_msg("Invalid node n = " << n);
481 }
482}
483
484
485
486
487
488unsigned short int Prism18::second_order_adjacent_vertex (const unsigned int n,
489 const unsigned int v) const
490{
491 libmesh_assert_greater_equal (n, this->n_vertices());
492 libmesh_assert_less (n, this->n_nodes());
493
494 switch (n)
495 {
496 /*
497 * These nodes are unique to \p Prism18,
498 * let our _remaining_... matrix handle
499 * this.
500 */
501 case 15:
502 case 16:
503 case 17:
504 {
505 libmesh_assert_less (v, 4);
506 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 default:
516 {
517 libmesh_assert_less (v, 2);
518 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
529const 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
538std::pair<unsigned short int, unsigned short int>
539Prism18::second_order_child_vertex (const unsigned int n) const
540{
541 libmesh_assert_greater_equal (n, this->n_vertices());
542 libmesh_assert_less (n, this->n_nodes());
543
544 return std::pair<unsigned short int, unsigned short int>
545 (_second_order_vertex_child_number[n],
546 _second_order_vertex_child_index[n]);
547}
548
549
550
551Real Prism18::volume () const
552{
553 // This specialization is good for Lagrange mappings only in general
554 if (this->mapping_type() != LAGRANGE_MAP)
555 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 x0 = point(0), x1 = point(1), x2 = point(2), x3 = point(3), x4 = point(4), x5 = point(5),
561 x6 = point(6), x7 = point(7), x8 = point(8), x9 = point(9), x10 = point(10), x11 = point(11),
562 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 const int n_components = 16;
566
567 // Terms are copied directly from a Python script.
568 Point dx_dxi[n_components] =
569 {
570 -x10 + 4*x15 - 3*x9,
571 3*x0/2 + x1/2 + 2*x12 - 3*x3/2 - x4/2 - 2*x6,
572 -3*x0/2 - x1/2 + x10 + 2*x12 - 4*x15 - 3*x3/2 - x4/2 + 2*x6 + 3*x9,
573 -4*x15 + 4*x16 - 4*x17 + 4*x9,
574 -2*x0 - 2*x12 + 2*x13 - 2*x14 + 2*x3 + 2*x6 - 2*x7 + 2*x8,
575 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 4*x10 - 8*x15 + 4*x9,
579 -2*x0 - 2*x1 - 4*x12 + 2*x3 + 2*x4 + 4*x6,
580 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 };
587
588 Point dx_deta[n_components] =
589 {
590 -x11 + 4*x17 - 3*x9,
591 3*x0/2 + 2*x14 + x2/2 - 3*x3/2 - x5/2 - 2*x8,
592 -3*x0/2 + x11 + 2*x14 - 4*x17 - x2/2 - 3*x3/2 - x5/2 + 2*x8 + 3*x9,
593 4*x11 - 8*x17 + 4*x9,
594 -2*x0 - 4*x14 - 2*x2 + 2*x3 + 2*x5 + 4*x8,
595 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 -4*x15 + 4*x16 - 4*x17 + 4*x9,
599 -2*x0 - 2*x12 + 2*x13 - 2*x14 + 2*x3 + 2*x6 - 2*x7 + 2*x8,
600 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 };
607
608 Point dx_dzeta[n_components] =
609 {
610 -x0/2 + x3/2,
611 x0 + x3 - 2*x9,
612 Point(0,0,0),
613 3*x0/2 + 2*x14 + x2/2 - 3*x3/2 - x5/2 - 2*x8,
614 -3*x0 + 2*x11 + 4*x14 - 8*x17 - x2 - 3*x3 - x5 + 4*x8 + 6*x9,
615 Point(0,0,0),
616 -x0 - 2*x14 - x2 + x3 + x5 + 2*x8,
617 2*x0 - 4*x11 - 4*x14 + 8*x17 + 2*x2 + 2*x3 + 2*x5 - 4*x8 - 4*x9,
618 3*x0/2 + x1/2 + 2*x12 - 3*x3/2 - x4/2 - 2*x6,
619 -3*x0 - x1 + 2*x10 + 4*x12 - 8*x15 - 3*x3 - x4 + 4*x6 + 6*x9,
620 Point(0,0,0),
621 -2*x0 - 2*x12 + 2*x13 - 2*x14 + 2*x3 + 2*x6 - 2*x7 + 2*x8,
622 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 -x0 - x1 - 2*x12 + x3 + x4 + 2*x6,
625 2*x0 + 2*x1 - 4*x10 - 4*x12 + 8*x15 + 2*x3 + 2*x4 - 4*x6 - 4*x9
626 };
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 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 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 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 Real vol = 0.;
707 for (int i=0; i<N2D; ++i)
708 for (int j=0; j<N1D; ++j)
709 {
710 // Compute dx_dxi, dx_deta, dx_dzeta at the current quadrature point.
711 Point dx_dxi_q, dx_deta_q, dx_dzeta_q;
712 for (int c=0; c<n_components; ++c)
713 {
714 Real coeff =
715 xi[i][exponents[c][0]]*
716 eta[i][exponents[c][1]]*
717 zeta[j][exponents[c][2]];
718
719 dx_dxi_q += coeff * dx_dxi[c];
720 dx_deta_q += coeff * dx_deta[c];
721 dx_dzeta_q += coeff * dx_dzeta[c];
722 }
723
724 // Compute scalar triple product, multiply by weight, and accumulate volume.
725 vol += w2D[i] * w1D[j] * triple_product(dx_dxi_q, dx_deta_q, dx_dzeta_q);
726 }
727
728 return vol;
729}
730
731
732
733
734#ifdef LIBMESH_ENABLE_AMR
735
736const 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
926void
927Prism18::permute(unsigned int perm_num)
928{
929 libmesh_assert_less (perm_num, 6);
930 const unsigned int side = perm_num % 2;
931 const unsigned int rotate = perm_num / 2;
932
933 for (unsigned int i = 0; i != rotate; ++i)
934 {
935 swap3nodes(0,1,2);
936 swap3nodes(3,4,5);
937 swap3nodes(6,7,8);
938 swap3nodes(9,10,11);
939 swap3nodes(12,13,14);
940 swap3nodes(15,16,17);
941 swap3neighbors(1,2,3);
942 }
943
944 switch (side) {
945 case 0:
946 break;
947 case 1:
948 swap2nodes(1,3);
949 swap2nodes(0,4);
950 swap2nodes(2,5);
951 swap2nodes(6,12);
952 swap2nodes(9,10);
953 swap2nodes(7,14);
954 swap2nodes(8,13);
955 swap2nodes(16,17);
956 swap2neighbors(0,4);
957 swap2neighbors(2,3);
958 break;
959 default:
960 libmesh_error();
961 }
962}
963
964
965void
966Prism18::flip(BoundaryInfo * boundary_info)
967{
968 libmesh_assert(boundary_info);
969
970 swap2nodes(0,1);
971 swap2nodes(3,4);
972 swap2nodes(7,8);
973 swap2nodes(9,10);
974 swap2nodes(13,14);
975 swap2nodes(16,17);
976 swap2neighbors(2,3);
977 swap2boundarysides(2,3,boundary_info);
978 swap2boundaryedges(0,1,boundary_info);
979 swap2boundaryedges(3,4,boundary_info);
980 swap2boundaryedges(7,8,boundary_info);
981}
982
983
984unsigned int Prism18::center_node_on_side(const unsigned short side) const
985{
986 libmesh_assert_less (side, Prism18::num_sides);
987 return (side >= 1 && side <= 3) ? side + 14 : invalid_uint;
988}
989
990
991ElemType
992Prism18::side_type (const unsigned int s) const
993{
994 libmesh_assert_less (s, 5);
995 if (s == 0 || s == 4)
996 return TRI6;
997 return QUAD9;
998}
999
1000
1001} // namespace libMesh
libMesh::BoundaryInfo
The BoundaryInfo class contains information relevant to boundary conditions including storing faces,...
Definition boundary_info.h:58
libMesh::Elem::affine_tol
static constexpr Real affine_tol
Default tolerance to use in has_affine_map().
Definition elem.h:2191
libMesh::Elem::point
const Point & point(const unsigned int i) const
Definition elem.h:2583
libMesh::Elem::swap2nodes
void swap2nodes(unsigned int n1, unsigned int n2)
Swaps two node_ptrs.
Definition elem.h:2222
libMesh::Elem::swap3neighbors
void swap3neighbors(unsigned int n1, unsigned int n2, unsigned int n3)
Swaps three neighbor_ptrs, "rotating" them.
Definition elem.h:2263
libMesh::Elem::_nodes
Node ** _nodes
Pointers to the nodes we are connected to.
Definition elem.h:2375
libMesh::Elem::swap2neighbors
void swap2neighbors(unsigned int n1, unsigned int n2)
Swaps two neighbor_ptrs.
Definition elem.h:2232
libMesh::Elem::swap2boundarysides
void swap2boundarysides(unsigned short s1, unsigned short s2, BoundaryInfo *boundary_info) const
Swaps two sides in boundary_info, if it is non-null.
Definition elem.C:3982
libMesh::Elem::mapping_type
ElemMappingType mapping_type() const
Definition elem.h:3258
libMesh::Elem::compute_key
static dof_id_type compute_key(dof_id_type n0)
Definition elem.h:3432
libMesh::Elem::swap3nodes
void swap3nodes(unsigned int n1, unsigned int n2, unsigned int n3)
Swaps three node_ptrs, "rotating" them.
Definition elem.h:2254
libMesh::Elem::volume
virtual Real volume() const
Definition elem.C:3865
libMesh::Elem::swap2boundaryedges
void swap2boundaryedges(unsigned short e1, unsigned short e2, BoundaryInfo *boundary_info) const
Swaps two edges in boundary_info, if it is non-null.
Definition elem.C:3998
libMesh::Elem::node_id
dof_id_type node_id(const unsigned int i) const
Definition elem.h:2605
libMesh::Point
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition point.h:41
libMesh::Prism18::is_vertex
virtual bool is_vertex(const unsigned int i) const override
Definition cell_prism18.C:42
libMesh::Prism18::_remaining_second_order_adjacent_vertices
static const unsigned short int _remaining_second_order_adjacent_vertices[3][4]
Matrix that tells which vertices define the location of mid-side (or second-order) nodes.
Definition cell_prism18.h:305
libMesh::Prism18::second_order_child_vertex
virtual std::pair< unsigned short int, unsigned short int > second_order_child_vertex(const unsigned int n) const override
Definition cell_prism18.C:539
libMesh::Prism18::connectivity
virtual void connectivity(const unsigned int sc, const IOPackage iop, std::vector< dof_id_type > &conn) const override
Definition cell_prism18.C:292
libMesh::Prism18::default_order
virtual Order default_order() const override
Definition cell_prism18.C:133
libMesh::Prism18::local_side_node
virtual unsigned int local_side_node(unsigned int side, unsigned int side_node) const override
Definition cell_prism18.C:171
libMesh::Prism18::is_face
virtual bool is_face(const unsigned int i) const override
Definition cell_prism18.C:58
libMesh::Prism18::_embedding_matrix
static const Real _embedding_matrix[num_children][num_nodes][num_nodes]
Matrix that computes new nodal locations/solution values from current nodes/solution.
Definition cell_prism18.h:291
libMesh::Prism18::local_edge_node
virtual unsigned int local_edge_node(unsigned int edge, unsigned int edge_node) const override
Definition cell_prism18.C:187
libMesh::Prism18::edge_nodes_map
static constexpr const unsigned int(& edge_nodes_map)[num_edges][nodes_per_edge]
Definition cell_prism18.h:248
libMesh::Prism18::permute
virtual void permute(unsigned int perm_num) override final
Permutes the element (by swapping node and neighbor pointers) according to the specified index.
Definition cell_prism18.C:927
libMesh::Prism18::is_edge
virtual bool is_edge(const unsigned int i) const override
Definition cell_prism18.C:49
libMesh::Prism18::n_nodes
virtual unsigned int n_nodes() const override
Definition cell_prism18.h:103
libMesh::Prism18::center_node_on_side
unsigned int center_node_on_side(const unsigned short side) const override final
Definition cell_prism18.C:984
libMesh::Prism18::volume
virtual Real volume() const override
A specialization for computing the volume of a Prism18.
Definition cell_prism18.C:551
libMesh::Prism18::has_affine_map
virtual bool has_affine_map() const override
Definition cell_prism18.C:100
libMesh::Prism18::key
virtual dof_id_type key() const
Don't hide Elem::key() defined in the base class.
Definition elem.C:678
libMesh::Prism18::n_sub_elem
virtual unsigned int n_sub_elem() const override
Definition cell_prism18.h:108
libMesh::Prism18::n_second_order_adjacent_vertices
virtual unsigned int n_second_order_adjacent_vertices(const unsigned int) const override
Definition cell_prism18.C:459
libMesh::Prism18::build_side_ptr
virtual std::unique_ptr< Elem > build_side_ptr(const unsigned int i) override
Builds a QUAD9 or TRI6 built coincident with face i.
Definition cell_prism18.C:198
libMesh::Prism18::nodes_on_side
virtual std::vector< unsigned int > nodes_on_side(const unsigned int s) const override
Definition cell_prism18.C:75
libMesh::Prism18::is_node_on_side
virtual bool is_node_on_side(const unsigned int n, const unsigned int s) const override
Definition cell_prism18.C:65
libMesh::Prism18::second_order_adjacent_vertex
virtual unsigned short int second_order_adjacent_vertex(const unsigned int n, const unsigned int v) const override
Definition cell_prism18.C:488
libMesh::Prism18::nodes_per_side
static const int nodes_per_side
Definition cell_prism18.h:238
libMesh::Prism18::is_node_on_edge
virtual bool is_node_on_edge(const unsigned int n, const unsigned int e) const override
Definition cell_prism18.C:89
libMesh::Prism18::nodes_per_edge
static const int nodes_per_edge
Definition cell_prism18.h:239
libMesh::Prism18::nodes_on_edge
virtual std::vector< unsigned int > nodes_on_edge(const unsigned int e) const override
Definition cell_prism18.C:83
libMesh::Prism18::flip
virtual void flip(BoundaryInfo *) override final
Flips the element (by swapping node and neighbor pointers) to have a mapping Jacobian of opposite sig...
Definition cell_prism18.C:966
libMesh::Prism18::num_nodes
static const int num_nodes
Geometric constants for Prism18.
Definition cell_prism18.h:237
libMesh::Prism18::side_type
ElemType side_type(const unsigned int s) const override final
Definition cell_prism18.C:992
libMesh::Prism18::side_nodes_map
static constexpr const unsigned int(& side_nodes_map)[num_sides][nodes_per_side]
Definition cell_prism18.h:253
libMesh::Prism18::build_edge_ptr
virtual std::unique_ptr< Elem > build_edge_ptr(const unsigned int i) override
Builds a EDGE3 or INFEDGE2 built coincident with edge i.
Definition cell_prism18.C:278
libMesh::Prism::n_vertices
virtual unsigned int n_vertices() const override final
Definition cell_prism.h:91
libMesh::Prism::n_edges
virtual unsigned int n_edges() const override final
Definition cell_prism.h:96
libMesh::Prism::num_children
static const int num_children
Definition cell_prism.h:75
libMesh::Prism::n_sides
virtual unsigned int n_sides() const override final
Definition cell_prism.h:86
libMesh::Prism::_second_order_vertex_child_number
static const unsigned short int _second_order_vertex_child_number[18]
Vector that names a child sharing each second order node.
Definition cell_prism.h:213
libMesh::Prism::key
virtual dof_id_type key() const
Don't hide Elem::key() defined in the base class.
Definition elem.C:678
libMesh::Prism::num_sides
static const int num_sides
Geometric constants for all Prisms.
Definition cell_prism.h:73
libMesh::Prism::_second_order_adjacent_vertices
static const unsigned short int _second_order_adjacent_vertices[9][2]
Matrix that tells which vertices define the location of mid-side (or second-order) nodes.
Definition cell_prism.h:208
libMesh::Prism::_second_order_vertex_child_index
static const unsigned short int _second_order_vertex_child_index[18]
Vector that names the child vertex index for each second order node.
Definition cell_prism.h:218
libMesh::TypeVector::relative_fuzzy_equals
bool relative_fuzzy_equals(const TypeVector< T > &rhs, Real tol=TOLERANCE) const
Definition type_vector.h:1052
libMesh
The libMesh namespace provides an interface to certain functionality in the library.
libMesh::index_range
auto index_range(const T &sizable)
Helper function that returns an IntRange<std::size_t> representing all the indices of the passed-in v...
Definition int_range.h:153
libMesh::IOPackage
IOPackage
libMesh interfaces with several different software packages for the purposes of creating,...
Definition enum_io_package.h:38
libMesh::triple_product
LIBMESH_DEVICE_INLINE T triple_product(const TypeVector< T > &a, const TypeVector< T > &b, const TypeVector< T > &c)
Definition type_vector.h:1100
libMesh::ElemType
ElemType
Defines an enum for geometric element types.
Definition enum_elem_type.h:33
libMesh::libmesh_assert
libmesh_assert(ctx)
libMesh::invalid_uint
const unsigned int invalid_uint
A number which is used quite often to represent an invalid or uninitialized value for an unsigned int...
Definition libmesh.h:303
libMesh::dof_id_type
uint8_t dof_id_type
Definition id_types.h:67
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition libmesh_common.h:148
Generated on Thu Nov 12 2026 20:15:22 for libMesh by doxygen 1.9.8