libMesh
Loading...
Searching...
No Matches
cell_prism15.C
Go to the documentation of this file.
1// The libMesh Finite Element Library.
2// Copyright (C) 2002-2026 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3
4// This library is free software; you can redistribute it and/or
5// modify it under the terms of the GNU Lesser General Public
6// License as published by the Free Software Foundation; either
7// version 2.1 of the License, or (at your option) any later version.
8
9// This library is distributed in the hope that it will be useful,
10// but WITHOUT ANY WARRANTY; without even the implied warranty of
11// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12// Lesser General Public License for more details.
13
14// You should have received a copy of the GNU Lesser General Public
15// License along with this library; if not, write to the Free Software
16// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17
18
19// Local includes
20#include "libmesh/cell_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
27namespace libMesh
28{
29
30
31
32// ------------------------------------------------------------
33// Prism15 class static member initializations
34const int Prism15::num_nodes;
35const int Prism15::nodes_per_side;
36const int Prism15::nodes_per_edge;
37
38const 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
47const 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
63bool Prism15::is_vertex(const unsigned int i) const
64{
65 if (i < 6)
66 return true;
67 return false;
68}
69
70bool Prism15::is_edge(const unsigned int i) const
71{
72 if (i < 6)
73 return false;
74 return true;
75}
76
77bool Prism15::is_face(const unsigned int) const
78{
79 return false;
80}
81
82bool Prism15::is_node_on_side(const unsigned int n,
83 const unsigned int s) const
84{
85 libmesh_assert_less (s, n_sides());
86 return std::find(std::begin(side_nodes_map[s]),
87 std::end(side_nodes_map[s]),
88 n) != std::end(side_nodes_map[s]);
89}
90
91std::vector<unsigned int>
92Prism15::nodes_on_side(const unsigned int s) const
93{
94 libmesh_assert_less(s, n_sides());
95 auto trim = (s > 0 && s < 4) ? 0 : 2;
96 return {std::begin(side_nodes_map[s]), std::end(side_nodes_map[s]) - trim};
97}
98
99std::vector<unsigned>
100Prism15::nodes_on_edge(const unsigned int e) const
101{
102 libmesh_assert_less(e, n_edges());
103 return {std::begin(edge_nodes_map[e]), std::end(edge_nodes_map[e])};
104}
105
106bool Prism15::is_node_on_edge(const unsigned int n,
107 const unsigned int e) const
108{
109 libmesh_assert_less (e, n_edges());
110 return std::find(std::begin(edge_nodes_map[e]),
111 std::end(edge_nodes_map[e]),
112 n) != std::end(edge_nodes_map[e]);
113}
114
115
116
117bool Prism15::has_affine_map() const
118{
119 // Make sure z edges are affine
120 Point v = this->point(3) - this->point(0);
121 if (!v.relative_fuzzy_equals(this->point(4) - this->point(1), affine_tol) ||
122 !v.relative_fuzzy_equals(this->point(5) - this->point(2), affine_tol))
123 return false;
124 // Make sure edges are straight
125 v /= 2;
126 if (!v.relative_fuzzy_equals(this->point(9) - this->point(0), affine_tol) ||
127 !v.relative_fuzzy_equals(this->point(10) - this->point(1), affine_tol) ||
128 !v.relative_fuzzy_equals(this->point(11) - this->point(2), affine_tol))
129 return false;
130 v = (this->point(1) - this->point(0))/2;
131 if (!v.relative_fuzzy_equals(this->point(6) - this->point(0), affine_tol) ||
132 !v.relative_fuzzy_equals(this->point(12) - this->point(3), affine_tol))
133 return false;
134 v = (this->point(2) - this->point(0))/2;
135 if (!v.relative_fuzzy_equals(this->point(8) - this->point(0), affine_tol) ||
136 !v.relative_fuzzy_equals(this->point(14) - this->point(3), affine_tol))
137 return false;
138 v = (this->point(2) - this->point(1))/2;
139 if (!v.relative_fuzzy_equals(this->point(7) - this->point(1), affine_tol) ||
140 !v.relative_fuzzy_equals(this->point(13) - this->point(4), affine_tol))
141 return false;
142 return true;
143}
144
145
146
147Order Prism15::default_order() const
148{
149 return SECOND;
150}
151
152
153
154unsigned int Prism15::local_side_node(unsigned int side,
155 unsigned int side_node) const
156{
157 libmesh_assert_less (side, this->n_sides());
158
159 // Never more than 8 nodes per side.
160 libmesh_assert_less(side_node, Prism15::nodes_per_side);
161
162 // Some sides have 6 nodes.
163 libmesh_assert(!(side==0 || side==4) || side_node < 6);
164
165 return Prism15::side_nodes_map[side][side_node];
166}
167
168
169
170unsigned int Prism15::local_edge_node(unsigned int edge,
171 unsigned int edge_node) const
172{
173 libmesh_assert_less(edge, this->n_edges());
174 libmesh_assert_less(edge_node, Prism15::nodes_per_edge);
175
176 return Prism15::edge_nodes_map[edge][edge_node];
177}
178
179
180
181std::unique_ptr<Elem> Prism15::build_side_ptr (const unsigned int i)
182{
183 libmesh_assert_less (i, this->n_sides());
184
185 std::unique_ptr<Elem> face;
186
187 switch (i)
188 {
189 case 0: // the triangular face at z=-1
190 case 4: // the triangular face at z=1
191 {
192 face = std::make_unique<Tri6>();
193 break;
194 }
195 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 face = std::make_unique<Quad8>();
200 break;
201 }
202 default:
203 libmesh_error_msg("Invalid side i = " << i);
204 }
205
206 // Set the nodes
207 for (auto n : face->node_index_range())
208 face->set_node(n, this->node_ptr(Prism15::side_nodes_map[i][n]));
209
210 face->set_interior_parent(this);
211 face->inherit_data_from(*this);
212
213 return face;
214}
215
216
217void Prism15::build_side_ptr (std::unique_ptr<Elem> & side,
218 const unsigned int i)
219{
220 libmesh_assert_less (i, this->n_sides());
221
222 switch (i)
223 {
224 case 0: // the triangular face at z=-1
225 case 4: // the triangular face at z=1
226 {
227 if (!side.get() || side->type() != TRI6)
228 {
229 side = this->build_side_ptr(i);
230 return;
231 }
232 break;
233 }
234
235 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 if (!side.get() || side->type() != QUAD8)
240 {
241 side = this->build_side_ptr(i);
242 return;
243 }
244 break;
245 }
246
247 default:
248 libmesh_error_msg("Invalid side i = " << i);
249 }
250
251 side->inherit_data_from(*this);
252
253 // Set the nodes
254 for (auto n : side->node_index_range())
255 side->set_node(n, this->node_ptr(Prism15::side_nodes_map[i][n]));
256}
257
258
259
260std::unique_ptr<Elem> Prism15::build_edge_ptr (const unsigned int i)
261{
262 return this->simple_build_edge_ptr<Edge3,Prism15>(i);
263}
264
265
266
267void Prism15::build_edge_ptr (std::unique_ptr<Elem> & edge, const unsigned int i)
268{
269 this->simple_build_edge_ptr<Prism15>(edge, i, EDGE3);
270}
271
272
273
274void Prism15::connectivity(const unsigned int libmesh_dbg_var(sc),
275 const IOPackage iop,
276 std::vector<dof_id_type> & conn) const
277{
278 libmesh_assert(_nodes);
279 libmesh_assert_less (sc, this->n_sub_elem());
280 libmesh_assert_not_equal_to (iop, INVALID_IO_PACKAGE);
281
282 switch (iop)
283 {
284 case TECPLOT:
285 {
286 conn.resize(8);
287 conn[0] = this->node_id(0)+1;
288 conn[1] = this->node_id(1)+1;
289 conn[2] = this->node_id(2)+1;
290 conn[3] = this->node_id(2)+1;
291 conn[4] = this->node_id(3)+1;
292 conn[5] = this->node_id(4)+1;
293 conn[6] = this->node_id(5)+1;
294 conn[7] = this->node_id(5)+1;
295 return;
296 }
297
298 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 conn.resize(15);
304 for (unsigned i=0; i<9; ++i)
305 conn[i] = this->node_id(i);
306
307 // top "ring" of mid-edge nodes
308 conn[9] = this->node_id(12);
309 conn[10] = this->node_id(13);
310 conn[11] = this->node_id(14);
311
312 // middle "ring" of mid-edge nodes
313 conn[12] = this->node_id(9);
314 conn[13] = this->node_id(10);
315 conn[14] = this->node_id(11);
316
317 return;
318 }
319
320 default:
321 libmesh_error_msg("Unsupported IO package " << iop);
322 }
323}
324
325
326
327
328unsigned short int Prism15::second_order_adjacent_vertex (const unsigned int n,
329 const unsigned int v) const
330{
331 libmesh_assert_greater_equal (n, this->n_vertices());
332 libmesh_assert_less (n, this->n_nodes());
333 libmesh_assert_less (v, 2);
334 return _second_order_adjacent_vertices[n-this->n_vertices()][v];
335}
336
337
338
339std::pair<unsigned short int, unsigned short int>
340Prism15::second_order_child_vertex (const unsigned int n) const
341{
342 libmesh_assert_greater_equal (n, this->n_vertices());
343 libmesh_assert_less (n, this->n_nodes());
344
345 return std::pair<unsigned short int, unsigned short int>
346 (_second_order_vertex_child_number[n],
347 _second_order_vertex_child_index[n]);
348}
349
350
351
352Real Prism15::volume () const
353{
354 // This specialization is good for Lagrange mappings only in general
355 if (this->mapping_type() != LAGRANGE_MAP)
356 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 x0 = point(0), x1 = point(1), x2 = point(2), x3 = point(3), x4 = point(4),
362 x5 = point(5), x6 = point(6), x7 = point(7), x8 = point(8), x9 = point(9),
363 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 -x0 - x1 + x10 + 2*x12 - x3 - x4 + 2*x6 - x9,
369 3*x0/2 + x1/2 + 2*x12 - 3*x3/2 - x4/2 - 2*x6,
370 -x0/2 + x1/2 - x10 - x3/2 + x4/2 + x9,
371 2*x0 - 2*x12 + 2*x13 - 2*x14 + 2*x3 - 2*x6 + 2*x7 - 2*x8,
372 -2*x0 - 2*x12 + 2*x13 - 2*x14 + 2*x3 + 2*x6 - 2*x7 + 2*x8,
373 Point(0,0,0),
374 2*x0 + 2*x1 - 4*x12 + 2*x3 + 2*x4 - 4*x6,
375 -2*x0 - 2*x1 - 4*x12 + 2*x3 + 2*x4 + 4*x6,
376 Point(0,0,0),
377 Point(0,0,0)
378 };
379
380 Point dx_deta[10] =
381 {
382 -x0 + x11 + 2*x14 - x2 - x3 - x5 + 2*x8 - x9,
383 3*x0/2 + 2*x14 + x2/2 - 3*x3/2 - x5/2 - 2*x8,
384 -x0/2 - x11 + x2/2 - x3/2 + x5/2 + x9,
385 2*x0 - 4*x14 + 2*x2 + 2*x3 + 2*x5 - 4*x8,
386 -2*x0 - 4*x14 - 2*x2 + 2*x3 + 2*x5 + 4*x8,
387 Point(0,0,0),
388 2*x0 - 2*x12 + 2*x13 - 2*x14 + 2*x3 - 2*x6 + 2*x7 - 2*x8,
389 -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 };
393
394 Point dx_dzeta[10] =
395 {
396 -x0/2 + x3/2,
397 x0 + x3 - 2*x9,
398 Point(0,0,0),
399 3*x0/2 + 2*x14 + x2/2 - 3*x3/2 - x5/2 - 2*x8,
400 -x0 - 2*x11 + x2 - x3 + x5 + 2*x9,
401 -x0 - 2*x14 - x2 + x3 + x5 + 2*x8,
402 3*x0/2 + x1/2 + 2*x12 - 3*x3/2 - x4/2 - 2*x6,
403 -x0 + x1 - 2*x10 - x3 + x4 + 2*x9,
404 -2*x0 - 2*x12 + 2*x13 - 2*x14 + 2*x3 + 2*x6 - 2*x7 + 2*x8,
405 -x0 - x1 - 2*x12 + x3 + x4 + 2*x6
406 };
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 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 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 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 Real vol = 0.;
481 for (int i=0; i<N2D; ++i)
482 for (int j=0; j<N1D; ++j)
483 {
484 // Compute dx_dxi, dx_deta, dx_dzeta at the current quadrature point.
485 Point dx_dxi_q, dx_deta_q, dx_dzeta_q;
486 for (int c=0; c<10; ++c)
487 {
488 Real coeff =
489 xi[i][exponents[c][0]]*
490 eta[i][exponents[c][1]]*
491 zeta[j][exponents[c][2]];
492
493 dx_dxi_q += coeff * dx_dxi[c];
494 dx_deta_q += coeff * dx_deta[c];
495 dx_dzeta_q += coeff * dx_dzeta[c];
496 }
497
498 // Compute scalar triple product, multiply by weight, and accumulate volume.
499 vol += w2D[i] * w1D[j] * triple_product(dx_dxi_q, dx_deta_q, dx_dzeta_q);
500 }
501
502 return vol;
503}
504
505
506
507#ifdef LIBMESH_ENABLE_AMR
508
509const 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
675void
676Prism15::permute(unsigned int perm_num)
677{
678 libmesh_assert_less (perm_num, 6);
679 const unsigned int side = perm_num % 2;
680 const unsigned int rotate = perm_num / 2;
681
682 for (unsigned int i = 0; i != rotate; ++i)
683 {
684 swap3nodes(0,1,2);
685 swap3nodes(3,4,5);
686 swap3nodes(6,7,8);
687 swap3nodes(9,10,11);
688 swap3nodes(12,13,14);
689 swap3neighbors(1,2,3);
690 }
691
692 switch (side) {
693 case 0:
694 break;
695 case 1:
696 swap2nodes(1,3);
697 swap2nodes(0,4);
698 swap2nodes(2,5);
699 swap2nodes(6,12);
700 swap2nodes(9,10);
701 swap2nodes(7,14);
702 swap2nodes(8,13);
703 swap2neighbors(0,4);
704 swap2neighbors(2,3);
705 break;
706 default:
707 libmesh_error();
708 }
709
710}
711
712
713void
714Prism15::flip(BoundaryInfo * boundary_info)
715{
716 libmesh_assert(boundary_info);
717
718 swap2nodes(0,1);
719 swap2nodes(3,4);
720 swap2nodes(7,8);
721 swap2nodes(9,10);
722 swap2nodes(13,14);
723 swap2neighbors(2,3);
724 swap2boundarysides(2,3,boundary_info);
725 swap2boundaryedges(0,1,boundary_info);
726 swap2boundaryedges(3,4,boundary_info);
727 swap2boundaryedges(7,8,boundary_info);
728}
729
730
731ElemType
732Prism15::side_type (const unsigned int s) const
733{
734 libmesh_assert_less (s, 5);
735 if (s == 0 || s == 4)
736 return TRI6;
737 return QUAD8;
738}
739
740} // namespace libMesh
libMesh::BoundaryInfo
The BoundaryInfo class contains information relevant to boundary conditions including storing faces,...
Definition boundary_info.h:58
libMesh::Elem::affine_tol
static constexpr Real affine_tol
Default tolerance to use in has_affine_map().
Definition elem.h:2070
libMesh::Elem::point
const Point & point(const unsigned int i) const
Definition elem.h:2462
libMesh::Elem::swap2nodes
void swap2nodes(unsigned int n1, unsigned int n2)
Swaps two node_ptrs.
Definition elem.h:2101
libMesh::Elem::swap3neighbors
void swap3neighbors(unsigned int n1, unsigned int n2, unsigned int n3)
Swaps three neighbor_ptrs, "rotating" them.
Definition elem.h:2142
libMesh::Elem::_nodes
Node ** _nodes
Pointers to the nodes we are connected to.
Definition elem.h:2254
libMesh::Elem::swap2neighbors
void swap2neighbors(unsigned int n1, unsigned int n2)
Swaps two neighbor_ptrs.
Definition elem.h:2111
libMesh::Elem::swap2boundarysides
void swap2boundarysides(unsigned short s1, unsigned short s2, BoundaryInfo *boundary_info) const
Swaps two sides in boundary_info, if it is non-null.
Definition elem.C:3579
libMesh::Elem::mapping_type
ElemMappingType mapping_type() const
Definition elem.h:3137
libMesh::Elem::swap3nodes
void swap3nodes(unsigned int n1, unsigned int n2, unsigned int n3)
Swaps three node_ptrs, "rotating" them.
Definition elem.h:2133
libMesh::Elem::volume
virtual Real volume() const
Definition elem.C:3462
libMesh::Elem::swap2boundaryedges
void swap2boundaryedges(unsigned short e1, unsigned short e2, BoundaryInfo *boundary_info) const
Swaps two edges in boundary_info, if it is non-null.
Definition elem.C:3595
libMesh::Elem::node_id
dof_id_type node_id(const unsigned int i) const
Definition elem.h:2484
libMesh::Point
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition point.h:40
libMesh::Prism15::is_edge
virtual bool is_edge(const unsigned int i) const override
Definition cell_prism15.C:70
libMesh::Prism15::nodes_per_side
static const int nodes_per_side
Definition cell_prism15.h:221
libMesh::Prism15::num_nodes
static const int num_nodes
Geometric constants for Prism15.
Definition cell_prism15.h:220
libMesh::Prism15::edge_nodes_map
static const unsigned int edge_nodes_map[num_edges][nodes_per_edge]
This maps the node of the edge to element node numbers.
Definition cell_prism15.h:234
libMesh::Prism15::second_order_child_vertex
virtual std::pair< unsigned short int, unsigned short int > second_order_child_vertex(const unsigned int n) const override
Definition cell_prism15.C:340
libMesh::Prism15::nodes_on_side
virtual std::vector< unsigned int > nodes_on_side(const unsigned int s) const override
Definition cell_prism15.C:92
libMesh::Prism15::build_edge_ptr
virtual std::unique_ptr< Elem > build_edge_ptr(const unsigned int i) override
Builds a EDGE3 or INFEDGE2 coincident with edge i.
Definition cell_prism15.C:260
libMesh::Prism15::is_node_on_side
virtual bool is_node_on_side(const unsigned int n, const unsigned int s) const override
Definition cell_prism15.C:82
libMesh::Prism15::connectivity
virtual void connectivity(const unsigned int sc, const IOPackage iop, std::vector< dof_id_type > &conn) const override
Definition cell_prism15.C:274
libMesh::Prism15::is_face
virtual bool is_face(const unsigned int i) const override
Definition cell_prism15.C:77
libMesh::Prism15::volume
virtual Real volume() const override
A specialization for computing the volume of a Prism15.
Definition cell_prism15.C:352
libMesh::Prism15::side_nodes_map
static const unsigned int side_nodes_map[num_sides][nodes_per_side]
This maps the node of the side to element node numbers.
Definition cell_prism15.h:228
libMesh::Prism15::flip
virtual void flip(BoundaryInfo *) override final
Flips the element (by swapping node and neighbor pointers) to have a mapping Jacobian of opposite sig...
Definition cell_prism15.C:714
libMesh::Prism15::local_side_node
virtual unsigned int local_side_node(unsigned int side, unsigned int side_node) const override
Definition cell_prism15.C:154
libMesh::Prism15::local_edge_node
virtual unsigned int local_edge_node(unsigned int edge, unsigned int edge_node) const override
Definition cell_prism15.C:170
libMesh::Prism15::default_order
virtual Order default_order() const override
Definition cell_prism15.C:147
libMesh::Prism15::nodes_per_edge
static const int nodes_per_edge
Definition cell_prism15.h:222
libMesh::Prism15::is_node_on_edge
virtual bool is_node_on_edge(const unsigned int n, const unsigned int e) const override
Definition cell_prism15.C:106
libMesh::Prism15::has_affine_map
virtual bool has_affine_map() const override
Definition cell_prism15.C:117
libMesh::Prism15::is_vertex
virtual bool is_vertex(const unsigned int i) const override
Definition cell_prism15.C:63
libMesh::Prism15::_embedding_matrix
static const Real _embedding_matrix[num_children][num_nodes][num_nodes]
Matrix that computes new nodal locations/solution values from current nodes/solution.
Definition cell_prism15.h:270
libMesh::Prism15::side_type
ElemType side_type(const unsigned int s) const override final
Definition cell_prism15.C:732
libMesh::Prism15::permute
virtual void permute(unsigned int perm_num) override final
Permutes the element (by swapping node and neighbor pointers) according to the specified index.
Definition cell_prism15.C:676
libMesh::Prism15::n_sub_elem
virtual unsigned int n_sub_elem() const override
Definition cell_prism15.h:106
libMesh::Prism15::second_order_adjacent_vertex
virtual unsigned short int second_order_adjacent_vertex(const unsigned int n, const unsigned int v) const override
Definition cell_prism15.C:328
libMesh::Prism15::n_nodes
virtual unsigned int n_nodes() const override
Definition cell_prism15.h:101
libMesh::Prism15::nodes_on_edge
virtual std::vector< unsigned int > nodes_on_edge(const unsigned int e) const override
Definition cell_prism15.C:100
libMesh::Prism15::build_side_ptr
virtual std::unique_ptr< Elem > build_side_ptr(const unsigned int i) override
Builds a QUAD8 or TRI6 built coincident with face i.
Definition cell_prism15.C:181
libMesh::Prism::n_vertices
virtual unsigned int n_vertices() const override final
Definition cell_prism.h:91
libMesh::Prism::n_edges
virtual unsigned int n_edges() const override final
Definition cell_prism.h:96
libMesh::Prism::num_children
static const int num_children
Definition cell_prism.h:75
libMesh::Prism::n_sides
virtual unsigned int n_sides() const override final
Definition cell_prism.h:86
libMesh::Prism::_second_order_vertex_child_number
static const unsigned short int _second_order_vertex_child_number[18]
Vector that names a child sharing each second order node.
Definition cell_prism.h:204
libMesh::Prism::num_edges
static const int num_edges
Definition cell_prism.h:74
libMesh::Prism::num_sides
static const int num_sides
Geometric constants for all Prisms.
Definition cell_prism.h:73
libMesh::Prism::_second_order_adjacent_vertices
static const unsigned short int _second_order_adjacent_vertices[9][2]
Matrix that tells which vertices define the location of mid-side (or second-order) nodes.
Definition cell_prism.h:199
libMesh::Prism::_second_order_vertex_child_index
static const unsigned short int _second_order_vertex_child_index[18]
Vector that names the child vertex index for each second order node.
Definition cell_prism.h:209
libMesh::TypeVector::relative_fuzzy_equals
bool relative_fuzzy_equals(const TypeVector< T > &rhs, Real tol=TOLERANCE) const
Definition type_vector.h:984
libMesh
The libMesh namespace provides an interface to certain functionality in the library.
libMesh::IOPackage
IOPackage
libMesh interfaces with several different software packages for the purposes of creating,...
Definition enum_io_package.h:38
libMesh::ElemType
ElemType
Defines an enum for geometric element types.
Definition enum_elem_type.h:33
libMesh::triple_product
T triple_product(const TypeVector< T > &a, const TypeVector< T > &b, const TypeVector< T > &c)
Definition type_vector.h:1032
libMesh::libmesh_assert
libmesh_assert(ctx)
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition libmesh_common.h:144
Generated on Fri Aug 21 2026 13:32:29 for libMesh by doxygen 1.9.8