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_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 2544 : void accumulate_qp(Real /*xi*/, Real /*eta*/, Real /*zeta*/, Real JxW)
48 : {
49 36840 : vol += JxW;
50 2544 : };
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_elems_map[Prism6::num_sides][Prism6::nodes_per_side] =
114 : {
115 : {0, 1, 2, 3}, // Side 0
116 : {0, 1, 4, 5}, // Side 1
117 : {1, 2, 5, 6}, // Side 2
118 : {0, 2, 4, 6}, // Side 3
119 : {4, 5, 6, 7} // Side 4
120 : };
121 :
122 : // ------------------------------------------------------------
123 : // Prism6 class member functions
124 :
125 3627814 : bool Prism6::is_vertex(const unsigned int libmesh_dbg_var(n)) const
126 : {
127 158268 : libmesh_assert_not_equal_to (n, invalid_uint);
128 3627814 : return true;
129 : }
130 :
131 0 : bool Prism6::is_edge(const unsigned int) const
132 : {
133 0 : return false;
134 : }
135 :
136 :
137 0 : bool Prism6::is_face(const unsigned int) const
138 : {
139 0 : return false;
140 : }
141 :
142 129810 : bool Prism6::is_node_on_side(const unsigned int n,
143 : const unsigned int s) const
144 : {
145 10700 : libmesh_assert_less (s, n_sides());
146 10700 : return std::find(std::begin(side_nodes_map[s]),
147 129810 : std::end(side_nodes_map[s]),
148 129810 : n) != std::end(side_nodes_map[s]);
149 : }
150 :
151 : std::vector<unsigned>
152 53945216 : Prism6::nodes_on_side(const unsigned int s) const
153 : {
154 4819428 : libmesh_assert_less(s, n_sides());
155 53945216 : auto trim = (s > 0 && s < 4) ? 0 : 1;
156 53945216 : return {std::begin(side_nodes_map[s]), std::end(side_nodes_map[s]) - trim};
157 : }
158 :
159 : std::vector<unsigned>
160 11007 : Prism6::nodes_on_edge(const unsigned int e) const
161 : {
162 882 : libmesh_assert_less(e, n_edges());
163 11889 : return {std::begin(edge_nodes_map[e]), std::end(edge_nodes_map[e])};
164 : }
165 :
166 1759482 : bool Prism6::is_node_on_edge(const unsigned int n,
167 : const unsigned int e) const
168 : {
169 65196 : libmesh_assert_less (e, n_edges());
170 65196 : return std::find(std::begin(edge_nodes_map[e]),
171 1759482 : std::end(edge_nodes_map[e]),
172 1759482 : n) != std::end(edge_nodes_map[e]);
173 : }
174 :
175 :
176 :
177 592170 : bool Prism6::has_affine_map() const
178 : {
179 : // Make sure z edges are affine
180 115484 : Point v = this->point(3) - this->point(0);
181 697030 : if (!v.relative_fuzzy_equals(this->point(4) - this->point(1), affine_tol) ||
182 636018 : !v.relative_fuzzy_equals(this->point(5) - this->point(2), affine_tol))
183 66456 : return false;
184 52298 : return true;
185 : }
186 :
187 :
188 :
189 12303537 : Order Prism6::default_order() const
190 : {
191 12303537 : return FIRST;
192 : }
193 :
194 :
195 :
196 175186 : std::unique_ptr<Elem> Prism6::build_side_ptr (const unsigned int i)
197 : {
198 13965 : libmesh_assert_less (i, this->n_sides());
199 :
200 175186 : std::unique_ptr<Elem> face;
201 :
202 175186 : switch (i)
203 : {
204 68394 : case 0: // the triangular face at z=-1
205 : case 4: // the triangular face at z=1
206 : {
207 68394 : face = std::make_unique<Tri3>();
208 68394 : break;
209 : }
210 106792 : case 1: // the quad face at y=0
211 : case 2: // the other quad face
212 : case 3: // the quad face at x=0
213 : {
214 106792 : face = std::make_unique<Quad4>();
215 106792 : break;
216 : }
217 0 : default:
218 0 : libmesh_error_msg("Invalid side i = " << i);
219 : }
220 :
221 : // Set the nodes
222 807536 : for (auto n : face->node_index_range())
223 682762 : face->set_node(n, this->node_ptr(Prism6::side_nodes_map[i][n]));
224 :
225 175186 : face->set_interior_parent(this);
226 161221 : face->inherit_data_from(*this);
227 :
228 175186 : return face;
229 0 : }
230 :
231 :
232 :
233 7835 : void Prism6::build_side_ptr (std::unique_ptr<Elem> & side,
234 : const unsigned int i)
235 : {
236 7835 : this->side_ptr(side, i);
237 7835 : side->set_interior_parent(this);
238 6993 : side->inherit_data_from(*this);
239 7835 : }
240 :
241 :
242 :
243 1612975 : std::unique_ptr<Elem> Prism6::build_edge_ptr (const unsigned int i)
244 : {
245 1612975 : return this->simple_build_edge_ptr<Edge2,Prism6>(i);
246 : }
247 :
248 :
249 :
250 0 : void Prism6::build_edge_ptr (std::unique_ptr<Elem> & edge, const unsigned int i)
251 : {
252 0 : this->simple_build_edge_ptr<Prism6>(edge, i, EDGE2);
253 0 : }
254 :
255 :
256 :
257 23040 : void Prism6::connectivity(const unsigned int libmesh_dbg_var(sc),
258 : const IOPackage iop,
259 : std::vector<dof_id_type> & conn) const
260 : {
261 1920 : libmesh_assert(_nodes);
262 1920 : libmesh_assert_less (sc, this->n_sub_elem());
263 1920 : libmesh_assert_not_equal_to (iop, INVALID_IO_PACKAGE);
264 :
265 23040 : switch (iop)
266 : {
267 23040 : case TECPLOT:
268 : {
269 23040 : conn.resize(8);
270 24960 : conn[0] = this->node_id(0)+1;
271 23040 : conn[1] = this->node_id(1)+1;
272 23040 : conn[2] = this->node_id(2)+1;
273 23040 : conn[3] = this->node_id(2)+1;
274 23040 : conn[4] = this->node_id(3)+1;
275 23040 : conn[5] = this->node_id(4)+1;
276 23040 : conn[6] = this->node_id(5)+1;
277 23040 : conn[7] = this->node_id(5)+1;
278 23040 : return;
279 : }
280 :
281 0 : case VTK:
282 : {
283 0 : conn.resize(6);
284 0 : conn[0] = this->node_id(0);
285 0 : conn[1] = this->node_id(2);
286 0 : conn[2] = this->node_id(1);
287 0 : conn[3] = this->node_id(3);
288 0 : conn[4] = this->node_id(5);
289 0 : conn[5] = this->node_id(4);
290 0 : return;
291 : }
292 :
293 0 : default:
294 0 : libmesh_error_msg("Unsupported IO package " << iop);
295 : }
296 : }
297 :
298 :
299 :
300 : #ifdef LIBMESH_ENABLE_AMR
301 :
302 : const Real Prism6::_embedding_matrix[Prism6::num_children][Prism6::num_nodes][Prism6::num_nodes] =
303 : {
304 : // embedding matrix for child 0
305 : {
306 : // 0 1 2 3 4 5
307 : { 1.0, 0.0, 0.0, 0.0, 0.0, 0.0}, // 0
308 : { 0.5, 0.5, 0.0, 0.0, 0.0, 0.0}, // 1
309 : { 0.5, 0.0, 0.5, 0.0, 0.0, 0.0}, // 2
310 : { 0.5, 0.0, 0.0, 0.5, 0.0, 0.0}, // 3
311 : { .25, .25, 0.0, .25, .25, 0.0}, // 4
312 : { .25, 0.0, .25, .25, 0.0, .25} // 5
313 : },
314 :
315 : // embedding matrix for child 1
316 : {
317 : // 0 1 2 3 4 5
318 : { 0.5, 0.5, 0.0, 0.0, 0.0, 0.0}, // 0
319 : { 0.0, 1.0, 0.0, 0.0, 0.0, 0.0}, // 1
320 : { 0.0, 0.5, 0.5, 0.0, 0.0, 0.0}, // 2
321 : { .25, .25, 0.0, .25, .25, 0.0}, // 3
322 : { 0.0, 0.5, 0.0, 0.0, 0.5, 0.0}, // 4
323 : { 0.0, .25, .25, 0.0, .25, .25} // 5
324 : },
325 :
326 : // embedding matrix for child 2
327 : {
328 : // 0 1 2 3 4 5
329 : { 0.5, 0.0, 0.5, 0.0, 0.0, 0.0}, // 0
330 : { 0.0, 0.5, 0.5, 0.0, 0.0, 0.0}, // 1
331 : { 0.0, 0.0, 1.0, 0.0, 0.0, 0.0}, // 2
332 : { .25, 0.0, .25, .25, 0.0, .25}, // 3
333 : { 0.0, .25, .25, 0.0, .25, .25}, // 4
334 : { 0.0, 0.0, 0.5, 0.0, 0.0, 0.5} // 5
335 : },
336 :
337 : // embedding matrix for child 3
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, 0.5, 0.5, 0.0, 0.0, 0.0}, // 1
342 : { 0.5, 0.0, 0.5, 0.0, 0.0, 0.0}, // 2
343 : { .25, .25, 0.0, .25, .25, 0.0}, // 3
344 : { 0.0, .25, .25, 0.0, .25, .25}, // 4
345 : { .25, 0.0, .25, .25, 0.0, .25} // 5
346 : },
347 :
348 : // embedding matrix for child 4
349 : {
350 : // 0 1 2 3 4 5
351 : { 0.5, 0.0, 0.0, 0.5, 0.0, 0.0}, // 0
352 : { .25, .25, 0.0, .25, .25, 0.0}, // 1
353 : { .25, 0.0, .25, .25, 0.0, .25}, // 2
354 : { 0.0, 0.0, 0.0, 1.0, 0.0, 0.0}, // 3
355 : { 0.0, 0.0, 0.0, 0.5, 0.5, 0.0}, // 4
356 : { 0.0, 0.0, 0.0, 0.5, 0.0, 0.5} // 5
357 : },
358 :
359 : // embedding matrix for child 5
360 : {
361 : // 0 1 2 3 4 5
362 : { .25, .25, 0.0, .25, .25, 0.0}, // 0
363 : { 0.0, 0.5, 0.0, 0.0, 0.5, 0.0}, // 1
364 : { 0.0, .25, .25, 0.0, .25, .25}, // 2
365 : { 0.0, 0.0, 0.0, 0.5, 0.5, 0.0}, // 3
366 : { 0.0, 0.0, 0.0, 0.0, 1.0, 0.0}, // 4
367 : { 0.0, 0.0, 0.0, 0.0, 0.5, 0.5} // 5
368 : },
369 :
370 : // embedding matrix for child 6
371 : {
372 : // 0 1 2 3 4 5
373 : { .25, 0.0, .25, .25, 0.0, .25}, // 0
374 : { 0.0, .25, .25, 0.0, .25, .25}, // 1
375 : { 0.0, 0.0, 0.5, 0.0, 0.0, 0.5}, // 2
376 : { 0.0, 0.0, 0.0, 0.5, 0.0, 0.5}, // 3
377 : { 0.0, 0.0, 0.0, 0.0, 0.5, 0.5}, // 4
378 : { 0.0, 0.0, 0.0, 0.0, 0.0, 1.0} // 5
379 : },
380 :
381 : // embedding matrix for child 7
382 : {
383 : // 0 1 2 3 4 5
384 : { .25, .25, 0.0, .25, .25, 0.0}, // 0
385 : { 0.0, .25, .25, 0.0, .25, .25}, // 1
386 : { .25, 0.0, .25, .25, 0.0, .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, 0.5, 0.5}, // 4
389 : { 0.0, 0.0, 0.0, 0.5, 0.0, 0.5} // 5
390 : }
391 : };
392 :
393 : #endif
394 :
395 :
396 :
397 2480 : Point Prism6::true_centroid () const
398 : {
399 2480 : NodalVolumeIntegral nvi;
400 2480 : cell_integral(this, nvi);
401 :
402 : // Compute centroid
403 144 : Point cp;
404 144 : Real vol = 0.;
405 :
406 17360 : for (int i=0; i<6; ++i)
407 : {
408 15744 : cp += this->point(i) * nvi.V[i];
409 14880 : vol += nvi.V[i];
410 : }
411 :
412 2480 : return cp / vol;
413 : }
414 :
415 4605 : Real Prism6::volume () const
416 : {
417 : // This specialization is good for Lagrange mappings only in general
418 4605 : if (this->mapping_type() != LAGRANGE_MAP)
419 0 : return this->Elem::volume();
420 :
421 4605 : VolumeIntegral vi;
422 4605 : cell_integral(this, vi);
423 4605 : return vi.vol;
424 : }
425 :
426 : BoundingBox
427 197247 : Prism6::loose_bounding_box () const
428 : {
429 197247 : return Elem::loose_bounding_box();
430 : }
431 :
432 :
433 : void
434 8208 : Prism6::permute(unsigned int perm_num)
435 : {
436 720 : libmesh_assert_less (perm_num, 6);
437 8208 : const unsigned int side = perm_num % 2;
438 8208 : const unsigned int rotate = perm_num / 2;
439 :
440 18792 : for (unsigned int i = 0; i != rotate; ++i)
441 : {
442 10584 : swap3nodes(0,1,2);
443 9648 : swap3nodes(3,4,5);
444 9648 : swap3neighbors(1,2,3);
445 : }
446 :
447 8208 : switch (side) {
448 144 : case 0:
449 144 : break;
450 6480 : case 1:
451 6480 : swap2nodes(1,3);
452 6480 : swap2nodes(0,4);
453 6480 : swap2nodes(2,5);
454 576 : swap2neighbors(0,4);
455 576 : swap2neighbors(2,3);
456 576 : break;
457 0 : default:
458 0 : libmesh_error();
459 : }
460 :
461 8208 : }
462 :
463 :
464 : void
465 576 : Prism6::flip(BoundaryInfo * boundary_info)
466 : {
467 48 : libmesh_assert(boundary_info);
468 :
469 576 : swap2nodes(0,1);
470 576 : swap2nodes(3,4);
471 48 : swap2neighbors(2,3);
472 576 : swap2boundarysides(2,3,boundary_info);
473 576 : swap2boundaryedges(0,1,boundary_info);
474 576 : swap2boundaryedges(3,4,boundary_info);
475 576 : swap2boundaryedges(7,8,boundary_info);
476 576 : }
477 :
478 :
479 : ElemType
480 107663544 : Prism6::side_type (const unsigned int s) const
481 : {
482 9621248 : libmesh_assert_less (s, 5);
483 107663544 : if (s == 0 || s == 4)
484 18124896 : return TRI3;
485 8001464 : return QUAD4;
486 : }
487 :
488 : Point
489 355 : Prism6::side_vertex_average_normal(const unsigned int s) const
490 : {
491 10 : libmesh_assert_less (s, 5);
492 10 : libmesh_assert_equal_to(this->mapping_type(), LAGRANGE_MAP);
493 :
494 355 : switch (s)
495 : {
496 142 : case 0: // the triangular face at z=-1
497 : case 4: // the triangular face at z=1
498 : {
499 150 : const Point n1 = this->point(side_nodes_map[s][1]) -
500 142 : this->point(side_nodes_map[s][0]);
501 150 : const Point n2 = this->point(side_nodes_map[s][2]) -
502 4 : this->point(side_nodes_map[s][1]);
503 138 : const Point pointing_out = n1.cross(n2);
504 142 : return pointing_out.unit();
505 : }
506 6 : case 1: // the quad face at y=0
507 : case 2: // the other quad face
508 : case 3: // the quad face at x=0
509 : {
510 : // At the side vertex average, things simplify a bit
511 : // We get the side "plane" normal at all vertices, then average them
512 6 : Point normal;
513 225 : Point current_edge = this->point(side_nodes_map[s][1]) -
514 213 : this->point(side_nodes_map[s][0]);
515 1065 : for (auto i : make_range(4))
516 : {
517 900 : const Point next_edge = this->point(side_nodes_map[s][(i + 2) % 4]) -
518 852 : this->point(side_nodes_map[s][(i + 1) % 4]);
519 828 : const Point normal_at_vertex = current_edge.cross(next_edge);
520 24 : normal += normal_at_vertex;
521 828 : current_edge = next_edge;
522 : }
523 6 : normal /= 4;
524 213 : return normal.unit();
525 : }
526 0 : default:
527 0 : libmesh_error_msg("Invalid side s = " << s);
528 : }
529 : }
530 :
531 : } // namespace libMesh
532 :
533 : namespace // anonymous
534 : {
535 :
536 : template <typename T>
537 7085 : void cell_integral(const Elem * elem, T & t)
538 : {
539 : // Conveient references to our points.
540 : const Point
541 920 : &x0 = elem->point(0), &x1 = elem->point(1), &x2 = elem->point(2),
542 462 : &x3 = elem->point(3), &x4 = elem->point(4), &x5 = elem->point(5);
543 :
544 : // constant and zeta terms only. These are copied directly from a
545 : // Python script.
546 924 : Point dx_dxi[2] =
547 : {
548 462 : -x0/2 + x1/2 - x3/2 + x4/2, // constant
549 462 : x0/2 - x1/2 - x3/2 + x4/2, // zeta
550 : };
551 :
552 : // constant and zeta terms only. These are copied directly from a
553 : // Python script.
554 924 : Point dx_deta[2] =
555 : {
556 462 : -x0/2 + x2/2 - x3/2 + x5/2, // constant
557 462 : x0/2 - x2/2 - x3/2 + x5/2, // zeta
558 : };
559 :
560 : // Constant, xi, and eta terms
561 1386 : Point dx_dzeta[3] =
562 : {
563 462 : -x0/2 + x3/2, // constant
564 462 : x0/2 - x2/2 - x3/2 + x5/2, // eta
565 462 : x0/2 - x1/2 - x3/2 + x4/2 // xi
566 : };
567 :
568 : // The quadrature rule for the Prism6 is a tensor product between a
569 : // four-point TRI3 rule (in xi, eta) and a two-point EDGE2 rule (in
570 : // zeta) which is capable of integrating cubics exactly.
571 :
572 : // Number of points in the 2D quadrature rule.
573 462 : const int N2D = 4;
574 :
575 : static const Real w2D[N2D] =
576 : {
577 : 1.5902069087198858469718450103758e-01_R,
578 : 9.0979309128011415302815498962418e-02_R,
579 : 1.5902069087198858469718450103758e-01_R,
580 : 9.0979309128011415302815498962418e-02_R
581 : };
582 :
583 : static const Real xi[N2D] =
584 : {
585 : 1.5505102572168219018027159252941e-01_R,
586 : 6.4494897427831780981972840747059e-01_R,
587 : 1.5505102572168219018027159252941e-01_R,
588 : 6.4494897427831780981972840747059e-01_R
589 : };
590 :
591 : static const Real eta[N2D] =
592 : {
593 : 1.7855872826361642311703513337422e-01_R,
594 : 7.5031110222608118177475598324603e-02_R,
595 : 6.6639024601470138670269327409637e-01_R,
596 : 2.8001991549907407200279599420481e-01_R
597 : };
598 :
599 : // Number of points in the 1D quadrature rule. The weights of the
600 : // 1D quadrature rule are equal to 1.
601 462 : const int N1D = 2;
602 :
603 : // Points of the 1D quadrature rule
604 : static const Real zeta[N1D] =
605 : {
606 : -std::sqrt(Real(3))/3,
607 : std::sqrt(Real(3))/3.
608 : };
609 :
610 35425 : for (int i=0; i<N2D; ++i)
611 : {
612 : // dx_dzeta depends only on the 2D quadrature rule points.
613 5512 : Point dx_dzeta_q = dx_dzeta[0] + eta[i]*dx_dzeta[1] + xi[i]*dx_dzeta[2];
614 :
615 85020 : for (int j=0; j<N1D; ++j)
616 : {
617 : // dx_dxi and dx_deta only depend on the 1D quadrature rule points.
618 : Point
619 7360 : dx_dxi_q = dx_dxi[0] + zeta[j]*dx_dxi[1],
620 3696 : dx_deta_q = dx_deta[0] + zeta[j]*dx_deta[1];
621 :
622 : // Compute t-specific contributions at current qp
623 60344 : t.accumulate_qp(xi[i], eta[i], zeta[j], w2D[i] * triple_product(dx_dxi_q, dx_deta_q, dx_dzeta_q));
624 : }
625 : }
626 7085 : }
627 :
628 : }
|
