Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 :
19 : // Local includes
20 : #include "libmesh/cell_prism6.h"
21 : #include "libmesh/edge_edge2.h"
22 : #include "libmesh/face_quad4.h"
23 : #include "libmesh/face_tri3.h"
24 : #include "libmesh/enum_io_package.h"
25 : #include "libmesh/enum_order.h"
26 : #include "libmesh/fe_lagrange_shape_1D.h"
27 :
28 : // C++ includes
29 : #include <array>
30 :
31 : // Helper functions for computing volume() and true_centroid()
32 : namespace {
33 :
34 : using namespace libMesh;
35 :
36 : // Templated numerical quadrature implementation function used by both
37 : // Prism6::volume() and Prism6::true_centroid() routines.
38 : template <typename T>
39 : void cell_integral(const Elem * elem, T & t);
40 :
41 : /**
42 : * Object passed to cell_integral() routine for computing the cell volume.
43 : */
44 : struct VolumeIntegral
45 : {
46 : // API
47 256 : void accumulate_qp(Real /*xi*/, Real /*eta*/, Real /*zeta*/, Real JxW)
48 : {
49 9088 : vol += JxW;
50 256 : };
51 :
52 : Real vol = 0.;
53 : };
54 :
55 : /**
56 : * Object passed to cell_integral() routine for computing the nodal
57 : * volumes (one value per node)
58 : */
59 : struct NodalVolumeIntegral
60 : {
61 : // API
62 19840 : void accumulate_qp(Real xi, Real eta, Real zeta, Real JxW)
63 : {
64 : // Copied from fe_lagrange_shape_3D.C
65 : static const unsigned int i0[] = {0, 0, 0, 1, 1, 1};
66 : static const unsigned int i1[] = {0, 1, 2, 0, 1, 2};
67 :
68 : // Copied from fe_lagrange_shape_2D.C
69 119040 : auto tri3_phi = [](int i, Real x, Real y)
70 : {
71 119040 : switch(i)
72 : {
73 39680 : case 0:
74 39680 : return 1 - x - y;
75 :
76 2304 : case 1:
77 2304 : return x;
78 :
79 39680 : case 2:
80 37376 : return y;
81 :
82 0 : default:
83 0 : libmesh_error_msg("Invalid shape function index i = " << i);
84 : }
85 : };
86 :
87 138880 : for (int i=0; i<6; ++i)
88 : {
89 : Real phi =
90 119040 : tri3_phi(i1[i], xi, eta) * fe_lagrange_1D_linear_shape(i0[i], zeta);
91 :
92 119040 : V[i] += phi * JxW;
93 : }
94 19840 : };
95 :
96 : // nodal volumes
97 : std::array<Real, 6> V {};
98 : };
99 :
100 : }
101 :
102 : namespace libMesh
103 : {
104 :
105 :
106 :
107 : // ------------------------------------------------------------
108 : // Prism6 class static member initializations
109 : const int Prism6::num_nodes;
110 : const int Prism6::nodes_per_side;
111 : const int Prism6::nodes_per_edge;
112 :
113 : const unsigned int Prism6::side_nodes_map[Prism6::num_sides][Prism6::nodes_per_side] =
114 : {
115 : {0, 2, 1, 99}, // Side 0
116 : {0, 1, 4, 3}, // Side 1
117 : {1, 2, 5, 4}, // Side 2
118 : {2, 0, 3, 5}, // Side 3
119 : {3, 4, 5, 99} // Side 4
120 : };
121 :
122 : const unsigned int Prism6::side_elems_map[Prism6::num_sides][Prism6::nodes_per_side] =
123 : {
124 : {0, 1, 2, 3}, // Side 0
125 : {0, 1, 4, 5}, // Side 1
126 : {1, 2, 5, 6}, // Side 2
127 : {0, 2, 4, 6}, // Side 3
128 : {4, 5, 6, 7} // Side 4
129 : };
130 :
131 : const unsigned int Prism6::edge_nodes_map[Prism6::num_edges][Prism6::nodes_per_edge] =
132 : {
133 : {0, 1}, // Edge 0
134 : {1, 2}, // Edge 1
135 : {0, 2}, // Edge 2
136 : {0, 3}, // Edge 3
137 : {1, 4}, // Edge 4
138 : {2, 5}, // Edge 5
139 : {3, 4}, // Edge 6
140 : {4, 5}, // Edge 7
141 : {3, 5} // Edge 8
142 : };
143 :
144 : // ------------------------------------------------------------
145 : // Prism6 class member functions
146 :
147 1071502 : bool Prism6::is_vertex(const unsigned int libmesh_dbg_var(n)) const
148 : {
149 61836 : libmesh_assert_not_equal_to (n, invalid_uint);
150 1071502 : return true;
151 : }
152 :
153 0 : bool Prism6::is_edge(const unsigned int) const
154 : {
155 0 : return false;
156 : }
157 :
158 :
159 0 : bool Prism6::is_face(const unsigned int) const
160 : {
161 0 : return false;
162 : }
163 :
164 20706 : bool Prism6::is_node_on_side(const unsigned int n,
165 : const unsigned int s) const
166 : {
167 1608 : libmesh_assert_less (s, n_sides());
168 1608 : return std::find(std::begin(side_nodes_map[s]),
169 20706 : std::end(side_nodes_map[s]),
170 20706 : n) != std::end(side_nodes_map[s]);
171 : }
172 :
173 : std::vector<unsigned>
174 40170851 : Prism6::nodes_on_side(const unsigned int s) const
175 : {
176 3583908 : libmesh_assert_less(s, n_sides());
177 40170851 : auto trim = (s > 0 && s < 4) ? 0 : 1;
178 40170851 : return {std::begin(side_nodes_map[s]), std::end(side_nodes_map[s]) - trim};
179 : }
180 :
181 : std::vector<unsigned>
182 9279 : Prism6::nodes_on_edge(const unsigned int e) const
183 : {
184 738 : libmesh_assert_less(e, n_edges());
185 10017 : return {std::begin(edge_nodes_map[e]), std::end(edge_nodes_map[e])};
186 : }
187 :
188 7290 : bool Prism6::is_node_on_edge(const unsigned int n,
189 : const unsigned int e) const
190 : {
191 396 : libmesh_assert_less (e, n_edges());
192 396 : return std::find(std::begin(edge_nodes_map[e]),
193 7290 : std::end(edge_nodes_map[e]),
194 7290 : n) != std::end(edge_nodes_map[e]);
195 : }
196 :
197 :
198 :
199 475377 : bool Prism6::has_affine_map() const
200 : {
201 : // Make sure z edges are affine
202 95384 : Point v = this->point(3) - this->point(0);
203 570681 : if (!v.relative_fuzzy_equals(this->point(4) - this->point(1), affine_tol) ||
204 568893 : !v.relative_fuzzy_equals(this->point(5) - this->point(2), affine_tol))
205 1848 : return false;
206 47632 : return true;
207 : }
208 :
209 :
210 :
211 11892466 : Order Prism6::default_order() const
212 : {
213 11892466 : return FIRST;
214 : }
215 :
216 :
217 :
218 175186 : std::unique_ptr<Elem> Prism6::build_side_ptr (const unsigned int i)
219 : {
220 13965 : libmesh_assert_less (i, this->n_sides());
221 :
222 175186 : std::unique_ptr<Elem> face;
223 :
224 175186 : switch (i)
225 : {
226 68394 : case 0: // the triangular face at z=-1
227 : case 4: // the triangular face at z=1
228 : {
229 68394 : face = std::make_unique<Tri3>();
230 68394 : break;
231 : }
232 106792 : case 1: // the quad face at y=0
233 : case 2: // the other quad face
234 : case 3: // the quad face at x=0
235 : {
236 106792 : face = std::make_unique<Quad4>();
237 106792 : break;
238 : }
239 0 : default:
240 0 : libmesh_error_msg("Invalid side i = " << i);
241 : }
242 :
243 : // Set the nodes
244 807536 : for (auto n : face->node_index_range())
245 682762 : face->set_node(n, this->node_ptr(Prism6::side_nodes_map[i][n]));
246 :
247 175186 : face->set_interior_parent(this);
248 161221 : face->inherit_data_from(*this);
249 :
250 175186 : return face;
251 0 : }
252 :
253 :
254 :
255 7835 : void Prism6::build_side_ptr (std::unique_ptr<Elem> & side,
256 : const unsigned int i)
257 : {
258 7835 : this->side_ptr(side, i);
259 7835 : side->set_interior_parent(this);
260 6993 : side->inherit_data_from(*this);
261 7835 : }
262 :
263 :
264 :
265 1612815 : std::unique_ptr<Elem> Prism6::build_edge_ptr (const unsigned int i)
266 : {
267 1612815 : return this->simple_build_edge_ptr<Edge2,Prism6>(i);
268 : }
269 :
270 :
271 :
272 0 : void Prism6::build_edge_ptr (std::unique_ptr<Elem> & edge, const unsigned int i)
273 : {
274 0 : this->simple_build_edge_ptr<Prism6>(edge, i, EDGE2);
275 0 : }
276 :
277 :
278 :
279 23040 : void Prism6::connectivity(const unsigned int libmesh_dbg_var(sc),
280 : const IOPackage iop,
281 : std::vector<dof_id_type> & conn) const
282 : {
283 1920 : libmesh_assert(_nodes);
284 1920 : libmesh_assert_less (sc, this->n_sub_elem());
285 1920 : libmesh_assert_not_equal_to (iop, INVALID_IO_PACKAGE);
286 :
287 23040 : switch (iop)
288 : {
289 23040 : case TECPLOT:
290 : {
291 23040 : conn.resize(8);
292 24960 : conn[0] = this->node_id(0)+1;
293 23040 : conn[1] = this->node_id(1)+1;
294 23040 : conn[2] = this->node_id(2)+1;
295 23040 : conn[3] = this->node_id(2)+1;
296 23040 : conn[4] = this->node_id(3)+1;
297 23040 : conn[5] = this->node_id(4)+1;
298 23040 : conn[6] = this->node_id(5)+1;
299 23040 : conn[7] = this->node_id(5)+1;
300 23040 : return;
301 : }
302 :
303 0 : case VTK:
304 : {
305 0 : conn.resize(6);
306 0 : conn[0] = this->node_id(0);
307 0 : conn[1] = this->node_id(2);
308 0 : conn[2] = this->node_id(1);
309 0 : conn[3] = this->node_id(3);
310 0 : conn[4] = this->node_id(5);
311 0 : conn[5] = this->node_id(4);
312 0 : return;
313 : }
314 :
315 0 : default:
316 0 : libmesh_error_msg("Unsupported IO package " << iop);
317 : }
318 : }
319 :
320 :
321 :
322 : #ifdef LIBMESH_ENABLE_AMR
323 :
324 : const Real Prism6::_embedding_matrix[Prism6::num_children][Prism6::num_nodes][Prism6::num_nodes] =
325 : {
326 : // embedding matrix for child 0
327 : {
328 : // 0 1 2 3 4 5
329 : { 1.0, 0.0, 0.0, 0.0, 0.0, 0.0}, // 0
330 : { 0.5, 0.5, 0.0, 0.0, 0.0, 0.0}, // 1
331 : { 0.5, 0.0, 0.5, 0.0, 0.0, 0.0}, // 2
332 : { 0.5, 0.0, 0.0, 0.5, 0.0, 0.0}, // 3
333 : { .25, .25, 0.0, .25, .25, 0.0}, // 4
334 : { .25, 0.0, .25, .25, 0.0, .25} // 5
335 : },
336 :
337 : // embedding matrix for child 1
338 : {
339 : // 0 1 2 3 4 5
340 : { 0.5, 0.5, 0.0, 0.0, 0.0, 0.0}, // 0
341 : { 0.0, 1.0, 0.0, 0.0, 0.0, 0.0}, // 1
342 : { 0.0, 0.5, 0.5, 0.0, 0.0, 0.0}, // 2
343 : { .25, .25, 0.0, .25, .25, 0.0}, // 3
344 : { 0.0, 0.5, 0.0, 0.0, 0.5, 0.0}, // 4
345 : { 0.0, .25, .25, 0.0, .25, .25} // 5
346 : },
347 :
348 : // embedding matrix for child 2
349 : {
350 : // 0 1 2 3 4 5
351 : { 0.5, 0.0, 0.5, 0.0, 0.0, 0.0}, // 0
352 : { 0.0, 0.5, 0.5, 0.0, 0.0, 0.0}, // 1
353 : { 0.0, 0.0, 1.0, 0.0, 0.0, 0.0}, // 2
354 : { .25, 0.0, .25, .25, 0.0, .25}, // 3
355 : { 0.0, .25, .25, 0.0, .25, .25}, // 4
356 : { 0.0, 0.0, 0.5, 0.0, 0.0, 0.5} // 5
357 : },
358 :
359 : // embedding matrix for child 3
360 : {
361 : // 0 1 2 3 4 5
362 : { 0.5, 0.5, 0.0, 0.0, 0.0, 0.0}, // 0
363 : { 0.0, 0.5, 0.5, 0.0, 0.0, 0.0}, // 1
364 : { 0.5, 0.0, 0.5, 0.0, 0.0, 0.0}, // 2
365 : { .25, .25, 0.0, .25, .25, 0.0}, // 3
366 : { 0.0, .25, .25, 0.0, .25, .25}, // 4
367 : { .25, 0.0, .25, .25, 0.0, .25} // 5
368 : },
369 :
370 : // embedding matrix for child 4
371 : {
372 : // 0 1 2 3 4 5
373 : { 0.5, 0.0, 0.0, 0.5, 0.0, 0.0}, // 0
374 : { .25, .25, 0.0, .25, .25, 0.0}, // 1
375 : { .25, 0.0, .25, .25, 0.0, .25}, // 2
376 : { 0.0, 0.0, 0.0, 1.0, 0.0, 0.0}, // 3
377 : { 0.0, 0.0, 0.0, 0.5, 0.5, 0.0}, // 4
378 : { 0.0, 0.0, 0.0, 0.5, 0.0, 0.5} // 5
379 : },
380 :
381 : // embedding matrix for child 5
382 : {
383 : // 0 1 2 3 4 5
384 : { .25, .25, 0.0, .25, .25, 0.0}, // 0
385 : { 0.0, 0.5, 0.0, 0.0, 0.5, 0.0}, // 1
386 : { 0.0, .25, .25, 0.0, .25, .25}, // 2
387 : { 0.0, 0.0, 0.0, 0.5, 0.5, 0.0}, // 3
388 : { 0.0, 0.0, 0.0, 0.0, 1.0, 0.0}, // 4
389 : { 0.0, 0.0, 0.0, 0.0, 0.5, 0.5} // 5
390 : },
391 :
392 : // embedding matrix for child 6
393 : {
394 : // 0 1 2 3 4 5
395 : { .25, 0.0, .25, .25, 0.0, .25}, // 0
396 : { 0.0, .25, .25, 0.0, .25, .25}, // 1
397 : { 0.0, 0.0, 0.5, 0.0, 0.0, 0.5}, // 2
398 : { 0.0, 0.0, 0.0, 0.5, 0.0, 0.5}, // 3
399 : { 0.0, 0.0, 0.0, 0.0, 0.5, 0.5}, // 4
400 : { 0.0, 0.0, 0.0, 0.0, 0.0, 1.0} // 5
401 : },
402 :
403 : // embedding matrix for child 7
404 : {
405 : // 0 1 2 3 4 5
406 : { .25, .25, 0.0, .25, .25, 0.0}, // 0
407 : { 0.0, .25, .25, 0.0, .25, .25}, // 1
408 : { .25, 0.0, .25, .25, 0.0, .25}, // 2
409 : { 0.0, 0.0, 0.0, 0.5, 0.5, 0.0}, // 3
410 : { 0.0, 0.0, 0.0, 0.0, 0.5, 0.5}, // 4
411 : { 0.0, 0.0, 0.0, 0.5, 0.0, 0.5} // 5
412 : }
413 : };
414 :
415 : #endif
416 :
417 :
418 :
419 2480 : Point Prism6::true_centroid () const
420 : {
421 2480 : NodalVolumeIntegral nvi;
422 2480 : cell_integral(this, nvi);
423 :
424 : // Compute centroid
425 144 : Point cp;
426 144 : Real vol = 0.;
427 :
428 17360 : for (int i=0; i<6; ++i)
429 : {
430 15744 : cp += this->point(i) * nvi.V[i];
431 14880 : vol += nvi.V[i];
432 : }
433 :
434 2480 : return cp / vol;
435 : }
436 :
437 1136 : Real Prism6::volume () const
438 : {
439 1136 : VolumeIntegral vi;
440 1136 : cell_integral(this, vi);
441 1136 : return vi.vol;
442 : }
443 :
444 : BoundingBox
445 193721 : Prism6::loose_bounding_box () const
446 : {
447 193721 : return Elem::loose_bounding_box();
448 : }
449 :
450 :
451 : void
452 4716 : Prism6::permute(unsigned int perm_num)
453 : {
454 420 : libmesh_assert_less (perm_num, 6);
455 4716 : const unsigned int side = perm_num % 2;
456 4716 : const unsigned int rotate = perm_num / 2;
457 :
458 11214 : for (unsigned int i = 0; i != rotate; ++i)
459 : {
460 6498 : swap3nodes(0,1,2);
461 5916 : swap3nodes(3,4,5);
462 5916 : swap3neighbors(1,2,3);
463 : }
464 :
465 4716 : switch (side) {
466 48 : case 0:
467 48 : break;
468 4140 : case 1:
469 4140 : swap2nodes(1,3);
470 4140 : swap2nodes(0,4);
471 4140 : swap2nodes(2,5);
472 372 : swap2neighbors(0,4);
473 372 : swap2neighbors(2,3);
474 372 : break;
475 0 : default:
476 0 : libmesh_error();
477 : }
478 :
479 4716 : }
480 :
481 :
482 : void
483 576 : Prism6::flip(BoundaryInfo * boundary_info)
484 : {
485 48 : libmesh_assert(boundary_info);
486 :
487 576 : swap2nodes(0,1);
488 576 : swap2nodes(3,4);
489 48 : swap2neighbors(2,3);
490 576 : swap2boundarysides(2,3,boundary_info);
491 576 : swap2boundaryedges(0,1,boundary_info);
492 576 : swap2boundaryedges(3,4,boundary_info);
493 576 : swap2boundaryedges(7,8,boundary_info);
494 576 : }
495 :
496 :
497 : ElemType
498 2880 : Prism6::side_type (const unsigned int s) const
499 : {
500 240 : libmesh_assert_less (s, 5);
501 2880 : if (s == 0 || s == 4)
502 1152 : return TRI3;
503 144 : return QUAD4;
504 : }
505 :
506 : } // namespace libMesh
507 :
508 : namespace // anonymous
509 : {
510 :
511 : template <typename T>
512 3616 : void cell_integral(const Elem * elem, T & t)
513 : {
514 : // Conveient references to our points.
515 : const Point
516 352 : &x0 = elem->point(0), &x1 = elem->point(1), &x2 = elem->point(2),
517 176 : &x3 = elem->point(3), &x4 = elem->point(4), &x5 = elem->point(5);
518 :
519 : // constant and zeta terms only. These are copied directly from a
520 : // Python script.
521 352 : Point dx_dxi[2] =
522 : {
523 176 : -x0/2 + x1/2 - x3/2 + x4/2, // constant
524 176 : x0/2 - x1/2 - x3/2 + x4/2, // zeta
525 : };
526 :
527 : // constant and zeta terms only. These are copied directly from a
528 : // Python script.
529 352 : Point dx_deta[2] =
530 : {
531 176 : -x0/2 + x2/2 - x3/2 + x5/2, // constant
532 176 : x0/2 - x2/2 - x3/2 + x5/2, // zeta
533 : };
534 :
535 : // Constant, xi, and eta terms
536 528 : Point dx_dzeta[3] =
537 : {
538 176 : -x0/2 + x3/2, // constant
539 176 : x0/2 - x2/2 - x3/2 + x5/2, // eta
540 176 : x0/2 - x1/2 - x3/2 + x4/2 // xi
541 : };
542 :
543 : // The quadrature rule for the Prism6 is a tensor product between a
544 : // four-point TRI3 rule (in xi, eta) and a two-point EDGE2 rule (in
545 : // zeta) which is capable of integrating cubics exactly.
546 :
547 : // Number of points in the 2D quadrature rule.
548 176 : const int N2D = 4;
549 :
550 : static const Real w2D[N2D] =
551 : {
552 : 1.5902069087198858469718450103758e-01_R,
553 : 9.0979309128011415302815498962418e-02_R,
554 : 1.5902069087198858469718450103758e-01_R,
555 : 9.0979309128011415302815498962418e-02_R
556 : };
557 :
558 : static const Real xi[N2D] =
559 : {
560 : 1.5505102572168219018027159252941e-01_R,
561 : 6.4494897427831780981972840747059e-01_R,
562 : 1.5505102572168219018027159252941e-01_R,
563 : 6.4494897427831780981972840747059e-01_R
564 : };
565 :
566 : static const Real eta[N2D] =
567 : {
568 : 1.7855872826361642311703513337422e-01_R,
569 : 7.5031110222608118177475598324603e-02_R,
570 : 6.6639024601470138670269327409637e-01_R,
571 : 2.8001991549907407200279599420481e-01_R
572 : };
573 :
574 : // Number of points in the 1D quadrature rule. The weights of the
575 : // 1D quadrature rule are equal to 1.
576 176 : const int N1D = 2;
577 :
578 : // Points of the 1D quadrature rule
579 : static const Real zeta[N1D] =
580 : {
581 : -std::sqrt(Real(3))/3,
582 : std::sqrt(Real(3))/3.
583 : };
584 :
585 18080 : for (int i=0; i<N2D; ++i)
586 : {
587 : // dx_dzeta depends only on the 2D quadrature rule points.
588 2112 : Point dx_dzeta_q = dx_dzeta[0] + eta[i]*dx_dzeta[1] + xi[i]*dx_dzeta[2];
589 :
590 43392 : for (int j=0; j<N1D; ++j)
591 : {
592 : // dx_dxi and dx_deta only depend on the 1D quadrature rule points.
593 : Point
594 2816 : dx_dxi_q = dx_dxi[0] + zeta[j]*dx_dxi[1],
595 1408 : dx_deta_q = dx_deta[0] + zeta[j]*dx_deta[1];
596 :
597 : // Compute t-specific contributions at current qp
598 30336 : t.accumulate_qp(xi[i], eta[i], zeta[j], w2D[i] * triple_product(dx_dxi_q, dx_deta_q, dx_dzeta_q));
599 : }
600 : }
601 3616 : }
602 :
603 : }
|
