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
38// ------------------------------------------------------------
39// Prism15 class member functions
40
41bool Prism15::is_vertex(const unsigned int i) const
42{
43 if (i < 6)
44 return true;
45 return false;
46}
47
48bool Prism15::is_edge(const unsigned int i) const
49{
50 if (i < 6)
51 return false;
52 return true;
53}
54
55bool Prism15::is_face(const unsigned int) const
56{
57 return false;
58}
59
60bool Prism15::is_node_on_side(const unsigned int n,
61 const unsigned int s) const
62{
63 libmesh_assert_less (s, n_sides());
64 return std::find(std::begin(side_nodes_map[s]),
65 std::end(side_nodes_map[s]),
66 n) != std::end(side_nodes_map[s]);
67}
68
69std::vector<unsigned int>
70Prism15::nodes_on_side(const unsigned int s) const
71{
72 libmesh_assert_less(s, n_sides());
73 auto trim = (s > 0 && s < 4) ? 0 : 2;
74 return {std::begin(side_nodes_map[s]), std::end(side_nodes_map[s]) - trim};
75}
76
77std::vector<unsigned>
78Prism15::nodes_on_edge(const unsigned int e) const
79{
80 libmesh_assert_less(e, n_edges());
81 return {std::begin(edge_nodes_map[e]), std::end(edge_nodes_map[e])};
82}
83
84bool Prism15::is_node_on_edge(const unsigned int n,
85 const unsigned int e) const
86{
87 libmesh_assert_less (e, n_edges());
88 return std::find(std::begin(edge_nodes_map[e]),
89 std::end(edge_nodes_map[e]),
90 n) != std::end(edge_nodes_map[e]);
91}
92
93
94
95bool Prism15::has_affine_map() const
96{
97 // Make sure z edges are affine
98 Point v = this->point(3) - this->point(0);
99 if (!v.relative_fuzzy_equals(this->point(4) - this->point(1), affine_tol) ||
100 !v.relative_fuzzy_equals(this->point(5) - this->point(2), affine_tol))
101 return false;
102 // Make sure edges are straight
103 v /= 2;
104 if (!v.relative_fuzzy_equals(this->point(9) - this->point(0), affine_tol) ||
105 !v.relative_fuzzy_equals(this->point(10) - this->point(1), affine_tol) ||
106 !v.relative_fuzzy_equals(this->point(11) - this->point(2), affine_tol))
107 return false;
108 v = (this->point(1) - this->point(0))/2;
109 if (!v.relative_fuzzy_equals(this->point(6) - this->point(0), affine_tol) ||
110 !v.relative_fuzzy_equals(this->point(12) - this->point(3), affine_tol))
111 return false;
112 v = (this->point(2) - this->point(0))/2;
113 if (!v.relative_fuzzy_equals(this->point(8) - this->point(0), affine_tol) ||
114 !v.relative_fuzzy_equals(this->point(14) - this->point(3), affine_tol))
115 return false;
116 v = (this->point(2) - this->point(1))/2;
117 if (!v.relative_fuzzy_equals(this->point(7) - this->point(1), affine_tol) ||
118 !v.relative_fuzzy_equals(this->point(13) - this->point(4), affine_tol))
119 return false;
120 return true;
121}
122
123
124
125Order Prism15::default_order() const
126{
127 return SECOND;
128}
129
130
131
132unsigned int Prism15::local_side_node(unsigned int side,
133 unsigned int side_node) const
134{
135 libmesh_assert_less (side, this->n_sides());
136
137 // Never more than 8 nodes per side.
138 libmesh_assert_less(side_node, Prism15::nodes_per_side);
139
140 // Some sides have 6 nodes.
141 libmesh_assert(!(side==0 || side==4) || side_node < 6);
142
143 return Prism15::side_nodes_map[side][side_node];
144}
145
146
147
148unsigned int Prism15::local_edge_node(unsigned int edge,
149 unsigned int edge_node) const
150{
151 libmesh_assert_less(edge, this->n_edges());
152 libmesh_assert_less(edge_node, Prism15::nodes_per_edge);
153
154 return Prism15::edge_nodes_map[edge][edge_node];
155}
156
157
158
159std::unique_ptr<Elem> Prism15::build_side_ptr (const unsigned int i)
160{
161 libmesh_assert_less (i, this->n_sides());
162
163 std::unique_ptr<Elem> face;
164
165 switch (i)
166 {
167 case 0: // the triangular face at z=-1
168 case 4: // the triangular face at z=1
169 {
170 face = std::make_unique<Tri6>();
171 break;
172 }
173 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 face = std::make_unique<Quad8>();
178 break;
179 }
180 default:
181 libmesh_error_msg("Invalid side i = " << i);
182 }
183
184 // Set the nodes
185 for (auto n : face->node_index_range())
186 face->set_node(n, this->node_ptr(Prism15::side_nodes_map[i][n]));
187
188 face->set_interior_parent(this);
189 face->inherit_data_from(*this);
190
191 return face;
192}
193
194
195void Prism15::build_side_ptr (std::unique_ptr<Elem> & side,
196 const unsigned int i)
197{
198 libmesh_assert_less (i, this->n_sides());
199
200 switch (i)
201 {
202 case 0: // the triangular face at z=-1
203 case 4: // the triangular face at z=1
204 {
205 if (!side.get() || side->type() != TRI6)
206 {
207 side = this->build_side_ptr(i);
208 return;
209 }
210 break;
211 }
212
213 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 if (!side.get() || side->type() != QUAD8)
218 {
219 side = this->build_side_ptr(i);
220 return;
221 }
222 break;
223 }
224
225 default:
226 libmesh_error_msg("Invalid side i = " << i);
227 }
228
229 side->inherit_data_from(*this);
230
231 // Set the nodes
232 for (auto n : side->node_index_range())
233 side->set_node(n, this->node_ptr(Prism15::side_nodes_map[i][n]));
234}
235
236
237
238std::unique_ptr<Elem> Prism15::build_edge_ptr (const unsigned int i)
239{
240 return this->simple_build_edge_ptr<Edge3,Prism15>(i);
241}
242
243
244
245void Prism15::build_edge_ptr (std::unique_ptr<Elem> & edge, const unsigned int i)
246{
247 this->simple_build_edge_ptr<Prism15>(edge, i, EDGE3);
248}
249
250
251
252void Prism15::connectivity(const unsigned int libmesh_dbg_var(sc),
253 const IOPackage iop,
254 std::vector<dof_id_type> & conn) const
255{
256 libmesh_assert(_nodes);
257 libmesh_assert_less (sc, this->n_sub_elem());
258 libmesh_assert_not_equal_to (iop, INVALID_IO_PACKAGE);
259
260 switch (iop)
261 {
262 case TECPLOT:
263 {
264 conn.resize(8);
265 conn[0] = this->node_id(0)+1;
266 conn[1] = this->node_id(1)+1;
267 conn[2] = this->node_id(2)+1;
268 conn[3] = this->node_id(2)+1;
269 conn[4] = this->node_id(3)+1;
270 conn[5] = this->node_id(4)+1;
271 conn[6] = this->node_id(5)+1;
272 conn[7] = this->node_id(5)+1;
273 return;
274 }
275
276 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 conn.resize(15);
282 for (unsigned i=0; i<9; ++i)
283 conn[i] = this->node_id(i);
284
285 // top "ring" of mid-edge nodes
286 conn[9] = this->node_id(12);
287 conn[10] = this->node_id(13);
288 conn[11] = this->node_id(14);
289
290 // middle "ring" of mid-edge nodes
291 conn[12] = this->node_id(9);
292 conn[13] = this->node_id(10);
293 conn[14] = this->node_id(11);
294
295 return;
296 }
297
298 default:
299 libmesh_error_msg("Unsupported IO package " << iop);
300 }
301}
302
303
304
305
306unsigned short int Prism15::second_order_adjacent_vertex (const unsigned int n,
307 const unsigned int v) const
308{
309 libmesh_assert_greater_equal (n, this->n_vertices());
310 libmesh_assert_less (n, this->n_nodes());
311 libmesh_assert_less (v, 2);
312 return _second_order_adjacent_vertices[n-this->n_vertices()][v];
313}
314
315
316
317std::pair<unsigned short int, unsigned short int>
318Prism15::second_order_child_vertex (const unsigned int n) const
319{
320 libmesh_assert_greater_equal (n, this->n_vertices());
321 libmesh_assert_less (n, this->n_nodes());
322
323 return std::pair<unsigned short int, unsigned short int>
324 (_second_order_vertex_child_number[n],
325 _second_order_vertex_child_index[n]);
326}
327
328
329
330Real Prism15::volume () const
331{
332 // This specialization is good for Lagrange mappings only in general
333 if (this->mapping_type() != LAGRANGE_MAP)
334 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 x0 = point(0), x1 = point(1), x2 = point(2), x3 = point(3), x4 = point(4),
340 x5 = point(5), x6 = point(6), x7 = point(7), x8 = point(8), x9 = point(9),
341 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 -x0 - x1 + x10 + 2*x12 - x3 - x4 + 2*x6 - x9,
347 3*x0/2 + x1/2 + 2*x12 - 3*x3/2 - x4/2 - 2*x6,
348 -x0/2 + x1/2 - x10 - x3/2 + x4/2 + x9,
349 2*x0 - 2*x12 + 2*x13 - 2*x14 + 2*x3 - 2*x6 + 2*x7 - 2*x8,
350 -2*x0 - 2*x12 + 2*x13 - 2*x14 + 2*x3 + 2*x6 - 2*x7 + 2*x8,
351 Point(0,0,0),
352 2*x0 + 2*x1 - 4*x12 + 2*x3 + 2*x4 - 4*x6,
353 -2*x0 - 2*x1 - 4*x12 + 2*x3 + 2*x4 + 4*x6,
354 Point(0,0,0),
355 Point(0,0,0)
356 };
357
358 Point dx_deta[10] =
359 {
360 -x0 + x11 + 2*x14 - x2 - x3 - x5 + 2*x8 - x9,
361 3*x0/2 + 2*x14 + x2/2 - 3*x3/2 - x5/2 - 2*x8,
362 -x0/2 - x11 + x2/2 - x3/2 + x5/2 + x9,
363 2*x0 - 4*x14 + 2*x2 + 2*x3 + 2*x5 - 4*x8,
364 -2*x0 - 4*x14 - 2*x2 + 2*x3 + 2*x5 + 4*x8,
365 Point(0,0,0),
366 2*x0 - 2*x12 + 2*x13 - 2*x14 + 2*x3 - 2*x6 + 2*x7 - 2*x8,
367 -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 };
371
372 Point dx_dzeta[10] =
373 {
374 -x0/2 + x3/2,
375 x0 + x3 - 2*x9,
376 Point(0,0,0),
377 3*x0/2 + 2*x14 + x2/2 - 3*x3/2 - x5/2 - 2*x8,
378 -x0 - 2*x11 + x2 - x3 + x5 + 2*x9,
379 -x0 - 2*x14 - x2 + x3 + x5 + 2*x8,
380 3*x0/2 + x1/2 + 2*x12 - 3*x3/2 - x4/2 - 2*x6,
381 -x0 + x1 - 2*x10 - x3 + x4 + 2*x9,
382 -2*x0 - 2*x12 + 2*x13 - 2*x14 + 2*x3 + 2*x6 - 2*x7 + 2*x8,
383 -x0 - x1 - 2*x12 + x3 + x4 + 2*x6
384 };
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 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 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 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 Real vol = 0.;
459 for (int i=0; i<N2D; ++i)
460 for (int j=0; j<N1D; ++j)
461 {
462 // Compute dx_dxi, dx_deta, dx_dzeta at the current quadrature point.
463 Point dx_dxi_q, dx_deta_q, dx_dzeta_q;
464 for (int c=0; c<10; ++c)
465 {
466 Real coeff =
467 xi[i][exponents[c][0]]*
468 eta[i][exponents[c][1]]*
469 zeta[j][exponents[c][2]];
470
471 dx_dxi_q += coeff * dx_dxi[c];
472 dx_deta_q += coeff * dx_deta[c];
473 dx_dzeta_q += coeff * dx_dzeta[c];
474 }
475
476 // Compute scalar triple product, multiply by weight, and accumulate volume.
477 vol += w2D[i] * w1D[j] * triple_product(dx_dxi_q, dx_deta_q, dx_dzeta_q);
478 }
479
480 return vol;
481}
482
483
484
485#ifdef LIBMESH_ENABLE_AMR
486
487const 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
653void
654Prism15::permute(unsigned int perm_num)
655{
656 libmesh_assert_less (perm_num, 6);
657 const unsigned int side = perm_num % 2;
658 const unsigned int rotate = perm_num / 2;
659
660 for (unsigned int i = 0; i != rotate; ++i)
661 {
662 swap3nodes(0,1,2);
663 swap3nodes(3,4,5);
664 swap3nodes(6,7,8);
665 swap3nodes(9,10,11);
666 swap3nodes(12,13,14);
667 swap3neighbors(1,2,3);
668 }
669
670 switch (side) {
671 case 0:
672 break;
673 case 1:
674 swap2nodes(1,3);
675 swap2nodes(0,4);
676 swap2nodes(2,5);
677 swap2nodes(6,12);
678 swap2nodes(9,10);
679 swap2nodes(7,14);
680 swap2nodes(8,13);
681 swap2neighbors(0,4);
682 swap2neighbors(2,3);
683 break;
684 default:
685 libmesh_error();
686 }
687
688}
689
690
691void
692Prism15::flip(BoundaryInfo * boundary_info)
693{
694 libmesh_assert(boundary_info);
695
696 swap2nodes(0,1);
697 swap2nodes(3,4);
698 swap2nodes(7,8);
699 swap2nodes(9,10);
700 swap2nodes(13,14);
701 swap2neighbors(2,3);
702 swap2boundarysides(2,3,boundary_info);
703 swap2boundaryedges(0,1,boundary_info);
704 swap2boundaryedges(3,4,boundary_info);
705 swap2boundaryedges(7,8,boundary_info);
706}
707
708
709ElemType
710Prism15::side_type (const unsigned int s) const
711{
712 libmesh_assert_less (s, 5);
713 if (s == 0 || s == 4)
714 return TRI6;
715 return QUAD8;
716}
717
718} // 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:2191
libMesh::Elem::point
const Point & point(const unsigned int i) const
Definition elem.h:2583
libMesh::Elem::swap2nodes
void swap2nodes(unsigned int n1, unsigned int n2)
Swaps two node_ptrs.
Definition elem.h:2222
libMesh::Elem::swap3neighbors
void swap3neighbors(unsigned int n1, unsigned int n2, unsigned int n3)
Swaps three neighbor_ptrs, "rotating" them.
Definition elem.h:2263
libMesh::Elem::_nodes
Node ** _nodes
Pointers to the nodes we are connected to.
Definition elem.h:2375
libMesh::Elem::swap2neighbors
void swap2neighbors(unsigned int n1, unsigned int n2)
Swaps two neighbor_ptrs.
Definition elem.h:2232
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:3982
libMesh::Elem::mapping_type
ElemMappingType mapping_type() const
Definition elem.h:3258
libMesh::Elem::swap3nodes
void swap3nodes(unsigned int n1, unsigned int n2, unsigned int n3)
Swaps three node_ptrs, "rotating" them.
Definition elem.h:2254
libMesh::Elem::volume
virtual Real volume() const
Definition elem.C:3865
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:3998
libMesh::Elem::node_id
dof_id_type node_id(const unsigned int i) const
Definition elem.h:2605
libMesh::Point
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition point.h:41
libMesh::Prism15::is_edge
virtual bool is_edge(const unsigned int i) const override
Definition cell_prism15.C:48
libMesh::Prism15::nodes_per_side
static const int nodes_per_side
Definition cell_prism15.h:223
libMesh::Prism15::num_nodes
static const int num_nodes
Geometric constants for Prism15.
Definition cell_prism15.h:222
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:318
libMesh::Prism15::nodes_on_side
virtual std::vector< unsigned int > nodes_on_side(const unsigned int s) const override
Definition cell_prism15.C:70
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:238
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:60
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:252
libMesh::Prism15::is_face
virtual bool is_face(const unsigned int i) const override
Definition cell_prism15.C:55
libMesh::Prism15::volume
virtual Real volume() const override
A specialization for computing the volume of a Prism15.
Definition cell_prism15.C:330
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:692
libMesh::Prism15::local_side_node
virtual unsigned int local_side_node(unsigned int side, unsigned int side_node) const override
Definition cell_prism15.C:132
libMesh::Prism15::local_edge_node
virtual unsigned int local_edge_node(unsigned int edge, unsigned int edge_node) const override
Definition cell_prism15.C:148
libMesh::Prism15::default_order
virtual Order default_order() const override
Definition cell_prism15.C:125
libMesh::Prism15::nodes_per_edge
static const int nodes_per_edge
Definition cell_prism15.h:224
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:84
libMesh::Prism15::has_affine_map
virtual bool has_affine_map() const override
Definition cell_prism15.C:95
libMesh::Prism15::is_vertex
virtual bool is_vertex(const unsigned int i) const override
Definition cell_prism15.C:41
libMesh::Prism15::side_nodes_map
static constexpr const unsigned int(& side_nodes_map)[num_sides][nodes_per_side]
Definition cell_prism15.h:238
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:274
libMesh::Prism15::side_type
ElemType side_type(const unsigned int s) const override final
Definition cell_prism15.C:710
libMesh::Prism15::edge_nodes_map
static constexpr const unsigned int(& edge_nodes_map)[num_edges][nodes_per_edge]
Definition cell_prism15.h:233
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:654
libMesh::Prism15::n_sub_elem
virtual unsigned int n_sub_elem() const override
Definition cell_prism15.h:108
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:306
libMesh::Prism15::n_nodes
virtual unsigned int n_nodes() const override
Definition cell_prism15.h:103
libMesh::Prism15::nodes_on_edge
virtual std::vector< unsigned int > nodes_on_edge(const unsigned int e) const override
Definition cell_prism15.C:78
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:159
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:213
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:208
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:218
libMesh::TypeVector::relative_fuzzy_equals
bool relative_fuzzy_equals(const TypeVector< T > &rhs, Real tol=TOLERANCE) const
Definition type_vector.h:1052
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::triple_product
LIBMESH_DEVICE_INLINE T triple_product(const TypeVector< T > &a, const TypeVector< T > &b, const TypeVector< T > &c)
Definition type_vector.h:1100
libMesh::ElemType
ElemType
Defines an enum for geometric element types.
Definition enum_elem_type.h:33
libMesh::libmesh_assert
libmesh_assert(ctx)
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition libmesh_common.h:148
Generated on Thu Nov 12 2026 20:15:22 for libMesh by doxygen 1.9.8