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