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_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 : // ------------------------------------------------------------
39 : // Prism15 class member functions
40 :
41 1424426 : bool Prism15::is_vertex(const unsigned int i) const
42 : {
43 1424426 : if (i < 6)
44 567866 : return true;
45 68064 : return false;
46 : }
47 :
48 34560 : bool Prism15::is_edge(const unsigned int i) const
49 : {
50 34560 : if (i < 6)
51 0 : return false;
52 2880 : return true;
53 : }
54 :
55 0 : bool Prism15::is_face(const unsigned int) const
56 : {
57 0 : return false;
58 : }
59 :
60 44085 : bool Prism15::is_node_on_side(const unsigned int n,
61 : const unsigned int s) const
62 : {
63 3380 : libmesh_assert_less (s, n_sides());
64 3380 : return std::find(std::begin(side_nodes_map[s]),
65 44085 : std::end(side_nodes_map[s]),
66 44085 : n) != std::end(side_nodes_map[s]);
67 : }
68 :
69 : std::vector<unsigned int>
70 4098817 : Prism15::nodes_on_side(const unsigned int s) const
71 : {
72 365336 : libmesh_assert_less(s, n_sides());
73 4098817 : auto trim = (s > 0 && s < 4) ? 0 : 2;
74 4098817 : return {std::begin(side_nodes_map[s]), std::end(side_nodes_map[s]) - trim};
75 : }
76 :
77 : std::vector<unsigned>
78 28287 : Prism15::nodes_on_edge(const unsigned int e) const
79 : {
80 2322 : libmesh_assert_less(e, n_edges());
81 30609 : return {std::begin(edge_nodes_map[e]), std::end(edge_nodes_map[e])};
82 : }
83 :
84 98001 : bool Prism15::is_node_on_edge(const unsigned int n,
85 : const unsigned int e) const
86 : {
87 7638 : libmesh_assert_less (e, n_edges());
88 7638 : return std::find(std::begin(edge_nodes_map[e]),
89 98001 : std::end(edge_nodes_map[e]),
90 98001 : n) != std::end(edge_nodes_map[e]);
91 : }
92 :
93 :
94 :
95 505033 : bool Prism15::has_affine_map() const
96 : {
97 : // Make sure z edges are affine
98 87952 : Point v = this->point(3) - this->point(0);
99 592473 : if (!v.relative_fuzzy_equals(this->point(4) - this->point(1), affine_tol) ||
100 589657 : !v.relative_fuzzy_equals(this->point(5) - this->point(2), affine_tol))
101 3072 : return false;
102 : // Make sure edges are straight
103 43720 : v /= 2;
104 545241 : if (!v.relative_fuzzy_equals(this->point(9) - this->point(0), affine_tol) ||
105 1003346 : !v.relative_fuzzy_equals(this->point(10) - this->point(1), affine_tol) ||
106 545025 : !v.relative_fuzzy_equals(this->point(11) - this->point(2), affine_tol))
107 672 : return false;
108 544953 : v = (this->point(1) - this->point(0))/2;
109 588609 : if (!v.relative_fuzzy_equals(this->point(6) - this->point(0), affine_tol) ||
110 544905 : !v.relative_fuzzy_equals(this->point(12) - this->point(3), affine_tol))
111 48 : return false;
112 544901 : v = (this->point(2) - this->point(0))/2;
113 588553 : if (!v.relative_fuzzy_equals(this->point(8) - this->point(0), affine_tol) ||
114 544853 : !v.relative_fuzzy_equals(this->point(14) - this->point(3), affine_tol))
115 48 : return false;
116 544849 : v = (this->point(2) - this->point(1))/2;
117 588505 : if (!v.relative_fuzzy_equals(this->point(7) - this->point(1), affine_tol) ||
118 588505 : !v.relative_fuzzy_equals(this->point(13) - this->point(4), affine_tol))
119 0 : return false;
120 43656 : return true;
121 : }
122 :
123 :
124 :
125 15188907 : Order Prism15::default_order() const
126 : {
127 15188907 : return SECOND;
128 : }
129 :
130 :
131 :
132 24046985 : unsigned int Prism15::local_side_node(unsigned int side,
133 : unsigned int side_node) const
134 : {
135 2146890 : libmesh_assert_less (side, this->n_sides());
136 :
137 : // Never more than 8 nodes per side.
138 2146890 : libmesh_assert_less(side_node, Prism15::nodes_per_side);
139 :
140 : // Some sides have 6 nodes.
141 2146890 : libmesh_assert(!(side==0 || side==4) || side_node < 6);
142 :
143 24046985 : return Prism15::side_nodes_map[side][side_node];
144 : }
145 :
146 :
147 :
148 8337357 : unsigned int Prism15::local_edge_node(unsigned int edge,
149 : unsigned int edge_node) const
150 : {
151 742482 : libmesh_assert_less(edge, this->n_edges());
152 742482 : libmesh_assert_less(edge_node, Prism15::nodes_per_edge);
153 :
154 8337357 : return Prism15::edge_nodes_map[edge][edge_node];
155 : }
156 :
157 :
158 :
159 512691 : std::unique_ptr<Elem> Prism15::build_side_ptr (const unsigned int i)
160 : {
161 20070 : libmesh_assert_less (i, this->n_sides());
162 :
163 512691 : std::unique_ptr<Elem> face;
164 :
165 512691 : switch (i)
166 : {
167 300432 : case 0: // the triangular face at z=-1
168 : case 4: // the triangular face at z=1
169 : {
170 300432 : face = std::make_unique<Tri6>();
171 300432 : break;
172 : }
173 212259 : case 1: // the quad face at y=0
174 : case 2: // the other quad face
175 : case 3: // the quad face at x=0
176 : {
177 212259 : face = std::make_unique<Quad8>();
178 212259 : break;
179 : }
180 0 : default:
181 0 : libmesh_error_msg("Invalid side i = " << i);
182 : }
183 :
184 : // Set the nodes
185 4013355 : for (auto n : face->node_index_range())
186 3640812 : face->set_node(n, this->node_ptr(Prism15::side_nodes_map[i][n]));
187 :
188 512691 : face->set_interior_parent(this);
189 492621 : face->inherit_data_from(*this);
190 :
191 512691 : return face;
192 0 : }
193 :
194 :
195 479459 : void Prism15::build_side_ptr (std::unique_ptr<Elem> & side,
196 : const unsigned int i)
197 : {
198 11066 : libmesh_assert_less (i, this->n_sides());
199 :
200 479459 : switch (i)
201 : {
202 2804 : case 0: // the triangular face at z=-1
203 : case 4: // the triangular face at z=1
204 : {
205 120350 : if (!side.get() || side->type() != TRI6)
206 : {
207 233820 : side = this->build_side_ptr(i);
208 119652 : return;
209 : }
210 62 : break;
211 : }
212 :
213 8262 : case 1: // the quad face at y=0
214 : case 2: // the other quad face
215 : case 3: // the quad face at x=0
216 : {
217 359109 : if (!side.get() || side->type() != QUAD8)
218 : {
219 233330 : side = this->build_side_ptr(i);
220 119389 : return;
221 : }
222 5538 : break;
223 : }
224 :
225 0 : default:
226 0 : libmesh_error_msg("Invalid side i = " << i);
227 : }
228 :
229 234818 : side->inherit_data_from(*this);
230 :
231 : // Set the nodes
232 2162366 : for (auto n : side->node_index_range())
233 1966624 : side->set_node(n, this->node_ptr(Prism15::side_nodes_map[i][n]));
234 : }
235 :
236 :
237 :
238 755441 : std::unique_ptr<Elem> Prism15::build_edge_ptr (const unsigned int i)
239 : {
240 755441 : return this->simple_build_edge_ptr<Edge3,Prism15>(i);
241 : }
242 :
243 :
244 :
245 0 : void Prism15::build_edge_ptr (std::unique_ptr<Elem> & edge, const unsigned int i)
246 : {
247 0 : this->simple_build_edge_ptr<Prism15>(edge, i, EDGE3);
248 0 : }
249 :
250 :
251 :
252 0 : void Prism15::connectivity(const unsigned int libmesh_dbg_var(sc),
253 : const IOPackage iop,
254 : std::vector<dof_id_type> & conn) const
255 : {
256 0 : libmesh_assert(_nodes);
257 0 : libmesh_assert_less (sc, this->n_sub_elem());
258 0 : libmesh_assert_not_equal_to (iop, INVALID_IO_PACKAGE);
259 :
260 0 : switch (iop)
261 : {
262 0 : case TECPLOT:
263 : {
264 0 : conn.resize(8);
265 0 : conn[0] = this->node_id(0)+1;
266 0 : conn[1] = this->node_id(1)+1;
267 0 : conn[2] = this->node_id(2)+1;
268 0 : conn[3] = this->node_id(2)+1;
269 0 : conn[4] = this->node_id(3)+1;
270 0 : conn[5] = this->node_id(4)+1;
271 0 : conn[6] = this->node_id(5)+1;
272 0 : conn[7] = this->node_id(5)+1;
273 0 : return;
274 : }
275 :
276 0 : case VTK:
277 : {
278 : // VTK's VTK_QUADRATIC_WEDGE first 9 nodes match, then their
279 : // middle and top layers of mid-edge nodes are reversed from
280 : // LibMesh's.
281 0 : conn.resize(15);
282 0 : for (unsigned i=0; i<9; ++i)
283 0 : conn[i] = this->node_id(i);
284 :
285 : // top "ring" of mid-edge nodes
286 0 : conn[9] = this->node_id(12);
287 0 : conn[10] = this->node_id(13);
288 0 : conn[11] = this->node_id(14);
289 :
290 : // middle "ring" of mid-edge nodes
291 0 : conn[12] = this->node_id(9);
292 0 : conn[13] = this->node_id(10);
293 0 : conn[14] = this->node_id(11);
294 :
295 0 : return;
296 : }
297 :
298 0 : default:
299 0 : libmesh_error_msg("Unsupported IO package " << iop);
300 : }
301 : }
302 :
303 :
304 :
305 :
306 0 : unsigned short int Prism15::second_order_adjacent_vertex (const unsigned int n,
307 : const unsigned int v) const
308 : {
309 0 : libmesh_assert_greater_equal (n, this->n_vertices());
310 0 : libmesh_assert_less (n, this->n_nodes());
311 0 : libmesh_assert_less (v, 2);
312 0 : return _second_order_adjacent_vertices[n-this->n_vertices()][v];
313 : }
314 :
315 :
316 :
317 : std::pair<unsigned short int, unsigned short int>
318 0 : Prism15::second_order_child_vertex (const unsigned int n) const
319 : {
320 0 : libmesh_assert_greater_equal (n, this->n_vertices());
321 0 : libmesh_assert_less (n, this->n_nodes());
322 :
323 0 : return std::pair<unsigned short int, unsigned short int>
324 0 : (_second_order_vertex_child_number[n],
325 0 : _second_order_vertex_child_index[n]);
326 : }
327 :
328 :
329 :
330 264 : Real Prism15::volume () const
331 : {
332 : // This specialization is good for Lagrange mappings only in general
333 264 : if (this->mapping_type() != LAGRANGE_MAP)
334 0 : return this->Elem::volume();
335 :
336 : // Make copies of our points. It makes the subsequent calculations a bit
337 : // shorter and avoids dereferencing the same pointer multiple times.
338 : Point
339 354 : x0 = point(0), x1 = point(1), x2 = point(2), x3 = point(3), x4 = point(4),
340 336 : x5 = point(5), x6 = point(6), x7 = point(7), x8 = point(8), x9 = point(9),
341 336 : x10 = point(10), x11 = point(11), x12 = point(12), x13 = point(13), x14 = point(14);
342 :
343 : // Terms are copied directly from a Python script.
344 : Point dx_dxi[10] =
345 : {
346 20 : -x0 - x1 + x10 + 2*x12 - x3 - x4 + 2*x6 - x9,
347 20 : 3*x0/2 + x1/2 + 2*x12 - 3*x3/2 - x4/2 - 2*x6,
348 20 : -x0/2 + x1/2 - x10 - x3/2 + x4/2 + x9,
349 20 : 2*x0 - 2*x12 + 2*x13 - 2*x14 + 2*x3 - 2*x6 + 2*x7 - 2*x8,
350 20 : -2*x0 - 2*x12 + 2*x13 - 2*x14 + 2*x3 + 2*x6 - 2*x7 + 2*x8,
351 : Point(0,0,0),
352 20 : 2*x0 + 2*x1 - 4*x12 + 2*x3 + 2*x4 - 4*x6,
353 20 : -2*x0 - 2*x1 - 4*x12 + 2*x3 + 2*x4 + 4*x6,
354 : Point(0,0,0),
355 : Point(0,0,0)
356 160 : };
357 :
358 : Point dx_deta[10] =
359 : {
360 20 : -x0 + x11 + 2*x14 - x2 - x3 - x5 + 2*x8 - x9,
361 20 : 3*x0/2 + 2*x14 + x2/2 - 3*x3/2 - x5/2 - 2*x8,
362 20 : -x0/2 - x11 + x2/2 - x3/2 + x5/2 + x9,
363 20 : 2*x0 - 4*x14 + 2*x2 + 2*x3 + 2*x5 - 4*x8,
364 20 : -2*x0 - 4*x14 - 2*x2 + 2*x3 + 2*x5 + 4*x8,
365 : Point(0,0,0),
366 20 : 2*x0 - 2*x12 + 2*x13 - 2*x14 + 2*x3 - 2*x6 + 2*x7 - 2*x8,
367 20 : -2*x0 - 2*x12 + 2*x13 - 2*x14 + 2*x3 + 2*x6 - 2*x7 + 2*x8,
368 : Point(0,0,0),
369 : Point(0,0,0)
370 160 : };
371 :
372 : Point dx_dzeta[10] =
373 : {
374 20 : -x0/2 + x3/2,
375 20 : x0 + x3 - 2*x9,
376 : Point(0,0,0),
377 20 : 3*x0/2 + 2*x14 + x2/2 - 3*x3/2 - x5/2 - 2*x8,
378 20 : -x0 - 2*x11 + x2 - x3 + x5 + 2*x9,
379 20 : -x0 - 2*x14 - x2 + x3 + x5 + 2*x8,
380 20 : 3*x0/2 + x1/2 + 2*x12 - 3*x3/2 - x4/2 - 2*x6,
381 20 : -x0 + x1 - 2*x10 - x3 + x4 + 2*x9,
382 20 : -2*x0 - 2*x12 + 2*x13 - 2*x14 + 2*x3 + 2*x6 - 2*x7 + 2*x8,
383 20 : -x0 - x1 - 2*x12 + x3 + x4 + 2*x6
384 200 : };
385 :
386 : // The quadrature rule for the Prism15 is a tensor product between a
387 : // FOURTH-order TRI3 rule (in xi, eta) and a FIFTH-order EDGE2 rule
388 : // in zeta.
389 :
390 : // Number of points in the 2D quadrature rule.
391 20 : const int N2D = 6;
392 :
393 : // Parameters of the 2D rule
394 : static const Real
395 : w1 = 1.1169079483900573284750350421656140e-01_R,
396 : w2 = 5.4975871827660933819163162450105264e-02_R,
397 : a1 = 4.4594849091596488631832925388305199e-01_R,
398 : a2 = 9.1576213509770743459571463402201508e-02_R;
399 :
400 : // Points and weights of the 2D rule
401 : static const Real w2D[N2D] = {w1, w1, w1, w2, w2, w2};
402 :
403 : // Quadrature point locations raised to powers. xi[0][2] is
404 : // quadrature point 0, squared, xi[1][1] is quadrature point 1 to the
405 : // first power, etc. This lets us avoid calling std::pow inside the
406 : // loops below.
407 : static const Real xi[N2D][3] =
408 : {
409 : // ^0 ^1 ^2
410 : { 1., a1, a1*a1},
411 : { 1., 1-2*a1, (1-2*a1)*(1-2*a1)},
412 : { 1., a1, a1*a1},
413 : { 1., a2, a2*a2},
414 : { 1., 1-2*a2, (1-2*a2)*(1-2*a2)},
415 : { 1., a2, a2*a2}
416 : };
417 :
418 : static const Real eta[N2D][3] =
419 : {
420 : // ^0 ^1 ^2
421 : { 1., a1, a1*a1},
422 : { 1., a1, a1*a1},
423 : { 1., 1-2*a1, (1-2*a1)*(1-2*a1)},
424 : { 1., a2, a2*a2},
425 : { 1., a2, a2*a2},
426 : { 1., 1-2*a2, (1-2*a2)*(1-2*a2)}
427 : };
428 :
429 : // Number of points in the 1D quadrature rule.
430 20 : const int N1D = 3;
431 :
432 : // Points and weights of the 1D quadrature rule.
433 : static const Real w1D[N1D] = {5./9, 8./9, 5./9};
434 :
435 264 : const Real zeta[N1D][3] =
436 : {
437 : //^0 ^1 ^2
438 : { 1., -std::sqrt(15)/5., 15./25},
439 : { 1., 0., 0.},
440 : { 1., std::sqrt(15)/5., 15./25}
441 : };
442 :
443 : // The integer exponents for each term.
444 : static const int exponents[10][3] =
445 : {
446 : {0, 0, 0},
447 : {0, 0, 1},
448 : {0, 0, 2},
449 : {0, 1, 0},
450 : {0, 1, 1},
451 : {0, 2, 0},
452 : {1, 0, 0},
453 : {1, 0, 1},
454 : {1, 1, 0},
455 : {2, 0, 0}
456 : };
457 :
458 20 : Real vol = 0.;
459 1848 : for (int i=0; i<N2D; ++i)
460 6336 : for (int j=0; j<N1D; ++j)
461 : {
462 : // Compute dx_dxi, dx_deta, dx_dzeta at the current quadrature point.
463 360 : Point dx_dxi_q, dx_deta_q, dx_dzeta_q;
464 52272 : for (int c=0; c<10; ++c)
465 : {
466 3600 : Real coeff =
467 54360 : xi[i][exponents[c][0]]*
468 54360 : eta[i][exponents[c][1]]*
469 47520 : zeta[j][exponents[c][2]];
470 :
471 3600 : dx_dxi_q += coeff * dx_dxi[c];
472 3600 : dx_deta_q += coeff * dx_deta[c];
473 3600 : dx_dzeta_q += coeff * dx_dzeta[c];
474 : }
475 :
476 : // Compute scalar triple product, multiply by weight, and accumulate volume.
477 5076 : vol += w2D[i] * w1D[j] * triple_product(dx_dxi_q, dx_deta_q, dx_dzeta_q);
478 : }
479 :
480 20 : return vol;
481 : }
482 :
483 :
484 :
485 : #ifdef LIBMESH_ENABLE_AMR
486 :
487 : const Real Prism15::_embedding_matrix[Prism15::num_children][Prism15::num_nodes][Prism15::num_nodes] =
488 : {
489 : // Embedding matrix for child 0
490 : {
491 : // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
492 : { 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, // 0
493 : { 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0 }, // 1
494 : { 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0 }, // 2
495 : { 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0 }, // 3
496 : { -0.25, -0.25, 0, -0.25, -0.25, 0, 0.5, 0, 0, 0.5, 0.5, 0, 0.5, 0, 0 }, // 4
497 : { -0.25, 0, -0.25, -0.25, 0, -0.25, 0, 0, 0.5, 0.5, 0, 0.5, 0, 0, 0.5 }, // 5
498 : { 0.375, -0.125, 0, 0, 0, 0, 0.75, 0, 0, 0, 0, 0, 0, 0, 0 }, // 6
499 : { 0, -0.125, -0.125, 0, 0, 0, 0.5, 0.25, 0.5, 0, 0, 0, 0, 0, 0 }, // 7
500 : { 0.375, 0, -0.125, 0, 0, 0, 0, 0, 0.75, 0, 0, 0, 0, 0, 0 }, // 8
501 : { 0.375, 0, 0, -0.125, 0, 0, 0, 0, 0, 0.75, 0, 0, 0, 0, 0 }, // 9
502 : { -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.75, 0, 0, 0.375, 0.375, 0, 0.25, 0, 0 }, // 10
503 : { -0.1875, 0, -0.1875, -0.1875, 0, -0.1875, 0, 0, 0.75, 0.375, 0, 0.375, 0, 0, 0.25 }, // 11
504 : { -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.375, 0, 0, 0.75, 0.25, 0, 0.375, 0, 0 }, // 12
505 : { -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
506 : { -0.1875, 0, -0.1875, -0.1875, 0, -0.1875, 0, 0, 0.375, 0.75, 0, 0.25, 0, 0, 0.375 } // 14
507 : },
508 :
509 : // Embedding matrix for child 1
510 : {
511 : // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
512 : { 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0 }, // 0
513 : { 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, // 1
514 : { 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0 }, // 2
515 : { -0.25, -0.25, 0, -0.25, -0.25, 0, 0.5, 0, 0, 0.5, 0.5, 0, 0.5, 0, 0 }, // 3
516 : { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0 }, // 4
517 : { 0, -0.25, -0.25, 0, -0.25, -0.25, 0, 0.5, 0, 0, 0.5, 0.5, 0, 0.5, 0 }, // 5
518 : { -0.125, 0.375, 0, 0, 0, 0, 0.75, 0, 0, 0, 0, 0, 0, 0, 0 }, // 6
519 : { 0, 0.375, -0.125, 0, 0, 0, 0, 0.75, 0, 0, 0, 0, 0, 0, 0 }, // 7
520 : { -0.125, 0, -0.125, 0, 0, 0, 0.5, 0.5, 0.25, 0, 0, 0, 0, 0, 0 }, // 8
521 : { -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.75, 0, 0, 0.375, 0.375, 0, 0.25, 0, 0 }, // 9
522 : { 0, 0.375, 0, 0, -0.125, 0, 0, 0, 0, 0, 0.75, 0, 0, 0, 0 }, // 10
523 : { 0, -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.75, 0, 0, 0.375, 0.375, 0, 0.25, 0 }, // 11
524 : { -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.375, 0, 0, 0.25, 0.75, 0, 0.375, 0, 0 }, // 12
525 : { 0, -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.375, 0, 0, 0.75, 0.25, 0, 0.375, 0 }, // 13
526 : { -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
527 : },
528 :
529 : // Embedding matrix for child 2
530 : {
531 : // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
532 : { 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0 }, // 0
533 : { 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0 }, // 1
534 : { 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, // 2
535 : { -0.25, 0, -0.25, -0.25, 0, -0.25, 0, 0, 0.5, 0.5, 0, 0.5, 0, 0, 0.5 }, // 3
536 : { 0, -0.25, -0.25, 0, -0.25, -0.25, 0, 0.5, 0, 0, 0.5, 0.5, 0, 0.5, 0 }, // 4
537 : { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0 }, // 5
538 : { -0.125, -0.125, 0, 0, 0, 0, 0.25, 0.5, 0.5, 0, 0, 0, 0, 0, 0 }, // 6
539 : { 0, -0.125, 0.375, 0, 0, 0, 0, 0.75, 0, 0, 0, 0, 0, 0, 0 }, // 7
540 : { -0.125, 0, 0.375, 0, 0, 0, 0, 0, 0.75, 0, 0, 0, 0, 0, 0 }, // 8
541 : { -0.1875, 0, -0.1875, -0.1875, 0, -0.1875, 0, 0, 0.75, 0.375, 0, 0.375, 0, 0, 0.25 }, // 9
542 : { 0, -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.75, 0, 0, 0.375, 0.375, 0, 0.25, 0 }, // 10
543 : { 0, 0, 0.375, 0, 0, -0.125, 0, 0, 0, 0, 0, 0.75, 0, 0, 0 }, // 11
544 : { -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
545 : { 0, -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.375, 0, 0, 0.25, 0.75, 0, 0.375, 0 }, // 13
546 : { -0.1875, 0, -0.1875, -0.1875, 0, -0.1875, 0, 0, 0.375, 0.25, 0, 0.75, 0, 0, 0.375 } // 14
547 : },
548 :
549 : // Embedding matrix for child 3
550 : {
551 : // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
552 : { 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0 }, // 0
553 : { 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0 }, // 1
554 : { 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0 }, // 2
555 : { -0.25, -0.25, 0, -0.25, -0.25, 0, 0.5, 0, 0, 0.5, 0.5, 0, 0.5, 0, 0 }, // 3
556 : { 0, -0.25, -0.25, 0, -0.25, -0.25, 0, 0.5, 0, 0, 0.5, 0.5, 0, 0.5, 0 }, // 4
557 : { -0.25, 0, -0.25, -0.25, 0, -0.25, 0, 0, 0.5, 0.5, 0, 0.5, 0, 0, 0.5 }, // 5
558 : { -0.125, 0, -0.125, 0, 0, 0, 0.5, 0.5, 0.25, 0, 0, 0, 0, 0, 0 }, // 6
559 : { -0.125, -0.125, 0, 0, 0, 0, 0.25, 0.5, 0.5, 0, 0, 0, 0, 0, 0 }, // 7
560 : { 0, -0.125, -0.125, 0, 0, 0, 0.5, 0.25, 0.5, 0, 0, 0, 0, 0, 0 }, // 8
561 : { -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.75, 0, 0, 0.375, 0.375, 0, 0.25, 0, 0 }, // 9
562 : { 0, -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.75, 0, 0, 0.375, 0.375, 0, 0.25, 0 }, // 10
563 : { -0.1875, 0, -0.1875, -0.1875, 0, -0.1875, 0, 0, 0.75, 0.375, 0, 0.375, 0, 0, 0.25 }, // 11
564 : { -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
565 : { -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
566 : { -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
567 : },
568 :
569 : // Embedding matrix for child 4
570 : {
571 : // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
572 : { 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0 }, // 0
573 : { -0.25, -0.25, 0, -0.25, -0.25, 0, 0.5, 0, 0, 0.5, 0.5, 0, 0.5, 0, 0 }, // 1
574 : { -0.25, 0, -0.25, -0.25, 0, -0.25, 0, 0, 0.5, 0.5, 0, 0.5, 0, 0, 0.5 }, // 2
575 : { 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, // 3
576 : { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0 }, // 4
577 : { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1 }, // 5
578 : { -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.375, 0, 0, 0.75, 0.25, 0, 0.375, 0, 0 }, // 6
579 : { -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
580 : { -0.1875, 0, -0.1875, -0.1875, 0, -0.1875, 0, 0, 0.375, 0.75, 0, 0.25, 0, 0, 0.375 }, // 8
581 : { -0.125, 0, 0, 0.375, 0, 0, 0, 0, 0, 0.75, 0, 0, 0, 0, 0 }, // 9
582 : { -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.25, 0, 0, 0.375, 0.375, 0, 0.75, 0, 0 }, // 10
583 : { -0.1875, 0, -0.1875, -0.1875, 0, -0.1875, 0, 0, 0.25, 0.375, 0, 0.375, 0, 0, 0.75 }, // 11
584 : { 0, 0, 0, 0.375, -0.125, 0, 0, 0, 0, 0, 0, 0, 0.75, 0, 0 }, // 12
585 : { 0, 0, 0, 0, -0.125, -0.125, 0, 0, 0, 0, 0, 0, 0.5, 0.25, 0.5 }, // 13
586 : { 0, 0, 0, 0.375, 0, -0.125, 0, 0, 0, 0, 0, 0, 0, 0, 0.75 } // 14
587 : },
588 :
589 : // Embedding matrix for child 5
590 : {
591 : // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
592 : { -0.25, -0.25, 0, -0.25, -0.25, 0, 0.5, 0, 0, 0.5, 0.5, 0, 0.5, 0, 0 }, // 0
593 : { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0 }, // 1
594 : { 0, -0.25, -0.25, 0, -0.25, -0.25, 0, 0.5, 0, 0, 0.5, 0.5, 0, 0.5, 0 }, // 2
595 : { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0 }, // 3
596 : { 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, // 4
597 : { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0 }, // 5
598 : { -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.375, 0, 0, 0.25, 0.75, 0, 0.375, 0, 0 }, // 6
599 : { 0, -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.375, 0, 0, 0.75, 0.25, 0, 0.375, 0 }, // 7
600 : { -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
601 : { -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.25, 0, 0, 0.375, 0.375, 0, 0.75, 0, 0 }, // 9
602 : { 0, -0.125, 0, 0, 0.375, 0, 0, 0, 0, 0, 0.75, 0, 0, 0, 0 }, // 10
603 : { 0, -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.25, 0, 0, 0.375, 0.375, 0, 0.75, 0 }, // 11
604 : { 0, 0, 0, -0.125, 0.375, 0, 0, 0, 0, 0, 0, 0, 0.75, 0, 0 }, // 12
605 : { 0, 0, 0, 0, 0.375, -0.125, 0, 0, 0, 0, 0, 0, 0, 0.75, 0 }, // 13
606 : { 0, 0, 0, -0.125, 0, -0.125, 0, 0, 0, 0, 0, 0, 0.5, 0.5, 0.25 } // 14
607 : },
608 :
609 : // Embedding matrix for child 6
610 : {
611 : // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
612 : { -0.25, 0, -0.25, -0.25, 0, -0.25, 0, 0, 0.5, 0.5, 0, 0.5, 0, 0, 0.5 }, // 0
613 : { 0, -0.25, -0.25, 0, -0.25, -0.25, 0, 0.5, 0, 0, 0.5, 0.5, 0, 0.5, 0 }, // 1
614 : { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0 }, // 2
615 : { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1 }, // 3
616 : { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0 }, // 4
617 : { 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, // 5
618 : { -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
619 : { 0, -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.375, 0, 0, 0.25, 0.75, 0, 0.375, 0 }, // 7
620 : { -0.1875, 0, -0.1875, -0.1875, 0, -0.1875, 0, 0, 0.375, 0.25, 0, 0.75, 0, 0, 0.375 }, // 8
621 : { -0.1875, 0, -0.1875, -0.1875, 0, -0.1875, 0, 0, 0.25, 0.375, 0, 0.375, 0, 0, 0.75 }, // 9
622 : { 0, -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.25, 0, 0, 0.375, 0.375, 0, 0.75, 0 }, // 10
623 : { 0, 0, -0.125, 0, 0, 0.375, 0, 0, 0, 0, 0, 0.75, 0, 0, 0 }, // 11
624 : { 0, 0, 0, -0.125, -0.125, 0, 0, 0, 0, 0, 0, 0, 0.25, 0.5, 0.5 }, // 12
625 : { 0, 0, 0, 0, -0.125, 0.375, 0, 0, 0, 0, 0, 0, 0, 0.75, 0 }, // 13
626 : { 0, 0, 0, -0.125, 0, 0.375, 0, 0, 0, 0, 0, 0, 0, 0, 0.75 } // 14
627 : },
628 :
629 : // Embedding matrix for child 7
630 : {
631 : // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
632 : { -0.25, -0.25, 0, -0.25, -0.25, 0, 0.5, 0, 0, 0.5, 0.5, 0, 0.5, 0, 0 }, // 0
633 : { 0, -0.25, -0.25, 0, -0.25, -0.25, 0, 0.5, 0, 0, 0.5, 0.5, 0, 0.5, 0 }, // 1
634 : { -0.25, 0, -0.25, -0.25, 0, -0.25, 0, 0, 0.5, 0.5, 0, 0.5, 0, 0, 0.5 }, // 2
635 : { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0 }, // 3
636 : { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0 }, // 4
637 : { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1 }, // 5
638 : { -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
639 : { -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
640 : { -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
641 : { -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.25, 0, 0, 0.375, 0.375, 0, 0.75, 0, 0 }, // 9
642 : { 0, -0.1875, -0.1875, 0, -0.1875, -0.1875, 0, 0.25, 0, 0, 0.375, 0.375, 0, 0.75, 0 }, // 10
643 : { -0.1875, 0, -0.1875, -0.1875, 0, -0.1875, 0, 0, 0.25, 0.375, 0, 0.375, 0, 0, 0.75 }, // 11
644 : { 0, 0, 0, -0.125, 0, -0.125, 0, 0, 0, 0, 0, 0, 0.5, 0.5, 0.25 }, // 12
645 : { 0, 0, 0, -0.125, -0.125, 0, 0, 0, 0, 0, 0, 0, 0.25, 0.5, 0.5 }, // 13
646 : { 0, 0, 0, 0, -0.125, -0.125, 0, 0, 0, 0, 0, 0, 0.5, 0.25, 0.5 } // 14
647 : }
648 : };
649 :
650 : #endif
651 :
652 :
653 : void
654 8736 : Prism15::permute(unsigned int perm_num)
655 : {
656 768 : libmesh_assert_less (perm_num, 6);
657 8736 : const unsigned int side = perm_num % 2;
658 8736 : const unsigned int rotate = perm_num / 2;
659 :
660 20112 : for (unsigned int i = 0; i != rotate; ++i)
661 : {
662 11376 : swap3nodes(0,1,2);
663 10368 : swap3nodes(3,4,5);
664 10368 : swap3nodes(6,7,8);
665 10368 : swap3nodes(9,10,11);
666 10368 : swap3nodes(12,13,14);
667 10368 : swap3neighbors(1,2,3);
668 : }
669 :
670 8736 : switch (side) {
671 144 : case 0:
672 144 : break;
673 7008 : case 1:
674 7008 : swap2nodes(1,3);
675 7008 : swap2nodes(0,4);
676 7008 : swap2nodes(2,5);
677 7008 : swap2nodes(6,12);
678 7008 : swap2nodes(9,10);
679 7008 : swap2nodes(7,14);
680 7008 : swap2nodes(8,13);
681 624 : swap2neighbors(0,4);
682 624 : swap2neighbors(2,3);
683 624 : break;
684 0 : default:
685 0 : libmesh_error();
686 : }
687 :
688 8736 : }
689 :
690 :
691 : void
692 576 : Prism15::flip(BoundaryInfo * boundary_info)
693 : {
694 48 : libmesh_assert(boundary_info);
695 :
696 576 : swap2nodes(0,1);
697 576 : swap2nodes(3,4);
698 576 : swap2nodes(7,8);
699 576 : swap2nodes(9,10);
700 576 : swap2nodes(13,14);
701 48 : swap2neighbors(2,3);
702 576 : swap2boundarysides(2,3,boundary_info);
703 576 : swap2boundaryedges(0,1,boundary_info);
704 576 : swap2boundaryedges(3,4,boundary_info);
705 576 : swap2boundaryedges(7,8,boundary_info);
706 576 : }
707 :
708 :
709 : ElemType
710 8036172 : Prism15::side_type (const unsigned int s) const
711 : {
712 717416 : libmesh_assert_less (s, 5);
713 8036172 : if (s == 0 || s == 4)
714 10056 : return TRI6;
715 716576 : return QUAD8;
716 : }
717 :
718 : } // namespace libMesh
|
