Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 :
19 : // Local includes
20 : #include "libmesh/cell_prism15.h"
21 : #include "libmesh/edge_edge3.h"
22 : #include "libmesh/face_quad8.h"
23 : #include "libmesh/face_tri6.h"
24 : #include "libmesh/enum_io_package.h"
25 : #include "libmesh/enum_order.h"
26 :
27 : namespace libMesh
28 : {
29 :
30 :
31 :
32 : // ------------------------------------------------------------
33 : // Prism15 class static member initializations
34 : const int Prism15::num_nodes;
35 : const int Prism15::nodes_per_side;
36 : const int Prism15::nodes_per_edge;
37 :
38 : const unsigned int Prism15::side_nodes_map[Prism15::num_sides][Prism15::nodes_per_side] =
39 : {
40 : {0, 2, 1, 8, 7, 6, 99, 99}, // Side 0
41 : {0, 1, 4, 3, 6, 10, 12, 9}, // Side 1
42 : {1, 2, 5, 4, 7, 11, 13, 10}, // Side 2
43 : {2, 0, 3, 5, 8, 9, 14, 11}, // Side 3
44 : {3, 4, 5, 12, 13, 14, 99, 99} // Side 4
45 : };
46 :
47 : const unsigned int Prism15::edge_nodes_map[Prism15::num_edges][Prism15::nodes_per_edge] =
48 : {
49 : {0, 1, 6}, // Edge 0
50 : {1, 2, 7}, // Edge 1
51 : {0, 2, 8}, // Edge 2
52 : {0, 3, 9}, // Edge 3
53 : {1, 4, 10}, // Edge 4
54 : {2, 5, 11}, // Edge 5
55 : {3, 4, 12}, // Edge 6
56 : {4, 5, 13}, // Edge 7
57 : {3, 5, 14} // Edge 8
58 : };
59 :
60 : // ------------------------------------------------------------
61 : // Prism15 class member functions
62 :
63 1188296 : bool Prism15::is_vertex(const unsigned int i) const
64 : {
65 1188296 : if (i < 6)
66 453326 : return true;
67 58482 : return false;
68 : }
69 :
70 10368 : bool Prism15::is_edge(const unsigned int i) const
71 : {
72 10368 : if (i < 6)
73 0 : return false;
74 864 : return true;
75 : }
76 :
77 0 : bool Prism15::is_face(const unsigned int) const
78 : {
79 0 : return false;
80 : }
81 :
82 22605 : bool Prism15::is_node_on_side(const unsigned int n,
83 : const unsigned int s) const
84 : {
85 1590 : libmesh_assert_less (s, n_sides());
86 1590 : return std::find(std::begin(side_nodes_map[s]),
87 22605 : std::end(side_nodes_map[s]),
88 22605 : n) != std::end(side_nodes_map[s]);
89 : }
90 :
91 : std::vector<unsigned int>
92 3572659 : Prism15::nodes_on_side(const unsigned int s) const
93 : {
94 318224 : libmesh_assert_less(s, n_sides());
95 3572659 : auto trim = (s > 0 && s < 4) ? 0 : 2;
96 3572659 : return {std::begin(side_nodes_map[s]), std::end(side_nodes_map[s]) - trim};
97 : }
98 :
99 : std::vector<unsigned>
100 11007 : Prism15::nodes_on_edge(const unsigned int e) const
101 : {
102 882 : libmesh_assert_less(e, n_edges());
103 11889 : return {std::begin(edge_nodes_map[e]), std::end(edge_nodes_map[e])};
104 : }
105 :
106 14769 : bool Prism15::is_node_on_edge(const unsigned int n,
107 : const unsigned int e) const
108 : {
109 702 : libmesh_assert_less (e, n_edges());
110 702 : return std::find(std::begin(edge_nodes_map[e]),
111 14769 : std::end(edge_nodes_map[e]),
112 14769 : n) != std::end(edge_nodes_map[e]);
113 : }
114 :
115 :
116 :
117 452288 : bool Prism15::has_affine_map() const
118 : {
119 : // Make sure z edges are affine
120 78488 : Point v = this->point(3) - this->point(0);
121 530776 : if (!v.relative_fuzzy_equals(this->point(4) - this->point(1), affine_tol) ||
122 530776 : !v.relative_fuzzy_equals(this->point(5) - this->point(2), affine_tol))
123 0 : return false;
124 : // Make sure edges are straight
125 39244 : v /= 2;
126 491532 : if (!v.relative_fuzzy_equals(this->point(9) - this->point(0), affine_tol) ||
127 904576 : !v.relative_fuzzy_equals(this->point(10) - this->point(1), affine_tol) ||
128 491532 : !v.relative_fuzzy_equals(this->point(11) - this->point(2), affine_tol))
129 0 : return false;
130 491532 : v = (this->point(1) - this->point(0))/2;
131 530776 : if (!v.relative_fuzzy_equals(this->point(6) - this->point(0), affine_tol) ||
132 491532 : !v.relative_fuzzy_equals(this->point(12) - this->point(3), affine_tol))
133 0 : return false;
134 491532 : v = (this->point(2) - this->point(0))/2;
135 530776 : if (!v.relative_fuzzy_equals(this->point(8) - this->point(0), affine_tol) ||
136 491532 : !v.relative_fuzzy_equals(this->point(14) - this->point(3), affine_tol))
137 0 : return false;
138 491532 : v = (this->point(2) - this->point(1))/2;
139 530776 : if (!v.relative_fuzzy_equals(this->point(7) - this->point(1), affine_tol) ||
140 530776 : !v.relative_fuzzy_equals(this->point(13) - this->point(4), affine_tol))
141 0 : return false;
142 39244 : return true;
143 : }
144 :
145 :
146 :
147 14044221 : Order Prism15::default_order() const
148 : {
149 14044221 : return SECOND;
150 : }
151 :
152 :
153 :
154 1349 : unsigned int Prism15::local_side_node(unsigned int side,
155 : unsigned int side_node) const
156 : {
157 38 : libmesh_assert_less (side, this->n_sides());
158 :
159 : // Never more than 8 nodes per side.
160 38 : libmesh_assert_less(side_node, Prism15::nodes_per_side);
161 :
162 : // Some sides have 6 nodes.
163 38 : libmesh_assert(!(side==0 || side==4) || side_node < 6);
164 :
165 1349 : return Prism15::side_nodes_map[side][side_node];
166 : }
167 :
168 :
169 :
170 7170813 : unsigned int Prism15::local_edge_node(unsigned int edge,
171 : unsigned int edge_node) const
172 : {
173 638514 : libmesh_assert_less(edge, this->n_edges());
174 638514 : libmesh_assert_less(edge_node, Prism15::nodes_per_edge);
175 :
176 7170813 : return Prism15::edge_nodes_map[edge][edge_node];
177 : }
178 :
179 :
180 :
181 435517 : std::unique_ptr<Elem> Prism15::build_side_ptr (const unsigned int i)
182 : {
183 21814 : libmesh_assert_less (i, this->n_sides());
184 :
185 435517 : std::unique_ptr<Elem> face;
186 :
187 435517 : switch (i)
188 : {
189 249172 : case 0: // the triangular face at z=-1
190 : case 4: // the triangular face at z=1
191 : {
192 249172 : face = std::make_unique<Tri6>();
193 249172 : break;
194 : }
195 186345 : case 1: // the quad face at y=0
196 : case 2: // the other quad face
197 : case 3: // the quad face at x=0
198 : {
199 186345 : face = std::make_unique<Quad8>();
200 186345 : break;
201 : }
202 0 : default:
203 0 : libmesh_error_msg("Invalid side i = " << i);
204 : }
205 :
206 : // Set the nodes
207 3421309 : for (auto n : face->node_index_range())
208 3137556 : face->set_node(n, this->node_ptr(Prism15::side_nodes_map[i][n]));
209 :
210 435517 : face->set_interior_parent(this);
211 413703 : face->inherit_data_from(*this);
212 :
213 435517 : return face;
214 0 : }
215 :
216 :
217 377507 : void Prism15::build_side_ptr (std::unique_ptr<Elem> & side,
218 : const unsigned int i)
219 : {
220 13418 : libmesh_assert_less (i, this->n_sides());
221 :
222 377507 : switch (i)
223 : {
224 3392 : case 0: // the triangular face at z=-1
225 : case 4: // the triangular face at z=1
226 : {
227 94862 : if (!side.get() || side->type() != TRI6)
228 : {
229 181668 : side = this->build_side_ptr(i);
230 94164 : return;
231 : }
232 62 : break;
233 : }
234 :
235 10026 : case 1: // the quad face at y=0
236 : case 2: // the other quad face
237 : case 3: // the quad face at x=0
238 : {
239 282645 : if (!side.get() || side->type() != QUAD8)
240 : {
241 181178 : side = this->build_side_ptr(i);
242 93901 : return;
243 : }
244 6714 : break;
245 : }
246 :
247 0 : default:
248 0 : libmesh_error_msg("Invalid side i = " << i);
249 : }
250 :
251 182666 : side->inherit_data_from(*this);
252 :
253 : // Set the nodes
254 1703582 : for (auto n : side->node_index_range())
255 1568224 : side->set_node(n, this->node_ptr(Prism15::side_nodes_map[i][n]));
256 : }
257 :
258 :
259 :
260 753713 : std::unique_ptr<Elem> Prism15::build_edge_ptr (const unsigned int i)
261 : {
262 753713 : return this->simple_build_edge_ptr<Edge3,Prism15>(i);
263 : }
264 :
265 :
266 :
267 0 : void Prism15::build_edge_ptr (std::unique_ptr<Elem> & edge, const unsigned int i)
268 : {
269 0 : this->simple_build_edge_ptr<Prism15>(edge, i, EDGE3);
270 0 : }
271 :
272 :
273 :
274 0 : void Prism15::connectivity(const unsigned int libmesh_dbg_var(sc),
275 : const IOPackage iop,
276 : std::vector<dof_id_type> & conn) const
277 : {
278 0 : libmesh_assert(_nodes);
279 0 : libmesh_assert_less (sc, this->n_sub_elem());
280 0 : libmesh_assert_not_equal_to (iop, INVALID_IO_PACKAGE);
281 :
282 0 : switch (iop)
283 : {
284 0 : case TECPLOT:
285 : {
286 0 : conn.resize(8);
287 0 : conn[0] = this->node_id(0)+1;
288 0 : conn[1] = this->node_id(1)+1;
289 0 : conn[2] = this->node_id(2)+1;
290 0 : conn[3] = this->node_id(2)+1;
291 0 : conn[4] = this->node_id(3)+1;
292 0 : conn[5] = this->node_id(4)+1;
293 0 : conn[6] = this->node_id(5)+1;
294 0 : conn[7] = this->node_id(5)+1;
295 0 : return;
296 : }
297 :
298 0 : case VTK:
299 : {
300 : // VTK's VTK_QUADRATIC_WEDGE first 9 nodes match, then their
301 : // middle and top layers of mid-edge nodes are reversed from
302 : // LibMesh's.
303 0 : conn.resize(15);
304 0 : for (unsigned i=0; i<9; ++i)
305 0 : conn[i] = this->node_id(i);
306 :
307 : // top "ring" of mid-edge nodes
308 0 : conn[9] = this->node_id(12);
309 0 : conn[10] = this->node_id(13);
310 0 : conn[11] = this->node_id(14);
311 :
312 : // middle "ring" of mid-edge nodes
313 0 : conn[12] = this->node_id(9);
314 0 : conn[13] = this->node_id(10);
315 0 : conn[14] = this->node_id(11);
316 :
317 0 : return;
318 : }
319 :
320 0 : default:
321 0 : libmesh_error_msg("Unsupported IO package " << iop);
322 : }
323 : }
324 :
325 :
326 :
327 :
328 0 : unsigned short int Prism15::second_order_adjacent_vertex (const unsigned int n,
329 : const unsigned int v) const
330 : {
331 0 : libmesh_assert_greater_equal (n, this->n_vertices());
332 0 : libmesh_assert_less (n, this->n_nodes());
333 0 : libmesh_assert_less (v, 2);
334 0 : return _second_order_adjacent_vertices[n-this->n_vertices()][v];
335 : }
336 :
337 :
338 :
339 : std::pair<unsigned short int, unsigned short int>
340 0 : Prism15::second_order_child_vertex (const unsigned int n) const
341 : {
342 0 : libmesh_assert_greater_equal (n, this->n_vertices());
343 0 : libmesh_assert_less (n, this->n_nodes());
344 :
345 0 : return std::pair<unsigned short int, unsigned short int>
346 0 : (_second_order_vertex_child_number[n],
347 0 : _second_order_vertex_child_index[n]);
348 : }
349 :
350 :
351 :
352 0 : Real Prism15::volume () const
353 : {
354 : // This specialization is good for Lagrange mappings only
355 0 : if (this->mapping_type() != LAGRANGE_MAP)
356 0 : return this->Elem::volume();
357 :
358 : // Make copies of our points. It makes the subsequent calculations a bit
359 : // shorter and avoids dereferencing the same pointer multiple times.
360 : Point
361 0 : x0 = point(0), x1 = point(1), x2 = point(2), x3 = point(3), x4 = point(4),
362 0 : x5 = point(5), x6 = point(6), x7 = point(7), x8 = point(8), x9 = point(9),
363 0 : x10 = point(10), x11 = point(11), x12 = point(12), x13 = point(13), x14 = point(14);
364 :
365 : // Terms are copied directly from a Python script.
366 : Point dx_dxi[10] =
367 : {
368 0 : -x0 - x1 + x10 + 2*x12 - x3 - x4 + 2*x6 - x9,
369 0 : 3*x0/2 + x1/2 + 2*x12 - 3*x3/2 - x4/2 - 2*x6,
370 0 : -x0/2 + x1/2 - x10 - x3/2 + x4/2 + x9,
371 0 : 2*x0 - 2*x12 + 2*x13 - 2*x14 + 2*x3 - 2*x6 + 2*x7 - 2*x8,
372 0 : -2*x0 - 2*x12 + 2*x13 - 2*x14 + 2*x3 + 2*x6 - 2*x7 + 2*x8,
373 : Point(0,0,0),
374 0 : 2*x0 + 2*x1 - 4*x12 + 2*x3 + 2*x4 - 4*x6,
375 0 : -2*x0 - 2*x1 - 4*x12 + 2*x3 + 2*x4 + 4*x6,
376 : Point(0,0,0),
377 : Point(0,0,0)
378 0 : };
379 :
380 : Point dx_deta[10] =
381 : {
382 0 : -x0 + x11 + 2*x14 - x2 - x3 - x5 + 2*x8 - x9,
383 0 : 3*x0/2 + 2*x14 + x2/2 - 3*x3/2 - x5/2 - 2*x8,
384 0 : -x0/2 - x11 + x2/2 - x3/2 + x5/2 + x9,
385 0 : 2*x0 - 4*x14 + 2*x2 + 2*x3 + 2*x5 - 4*x8,
386 0 : -2*x0 - 4*x14 - 2*x2 + 2*x3 + 2*x5 + 4*x8,
387 : Point(0,0,0),
388 0 : 2*x0 - 2*x12 + 2*x13 - 2*x14 + 2*x3 - 2*x6 + 2*x7 - 2*x8,
389 0 : -2*x0 - 2*x12 + 2*x13 - 2*x14 + 2*x3 + 2*x6 - 2*x7 + 2*x8,
390 : Point(0,0,0),
391 : Point(0,0,0)
392 0 : };
393 :
394 : Point dx_dzeta[10] =
395 : {
396 0 : -x0/2 + x3/2,
397 0 : x0 + x3 - 2*x9,
398 : Point(0,0,0),
399 0 : 3*x0/2 + 2*x14 + x2/2 - 3*x3/2 - x5/2 - 2*x8,
400 0 : -x0 - 2*x11 + x2 - x3 + x5 + 2*x9,
401 0 : -x0 - 2*x14 - x2 + x3 + x5 + 2*x8,
402 0 : 3*x0/2 + x1/2 + 2*x12 - 3*x3/2 - x4/2 - 2*x6,
403 0 : -x0 + x1 - 2*x10 - x3 + x4 + 2*x9,
404 0 : -2*x0 - 2*x12 + 2*x13 - 2*x14 + 2*x3 + 2*x6 - 2*x7 + 2*x8,
405 0 : -x0 - x1 - 2*x12 + x3 + x4 + 2*x6
406 0 : };
407 :
408 : // The quadrature rule for the Prism15 is a tensor product between a
409 : // FOURTH-order TRI3 rule (in xi, eta) and a FIFTH-order EDGE2 rule
410 : // in zeta.
411 :
412 : // Number of points in the 2D quadrature rule.
413 0 : const int N2D = 6;
414 :
415 : // Parameters of the 2D rule
416 : static const Real
417 : w1 = 1.1169079483900573284750350421656140e-01_R,
418 : w2 = 5.4975871827660933819163162450105264e-02_R,
419 : a1 = 4.4594849091596488631832925388305199e-01_R,
420 : a2 = 9.1576213509770743459571463402201508e-02_R;
421 :
422 : // Points and weights of the 2D rule
423 : static const Real w2D[N2D] = {w1, w1, w1, w2, w2, w2};
424 :
425 : // Quadrature point locations raised to powers. xi[0][2] is
426 : // quadrature point 0, squared, xi[1][1] is quadrature point 1 to the
427 : // first power, etc. This lets us avoid calling std::pow inside the
428 : // loops below.
429 : static const Real xi[N2D][3] =
430 : {
431 : // ^0 ^1 ^2
432 : { 1., a1, a1*a1},
433 : { 1., 1-2*a1, (1-2*a1)*(1-2*a1)},
434 : { 1., a1, a1*a1},
435 : { 1., a2, a2*a2},
436 : { 1., 1-2*a2, (1-2*a2)*(1-2*a2)},
437 : { 1., a2, a2*a2}
438 : };
439 :
440 : static const Real eta[N2D][3] =
441 : {
442 : // ^0 ^1 ^2
443 : { 1., a1, a1*a1},
444 : { 1., a1, a1*a1},
445 : { 1., 1-2*a1, (1-2*a1)*(1-2*a1)},
446 : { 1., a2, a2*a2},
447 : { 1., a2, a2*a2},
448 : { 1., 1-2*a2, (1-2*a2)*(1-2*a2)}
449 : };
450 :
451 : // Number of points in the 1D quadrature rule.
452 0 : const int N1D = 3;
453 :
454 : // Points and weights of the 1D quadrature rule.
455 : static const Real w1D[N1D] = {5./9, 8./9, 5./9};
456 :
457 0 : const Real zeta[N1D][3] =
458 : {
459 : //^0 ^1 ^2
460 : { 1., -std::sqrt(15)/5., 15./25},
461 : { 1., 0., 0.},
462 : { 1., std::sqrt(15)/5., 15./25}
463 : };
464 :
465 : // The integer exponents for each term.
466 : static const int exponents[10][3] =
467 : {
468 : {0, 0, 0},
469 : {0, 0, 1},
470 : {0, 0, 2},
471 : {0, 1, 0},
472 : {0, 1, 1},
473 : {0, 2, 0},
474 : {1, 0, 0},
475 : {1, 0, 1},
476 : {1, 1, 0},
477 : {2, 0, 0}
478 : };
479 :
480 0 : Real vol = 0.;
481 0 : for (int i=0; i<N2D; ++i)
482 0 : for (int j=0; j<N1D; ++j)
483 : {
484 : // Compute dx_dxi, dx_deta, dx_dzeta at the current quadrature point.
485 0 : Point dx_dxi_q, dx_deta_q, dx_dzeta_q;
486 0 : for (int c=0; c<10; ++c)
487 : {
488 0 : Real coeff =
489 0 : xi[i][exponents[c][0]]*
490 0 : eta[i][exponents[c][1]]*
491 0 : zeta[j][exponents[c][2]];
492 :
493 0 : dx_dxi_q += coeff * dx_dxi[c];
494 0 : dx_deta_q += coeff * dx_deta[c];
495 0 : dx_dzeta_q += coeff * dx_dzeta[c];
496 : }
497 :
498 : // Compute scalar triple product, multiply by weight, and accumulate volume.
499 0 : vol += w2D[i] * w1D[j] * triple_product(dx_dxi_q, dx_deta_q, dx_dzeta_q);
500 : }
501 :
502 0 : return vol;
503 : }
504 :
505 :
506 :
507 : #ifdef LIBMESH_ENABLE_AMR
508 :
509 : const Real Prism15::_embedding_matrix[Prism15::num_children][Prism15::num_nodes][Prism15::num_nodes] =
510 : {
511 : // Embedding matrix for child 0
512 : {
513 : // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
514 : { 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, // 0
515 : { 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0 }, // 1
516 : { 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0 }, // 2
517 : { 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0 }, // 3
518 : { -0.25, -0.25, 0, -0.25, -0.25, 0, 0.5, 0, 0, 0.5, 0.5, 0, 0.5, 0, 0 }, // 4
519 : { -0.25, 0, -0.25, -0.25, 0, -0.25, 0, 0, 0.5, 0.5, 0, 0.5, 0, 0, 0.5 }, // 5
520 : { 0.375, -0.125, 0, 0, 0, 0, 0.75, 0, 0, 0, 0, 0, 0, 0, 0 }, // 6
521 : { 0, -0.125, -0.125, 0, 0, 0, 0.5, 0.25, 0.5, 0, 0, 0, 0, 0, 0 }, // 7
522 : { 0.375, 0, -0.125, 0, 0, 0, 0, 0, 0.75, 0, 0, 0, 0, 0, 0 }, // 8
523 : { 0.375, 0, 0, -0.125, 0, 0, 0, 0, 0, 0.75, 0, 0, 0, 0, 0 }, // 9
524 : { -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.75, 0, 0, 0.375, 0.375, 0, 0.25, 0, 0 }, // 10
525 : { -0.1875, 0, -0.1875, -0.1875, 0, -0.1875, 0, 0, 0.75, 0.375, 0, 0.375, 0, 0, 0.25 }, // 11
526 : { -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.375, 0, 0, 0.75, 0.25, 0, 0.375, 0, 0 }, // 12
527 : { -0.25, -0.1875, -0.1875, -0.25, -0.1875, -0.1875, 0.25, 0.125, 0.25, 0.5, 0.25, 0.25, 0.25, 0.125, 0.25 }, // 13
528 : { -0.1875, 0, -0.1875, -0.1875, 0, -0.1875, 0, 0, 0.375, 0.75, 0, 0.25, 0, 0, 0.375 } // 14
529 : },
530 :
531 : // Embedding matrix for child 1
532 : {
533 : // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
534 : { 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0 }, // 0
535 : { 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, // 1
536 : { 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0 }, // 2
537 : { -0.25, -0.25, 0, -0.25, -0.25, 0, 0.5, 0, 0, 0.5, 0.5, 0, 0.5, 0, 0 }, // 3
538 : { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0 }, // 4
539 : { 0, -0.25, -0.25, 0, -0.25, -0.25, 0, 0.5, 0, 0, 0.5, 0.5, 0, 0.5, 0 }, // 5
540 : { -0.125, 0.375, 0, 0, 0, 0, 0.75, 0, 0, 0, 0, 0, 0, 0, 0 }, // 6
541 : { 0, 0.375, -0.125, 0, 0, 0, 0, 0.75, 0, 0, 0, 0, 0, 0, 0 }, // 7
542 : { -0.125, 0, -0.125, 0, 0, 0, 0.5, 0.5, 0.25, 0, 0, 0, 0, 0, 0 }, // 8
543 : { -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.75, 0, 0, 0.375, 0.375, 0, 0.25, 0, 0 }, // 9
544 : { 0, 0.375, 0, 0, -0.125, 0, 0, 0, 0, 0, 0.75, 0, 0, 0, 0 }, // 10
545 : { 0, -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.75, 0, 0, 0.375, 0.375, 0, 0.25, 0 }, // 11
546 : { -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.375, 0, 0, 0.25, 0.75, 0, 0.375, 0, 0 }, // 12
547 : { 0, -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.375, 0, 0, 0.75, 0.25, 0, 0.375, 0 }, // 13
548 : { -0.1875, -0.25, -0.1875, -0.1875, -0.25, -0.1875, 0.25, 0.25, 0.125, 0.25, 0.5, 0.25, 0.25, 0.25, 0.125 } // 14
549 : },
550 :
551 : // Embedding matrix for child 2
552 : {
553 : // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
554 : { 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0 }, // 0
555 : { 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0 }, // 1
556 : { 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, // 2
557 : { -0.25, 0, -0.25, -0.25, 0, -0.25, 0, 0, 0.5, 0.5, 0, 0.5, 0, 0, 0.5 }, // 3
558 : { 0, -0.25, -0.25, 0, -0.25, -0.25, 0, 0.5, 0, 0, 0.5, 0.5, 0, 0.5, 0 }, // 4
559 : { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0 }, // 5
560 : { -0.125, -0.125, 0, 0, 0, 0, 0.25, 0.5, 0.5, 0, 0, 0, 0, 0, 0 }, // 6
561 : { 0, -0.125, 0.375, 0, 0, 0, 0, 0.75, 0, 0, 0, 0, 0, 0, 0 }, // 7
562 : { -0.125, 0, 0.375, 0, 0, 0, 0, 0, 0.75, 0, 0, 0, 0, 0, 0 }, // 8
563 : { -0.1875, 0, -0.1875, -0.1875, 0, -0.1875, 0, 0, 0.75, 0.375, 0, 0.375, 0, 0, 0.25 }, // 9
564 : { 0, -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.75, 0, 0, 0.375, 0.375, 0, 0.25, 0 }, // 10
565 : { 0, 0, 0.375, 0, 0, -0.125, 0, 0, 0, 0, 0, 0.75, 0, 0, 0 }, // 11
566 : { -0.1875, -0.1875, -0.25, -0.1875, -0.1875, -0.25, 0.125, 0.25, 0.25, 0.25, 0.25, 0.5, 0.125, 0.25, 0.25 }, // 12
567 : { 0, -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.375, 0, 0, 0.25, 0.75, 0, 0.375, 0 }, // 13
568 : { -0.1875, 0, -0.1875, -0.1875, 0, -0.1875, 0, 0, 0.375, 0.25, 0, 0.75, 0, 0, 0.375 } // 14
569 : },
570 :
571 : // Embedding matrix for child 3
572 : {
573 : // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
574 : { 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0 }, // 0
575 : { 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0 }, // 1
576 : { 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0 }, // 2
577 : { -0.25, -0.25, 0, -0.25, -0.25, 0, 0.5, 0, 0, 0.5, 0.5, 0, 0.5, 0, 0 }, // 3
578 : { 0, -0.25, -0.25, 0, -0.25, -0.25, 0, 0.5, 0, 0, 0.5, 0.5, 0, 0.5, 0 }, // 4
579 : { -0.25, 0, -0.25, -0.25, 0, -0.25, 0, 0, 0.5, 0.5, 0, 0.5, 0, 0, 0.5 }, // 5
580 : { -0.125, 0, -0.125, 0, 0, 0, 0.5, 0.5, 0.25, 0, 0, 0, 0, 0, 0 }, // 6
581 : { -0.125, -0.125, 0, 0, 0, 0, 0.25, 0.5, 0.5, 0, 0, 0, 0, 0, 0 }, // 7
582 : { 0, -0.125, -0.125, 0, 0, 0, 0.5, 0.25, 0.5, 0, 0, 0, 0, 0, 0 }, // 8
583 : { -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.75, 0, 0, 0.375, 0.375, 0, 0.25, 0, 0 }, // 9
584 : { 0, -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.75, 0, 0, 0.375, 0.375, 0, 0.25, 0 }, // 10
585 : { -0.1875, 0, -0.1875, -0.1875, 0, -0.1875, 0, 0, 0.75, 0.375, 0, 0.375, 0, 0, 0.25 }, // 11
586 : { -0.1875, -0.25, -0.1875, -0.1875, -0.25, -0.1875, 0.25, 0.25, 0.125, 0.25, 0.5, 0.25, 0.25, 0.25, 0.125 }, // 12
587 : { -0.1875, -0.1875, -0.25, -0.1875, -0.1875, -0.25, 0.125, 0.25, 0.25, 0.25, 0.25, 0.5, 0.125, 0.25, 0.25 }, // 13
588 : { -0.25, -0.1875, -0.1875, -0.25, -0.1875, -0.1875, 0.25, 0.125, 0.25, 0.5, 0.25, 0.25, 0.25, 0.125, 0.25 } // 14
589 : },
590 :
591 : // Embedding matrix for child 4
592 : {
593 : // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
594 : { 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0 }, // 0
595 : { -0.25, -0.25, 0, -0.25, -0.25, 0, 0.5, 0, 0, 0.5, 0.5, 0, 0.5, 0, 0 }, // 1
596 : { -0.25, 0, -0.25, -0.25, 0, -0.25, 0, 0, 0.5, 0.5, 0, 0.5, 0, 0, 0.5 }, // 2
597 : { 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, // 3
598 : { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0 }, // 4
599 : { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1 }, // 5
600 : { -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.375, 0, 0, 0.75, 0.25, 0, 0.375, 0, 0 }, // 6
601 : { -0.25, -0.1875, -0.1875, -0.25, -0.1875, -0.1875, 0.25, 0.125, 0.25, 0.5, 0.25, 0.25, 0.25, 0.125, 0.25 }, // 7
602 : { -0.1875, 0, -0.1875, -0.1875, 0, -0.1875, 0, 0, 0.375, 0.75, 0, 0.25, 0, 0, 0.375 }, // 8
603 : { -0.125, 0, 0, 0.375, 0, 0, 0, 0, 0, 0.75, 0, 0, 0, 0, 0 }, // 9
604 : { -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.25, 0, 0, 0.375, 0.375, 0, 0.75, 0, 0 }, // 10
605 : { -0.1875, 0, -0.1875, -0.1875, 0, -0.1875, 0, 0, 0.25, 0.375, 0, 0.375, 0, 0, 0.75 }, // 11
606 : { 0, 0, 0, 0.375, -0.125, 0, 0, 0, 0, 0, 0, 0, 0.75, 0, 0 }, // 12
607 : { 0, 0, 0, 0, -0.125, -0.125, 0, 0, 0, 0, 0, 0, 0.5, 0.25, 0.5 }, // 13
608 : { 0, 0, 0, 0.375, 0, -0.125, 0, 0, 0, 0, 0, 0, 0, 0, 0.75 } // 14
609 : },
610 :
611 : // Embedding matrix for child 5
612 : {
613 : // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
614 : { -0.25, -0.25, 0, -0.25, -0.25, 0, 0.5, 0, 0, 0.5, 0.5, 0, 0.5, 0, 0 }, // 0
615 : { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0 }, // 1
616 : { 0, -0.25, -0.25, 0, -0.25, -0.25, 0, 0.5, 0, 0, 0.5, 0.5, 0, 0.5, 0 }, // 2
617 : { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0 }, // 3
618 : { 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, // 4
619 : { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0 }, // 5
620 : { -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.375, 0, 0, 0.25, 0.75, 0, 0.375, 0, 0 }, // 6
621 : { 0, -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.375, 0, 0, 0.75, 0.25, 0, 0.375, 0 }, // 7
622 : { -0.1875, -0.25, -0.1875, -0.1875, -0.25, -0.1875, 0.25, 0.25, 0.125, 0.25, 0.5, 0.25, 0.25, 0.25, 0.125 }, // 8
623 : { -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.25, 0, 0, 0.375, 0.375, 0, 0.75, 0, 0 }, // 9
624 : { 0, -0.125, 0, 0, 0.375, 0, 0, 0, 0, 0, 0.75, 0, 0, 0, 0 }, // 10
625 : { 0, -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.25, 0, 0, 0.375, 0.375, 0, 0.75, 0 }, // 11
626 : { 0, 0, 0, -0.125, 0.375, 0, 0, 0, 0, 0, 0, 0, 0.75, 0, 0 }, // 12
627 : { 0, 0, 0, 0, 0.375, -0.125, 0, 0, 0, 0, 0, 0, 0, 0.75, 0 }, // 13
628 : { 0, 0, 0, -0.125, 0, -0.125, 0, 0, 0, 0, 0, 0, 0.5, 0.5, 0.25 } // 14
629 : },
630 :
631 : // Embedding matrix for child 6
632 : {
633 : // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
634 : { -0.25, 0, -0.25, -0.25, 0, -0.25, 0, 0, 0.5, 0.5, 0, 0.5, 0, 0, 0.5 }, // 0
635 : { 0, -0.25, -0.25, 0, -0.25, -0.25, 0, 0.5, 0, 0, 0.5, 0.5, 0, 0.5, 0 }, // 1
636 : { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0 }, // 2
637 : { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1 }, // 3
638 : { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0 }, // 4
639 : { 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, // 5
640 : { -0.1875, -0.1875, -0.25, -0.1875, -0.1875, -0.25, 0.125, 0.25, 0.25, 0.25, 0.25, 0.5, 0.125, 0.25, 0.25 }, // 6
641 : { 0, -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.375, 0, 0, 0.25, 0.75, 0, 0.375, 0 }, // 7
642 : { -0.1875, 0, -0.1875, -0.1875, 0, -0.1875, 0, 0, 0.375, 0.25, 0, 0.75, 0, 0, 0.375 }, // 8
643 : { -0.1875, 0, -0.1875, -0.1875, 0, -0.1875, 0, 0, 0.25, 0.375, 0, 0.375, 0, 0, 0.75 }, // 9
644 : { 0, -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.25, 0, 0, 0.375, 0.375, 0, 0.75, 0 }, // 10
645 : { 0, 0, -0.125, 0, 0, 0.375, 0, 0, 0, 0, 0, 0.75, 0, 0, 0 }, // 11
646 : { 0, 0, 0, -0.125, -0.125, 0, 0, 0, 0, 0, 0, 0, 0.25, 0.5, 0.5 }, // 12
647 : { 0, 0, 0, 0, -0.125, 0.375, 0, 0, 0, 0, 0, 0, 0, 0.75, 0 }, // 13
648 : { 0, 0, 0, -0.125, 0, 0.375, 0, 0, 0, 0, 0, 0, 0, 0, 0.75 } // 14
649 : },
650 :
651 : // Embedding matrix for child 7
652 : {
653 : // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
654 : { -0.25, -0.25, 0, -0.25, -0.25, 0, 0.5, 0, 0, 0.5, 0.5, 0, 0.5, 0, 0 }, // 0
655 : { 0, -0.25, -0.25, 0, -0.25, -0.25, 0, 0.5, 0, 0, 0.5, 0.5, 0, 0.5, 0 }, // 1
656 : { -0.25, 0, -0.25, -0.25, 0, -0.25, 0, 0, 0.5, 0.5, 0, 0.5, 0, 0, 0.5 }, // 2
657 : { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0 }, // 3
658 : { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0 }, // 4
659 : { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1 }, // 5
660 : { -0.1875, -0.25, -0.1875, -0.1875, -0.25, -0.1875, 0.25, 0.25, 0.125, 0.25, 0.5, 0.25, 0.25, 0.25, 0.125 }, // 6
661 : { -0.1875, -0.1875, -0.25, -0.1875, -0.1875, -0.25, 0.125, 0.25, 0.25, 0.25, 0.25, 0.5, 0.125, 0.25, 0.25 }, // 7
662 : { -0.25, -0.1875, -0.1875, -0.25, -0.1875, -0.1875, 0.25, 0.125, 0.25, 0.5, 0.25, 0.25, 0.25, 0.125, 0.25 }, // 8
663 : { -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.25, 0, 0, 0.375, 0.375, 0, 0.75, 0, 0 }, // 9
664 : { 0, -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.25, 0, 0, 0.375, 0.375, 0, 0.75, 0 }, // 10
665 : { -0.1875, 0, -0.1875, -0.1875, 0, -0.1875, 0, 0, 0.25, 0.375, 0, 0.375, 0, 0, 0.75 }, // 11
666 : { 0, 0, 0, -0.125, 0, -0.125, 0, 0, 0, 0, 0, 0, 0.5, 0.5, 0.25 }, // 12
667 : { 0, 0, 0, -0.125, -0.125, 0, 0, 0, 0, 0, 0, 0, 0.25, 0.5, 0.5 }, // 13
668 : { 0, 0, 0, 0, -0.125, -0.125, 0, 0, 0, 0, 0, 0, 0.5, 0.25, 0.5 } // 14
669 : }
670 : };
671 :
672 : #endif
673 :
674 :
675 : void
676 5112 : Prism15::permute(unsigned int perm_num)
677 : {
678 456 : libmesh_assert_less (perm_num, 6);
679 5112 : const unsigned int side = perm_num % 2;
680 5112 : const unsigned int rotate = perm_num / 2;
681 :
682 12204 : for (unsigned int i = 0; i != rotate; ++i)
683 : {
684 7092 : swap3nodes(0,1,2);
685 6456 : swap3nodes(3,4,5);
686 6456 : swap3nodes(6,7,8);
687 6456 : swap3nodes(9,10,11);
688 6456 : swap3nodes(12,13,14);
689 6456 : swap3neighbors(1,2,3);
690 : }
691 :
692 5112 : switch (side) {
693 48 : case 0:
694 48 : break;
695 4536 : case 1:
696 4536 : swap2nodes(1,3);
697 4536 : swap2nodes(0,4);
698 4536 : swap2nodes(2,5);
699 4536 : swap2nodes(6,12);
700 4536 : swap2nodes(9,10);
701 4536 : swap2nodes(7,14);
702 4536 : swap2nodes(8,13);
703 408 : swap2neighbors(0,4);
704 408 : swap2neighbors(2,3);
705 408 : break;
706 0 : default:
707 0 : libmesh_error();
708 : }
709 :
710 5112 : }
711 :
712 :
713 : void
714 576 : Prism15::flip(BoundaryInfo * boundary_info)
715 : {
716 48 : libmesh_assert(boundary_info);
717 :
718 576 : swap2nodes(0,1);
719 576 : swap2nodes(3,4);
720 576 : swap2nodes(7,8);
721 576 : swap2nodes(9,10);
722 576 : swap2nodes(13,14);
723 48 : swap2neighbors(2,3);
724 576 : swap2boundarysides(2,3,boundary_info);
725 576 : swap2boundaryedges(0,1,boundary_info);
726 576 : swap2boundaryedges(3,4,boundary_info);
727 576 : swap2boundaryedges(7,8,boundary_info);
728 576 : }
729 :
730 :
731 : ElemType
732 2880 : Prism15::side_type (const unsigned int s) const
733 : {
734 240 : libmesh_assert_less (s, 5);
735 2880 : if (s == 0 || s == 4)
736 1152 : return TRI6;
737 144 : return QUAD8;
738 : }
739 :
740 : } // namespace libMesh
|
