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_prism20.h"
21 : #include "libmesh/edge_edge3.h"
22 : #include "libmesh/face_quad9.h"
23 : #include "libmesh/face_tri7.h"
24 : #include "libmesh/enum_io_package.h"
25 : #include "libmesh/enum_order.h"
26 : #include "libmesh/int_range.h"
27 :
28 : namespace libMesh
29 : {
30 :
31 :
32 :
33 : // ------------------------------------------------------------
34 : // Prism20 class static member initializations
35 : const int Prism20::num_nodes;
36 : const int Prism20::nodes_per_side;
37 : const int Prism20::nodes_per_edge;
38 :
39 : const unsigned int Prism20::side_nodes_map[Prism20::num_sides][Prism20::nodes_per_side] =
40 : {
41 : {0, 2, 1, 8, 7, 6, 18, 99, 99}, // Side 0
42 : {0, 1, 4, 3, 6, 10, 12, 9, 15}, // Side 1
43 : {1, 2, 5, 4, 7, 11, 13, 10, 16}, // Side 2
44 : {2, 0, 3, 5, 8, 9, 14, 11, 17}, // Side 3
45 : {3, 4, 5, 12, 13, 14, 19, 99, 99} // Side 4
46 : };
47 :
48 : const unsigned int Prism20::edge_nodes_map[Prism20::num_edges][Prism20::nodes_per_edge] =
49 : {
50 : {0, 1, 6}, // Edge 0
51 : {1, 2, 7}, // Edge 1
52 : {0, 2, 8}, // Edge 2
53 : {0, 3, 9}, // Edge 3
54 : {1, 4, 10}, // Edge 4
55 : {2, 5, 11}, // Edge 5
56 : {3, 4, 12}, // Edge 6
57 : {4, 5, 13}, // Edge 7
58 : {3, 5, 14} // Edge 8
59 : };
60 :
61 : // ------------------------------------------------------------
62 : // Prism20 class member functions
63 :
64 1734190 : bool Prism20::is_vertex(const unsigned int i) const
65 : {
66 1734190 : if (i < 6)
67 520608 : return true;
68 70877 : return false;
69 : }
70 :
71 16128 : bool Prism20::is_edge(const unsigned int i) const
72 : {
73 16128 : if (i < 6)
74 0 : return false;
75 16128 : if (i > 14)
76 5760 : return false;
77 864 : return true;
78 : }
79 :
80 1200145 : bool Prism20::is_face(const unsigned int i) const
81 : {
82 1200145 : if (i > 14)
83 193345 : return true;
84 81450 : return false;
85 : }
86 :
87 670765 : bool Prism20::is_node_on_side(const unsigned int n,
88 : const unsigned int s) const
89 : {
90 51620 : libmesh_assert_less (s, n_sides());
91 51620 : return std::find(std::begin(side_nodes_map[s]),
92 670765 : std::end(side_nodes_map[s]),
93 670765 : n) != std::end(side_nodes_map[s]);
94 : }
95 :
96 : std::vector<unsigned>
97 107915869 : Prism20::nodes_on_side(const unsigned int s) const
98 : {
99 9622190 : libmesh_assert_less(s, n_sides());
100 107915869 : auto trim = (s > 0 && s < 4) ? 0 : 2;
101 107915869 : return {std::begin(side_nodes_map[s]), std::end(side_nodes_map[s]) - trim};
102 : }
103 :
104 : std::vector<unsigned>
105 11007 : Prism20::nodes_on_edge(const unsigned int e) const
106 : {
107 882 : libmesh_assert_less(e, n_edges());
108 11889 : return {std::begin(edge_nodes_map[e]), std::end(edge_nodes_map[e])};
109 : }
110 :
111 130284 : bool Prism20::is_node_on_edge(const unsigned int n,
112 : const unsigned int e) const
113 : {
114 10512 : libmesh_assert_less (e, n_edges());
115 10512 : return std::find(std::begin(edge_nodes_map[e]),
116 130284 : std::end(edge_nodes_map[e]),
117 130284 : n) != std::end(edge_nodes_map[e]);
118 : }
119 :
120 :
121 :
122 451408 : bool Prism20::has_affine_map() const
123 : {
124 : // Make sure z edges are affine
125 78708 : Point v = this->point(3) - this->point(0);
126 530116 : if (!v.relative_fuzzy_equals(this->point(4) - this->point(1), affine_tol) ||
127 530116 : !v.relative_fuzzy_equals(this->point(5) - this->point(2), affine_tol))
128 0 : return false;
129 :
130 : // Make sure edges are straight
131 39354 : v /= 2;
132 490762 : if (!v.relative_fuzzy_equals(this->point(9) - this->point(0), affine_tol) ||
133 490762 : !v.relative_fuzzy_equals(this->point(10) - this->point(1), affine_tol) ||
134 490762 : !v.relative_fuzzy_equals(this->point(11) - this->point(2), affine_tol) ||
135 490762 : !v.relative_fuzzy_equals(this->point(15) - this->point(6), affine_tol) ||
136 981524 : !v.relative_fuzzy_equals(this->point(16) - this->point(7), affine_tol) ||
137 490762 : !v.relative_fuzzy_equals(this->point(17) - this->point(8), affine_tol))
138 0 : return false;
139 490762 : v = (this->point(1) - this->point(0))/2;
140 530116 : if (!v.relative_fuzzy_equals(this->point(6) - this->point(0), affine_tol) ||
141 490762 : !v.relative_fuzzy_equals(this->point(12) - this->point(3), affine_tol))
142 0 : return false;
143 490762 : v = (this->point(2) - this->point(0))/2;
144 530116 : if (!v.relative_fuzzy_equals(this->point(8) - this->point(0), affine_tol) ||
145 490762 : !v.relative_fuzzy_equals(this->point(14) - this->point(3), affine_tol))
146 0 : return false;
147 490762 : v = (this->point(2) - this->point(1))/2;
148 530116 : if (!v.relative_fuzzy_equals(this->point(7) - this->point(1), affine_tol) ||
149 490762 : !v.relative_fuzzy_equals(this->point(13) - this->point(4), affine_tol))
150 0 : return false;
151 :
152 : // Make sure triangle face midpoints are centered
153 490762 : v = this->point(2) + this->point(1) - 2*this->point(0);
154 451408 : if (!v.relative_fuzzy_equals((this->point(18)-this->point(0))*3))
155 0 : return false;
156 490762 : v = this->point(5) + this->point(4) - 2*this->point(3);
157 451408 : if (!v.relative_fuzzy_equals((this->point(19)-this->point(3))*3))
158 0 : return false;
159 :
160 39354 : return true;
161 : }
162 :
163 :
164 :
165 14245173 : Order Prism20::default_order() const
166 : {
167 14245173 : return THIRD;
168 : }
169 :
170 0 : dof_id_type Prism20::key (const unsigned int s) const
171 : {
172 0 : libmesh_assert_less (s, this->n_sides());
173 :
174 0 : switch (s)
175 : {
176 0 : case 0: // the triangular face at z=0
177 : {
178 0 : return Elem::compute_key (this->node_id(18));
179 : }
180 0 : case 1: // the quad face at y=0
181 : {
182 0 : return Elem::compute_key (this->node_id(15));
183 : }
184 0 : case 2: // the other quad face
185 : {
186 0 : return Elem::compute_key (this->node_id(16));
187 : }
188 0 : case 3: // the quad face at x=0
189 : {
190 0 : return Elem::compute_key (this->node_id(17));
191 : }
192 0 : case 4: // the triangular face at z=1
193 : {
194 0 : return Elem::compute_key (this->node_id(19));
195 : }
196 0 : default:
197 0 : libmesh_error_msg("Invalid side " << s);
198 : }
199 : }
200 :
201 :
202 :
203 388736 : unsigned int Prism20::local_side_node(unsigned int side,
204 : unsigned int side_node) const
205 : {
206 30962 : libmesh_assert_less (side, this->n_sides());
207 :
208 : // Never more than 9 nodes per side.
209 30962 : libmesh_assert_less(side_node, Prism20::nodes_per_side);
210 :
211 : // Some sides have 7 nodes.
212 30962 : libmesh_assert(!(side==0 || side==4) || side_node < 7);
213 :
214 388736 : return Prism20::side_nodes_map[side][side_node];
215 : }
216 :
217 :
218 :
219 86875929 : unsigned int Prism20::local_edge_node(unsigned int edge,
220 : unsigned int edge_node) const
221 : {
222 7744896 : libmesh_assert_less(edge, this->n_edges());
223 7744896 : libmesh_assert_less(edge_node, Prism20::nodes_per_edge);
224 :
225 86875929 : return Prism20::edge_nodes_map[edge][edge_node];
226 : }
227 :
228 :
229 :
230 354400 : std::unique_ptr<Elem> Prism20::build_side_ptr (const unsigned int i)
231 : {
232 13676 : libmesh_assert_less (i, this->n_sides());
233 :
234 354400 : std::unique_ptr<Elem> face;
235 :
236 354400 : switch (i)
237 : {
238 222628 : case 0: // the triangular face at z=-1
239 : case 4: // the triangular face at z=1
240 : {
241 222628 : face = std::make_unique<Tri7>();
242 222628 : break;
243 : }
244 131772 : case 1: // the quad face at y=0
245 : case 2: // the other quad face
246 : case 3: // the quad face at x=0
247 : {
248 131772 : face = std::make_unique<Quad9>();
249 131772 : break;
250 : }
251 0 : default:
252 0 : libmesh_error_msg("Invalid side i = " << i);
253 : }
254 :
255 : // Set the nodes
256 3098744 : for (auto n : face->node_index_range())
257 2851356 : face->set_node(n, this->node_ptr(Prism20::side_nodes_map[i][n]));
258 :
259 354400 : face->set_interior_parent(this);
260 340724 : face->inherit_data_from(*this);
261 :
262 354400 : return face;
263 0 : }
264 :
265 :
266 :
267 405907 : void Prism20::build_side_ptr (std::unique_ptr<Elem> & side,
268 : const unsigned int i)
269 : {
270 13010 : libmesh_assert_less (i, this->n_sides());
271 :
272 405907 : switch (i)
273 : {
274 3294 : case 0: // the triangular face at z=-1
275 : case 4: // the triangular face at z=1
276 : {
277 102104 : if (!side.get() || side->type() != TRI7)
278 : {
279 196212 : side = this->build_side_ptr(i);
280 101336 : return;
281 : }
282 64 : break;
283 : }
284 :
285 9716 : case 1: // the quad face at y=0
286 : case 2: // the other quad face
287 : case 3: // the quad face at x=0
288 : {
289 303803 : if (!side.get() || side->type() != QUAD9)
290 : {
291 195584 : side = this->build_side_ptr(i);
292 101002 : return;
293 : }
294 6506 : break;
295 : }
296 :
297 0 : default:
298 0 : libmesh_error_msg("Invalid side i = " << i);
299 : }
300 :
301 196999 : side->inherit_data_from(*this);
302 :
303 : // Set the nodes
304 2034154 : for (auto n : side->node_index_range())
305 1889587 : side->set_node(n, this->node_ptr(Prism20::side_nodes_map[i][n]));
306 : }
307 :
308 :
309 :
310 7533 : std::unique_ptr<Elem> Prism20::build_edge_ptr (const unsigned int i)
311 : {
312 7533 : return this->simple_build_edge_ptr<Edge3,Prism20>(i);
313 : }
314 :
315 :
316 :
317 0 : void Prism20::build_edge_ptr (std::unique_ptr<Elem> & edge, const unsigned int i)
318 : {
319 0 : this->simple_build_edge_ptr<Prism20>(edge, i, EDGE3);
320 0 : }
321 :
322 :
323 :
324 0 : void Prism20::connectivity(const unsigned int /*sc*/,
325 : const IOPackage /*iop*/,
326 : std::vector<dof_id_type> & /*conn*/) const
327 : {
328 0 : libmesh_not_implemented(); // FIXME RHS
329 :
330 : /*
331 : libmesh_assert(_nodes);
332 : libmesh_assert_less (sc, this->n_sub_elem());
333 : libmesh_assert_not_equal_to (iop, INVALID_IO_PACKAGE);
334 :
335 : switch (iop)
336 : {
337 : case TECPLOT:
338 : {
339 : conn.resize(8);
340 : switch (sc)
341 : {
342 :
343 : case 0:
344 : {
345 : conn[0] = this->node_id(0)+1;
346 : conn[1] = this->node_id(6)+1;
347 : conn[2] = this->node_id(8)+1;
348 : conn[3] = this->node_id(8)+1;
349 : conn[4] = this->node_id(9)+1;
350 : conn[5] = this->node_id(15)+1;
351 : conn[6] = this->node_id(17)+1;
352 : conn[7] = this->node_id(17)+1;
353 :
354 : return;
355 : }
356 :
357 : case 1:
358 : {
359 : conn[0] = this->node_id(6)+1;
360 : conn[1] = this->node_id(1)+1;
361 : conn[2] = this->node_id(7)+1;
362 : conn[3] = this->node_id(7)+1;
363 : conn[4] = this->node_id(15)+1;
364 : conn[5] = this->node_id(10)+1;
365 : conn[6] = this->node_id(16)+1;
366 : conn[7] = this->node_id(16)+1;
367 :
368 : return;
369 : }
370 :
371 : case 2:
372 : {
373 : conn[0] = this->node_id(8)+1;
374 : conn[1] = this->node_id(7)+1;
375 : conn[2] = this->node_id(2)+1;
376 : conn[3] = this->node_id(2)+1;
377 : conn[4] = this->node_id(17)+1;
378 : conn[5] = this->node_id(16)+1;
379 : conn[6] = this->node_id(11)+1;
380 : conn[7] = this->node_id(11)+1;
381 :
382 : return;
383 : }
384 :
385 : case 3:
386 : {
387 : conn[0] = this->node_id(6)+1;
388 : conn[1] = this->node_id(7)+1;
389 : conn[2] = this->node_id(8)+1;
390 : conn[3] = this->node_id(8)+1;
391 : conn[4] = this->node_id(15)+1;
392 : conn[5] = this->node_id(16)+1;
393 : conn[6] = this->node_id(17)+1;
394 : conn[7] = this->node_id(17)+1;
395 :
396 : return;
397 : }
398 :
399 : case 4:
400 : {
401 : conn[0] = this->node_id(9)+1;
402 : conn[1] = this->node_id(15)+1;
403 : conn[2] = this->node_id(17)+1;
404 : conn[3] = this->node_id(17)+1;
405 : conn[4] = this->node_id(3)+1;
406 : conn[5] = this->node_id(12)+1;
407 : conn[6] = this->node_id(14)+1;
408 : conn[7] = this->node_id(14)+1;
409 :
410 : return;
411 : }
412 :
413 : case 5:
414 : {
415 : conn[0] = this->node_id(15)+1;
416 : conn[1] = this->node_id(10)+1;
417 : conn[2] = this->node_id(16)+1;
418 : conn[3] = this->node_id(16)+1;
419 : conn[4] = this->node_id(12)+1;
420 : conn[5] = this->node_id(4)+1;
421 : conn[6] = this->node_id(13)+1;
422 : conn[7] = this->node_id(13)+1;
423 :
424 : return;
425 : }
426 :
427 : case 6:
428 : {
429 : conn[0] = this->node_id(17)+1;
430 : conn[1] = this->node_id(16)+1;
431 : conn[2] = this->node_id(11)+1;
432 : conn[3] = this->node_id(11)+1;
433 : conn[4] = this->node_id(14)+1;
434 : conn[5] = this->node_id(13)+1;
435 : conn[6] = this->node_id(5)+1;
436 : conn[7] = this->node_id(5)+1;
437 :
438 : return;
439 : }
440 :
441 : case 7:
442 : {
443 : conn[0] = this->node_id(15)+1;
444 : conn[1] = this->node_id(16)+1;
445 : conn[2] = this->node_id(17)+1;
446 : conn[3] = this->node_id(17)+1;
447 : conn[4] = this->node_id(12)+1;
448 : conn[5] = this->node_id(13)+1;
449 : conn[6] = this->node_id(14)+1;
450 : conn[7] = this->node_id(14)+1;
451 :
452 : return;
453 : }
454 :
455 : default:
456 : libmesh_error_msg("Invalid sc = " << sc);
457 : }
458 :
459 : }
460 :
461 : case VTK:
462 : {
463 : // VTK now supports VTK_BIQUADRATIC_QUADRATIC_WEDGE directly
464 : const unsigned int conn_size = 18;
465 : conn.resize(conn_size);
466 :
467 : // VTK's VTK_BIQUADRATIC_QUADRATIC_WEDGE first 9 (vertex) and
468 : // last 3 (mid-face) nodes match. The middle and top layers
469 : // of mid-edge nodes are reversed from LibMesh's.
470 : for (auto i : index_range(conn))
471 : conn[i] = this->node_id(i);
472 :
473 : // top "ring" of mid-edge nodes
474 : conn[9] = this->node_id(12);
475 : conn[10] = this->node_id(13);
476 : conn[11] = this->node_id(14);
477 :
478 : // middle "ring" of mid-edge nodes
479 : conn[12] = this->node_id(9);
480 : conn[13] = this->node_id(10);
481 : conn[14] = this->node_id(11);
482 :
483 : return;
484 : }
485 :
486 : default:
487 : libmesh_error_msg("Unsupported IO package " << iop);
488 : }
489 : */
490 : }
491 :
492 :
493 :
494 :
495 0 : unsigned int Prism20::n_second_order_adjacent_vertices (const unsigned int n) const
496 : {
497 0 : switch (n)
498 : {
499 0 : case 6:
500 : case 7:
501 : case 8:
502 : case 9:
503 : case 10:
504 : case 11:
505 : case 12:
506 : case 13:
507 : case 14:
508 0 : return 2;
509 :
510 0 : case 15:
511 : case 16:
512 : case 17:
513 0 : return 4;
514 :
515 0 : case 18:
516 : case 19:
517 0 : return 3;
518 :
519 0 : default:
520 0 : libmesh_error_msg("Invalid node n = " << n);
521 : }
522 : }
523 :
524 :
525 :
526 :
527 :
528 0 : unsigned short int Prism20::second_order_adjacent_vertex (const unsigned int n,
529 : const unsigned int v) const
530 : {
531 0 : libmesh_assert_greater_equal (n, this->n_vertices());
532 0 : libmesh_assert_less (n, this->n_nodes());
533 :
534 0 : switch (n)
535 : {
536 : /*
537 : * These nodes are unique to \p Prism20,
538 : * let our _remaining_... matrix handle
539 : * this.
540 : */
541 0 : case 15:
542 : case 16:
543 : case 17:
544 : case 18:
545 : case 19:
546 : {
547 0 : libmesh_assert_less (v, 4);
548 0 : return _remaining_second_order_adjacent_vertices[n-15][v];
549 : }
550 :
551 : /*
552 : * All other second-order nodes (6,...,14) are
553 : * identical with Prism15 and are therefore
554 : * delegated to the _second_order matrix of
555 : * \p Prism
556 : */
557 0 : default:
558 : {
559 0 : libmesh_assert_less (v, 2);
560 0 : return _second_order_adjacent_vertices[n-this->n_vertices()][v];
561 : }
562 :
563 : }
564 :
565 : return static_cast<unsigned short int>(-1);
566 : }
567 :
568 :
569 :
570 : const unsigned short int Prism20::_remaining_second_order_adjacent_vertices[5][4] =
571 : {
572 : { 0, 1, 3, 4}, // vertices adjacent to node 15
573 : { 1, 2, 4, 5}, // vertices adjacent to node 16
574 : { 0, 2, 3, 5}, // vertices adjacent to node 17
575 : { 0, 1, 2, 99}, // vertices adjacent to node 18
576 : { 3, 4, 5, 99} // vertices adjacent to node 19
577 : };
578 :
579 :
580 :
581 : std::pair<unsigned short int, unsigned short int>
582 0 : Prism20::second_order_child_vertex (const unsigned int n) const
583 : {
584 0 : libmesh_assert_greater_equal (n, this->n_vertices());
585 0 : libmesh_assert_less (n, this->n_nodes());
586 :
587 0 : return std::pair<unsigned short int, unsigned short int>
588 0 : (_second_order_vertex_child_number[n],
589 0 : _second_order_vertex_child_index[n]);
590 : }
591 :
592 :
593 :
594 : void
595 28912 : Prism20::permute(unsigned int perm_num)
596 : {
597 1256 : libmesh_assert_less (perm_num, 6);
598 28912 : const unsigned int side = perm_num % 2;
599 28912 : const unsigned int rotate = perm_num / 2;
600 :
601 64579 : for (unsigned int i = 0; i != rotate; ++i)
602 : {
603 35667 : swap3nodes(0,1,2);
604 34081 : swap3nodes(3,4,5);
605 34081 : swap3nodes(6,7,8);
606 34081 : swap3nodes(9,10,11);
607 34081 : swap3nodes(12,13,14);
608 34081 : swap3nodes(15,16,17);
609 34081 : swap3neighbors(1,2,3);
610 : }
611 :
612 28912 : switch (side) {
613 398 : case 0:
614 398 : break;
615 17911 : case 1:
616 17911 : swap2nodes(1,3);
617 17911 : swap2nodes(0,4);
618 17911 : swap2nodes(2,5);
619 17911 : swap2nodes(6,12);
620 17911 : swap2nodes(9,10);
621 17911 : swap2nodes(7,14);
622 17911 : swap2nodes(8,13);
623 17911 : swap2nodes(16,17);
624 17911 : swap2nodes(18,19);
625 858 : swap2neighbors(0,4);
626 858 : swap2neighbors(2,3);
627 858 : break;
628 0 : default:
629 0 : libmesh_error();
630 : }
631 28912 : }
632 :
633 :
634 : void
635 576 : Prism20::flip(BoundaryInfo * boundary_info)
636 : {
637 48 : libmesh_assert(boundary_info);
638 :
639 576 : swap2nodes(0,1);
640 576 : swap2nodes(3,4);
641 576 : swap2nodes(7,8);
642 576 : swap2nodes(9,10);
643 576 : swap2nodes(13,14);
644 576 : swap2nodes(16,17);
645 48 : swap2neighbors(2,3);
646 576 : swap2boundarysides(2,3,boundary_info);
647 576 : swap2boundaryedges(0,1,boundary_info);
648 576 : swap2boundaryedges(3,4,boundary_info);
649 576 : swap2boundaryedges(7,8,boundary_info);
650 576 : }
651 :
652 :
653 960 : unsigned int Prism20::center_node_on_side(const unsigned short side) const
654 : {
655 80 : libmesh_assert_less (side, Prism20::num_sides);
656 960 : if (side >= 1 && side <= 3)
657 576 : return side + 14;
658 384 : if (side == 4)
659 192 : return 19;
660 16 : return 18;
661 : }
662 :
663 :
664 : ElemType
665 2880 : Prism20::side_type (const unsigned int s) const
666 : {
667 240 : libmesh_assert_less (s, 5);
668 2880 : if (s == 0 || s == 4)
669 1152 : return TRI7;
670 144 : return QUAD9;
671 : }
672 :
673 :
674 : } // namespace libMesh
|
