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