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