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