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
39const 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
48const 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
64bool Prism18::is_vertex(const unsigned int i) const
65{
66 if (i < 6)
67 return true;
68 return false;
69}
70
71bool Prism18::is_edge(const unsigned int i) const
72{
73 if (i < 6)
74 return false;
75 if (i > 14)
76 return false;
77 return true;
78}
79
80bool Prism18::is_face(const unsigned int i) const
81{
82 if (i > 14)
83 return true;
84 return false;
85}
86
87bool Prism18::is_node_on_side(const unsigned int n,
88 const unsigned int s) const
89{
90 libmesh_assert_less (s, n_sides());
91 return std::find(std::begin(side_nodes_map[s]),
92 std::end(side_nodes_map[s]),
93 n) != std::end(side_nodes_map[s]);
94}
95
96std::vector<unsigned>
97Prism18::nodes_on_side(const unsigned int s) const
98{
99 libmesh_assert_less(s, n_sides());
100 auto trim = (s > 0 && s < 4) ? 0 : 3;
101 return {std::begin(side_nodes_map[s]), std::end(side_nodes_map[s]) - trim};
102}
103
104std::vector<unsigned>
105Prism18::nodes_on_edge(const unsigned int e) const
106{
107 libmesh_assert_less(e, n_edges());
108 return {std::begin(edge_nodes_map[e]), std::end(edge_nodes_map[e])};
109}
110
111bool Prism18::is_node_on_edge(const unsigned int n,
112 const unsigned int e) const
113{
114 libmesh_assert_less (e, n_edges());
115 return std::find(std::begin(edge_nodes_map[e]),
116 std::end(edge_nodes_map[e]),
117 n) != std::end(edge_nodes_map[e]);
118}
119
120
121
122bool Prism18::has_affine_map() const
123{
124 // Make sure z edges are affine
125 Point v = this->point(3) - this->point(0);
126 if (!v.relative_fuzzy_equals(this->point(4) - this->point(1), affine_tol) ||
127 !v.relative_fuzzy_equals(this->point(5) - this->point(2), affine_tol))
128 return false;
129 // Make sure edges are straight
130 v /= 2;
131 if (!v.relative_fuzzy_equals(this->point(9) - this->point(0), affine_tol) ||
132 !v.relative_fuzzy_equals(this->point(10) - this->point(1), affine_tol) ||
133 !v.relative_fuzzy_equals(this->point(11) - this->point(2), affine_tol) ||
134 !v.relative_fuzzy_equals(this->point(15) - this->point(6), affine_tol) ||
135 !v.relative_fuzzy_equals(this->point(16) - this->point(7), affine_tol) ||
136 !v.relative_fuzzy_equals(this->point(17) - this->point(8), affine_tol))
137 return false;
138 v = (this->point(1) - this->point(0))/2;
139 if (!v.relative_fuzzy_equals(this->point(6) - this->point(0), affine_tol) ||
140 !v.relative_fuzzy_equals(this->point(12) - this->point(3), affine_tol))
141 return false;
142 v = (this->point(2) - this->point(0))/2;
143 if (!v.relative_fuzzy_equals(this->point(8) - this->point(0), affine_tol) ||
144 !v.relative_fuzzy_equals(this->point(14) - this->point(3), affine_tol))
145 return false;
146 v = (this->point(2) - this->point(1))/2;
147 if (!v.relative_fuzzy_equals(this->point(7) - this->point(1), affine_tol) ||
148 !v.relative_fuzzy_equals(this->point(13) - this->point(4), affine_tol))
149 return false;
150 return true;
151}
152
153
154
155Order Prism18::default_order() const
156{
157 return SECOND;
158}
159
160dof_id_type Prism18::key (const unsigned int s) const
161{
162 libmesh_assert_less (s, this->n_sides());
163
164 switch (s)
165 {
166 case 0: // the triangular face at z=0
167 {
168 return Prism::key(0);
169 }
170 case 1: // the quad face at y=0
171 {
172 return Elem::compute_key (this->node_id(15));
173 }
174 case 2: // the other quad face
175 {
176 return Elem::compute_key (this->node_id(16));
177 }
178 case 3: // the quad face at x=0
179 {
180 return Elem::compute_key (this->node_id(17));
181 }
182 case 4: // the triangular face at z=1
183 {
184 return Prism::key(4);
185 }
186 default:
187 libmesh_error_msg("Invalid side " << s);
188 }
189}
190
191
192
193unsigned int Prism18::local_side_node(unsigned int side,
194 unsigned int side_node) const
195{
196 libmesh_assert_less (side, this->n_sides());
197
198 // Never more than 9 nodes per side.
199 libmesh_assert_less(side_node, Prism18::nodes_per_side);
200
201 // Some sides have 6 nodes.
202 libmesh_assert(!(side==0 || side==4) || side_node < 6);
203
204 return Prism18::side_nodes_map[side][side_node];
205}
206
207
208
209unsigned int Prism18::local_edge_node(unsigned int edge,
210 unsigned int edge_node) const
211{
212 libmesh_assert_less(edge, this->n_edges());
213 libmesh_assert_less(edge_node, Prism18::nodes_per_edge);
214
215 return Prism18::edge_nodes_map[edge][edge_node];
216}
217
218
219
220std::unique_ptr<Elem> Prism18::build_side_ptr (const unsigned int i)
221{
222 libmesh_assert_less (i, this->n_sides());
223
224 std::unique_ptr<Elem> face;
225
226 switch (i)
227 {
228 case 0: // the triangular face at z=-1
229 case 4: // the triangular face at z=1
230 {
231 face = std::make_unique<Tri6>();
232 break;
233 }
234 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 face = std::make_unique<Quad9>();
239 break;
240 }
241 default:
242 libmesh_error_msg("Invalid side i = " << i);
243 }
244
245 // Set the nodes
246 for (auto n : face->node_index_range())
247 face->set_node(n, this->node_ptr(Prism18::side_nodes_map[i][n]));
248
249 face->set_interior_parent(this);
250 face->inherit_data_from(*this);
251
252 return face;
253}
254
255
256
257void Prism18::build_side_ptr (std::unique_ptr<Elem> & side,
258 const unsigned int i)
259{
260 libmesh_assert_less (i, this->n_sides());
261
262 switch (i)
263 {
264 case 0: // the triangular face at z=-1
265 case 4: // the triangular face at z=1
266 {
267 if (!side.get() || side->type() != TRI6)
268 {
269 side = this->build_side_ptr(i);
270 return;
271 }
272 break;
273 }
274
275 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 if (!side.get() || side->type() != QUAD9)
280 {
281 side = this->build_side_ptr(i);
282 return;
283 }
284 break;
285 }
286
287 default:
288 libmesh_error_msg("Invalid side i = " << i);
289 }
290
291 side->inherit_data_from(*this);
292
293 // Set the nodes
294 for (auto n : side->node_index_range())
295 side->set_node(n, this->node_ptr(Prism18::side_nodes_map[i][n]));
296}
297
298
299
300std::unique_ptr<Elem> Prism18::build_edge_ptr (const unsigned int i)
301{
302 return this->simple_build_edge_ptr<Edge3,Prism18>(i);
303}
304
305
306
307void Prism18::build_edge_ptr (std::unique_ptr<Elem> & edge, const unsigned int i)
308{
309 this->simple_build_edge_ptr<Prism18>(edge, i, EDGE3);
310}
311
312
313
314void Prism18::connectivity(const unsigned int sc,
315 const IOPackage iop,
316 std::vector<dof_id_type> & conn) const
317{
318 libmesh_assert(_nodes);
319 libmesh_assert_less (sc, this->n_sub_elem());
320 libmesh_assert_not_equal_to (iop, INVALID_IO_PACKAGE);
321
322 switch (iop)
323 {
324 case TECPLOT:
325 {
326 conn.resize(8);
327 switch (sc)
328 {
329
330 case 0:
331 {
332 conn[0] = this->node_id(0)+1;
333 conn[1] = this->node_id(6)+1;
334 conn[2] = this->node_id(8)+1;
335 conn[3] = this->node_id(8)+1;
336 conn[4] = this->node_id(9)+1;
337 conn[5] = this->node_id(15)+1;
338 conn[6] = this->node_id(17)+1;
339 conn[7] = this->node_id(17)+1;
340
341 return;
342 }
343
344 case 1:
345 {
346 conn[0] = this->node_id(6)+1;
347 conn[1] = this->node_id(1)+1;
348 conn[2] = this->node_id(7)+1;
349 conn[3] = this->node_id(7)+1;
350 conn[4] = this->node_id(15)+1;
351 conn[5] = this->node_id(10)+1;
352 conn[6] = this->node_id(16)+1;
353 conn[7] = this->node_id(16)+1;
354
355 return;
356 }
357
358 case 2:
359 {
360 conn[0] = this->node_id(8)+1;
361 conn[1] = this->node_id(7)+1;
362 conn[2] = this->node_id(2)+1;
363 conn[3] = this->node_id(2)+1;
364 conn[4] = this->node_id(17)+1;
365 conn[5] = this->node_id(16)+1;
366 conn[6] = this->node_id(11)+1;
367 conn[7] = this->node_id(11)+1;
368
369 return;
370 }
371
372 case 3:
373 {
374 conn[0] = this->node_id(6)+1;
375 conn[1] = this->node_id(7)+1;
376 conn[2] = this->node_id(8)+1;
377 conn[3] = this->node_id(8)+1;
378 conn[4] = this->node_id(15)+1;
379 conn[5] = this->node_id(16)+1;
380 conn[6] = this->node_id(17)+1;
381 conn[7] = this->node_id(17)+1;
382
383 return;
384 }
385
386 case 4:
387 {
388 conn[0] = this->node_id(9)+1;
389 conn[1] = this->node_id(15)+1;
390 conn[2] = this->node_id(17)+1;
391 conn[3] = this->node_id(17)+1;
392 conn[4] = this->node_id(3)+1;
393 conn[5] = this->node_id(12)+1;
394 conn[6] = this->node_id(14)+1;
395 conn[7] = this->node_id(14)+1;
396
397 return;
398 }
399
400 case 5:
401 {
402 conn[0] = this->node_id(15)+1;
403 conn[1] = this->node_id(10)+1;
404 conn[2] = this->node_id(16)+1;
405 conn[3] = this->node_id(16)+1;
406 conn[4] = this->node_id(12)+1;
407 conn[5] = this->node_id(4)+1;
408 conn[6] = this->node_id(13)+1;
409 conn[7] = this->node_id(13)+1;
410
411 return;
412 }
413
414 case 6:
415 {
416 conn[0] = this->node_id(17)+1;
417 conn[1] = this->node_id(16)+1;
418 conn[2] = this->node_id(11)+1;
419 conn[3] = this->node_id(11)+1;
420 conn[4] = this->node_id(14)+1;
421 conn[5] = this->node_id(13)+1;
422 conn[6] = this->node_id(5)+1;
423 conn[7] = this->node_id(5)+1;
424
425 return;
426 }
427
428 case 7:
429 {
430 conn[0] = this->node_id(15)+1;
431 conn[1] = this->node_id(16)+1;
432 conn[2] = this->node_id(17)+1;
433 conn[3] = this->node_id(17)+1;
434 conn[4] = this->node_id(12)+1;
435 conn[5] = this->node_id(13)+1;
436 conn[6] = this->node_id(14)+1;
437 conn[7] = this->node_id(14)+1;
438
439 return;
440 }
441
442 default:
443 libmesh_error_msg("Invalid sc = " << sc);
444 }
445
446 }
447
448 case VTK:
449 {
450 // VTK now supports VTK_BIQUADRATIC_QUADRATIC_WEDGE directly
451 const unsigned int conn_size = 18;
452 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 for (auto i : index_range(conn))
458 conn[i] = this->node_id(i);
459
460 // top "ring" of mid-edge nodes
461 conn[9] = this->node_id(12);
462 conn[10] = this->node_id(13);
463 conn[11] = this->node_id(14);
464
465 // middle "ring" of mid-edge nodes
466 conn[12] = this->node_id(9);
467 conn[13] = this->node_id(10);
468 conn[14] = this->node_id(11);
469
470 return;
471 }
472
473 default:
474 libmesh_error_msg("Unsupported IO package " << iop);
475 }
476}
477
478
479
480
481unsigned int Prism18::n_second_order_adjacent_vertices (const unsigned int n) const
482{
483 switch (n)
484 {
485 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 return 2;
495
496 case 15:
497 case 16:
498 case 17:
499 return 4;
500
501 default:
502 libmesh_error_msg("Invalid node n = " << n);
503 }
504}
505
506
507
508
509
510unsigned short int Prism18::second_order_adjacent_vertex (const unsigned int n,
511 const unsigned int v) const
512{
513 libmesh_assert_greater_equal (n, this->n_vertices());
514 libmesh_assert_less (n, this->n_nodes());
515
516 switch (n)
517 {
518 /*
519 * These nodes are unique to \p Prism18,
520 * let our _remaining_... matrix handle
521 * this.
522 */
523 case 15:
524 case 16:
525 case 17:
526 {
527 libmesh_assert_less (v, 4);
528 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 default:
538 {
539 libmesh_assert_less (v, 2);
540 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
551const 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
560std::pair<unsigned short int, unsigned short int>
561Prism18::second_order_child_vertex (const unsigned int n) const
562{
563 libmesh_assert_greater_equal (n, this->n_vertices());
564 libmesh_assert_less (n, this->n_nodes());
565
566 return std::pair<unsigned short int, unsigned short int>
567 (_second_order_vertex_child_number[n],
568 _second_order_vertex_child_index[n]);
569}
570
571
572
573Real Prism18::volume () const
574{
575 // This specialization is good for Lagrange mappings only in general
576 if (this->mapping_type() != LAGRANGE_MAP)
577 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 x0 = point(0), x1 = point(1), x2 = point(2), x3 = point(3), x4 = point(4), x5 = point(5),
583 x6 = point(6), x7 = point(7), x8 = point(8), x9 = point(9), x10 = point(10), x11 = point(11),
584 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 const int n_components = 16;
588
589 // Terms are copied directly from a Python script.
590 Point dx_dxi[n_components] =
591 {
592 -x10 + 4*x15 - 3*x9,
593 3*x0/2 + x1/2 + 2*x12 - 3*x3/2 - x4/2 - 2*x6,
594 -3*x0/2 - x1/2 + x10 + 2*x12 - 4*x15 - 3*x3/2 - x4/2 + 2*x6 + 3*x9,
595 -4*x15 + 4*x16 - 4*x17 + 4*x9,
596 -2*x0 - 2*x12 + 2*x13 - 2*x14 + 2*x3 + 2*x6 - 2*x7 + 2*x8,
597 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 4*x10 - 8*x15 + 4*x9,
601 -2*x0 - 2*x1 - 4*x12 + 2*x3 + 2*x4 + 4*x6,
602 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 };
609
610 Point dx_deta[n_components] =
611 {
612 -x11 + 4*x17 - 3*x9,
613 3*x0/2 + 2*x14 + x2/2 - 3*x3/2 - x5/2 - 2*x8,
614 -3*x0/2 + x11 + 2*x14 - 4*x17 - x2/2 - 3*x3/2 - x5/2 + 2*x8 + 3*x9,
615 4*x11 - 8*x17 + 4*x9,
616 -2*x0 - 4*x14 - 2*x2 + 2*x3 + 2*x5 + 4*x8,
617 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 -4*x15 + 4*x16 - 4*x17 + 4*x9,
621 -2*x0 - 2*x12 + 2*x13 - 2*x14 + 2*x3 + 2*x6 - 2*x7 + 2*x8,
622 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 };
629
630 Point dx_dzeta[n_components] =
631 {
632 -x0/2 + x3/2,
633 x0 + x3 - 2*x9,
634 Point(0,0,0),
635 3*x0/2 + 2*x14 + x2/2 - 3*x3/2 - x5/2 - 2*x8,
636 -3*x0 + 2*x11 + 4*x14 - 8*x17 - x2 - 3*x3 - x5 + 4*x8 + 6*x9,
637 Point(0,0,0),
638 -x0 - 2*x14 - x2 + x3 + x5 + 2*x8,
639 2*x0 - 4*x11 - 4*x14 + 8*x17 + 2*x2 + 2*x3 + 2*x5 - 4*x8 - 4*x9,
640 3*x0/2 + x1/2 + 2*x12 - 3*x3/2 - x4/2 - 2*x6,
641 -3*x0 - x1 + 2*x10 + 4*x12 - 8*x15 - 3*x3 - x4 + 4*x6 + 6*x9,
642 Point(0,0,0),
643 -2*x0 - 2*x12 + 2*x13 - 2*x14 + 2*x3 + 2*x6 - 2*x7 + 2*x8,
644 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 -x0 - x1 - 2*x12 + x3 + x4 + 2*x6,
647 2*x0 + 2*x1 - 4*x10 - 4*x12 + 8*x15 + 2*x3 + 2*x4 - 4*x6 - 4*x9
648 };
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 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 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 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 Real vol = 0.;
729 for (int i=0; i<N2D; ++i)
730 for (int j=0; j<N1D; ++j)
731 {
732 // Compute dx_dxi, dx_deta, dx_dzeta at the current quadrature point.
733 Point dx_dxi_q, dx_deta_q, dx_dzeta_q;
734 for (int c=0; c<n_components; ++c)
735 {
736 Real coeff =
737 xi[i][exponents[c][0]]*
738 eta[i][exponents[c][1]]*
739 zeta[j][exponents[c][2]];
740
741 dx_dxi_q += coeff * dx_dxi[c];
742 dx_deta_q += coeff * dx_deta[c];
743 dx_dzeta_q += coeff * dx_dzeta[c];
744 }
745
746 // Compute scalar triple product, multiply by weight, and accumulate volume.
747 vol += w2D[i] * w1D[j] * triple_product(dx_dxi_q, dx_deta_q, dx_dzeta_q);
748 }
749
750 return vol;
751}
752
753
754
755
756#ifdef LIBMESH_ENABLE_AMR
757
758const 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
948void
949Prism18::permute(unsigned int perm_num)
950{
951 libmesh_assert_less (perm_num, 6);
952 const unsigned int side = perm_num % 2;
953 const unsigned int rotate = perm_num / 2;
954
955 for (unsigned int i = 0; i != rotate; ++i)
956 {
957 swap3nodes(0,1,2);
958 swap3nodes(3,4,5);
959 swap3nodes(6,7,8);
960 swap3nodes(9,10,11);
961 swap3nodes(12,13,14);
962 swap3nodes(15,16,17);
963 swap3neighbors(1,2,3);
964 }
965
966 switch (side) {
967 case 0:
968 break;
969 case 1:
970 swap2nodes(1,3);
971 swap2nodes(0,4);
972 swap2nodes(2,5);
973 swap2nodes(6,12);
974 swap2nodes(9,10);
975 swap2nodes(7,14);
976 swap2nodes(8,13);
977 swap2nodes(16,17);
978 swap2neighbors(0,4);
979 swap2neighbors(2,3);
980 break;
981 default:
982 libmesh_error();
983 }
984}
985
986
987void
988Prism18::flip(BoundaryInfo * boundary_info)
989{
990 libmesh_assert(boundary_info);
991
992 swap2nodes(0,1);
993 swap2nodes(3,4);
994 swap2nodes(7,8);
995 swap2nodes(9,10);
996 swap2nodes(13,14);
997 swap2nodes(16,17);
998 swap2neighbors(2,3);
999 swap2boundarysides(2,3,boundary_info);
1000 swap2boundaryedges(0,1,boundary_info);
1001 swap2boundaryedges(3,4,boundary_info);
1002 swap2boundaryedges(7,8,boundary_info);
1003}
1004
1005
1006unsigned int Prism18::center_node_on_side(const unsigned short side) const
1007{
1008 libmesh_assert_less (side, Prism18::num_sides);
1009 return (side >= 1 && side <= 3) ? side + 14 : invalid_uint;
1010}
1011
1012
1013ElemType
1014Prism18::side_type (const unsigned int s) const
1015{
1016 libmesh_assert_less (s, 5);
1017 if (s == 0 || s == 4)
1018 return TRI6;
1019 return QUAD9;
1020}
1021
1022
1023} // 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:2070
libMesh::Elem::point
const Point & point(const unsigned int i) const
Definition elem.h:2462
libMesh::Elem::swap2nodes
void swap2nodes(unsigned int n1, unsigned int n2)
Swaps two node_ptrs.
Definition elem.h:2101
libMesh::Elem::swap3neighbors
void swap3neighbors(unsigned int n1, unsigned int n2, unsigned int n3)
Swaps three neighbor_ptrs, "rotating" them.
Definition elem.h:2142
libMesh::Elem::_nodes
Node ** _nodes
Pointers to the nodes we are connected to.
Definition elem.h:2254
libMesh::Elem::swap2neighbors
void swap2neighbors(unsigned int n1, unsigned int n2)
Swaps two neighbor_ptrs.
Definition elem.h:2111
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:3579
libMesh::Elem::mapping_type
ElemMappingType mapping_type() const
Definition elem.h:3137
libMesh::Elem::compute_key
static dof_id_type compute_key(dof_id_type n0)
Definition elem.h:3311
libMesh::Elem::swap3nodes
void swap3nodes(unsigned int n1, unsigned int n2, unsigned int n3)
Swaps three node_ptrs, "rotating" them.
Definition elem.h:2133
libMesh::Elem::volume
virtual Real volume() const
Definition elem.C:3462
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:3595
libMesh::Elem::node_id
dof_id_type node_id(const unsigned int i) const
Definition elem.h:2484
libMesh::Point
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition point.h:40
libMesh::Prism18::is_vertex
virtual bool is_vertex(const unsigned int i) const override
Definition cell_prism18.C:64
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:301
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:561
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:314
libMesh::Prism18::side_nodes_map
static const unsigned int side_nodes_map[num_sides][nodes_per_side]
This maps the node of the side to element node numbers.
Definition cell_prism18.h:243
libMesh::Prism18::default_order
virtual Order default_order() const override
Definition cell_prism18.C:155
libMesh::Prism18::local_side_node
virtual unsigned int local_side_node(unsigned int side, unsigned int side_node) const override
Definition cell_prism18.C:193
libMesh::Prism18::is_face
virtual bool is_face(const unsigned int i) const override
Definition cell_prism18.C:80
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:287
libMesh::Prism18::local_edge_node
virtual unsigned int local_edge_node(unsigned int edge, unsigned int edge_node) const override
Definition cell_prism18.C:209
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:949
libMesh::Prism18::is_edge
virtual bool is_edge(const unsigned int i) const override
Definition cell_prism18.C:71
libMesh::Prism18::n_nodes
virtual unsigned int n_nodes() const override
Definition cell_prism18.h:101
libMesh::Prism18::center_node_on_side
unsigned int center_node_on_side(const unsigned short side) const override final
Definition cell_prism18.C:1006
libMesh::Prism18::volume
virtual Real volume() const override
A specialization for computing the volume of a Prism18.
Definition cell_prism18.C:573
libMesh::Prism18::has_affine_map
virtual bool has_affine_map() const override
Definition cell_prism18.C:122
libMesh::Prism18::key
virtual dof_id_type key() const
Don't hide Elem::key() defined in the base class.
Definition elem.C:738
libMesh::Prism18::n_sub_elem
virtual unsigned int n_sub_elem() const override
Definition cell_prism18.h:106
libMesh::Prism18::n_second_order_adjacent_vertices
virtual unsigned int n_second_order_adjacent_vertices(const unsigned int) const override
Definition cell_prism18.C:481
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:220
libMesh::Prism18::nodes_on_side
virtual std::vector< unsigned int > nodes_on_side(const unsigned int s) const override
Definition cell_prism18.C:97
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:87
libMesh::Prism18::edge_nodes_map
static const unsigned int edge_nodes_map[num_edges][nodes_per_edge]
This maps the node of the edge to element node numbers.
Definition cell_prism18.h:249
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:510
libMesh::Prism18::nodes_per_side
static const int nodes_per_side
Definition cell_prism18.h:236
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:111
libMesh::Prism18::nodes_per_edge
static const int nodes_per_edge
Definition cell_prism18.h:237
libMesh::Prism18::nodes_on_edge
virtual std::vector< unsigned int > nodes_on_edge(const unsigned int e) const override
Definition cell_prism18.C:105
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:988
libMesh::Prism18::num_nodes
static const int num_nodes
Geometric constants for Prism18.
Definition cell_prism18.h:235
libMesh::Prism18::side_type
ElemType side_type(const unsigned int s) const override final
Definition cell_prism18.C:1014
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:300
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:204
libMesh::Prism::key
virtual dof_id_type key() const
Don't hide Elem::key() defined in the base class.
Definition elem.C:738
libMesh::Prism::num_edges
static const int num_edges
Definition cell_prism.h:74
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:199
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:209
libMesh::TypeVector::relative_fuzzy_equals
bool relative_fuzzy_equals(const TypeVector< T > &rhs, Real tol=TOLERANCE) const
Definition type_vector.h:984
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::ElemType
ElemType
Defines an enum for geometric element types.
Definition enum_elem_type.h:33
libMesh::triple_product
T triple_product(const TypeVector< T > &a, const TypeVector< T > &b, const TypeVector< T > &c)
Definition type_vector.h:1032
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:144
Generated on Fri Aug 21 2026 13:32:29 for libMesh by doxygen 1.9.8