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 : // Reference-element topology and node locations shared between the host
19 : // element classes and Kokkos device code: constexpr side/edge tables and
20 : // lookups, with second-order side rows and higher-order node coordinates
21 : // derived from the stored linear facts.
22 :
23 : #ifndef LIBMESH_FE_REFERENCE_ELEMENT_TRAITS_H
24 : #define LIBMESH_FE_REFERENCE_ELEMENT_TRAITS_H
25 :
26 : #include "libmesh/enum_elem_type.h"
27 : #include "libmesh/libmesh.h"
28 : #include "libmesh/libmesh_device.h"
29 : #include "libmesh/point.h"
30 :
31 : namespace libMesh
32 : {
33 :
34 : template <unsigned int N>
35 : struct ReferenceElementVector
36 : {
37 : unsigned int values[N];
38 :
39 : LIBMESH_DEVICE_INLINE constexpr unsigned int operator[](unsigned int i) const
40 : { return values[i]; }
41 : };
42 :
43 : template <unsigned int Rows, unsigned int Cols>
44 : struct ReferenceElementTable
45 : {
46 : using Row = unsigned int[Cols];
47 :
48 : Row values[Rows];
49 :
50 359263658 : LIBMESH_DEVICE_INLINE constexpr const Row & operator[](unsigned int i) const
51 550847348 : { return values[i]; }
52 :
53 : LIBMESH_DEVICE_INLINE constexpr unsigned int operator()(unsigned int i, unsigned int j) const
54 : { return values[i][j]; }
55 : };
56 :
57 :
58 : LIBMESH_DEVICE_INLINE constexpr ReferenceElementTable<5, 4>
59 : prism6_side_nodes()
60 : {
61 : return {{
62 : {0, 2, 1, 99},
63 : {0, 1, 4, 3},
64 : {1, 2, 5, 4},
65 : {2, 0, 3, 5},
66 : {3, 4, 5, 99}
67 : }};
68 : }
69 :
70 : LIBMESH_DEVICE_INLINE constexpr ReferenceElementTable<5, 4>
71 : pyramid5_side_nodes()
72 : {
73 : return {{
74 : {0, 1, 4, 99},
75 : {1, 2, 4, 99},
76 : {2, 3, 4, 99},
77 : {3, 0, 4, 99},
78 : {0, 3, 2, 1}
79 : }};
80 : }
81 :
82 : LIBMESH_DEVICE_INLINE constexpr ReferenceElementTable<3, 2>
83 : tri3_side_nodes()
84 : {
85 : return {{
86 : {0, 1},
87 : {1, 2},
88 : {2, 0}
89 : }};
90 : }
91 :
92 : LIBMESH_DEVICE_INLINE constexpr ReferenceElementTable<4, 2>
93 : quad4_side_nodes()
94 : {
95 : return {{
96 : {0, 1},
97 : {1, 2},
98 : {2, 3},
99 : {3, 0}
100 : }};
101 : }
102 :
103 : LIBMESH_DEVICE_INLINE constexpr ReferenceElementTable<4, 3>
104 : tet4_side_nodes()
105 : {
106 : return {{
107 : {0, 2, 1},
108 : {0, 1, 3},
109 : {1, 2, 3},
110 : {2, 0, 3}
111 : }};
112 : }
113 :
114 : LIBMESH_DEVICE_INLINE constexpr ReferenceElementTable<6, 4>
115 : hex8_side_nodes()
116 : {
117 : return {{
118 : {0, 3, 2, 1},
119 : {0, 1, 5, 4},
120 : {1, 2, 6, 5},
121 : {2, 3, 7, 6},
122 : {3, 0, 4, 7},
123 : {4, 5, 6, 7}
124 : }};
125 : }
126 :
127 : LIBMESH_DEVICE_INLINE constexpr ReferenceElementTable<6, 3>
128 : tet_edge_nodes()
129 : {
130 : return {{
131 : {0, 1, 4},
132 : {1, 2, 5},
133 : {0, 2, 6},
134 : {0, 3, 7},
135 : {1, 3, 8},
136 : {2, 3, 9}
137 : }};
138 : }
139 :
140 : LIBMESH_DEVICE_INLINE constexpr ReferenceElementTable<12, 3>
141 : hex_edge_nodes()
142 : {
143 : return {{
144 : {0, 1, 8},
145 : {1, 2, 9},
146 : {2, 3, 10},
147 : {0, 3, 11},
148 : {0, 4, 12},
149 : {1, 5, 13},
150 : {2, 6, 14},
151 : {3, 7, 15},
152 : {4, 5, 16},
153 : {5, 6, 17},
154 : {6, 7, 18},
155 : {4, 7, 19}
156 : }};
157 : }
158 :
159 : LIBMESH_DEVICE_INLINE constexpr ReferenceElementTable<9, 3>
160 : prism_edge_nodes()
161 : {
162 : return {{
163 : {0, 1, 6},
164 : {1, 2, 7},
165 : {0, 2, 8},
166 : {0, 3, 9},
167 : {1, 4, 10},
168 : {2, 5, 11},
169 : {3, 4, 12},
170 : {4, 5, 13},
171 : {3, 5, 14}
172 : }};
173 : }
174 :
175 : LIBMESH_DEVICE_INLINE constexpr ReferenceElementTable<8, 3>
176 : pyramid_edge_nodes()
177 : {
178 : return {{
179 : {0, 1, 5},
180 : {1, 2, 6},
181 : {2, 3, 7},
182 : {0, 3, 8},
183 : {0, 4, 9},
184 : {1, 4, 10},
185 : {2, 4, 11},
186 : {3, 4, 12}
187 : }};
188 : }
189 :
190 : LIBMESH_DEVICE_INLINE bool
191 : requires_side_specific_topology(ElemType parent)
192 : {
193 : switch (parent)
194 : {
195 : case PRISM6:
196 : case PRISM15:
197 : case PRISM18:
198 : case PRISM20:
199 : case PRISM21:
200 : case PYRAMID5:
201 : case PYRAMID13:
202 : case PYRAMID14:
203 : case PYRAMID18:
204 : return true;
205 : default:
206 : return false;
207 : }
208 : }
209 :
210 : LIBMESH_DEVICE_INLINE ElemType
211 : side_topology_or_invalid(ElemType parent,
212 : unsigned int side)
213 : {
214 : if (side > 4)
215 : return INVALID_ELEM;
216 :
217 : // Prism sides 0 and 4 are the triangles; pyramid side 4 is the quad base.
218 : // Every other supported element has the same topology on all its sides.
219 : const bool prism_tri = (side == 0 || side == 4);
220 : const bool pyramid_tri = (side != 4);
221 :
222 : switch (parent)
223 : {
224 : case EDGE2:
225 : case EDGE3:
226 : case EDGE4:
227 : return NODEELEM;
228 : case TRI3:
229 : case QUAD4:
230 : return EDGE2;
231 : case TRI6:
232 : case TRI7:
233 : case QUAD8:
234 : case QUAD9:
235 : return EDGE3;
236 : case TET4:
237 : return TRI3;
238 : case HEX8:
239 : return QUAD4;
240 : case TET10:
241 : return TRI6;
242 : case TET14:
243 : return TRI7;
244 : case HEX20:
245 : return QUAD8;
246 : case HEX27:
247 : return QUAD9;
248 : case PRISM6:
249 : return prism_tri ? TRI3 : QUAD4;
250 : case PRISM15:
251 : return prism_tri ? TRI6 : QUAD8;
252 : case PRISM18:
253 : return prism_tri ? TRI6 : QUAD9;
254 : case PRISM20:
255 : case PRISM21:
256 : return prism_tri ? TRI7 : QUAD9;
257 : case PYRAMID5:
258 : return pyramid_tri ? TRI3 : QUAD4;
259 : case PYRAMID13:
260 : return pyramid_tri ? TRI6 : QUAD8;
261 : case PYRAMID14:
262 : return pyramid_tri ? TRI6 : QUAD9;
263 : case PYRAMID18:
264 : return pyramid_tri ? TRI7 : QUAD9;
265 : default:
266 : return INVALID_ELEM;
267 : }
268 : }
269 :
270 : // Mid-side node counts are uniform per element type except for the
271 : // prisms and pyramids, whose triangular and quadrilateral faces differ.
272 : LIBMESH_DEVICE_INLINE constexpr unsigned int
273 : side_node_count_or_zero(ElemType parent,
274 : unsigned int side)
275 : {
276 : switch (parent)
277 : {
278 : case EDGE2:
279 : case EDGE3:
280 : case EDGE4:
281 : return side < 2 ? 1 : 0;
282 : case TRI3:
283 : case TRISHELL3:
284 : return side < 3 ? 2 : 0;
285 : case TRI6:
286 : case TRI7:
287 : return side < 3 ? 3 : 0;
288 : case QUAD4:
289 : case QUADSHELL4:
290 : return side < 4 ? 2 : 0;
291 : case QUAD8:
292 : case QUADSHELL8:
293 : case QUAD9:
294 : case QUADSHELL9:
295 : return side < 4 ? 3 : 0;
296 : case TET4:
297 : return side < 4 ? 3 : 0;
298 : case TET10:
299 : return side < 4 ? 6 : 0;
300 : case TET14:
301 : return side < 4 ? 7 : 0;
302 : case HEX8:
303 : return side < 6 ? 4 : 0;
304 : case HEX20:
305 : return side < 6 ? 8 : 0;
306 : case HEX27:
307 : return side < 6 ? 9 : 0;
308 : case PRISM6:
309 : return side < 5 ? ((side == 0 || side == 4) ? 3 : 4) : 0;
310 : case PRISM15:
311 : return side < 5 ? ((side == 0 || side == 4) ? 6 : 8) : 0;
312 : case PRISM18:
313 : return side < 5 ? ((side == 0 || side == 4) ? 6 : 9) : 0;
314 : case PRISM20:
315 : case PRISM21:
316 : return side < 5 ? ((side == 0 || side == 4) ? 7 : 9) : 0;
317 : case PYRAMID5:
318 : return side < 5 ? (side == 4 ? 4 : 3) : 0;
319 : case PYRAMID13:
320 : return side < 5 ? (side == 4 ? 8 : 6) : 0;
321 : case PYRAMID14:
322 : return side < 5 ? (side == 4 ? 9 : 6) : 0;
323 : case PYRAMID18:
324 : return side < 5 ? (side == 4 ? 9 : 7) : 0;
325 : default:
326 : return 0;
327 : }
328 : }
329 :
330 : // Every element type with edge nodes has 3 nodes (two vertices plus a
331 : // midpoint) on each of its edges; only the edge count varies. The edge
332 : // tables cover the second-order families only: the linear elements'
333 : // vertex-pair edges stay with their classes, which device code never
334 : // queries.
335 : LIBMESH_DEVICE_INLINE constexpr unsigned int
336 : edge_node_count_or_zero(ElemType parent,
337 : unsigned int edge)
338 : {
339 : switch (parent)
340 : {
341 : case TET10:
342 : case TET14:
343 : return edge < 6 ? 3 : 0;
344 : case HEX20:
345 : case HEX27:
346 : return edge < 12 ? 3 : 0;
347 : case PRISM15:
348 : case PRISM18:
349 : case PRISM20:
350 : case PRISM21:
351 : return edge < 9 ? 3 : 0;
352 : case PYRAMID13:
353 : case PYRAMID14:
354 : case PYRAMID18:
355 : return edge < 8 ? 3 : 0;
356 : default:
357 : return 0;
358 : }
359 : }
360 :
361 : LIBMESH_DEVICE_INLINE constexpr bool
362 : try_local_edge_node(ElemType parent,
363 : unsigned int edge,
364 : unsigned int edge_node,
365 : unsigned int & node)
366 : {
367 : const unsigned int count = edge_node_count_or_zero(parent, edge);
368 : if (!count || edge_node >= count)
369 : return false;
370 :
371 : switch (parent)
372 : {
373 : case TET10:
374 : case TET14:
375 : node = tet_edge_nodes()(edge, edge_node);
376 : return true;
377 : case HEX20:
378 : case HEX27:
379 : node = hex_edge_nodes()(edge, edge_node);
380 : return true;
381 : case PRISM15:
382 : case PRISM18:
383 : case PRISM20:
384 : case PRISM21:
385 : node = prism_edge_nodes()(edge, edge_node);
386 : return true;
387 : case PYRAMID13:
388 : case PYRAMID14:
389 : case PYRAMID18:
390 : node = pyramid_edge_nodes()(edge, edge_node);
391 : return true;
392 : default:
393 : return false;
394 : }
395 : }
396 :
397 : LIBMESH_DEVICE_INLINE constexpr ElemType
398 : linear_sibling_or_invalid(ElemType type)
399 : {
400 : switch (type)
401 : {
402 : case EDGE2:
403 : case EDGE3:
404 : case EDGE4:
405 : return EDGE2;
406 : case TRI3:
407 : case TRISHELL3:
408 : case TRI6:
409 : case TRI7:
410 : return TRI3;
411 : case QUAD4:
412 : case QUADSHELL4:
413 : case QUAD8:
414 : case QUADSHELL8:
415 : case QUAD9:
416 : case QUADSHELL9:
417 : return QUAD4;
418 : case TET4:
419 : case TET10:
420 : case TET14:
421 : return TET4;
422 : case HEX8:
423 : case HEX20:
424 : case HEX27:
425 : return HEX8;
426 : case PRISM6:
427 : case PRISM15:
428 : case PRISM18:
429 : case PRISM20:
430 : case PRISM21:
431 : return PRISM6;
432 : case PYRAMID5:
433 : case PYRAMID13:
434 : case PYRAMID14:
435 : case PYRAMID18:
436 : return PYRAMID5;
437 : default:
438 : return INVALID_ELEM;
439 : }
440 : }
441 :
442 : // The node id of a face's center node, for the element types that have them.
443 : LIBMESH_DEVICE_INLINE constexpr unsigned int
444 : face_center_node_or_invalid(ElemType type,
445 : unsigned int side)
446 : {
447 : switch (type)
448 : {
449 : case TET14:
450 : return side < 4 ? 10 + side : invalid_uint;
451 : case HEX27:
452 : return side < 6 ? 20 + side : invalid_uint;
453 : case PRISM18:
454 : return (side >= 1 && side <= 3) ? 14 + side : invalid_uint;
455 : case PRISM20:
456 : case PRISM21:
457 : return side == 0 ? 18 :
458 : side == 4 ? 19 :
459 : side <= 3 ? 14 + side : invalid_uint;
460 : case PYRAMID14:
461 : return side == 4 ? 13 : invalid_uint;
462 : case PYRAMID18:
463 : return side == 4 ? 13 :
464 : side < 4 ? 14 + side : invalid_uint;
465 : default:
466 : return invalid_uint;
467 : }
468 : }
469 :
470 : // Corner k of a linear element's side, from the stored linear tables.
471 : LIBMESH_DEVICE_INLINE constexpr bool
472 : try_linear_corner(ElemType linear,
473 : unsigned int side,
474 : unsigned int k,
475 : unsigned int & node)
476 : {
477 : switch (linear)
478 : {
479 : case EDGE2:
480 : node = side;
481 : return true;
482 : case TRI3:
483 : node = tri3_side_nodes()(side, k);
484 : return true;
485 : case QUAD4:
486 : node = quad4_side_nodes()(side, k);
487 : return true;
488 : case TET4:
489 : node = tet4_side_nodes()(side, k);
490 : return true;
491 : case HEX8:
492 : node = hex8_side_nodes()(side, k);
493 : return true;
494 : case PRISM6:
495 : node = prism6_side_nodes()(side, k);
496 : return true;
497 : case PYRAMID5:
498 : node = pyramid5_side_nodes()(side, k);
499 : return true;
500 : default:
501 : return false;
502 : }
503 : }
504 :
505 : // The mid-edge node of the edge joining vertices a and b, if any.
506 : LIBMESH_DEVICE_INLINE constexpr bool
507 : try_edge_mid_between(ElemType parent,
508 : unsigned int a,
509 : unsigned int b,
510 : unsigned int & node)
511 : {
512 : for (unsigned int e = 0; edge_node_count_or_zero(parent, e); ++e)
513 : {
514 : unsigned int v0 = 0, v1 = 0;
515 : if (try_local_edge_node(parent, e, 0, v0) &&
516 : try_local_edge_node(parent, e, 1, v1) &&
517 : ((v0 == a && v1 == b) || (v0 == b && v1 == a)))
518 : return try_local_edge_node(parent, e, 2, node);
519 : }
520 : return false;
521 : }
522 :
523 : // A second-order side row is fully determined by the linear sibling's
524 : // corner row plus the element's edge table: the corners come first, then
525 : // the midpoint of each consecutive corner pair (wrapping), then the face's
526 : // center node where one exists. Only the linear tables are stored; the
527 : // second-order rows are derived, so the side and edge topologies cannot
528 : // disagree.
529 : LIBMESH_DEVICE_INLINE constexpr bool
530 : derive_local_side_node(ElemType parent,
531 : unsigned int side,
532 : unsigned int side_node,
533 : unsigned int & node)
534 : {
535 : const unsigned int count = side_node_count_or_zero(parent, side);
536 : if (!count || side_node >= count)
537 : return false;
538 :
539 : const ElemType linear = linear_sibling_or_invalid(parent);
540 : const unsigned int corners = side_node_count_or_zero(linear, side);
541 :
542 : if (side_node < corners)
543 : return try_linear_corner(linear, side, side_node, node);
544 :
545 : // 2D second-order elements: the single mid-side node
546 : if (count == 3 && corners == 2)
547 : {
548 : node = (linear == TRI3 ? 3u : 4u) + side;
549 : return true;
550 : }
551 :
552 : // 3D mid-edge nodes
553 : if (side_node < 2 * corners)
554 : {
555 : unsigned int a = 0, b = 0;
556 : if (!try_linear_corner(linear, side, side_node - corners, a) ||
557 : !try_linear_corner(linear, side, (side_node - corners + 1) % corners, b))
558 : return false;
559 : return try_edge_mid_between(parent, a, b, node);
560 : }
561 :
562 : // face center
563 : node = face_center_node_or_invalid(parent, side);
564 : return node != invalid_uint;
565 : }
566 :
567 : // Materialize a type's full side-node table at compile time from the
568 : // derivation above, so runtime lookups are direct indexing while the
569 : // derivation remains the only authority (99 pads short rows, matching the
570 : // stored linear tables).
571 : template <unsigned int Rows, unsigned int Cols>
572 : LIBMESH_DEVICE_INLINE constexpr ReferenceElementTable<Rows, Cols>
573 : build_side_nodes(ElemType parent)
574 : {
575 : ReferenceElementTable<Rows, Cols> t {};
576 : for (unsigned int r = 0; r != Rows; ++r)
577 : for (unsigned int c = 0; c != Cols; ++c)
578 : {
579 : unsigned int n = 99;
580 : derive_local_side_node(parent, r, c, n);
581 : t.values[r][c] = n;
582 : }
583 : return t;
584 : }
585 :
586 : LIBMESH_DEVICE_INLINE bool
587 : try_local_side_node(ElemType parent,
588 : unsigned int side,
589 : unsigned int side_node,
590 : unsigned int & node)
591 : {
592 : const unsigned int count = side_node_count_or_zero(parent, side);
593 : if (!count || side_node >= count)
594 : return false;
595 :
596 : switch (parent)
597 : {
598 : case EDGE2:
599 : case EDGE3:
600 : case EDGE4:
601 : node = side;
602 : return true;
603 : case TRI3:
604 : case TRISHELL3:
605 : node = tri3_side_nodes()(side, side_node);
606 : return true;
607 : case QUAD4:
608 : case QUADSHELL4:
609 : node = quad4_side_nodes()(side, side_node);
610 : return true;
611 : case TET4:
612 : node = tet4_side_nodes()(side, side_node);
613 : return true;
614 : case HEX8:
615 : node = hex8_side_nodes()(side, side_node);
616 : return true;
617 : case PRISM6:
618 : node = prism6_side_nodes()(side, side_node);
619 : return true;
620 : case PYRAMID5:
621 : node = pyramid5_side_nodes()(side, side_node);
622 : return true;
623 : case TRI6:
624 : case TRI7:
625 : {
626 : constexpr auto t = build_side_nodes<3, 3>(TRI6);
627 : node = t(side, side_node);
628 : return true;
629 : }
630 : case QUAD8:
631 : case QUADSHELL8:
632 : case QUAD9:
633 : case QUADSHELL9:
634 : {
635 : constexpr auto t = build_side_nodes<4, 3>(QUAD8);
636 : node = t(side, side_node);
637 : return true;
638 : }
639 : case TET10:
640 : {
641 : constexpr auto t = build_side_nodes<4, 6>(TET10);
642 : node = t(side, side_node);
643 : return true;
644 : }
645 : case TET14:
646 : {
647 : constexpr auto t = build_side_nodes<4, 7>(TET14);
648 : node = t(side, side_node);
649 : return true;
650 : }
651 : case HEX20:
652 : {
653 : constexpr auto t = build_side_nodes<6, 8>(HEX20);
654 : node = t(side, side_node);
655 : return true;
656 : }
657 : case HEX27:
658 : {
659 : constexpr auto t = build_side_nodes<6, 9>(HEX27);
660 : node = t(side, side_node);
661 : return true;
662 : }
663 : case PRISM15:
664 : {
665 : constexpr auto t = build_side_nodes<5, 8>(PRISM15);
666 : node = t(side, side_node);
667 : return true;
668 : }
669 : case PRISM18:
670 : {
671 : constexpr auto t = build_side_nodes<5, 9>(PRISM18);
672 : node = t(side, side_node);
673 : return true;
674 : }
675 : case PRISM20:
676 : case PRISM21:
677 : {
678 : constexpr auto t = build_side_nodes<5, 9>(PRISM20);
679 : node = t(side, side_node);
680 : return true;
681 : }
682 : case PYRAMID13:
683 : {
684 : constexpr auto t = build_side_nodes<5, 8>(PYRAMID13);
685 : node = t(side, side_node);
686 : return true;
687 : }
688 : case PYRAMID14:
689 : {
690 : constexpr auto t = build_side_nodes<5, 9>(PYRAMID14);
691 : node = t(side, side_node);
692 : return true;
693 : }
694 : case PYRAMID18:
695 : {
696 : constexpr auto t = build_side_nodes<5, 9>(PYRAMID18);
697 : node = t(side, side_node);
698 : return true;
699 : }
700 : default:
701 : return false;
702 : }
703 : }
704 :
705 : // Reference-element vertex locations, per family (higher-order family
706 : // members share their linear sibling's vertices).
707 : LIBMESH_DEVICE_INLINE unsigned int
708 : reference_vertex_count(ElemType type)
709 : {
710 : switch (type)
711 : {
712 : case EDGE2:
713 : case EDGE3:
714 : case EDGE4:
715 : return 2;
716 : case TRI3:
717 : case TRISHELL3:
718 : case TRI6:
719 : case TRI7:
720 : return 3;
721 : case QUAD4:
722 : case QUADSHELL4:
723 : case QUAD8:
724 : case QUADSHELL8:
725 : case QUAD9:
726 : case QUADSHELL9:
727 : case TET4:
728 : case TET10:
729 : case TET14:
730 : return 4;
731 : case PYRAMID5:
732 : case PYRAMID13:
733 : case PYRAMID14:
734 : case PYRAMID18:
735 : return 5;
736 : case PRISM6:
737 : case PRISM15:
738 : case PRISM18:
739 : case PRISM20:
740 : case PRISM21:
741 : return 6;
742 : case HEX8:
743 : case HEX20:
744 : case HEX27:
745 : return 8;
746 : default:
747 : return 0;
748 : }
749 : }
750 :
751 : LIBMESH_DEVICE_INLINE bool
752 : reference_vertex(ElemType type,
753 : unsigned int v,
754 : Point & pt)
755 : {
756 : if (v >= reference_vertex_count(type))
757 : return false;
758 :
759 : switch (type)
760 : {
761 : case EDGE2:
762 : case EDGE3:
763 : case EDGE4:
764 : pt = Point(v == 0 ? -1.0 : 1.0);
765 : return true;
766 : case TRI3:
767 : case TRISHELL3:
768 : case TRI6:
769 : case TRI7:
770 : pt = Point(v == 1 ? 1.0 : 0.0, v == 2 ? 1.0 : 0.0);
771 : return true;
772 : case QUAD4:
773 : case QUADSHELL4:
774 : case QUAD8:
775 : case QUADSHELL8:
776 : case QUAD9:
777 : case QUADSHELL9:
778 : pt = Point((v == 1 || v == 2) ? 1.0 : -1.0,
779 : (v == 2 || v == 3) ? 1.0 : -1.0);
780 : return true;
781 : case TET4:
782 : case TET10:
783 : case TET14:
784 : pt = Point(v == 1 ? 1.0 : 0.0, v == 2 ? 1.0 : 0.0, v == 3 ? 1.0 : 0.0);
785 : return true;
786 : case PYRAMID5:
787 : case PYRAMID13:
788 : case PYRAMID14:
789 : case PYRAMID18:
790 : pt = v == 4 ? Point(0.0, 0.0, 1.0)
791 : : Point((v == 1 || v == 2) ? 1.0 : -1.0,
792 : (v == 2 || v == 3) ? 1.0 : -1.0,
793 : 0.0);
794 : return true;
795 : case PRISM6:
796 : case PRISM15:
797 : case PRISM18:
798 : case PRISM20:
799 : case PRISM21:
800 : pt = Point(v % 3 == 1 ? 1.0 : 0.0,
801 : v % 3 == 2 ? 1.0 : 0.0,
802 : v < 3 ? -1.0 : 1.0);
803 : return true;
804 : case HEX8:
805 : case HEX20:
806 : case HEX27:
807 : pt = Point((v % 4 == 1 || v % 4 == 2) ? 1.0 : -1.0,
808 : (v % 4 == 2 || v % 4 == 3) ? 1.0 : -1.0,
809 : v < 4 ? -1.0 : 1.0);
810 : return true;
811 : default:
812 : return false;
813 : }
814 : }
815 :
816 : // The node holding an element's vertex-centroid, if it has one.
817 : LIBMESH_DEVICE_INLINE unsigned int
818 : centroid_node_or_invalid(ElemType type)
819 : {
820 : switch (type)
821 : {
822 : case EDGE3:
823 : return 2;
824 : case TRI7:
825 : return 6;
826 : case QUAD9:
827 : case QUADSHELL9:
828 : return 8;
829 : case HEX27:
830 : return 26;
831 : case PRISM21:
832 : return 20;
833 : default:
834 : return invalid_uint;
835 : }
836 : }
837 :
838 : // Every higher-order reference node sits at the centroid of its
839 : // subentity's vertices -- mid-edge nodes at edge midpoints, face nodes at
840 : // face-corner centroids, interior nodes at the vertex centroid -- so only
841 : // the vertices are tabulated and the rest is derived through the same
842 : // side/edge topology tables everything else uses. The one exception is
843 : // the cubic EDGE4, whose two interior nodes trisect the edge.
844 : LIBMESH_DEVICE_INLINE bool
845 : try_reference_node(ElemType type,
846 : unsigned int node,
847 : Point & pt)
848 : {
849 : const unsigned int nv = reference_vertex_count(type);
850 : if (node < nv)
851 : return reference_vertex(type, node, pt);
852 :
853 : if (type == EDGE4)
854 : {
855 : if (node > 3)
856 : return false;
857 : pt = Point(node == 2 ? Real(-1) / 3 : Real(1) / 3);
858 : return true;
859 : }
860 :
861 : if (node == centroid_node_or_invalid(type))
862 : {
863 : Point sum;
864 : for (unsigned int v = 0; v != nv; ++v)
865 : {
866 : Point pv;
867 : reference_vertex(type, v, pv);
868 : sum += pv;
869 : }
870 : pt = sum / Real(nv);
871 : return true;
872 : }
873 :
874 : // Mid-edge nodes: the third entry of an edge-table row {v0, v1, mid}.
875 : for (unsigned int e = 0; edge_node_count_or_zero(type, e); ++e)
876 : {
877 : unsigned int mid, v0, v1;
878 : if (try_local_edge_node(type, e, 2, mid) && mid == node &&
879 : try_local_edge_node(type, e, 0, v0) &&
880 : try_local_edge_node(type, e, 1, v1))
881 : {
882 : Point p0, p1;
883 : reference_vertex(type, v0, p0);
884 : reference_vertex(type, v1, p1);
885 : pt = (p0 + p1) / 2;
886 : return true;
887 : }
888 : }
889 :
890 : // 2D mid-side nodes ({v0, v1, mid} side rows) and 3D face-center nodes
891 : // (the last entry of a 7- or 9-node side row, behind 3 or 4 corners).
892 : for (unsigned int s = 0; ; ++s)
893 : {
894 : const unsigned int count = side_node_count_or_zero(type, s);
895 : if (!count)
896 : break;
897 : const unsigned int corners =
898 : count == 3 ? 2 : count == 7 ? 3 : count == 9 ? 4 : 0;
899 : unsigned int last;
900 : if (!corners ||
901 : !try_local_side_node(type, s, count - 1, last) || last != node)
902 : continue;
903 : Point sum;
904 : for (unsigned int k = 0; k != corners; ++k)
905 : {
906 : unsigned int v;
907 : try_local_side_node(type, s, k, v);
908 : Point pv;
909 : reference_vertex(type, v, pv);
910 : sum += pv;
911 : }
912 : pt = sum / Real(corners);
913 : return true;
914 : }
915 :
916 : return false;
917 : }
918 : LIBMESH_DEVICE_INLINE bool
919 : try_refspace_node(ElemType type,
920 : unsigned int node,
921 : Point & pt)
922 : {
923 : switch (type)
924 : {
925 : case NODEELEM:
926 : if (!node)
927 : {
928 : pt = Point(0.0, 0.0, 0.0);
929 : return true;
930 : }
931 : return false;
932 :
933 : case TRISHELL3:
934 : return try_reference_node(TRI3, node, pt);
935 :
936 : case QUADSHELL4:
937 : return try_reference_node(QUAD4, node, pt);
938 :
939 : case QUADSHELL8:
940 : return try_reference_node(QUAD8, node, pt);
941 :
942 : case QUADSHELL9:
943 : return try_reference_node(QUAD9, node, pt);
944 :
945 : default:
946 : return try_reference_node(type, node, pt);
947 : }
948 : }
949 :
950 : LIBMESH_DEVICE_INLINE bool
951 : try_reference_side_node(ElemType parent,
952 : unsigned int side,
953 : unsigned int side_node,
954 : Point & pt)
955 : {
956 : unsigned int node = libMesh::invalid_uint;
957 : if (!try_local_side_node(parent, side, side_node, node))
958 : return false;
959 :
960 : return try_reference_node(parent, node, pt);
961 : }
962 :
963 : } // namespace libMesh
964 :
965 : #endif // LIBMESH_FE_REFERENCE_ELEMENT_TRAITS_H
|
