libMesh
Loading...
Searching...
No Matches
fe_hierarchic_shape_3D.C
Go to the documentation of this file.
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/fe.h"
21#include "libmesh/elem.h"
22#include "libmesh/number_lookups.h"
23#include "libmesh/enum_to_string.h"
24#include "libmesh/cell_tet4.h" // We need edge_nodes_map + side_nodes_map
25#include "libmesh/cell_prism6.h"
26#include "libmesh/face_tri3.h" // Faster to construct these on the stack
27#include "libmesh/face_quad4.h"
28
29// Anonymous namespace for functions shared by HIERARCHIC and
30// L2_HIERARCHIC implementations. Implementations appear at the bottom
31// of this file.
32namespace
33{
34using namespace libMesh;
35
36unsigned int cube_side(const Point & p);
37
38Point cube_side_point(unsigned int sidenum, const Point & interior_point);
39
40std::array<unsigned int, 4> oriented_prism_nodes(const Elem & elem,
41 unsigned int face_num);
42
43std::array<unsigned int, 3> oriented_tet_nodes(const Elem & elem,
44 unsigned int face_num);
45
46template <FEFamily T>
47Real fe_hierarchic_3D_shape(const Elem * elem,
48 const Order order,
49 const unsigned int i,
50 const Point & p,
51 const bool add_p_level);
52
53template <FEFamily T>
54Real fe_hierarchic_3D_shape_deriv(const Elem * elem,
55 const Order order,
56 const unsigned int i,
57 const unsigned int j,
58 const Point & p,
59 const bool add_p_level);
60
61#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
62
63template <FEFamily T>
64Real fe_hierarchic_3D_shape_second_deriv(const Elem * elem,
65 const Order order,
66 const unsigned int i,
67 const unsigned int j,
68 const Point & p,
69 const bool add_p_level);
70
71#endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
72
73#if LIBMESH_DIM > 2
74Point get_min_point(const Elem * elem,
75 unsigned int a,
76 unsigned int b,
77 unsigned int c,
78 unsigned int d)
79{
80 return std::min(std::min(elem->point(a),elem->point(b)),
81 std::min(elem->point(c),elem->point(d)));
82}
83
84// Remap non-face-nodes based on point ordering
85template <unsigned int N_nodes>
86unsigned int remap_node(unsigned int n,
87 const Elem & elem,
88 unsigned int nodebegin)
89{
90 std::array<const Point *, N_nodes> points;
91
92 for (auto i : IntRange<unsigned int>(0, N_nodes))
93 points[i] = &elem.point(nodebegin+i);
94
95 std::sort(points.begin(), points.end(),
96 [](const Point * a, const Point * b)
97 { return *a < *b; });
98
99 const Point * pn = points[n-nodebegin];
100
101 for (auto i : IntRange<unsigned int>(nodebegin, nodebegin+N_nodes))
102 if (pn == &elem.point(i))
103 return i;
104
105 libmesh_assert(false);
106 return libMesh::invalid_uint;
107}
108
109
110void cube_remap(unsigned int & side_i,
111 const Elem & side,
112 unsigned int totalorder,
113 Point & sidep)
114{
115 // "vertex" nodes are now decoupled from vertices, so we have
116 // to order them consistently otherwise
117 if (side_i < 4)
118 side_i = remap_node<4>(side_i, side, 0);
119
120 // And "edge" nodes are decoupled from edges, so we have to
121 // reorder them too!
122 else if (side_i < 4u*totalorder)
123 {
124 unsigned int side_node = (side_i - 4)/(totalorder-1)+4;
125 side_node = remap_node<4>(side_node, side, 4);
126 side_i = ((side_i - 4) % (totalorder - 1)) // old local edge_i
127 + 4 + (side_node-4)*(totalorder-1);
128 }
129
130 // Interior dofs in 2D don't care about where xi/eta point in
131 // physical space, but here we need them to match from both
132 // sides of a face!
133 else
134 {
135 unsigned int min_side_node = remap_node<4>(0, side, 0);
136 const bool flip = (side.point(min_side_node) < side.point((min_side_node+1)%4));
137
138 switch (min_side_node) {
139 case 0:
140 if (flip)
141 std::swap(sidep(0), sidep(1));
142 break;
143 case 1:
144 sidep(0) = -sidep(0);
145 if (!flip)
146 std::swap(sidep(0), sidep(1));
147 break;
148 case 2:
149 sidep(0) = -sidep(0);
150 sidep(1) = -sidep(1);
151 if (flip)
152 std::swap(sidep(0), sidep(1));
153 break;
154 case 3:
155 sidep(1) = -sidep(1);
156 if (!flip)
157 std::swap(sidep(0), sidep(1));
158 break;
159 default:
160 libmesh_error();
161 }
162 }
163}
164
165
166void cube_indices(const Elem * elem,
167 const unsigned int totalorder,
168 const unsigned int i,
169 Real & xi, Real & eta, Real & zeta,
170 unsigned int & i0,
171 unsigned int & i1,
172 unsigned int & i2)
173{
174 // The only way to make any sense of this
175 // is to look at the mgflo/mg2/mgf documentation
176 // and make the cut-out cube!
177 // Example i0 and i1 values for totalorder = 3:
178 // FIXME - these examples are incorrect now that we've got truly
179 // hierarchic basis functions
180 // Nodes 0 1 2 3 4 5 6 7 8 8 9 9 10 10 11 11 12 12 13 13 14 14 15 15 16 16 17 17 18 18 19 19 20 20 20 20 21 21 21 21 22 22 22 22 23 23 23 23 24 24 24 24 25 25 25 25 26 26 26 26 26 26 26 26
181 // DOFS 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 18 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 60 62 63
182 // static const unsigned int i0[] = {0, 1, 1, 0, 0, 1, 1, 0, 2, 3, 1, 1, 2, 3, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 2, 3, 1, 1, 2, 3, 0, 0, 2, 3, 2, 3, 2, 3, 2, 3, 1, 1, 1, 1, 2, 3, 2, 3, 0, 0, 0, 0, 2, 3, 2, 3, 2, 3, 2, 3, 2, 3, 2, 3};
183 // static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 2, 3, 1, 1, 2, 3, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 2, 3, 1, 1, 2, 3, 2, 2, 3, 3, 0, 0, 0, 0, 2, 3, 2, 3, 1, 1, 1, 1, 2, 3, 2, 3, 2, 2, 3, 3, 2, 2, 3, 3, 2, 2, 3, 3};
184 // static const unsigned int i2[] = {0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 2, 3, 2, 3, 2, 3, 2, 3, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 2, 2, 3, 3, 2, 2, 3, 3, 2, 2, 3, 3, 2, 2, 3, 3, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3};
185
186 // the number of DoFs per edge appears everywhere:
187 const unsigned int e = totalorder - 1u;
188
189 libmesh_assert_less (i, (totalorder+1u)*(totalorder+1u)*(totalorder+1u));
190
191 Real xi_saved = xi, eta_saved = eta, zeta_saved = zeta;
192
193 // Vertices:
194 if (i == 0)
195 {
196 i0 = 0;
197 i1 = 0;
198 i2 = 0;
199 }
200 else if (i == 1)
201 {
202 i0 = 1;
203 i1 = 0;
204 i2 = 0;
205 }
206 else if (i == 2)
207 {
208 i0 = 1;
209 i1 = 1;
210 i2 = 0;
211 }
212 else if (i == 3)
213 {
214 i0 = 0;
215 i1 = 1;
216 i2 = 0;
217 }
218 else if (i == 4)
219 {
220 i0 = 0;
221 i1 = 0;
222 i2 = 1;
223 }
224 else if (i == 5)
225 {
226 i0 = 1;
227 i1 = 0;
228 i2 = 1;
229 }
230 else if (i == 6)
231 {
232 i0 = 1;
233 i1 = 1;
234 i2 = 1;
235 }
236 else if (i == 7)
237 {
238 i0 = 0;
239 i1 = 1;
240 i2 = 1;
241 }
242 // Edge 0
243 else if (i < 8 + e)
244 {
245 i0 = i - 6;
246 i1 = 0;
247 i2 = 0;
248 if (elem->positive_edge_orientation(0))
249 xi = -xi_saved;
250 }
251 // Edge 1
252 else if (i < 8 + 2*e)
253 {
254 i0 = 1;
255 i1 = i - e - 6;
256 i2 = 0;
257 if (elem->positive_edge_orientation(1))
258 eta = -eta_saved;
259 }
260 // Edge 2
261 else if (i < 8 + 3*e)
262 {
263 i0 = i - 2*e - 6;
264 i1 = 1;
265 i2 = 0;
266 if (!elem->positive_edge_orientation(2))
267 xi = -xi_saved;
268 }
269 // Edge 3
270 else if (i < 8 + 4*e)
271 {
272 i0 = 0;
273 i1 = i - 3*e - 6;
274 i2 = 0;
275 if (elem->positive_edge_orientation(3))
276 eta = -eta_saved;
277 }
278 // Edge 4
279 else if (i < 8 + 5*e)
280 {
281 i0 = 0;
282 i1 = 0;
283 i2 = i - 4*e - 6;
284 if (elem->positive_edge_orientation(4))
285 zeta = -zeta_saved;
286 }
287 // Edge 5
288 else if (i < 8 + 6*e)
289 {
290 i0 = 1;
291 i1 = 0;
292 i2 = i - 5*e - 6;
293 if (elem->positive_edge_orientation(5))
294 zeta = -zeta_saved;
295 }
296 // Edge 6
297 else if (i < 8 + 7*e)
298 {
299 i0 = 1;
300 i1 = 1;
301 i2 = i - 6*e - 6;
302 if (elem->positive_edge_orientation(6))
303 zeta = -zeta_saved;
304 }
305 // Edge 7
306 else if (i < 8 + 8*e)
307 {
308 i0 = 0;
309 i1 = 1;
310 i2 = i - 7*e - 6;
311 if (elem->positive_edge_orientation(7))
312 zeta = -zeta_saved;
313 }
314 // Edge 8
315 else if (i < 8 + 9*e)
316 {
317 i0 = i - 8*e - 6;
318 i1 = 0;
319 i2 = 1;
320 if (elem->positive_edge_orientation(8))
321 xi = -xi_saved;
322 }
323 // Edge 9
324 else if (i < 8 + 10*e)
325 {
326 i0 = 1;
327 i1 = i - 9*e - 6;
328 i2 = 1;
329 if (elem->positive_edge_orientation(9))
330 eta = -eta_saved;
331 }
332 // Edge 10
333 else if (i < 8 + 11*e)
334 {
335 i0 = i - 10*e - 6;
336 i1 = 1;
337 i2 = 1;
338 if (!elem->positive_edge_orientation(10))
339 xi = -xi_saved;
340 }
341 // Edge 11
342 else if (i < 8 + 12*e)
343 {
344 i0 = 0;
345 i1 = i - 11*e - 6;
346 i2 = 1;
347 if (elem->positive_edge_orientation(11))
348 eta = -eta_saved;
349 }
350 // Face 0
351 else if (i < 8 + 12*e + e*e)
352 {
353 unsigned int basisnum = i - 8 - 12*e;
354 i0 = square_number_row[basisnum] + 2;
355 i1 = square_number_column[basisnum] + 2;
356 i2 = 0;
357 const Point min_point = get_min_point(elem, 1, 2, 0, 3);
358
359 if (elem->point(0) == min_point)
360 if (elem->positive_face_orientation(0))
361 {
362 // Case 1
363 xi = xi_saved;
364 eta = eta_saved;
365 }
366 else
367 {
368 // Case 2
369 xi = eta_saved;
370 eta = xi_saved;
371 }
372
373 else if (elem->point(3) == min_point)
374 if (elem->positive_face_orientation(0))
375 {
376 // Case 3
377 xi = -eta_saved;
378 eta = xi_saved;
379 }
380 else
381 {
382 // Case 4
383 xi = xi_saved;
384 eta = -eta_saved;
385 }
386
387 else if (elem->point(2) == min_point)
388 if (elem->positive_face_orientation(0))
389 {
390 // Case 5
391 xi = -xi_saved;
392 eta = -eta_saved;
393 }
394 else
395 {
396 // Case 6
397 xi = -eta_saved;
398 eta = -xi_saved;
399 }
400
401 else if (elem->point(1) == min_point)
402 {
403 if (elem->positive_face_orientation(0))
404 {
405 // Case 7
406 xi = eta_saved;
407 eta = -xi_saved;
408 }
409 else
410 {
411 // Case 8
412 xi = -xi_saved;
413 eta = eta_saved;
414 }
415 }
416 }
417 // Face 1
418 else if (i < 8 + 12*e + 2*e*e)
419 {
420 unsigned int basisnum = i - 8 - 12*e - e*e;
421 i0 = square_number_row[basisnum] + 2;
422 i1 = 0;
423 i2 = square_number_column[basisnum] + 2;
424 const Point min_point = get_min_point(elem, 0, 1, 5, 4);
425
426 if (elem->point(0) == min_point)
427 if (!elem->positive_face_orientation(1))
428 {
429 // Case 1
430 xi = xi_saved;
431 zeta = zeta_saved;
432 }
433 else
434 {
435 // Case 2
436 xi = zeta_saved;
437 zeta = xi_saved;
438 }
439
440 else if (elem->point(1) == min_point)
441 if (!elem->positive_face_orientation(1))
442 {
443 // Case 3
444 xi = zeta_saved;
445 zeta = -xi_saved;
446 }
447 else
448 {
449 // Case 4
450 xi = -xi_saved;
451 zeta = zeta_saved;
452 }
453
454 else if (elem->point(5) == min_point)
455 if (!elem->positive_face_orientation(1))
456 {
457 // Case 5
458 xi = -xi_saved;
459 zeta = -zeta_saved;
460 }
461 else
462 {
463 // Case 6
464 xi = -zeta_saved;
465 zeta = -xi_saved;
466 }
467
468 else if (elem->point(4) == min_point)
469 {
470 if (!elem->positive_face_orientation(1))
471 {
472 // Case 7
473 xi = -xi_saved;
474 zeta = zeta_saved;
475 }
476 else
477 {
478 // Case 8
479 xi = xi_saved;
480 zeta = -zeta_saved;
481 }
482 }
483 }
484 // Face 2
485 else if (i < 8 + 12*e + 3*e*e)
486 {
487 unsigned int basisnum = i - 8 - 12*e - 2*e*e;
488 i0 = 1;
489 i1 = square_number_row[basisnum] + 2;
490 i2 = square_number_column[basisnum] + 2;
491 const Point min_point = get_min_point(elem, 1, 2, 6, 5);
492
493 if (elem->point(1) == min_point)
494 if (!elem->positive_face_orientation(2))
495 {
496 // Case 1
497 eta = eta_saved;
498 zeta = zeta_saved;
499 }
500 else
501 {
502 // Case 2
503 eta = zeta_saved;
504 zeta = eta_saved;
505 }
506
507 else if (elem->point(2) == min_point)
508 if (!elem->positive_face_orientation(2))
509 {
510 // Case 3
511 eta = zeta_saved;
512 zeta = -eta_saved;
513 }
514 else
515 {
516 // Case 4
517 eta = -eta_saved;
518 zeta = zeta_saved;
519 }
520
521 else if (elem->point(6) == min_point)
522 if (!elem->positive_face_orientation(2))
523 {
524 // Case 5
525 eta = -eta_saved;
526 zeta = -zeta_saved;
527 }
528 else
529 {
530 // Case 6
531 eta = -zeta_saved;
532 zeta = -eta_saved;
533 }
534
535 else if (elem->point(5) == min_point)
536 {
537 if (!elem->positive_face_orientation(2))
538 {
539 // Case 7
540 eta = -zeta_saved;
541 zeta = eta_saved;
542 }
543 else
544 {
545 // Case 8
546 eta = eta_saved;
547 zeta = -zeta_saved;
548 }
549 }
550 }
551 // Face 3
552 else if (i < 8 + 12*e + 4*e*e)
553 {
554 unsigned int basisnum = i - 8 - 12*e - 3*e*e;
555 i0 = square_number_row[basisnum] + 2;
556 i1 = 1;
557 i2 = square_number_column[basisnum] + 2;
558 const Point min_point = get_min_point(elem, 2, 3, 7, 6);
559
560 if (elem->point(3) == min_point)
561 if (elem->positive_face_orientation(3))
562 {
563 // Case 1
564 xi = xi_saved;
565 zeta = zeta_saved;
566 }
567 else
568 {
569 // Case 2
570 xi = zeta_saved;
571 zeta = xi_saved;
572 }
573
574 else if (elem->point(7) == min_point)
575 if (elem->positive_face_orientation(3))
576 {
577 // Case 3
578 xi = -zeta_saved;
579 zeta = xi_saved;
580 }
581 else
582 {
583 // Case 4
584 xi = xi_saved;
585 zeta = -zeta_saved;
586 }
587
588 else if (elem->point(6) == min_point)
589 if (elem->positive_face_orientation(3))
590 {
591 // Case 5
592 xi = -xi_saved;
593 zeta = -zeta_saved;
594 }
595 else
596 {
597 // Case 6
598 xi = -zeta_saved;
599 zeta = -xi_saved;
600 }
601
602 else if (elem->point(2) == min_point)
603 {
604 if (elem->positive_face_orientation(3))
605 {
606 // Case 7
607 xi = zeta_saved;
608 zeta = -xi_saved;
609 }
610 else
611 {
612 // Case 8
613 xi = -xi_saved;
614 zeta = zeta_saved;
615 }
616 }
617 }
618 // Face 4
619 else if (i < 8 + 12*e + 5*e*e)
620 {
621 unsigned int basisnum = i - 8 - 12*e - 4*e*e;
622 i0 = 0;
623 i1 = square_number_row[basisnum] + 2;
624 i2 = square_number_column[basisnum] + 2;
625 const Point min_point = get_min_point(elem, 3, 0, 4, 7);
626
627 if (elem->point(0) == min_point)
628 if (elem->positive_face_orientation(4))
629 {
630 // Case 1
631 eta = eta_saved;
632 zeta = zeta_saved;
633 }
634 else
635 {
636 // Case 2
637 eta = zeta_saved;
638 zeta = eta_saved;
639 }
640
641 else if (elem->point(4) == min_point)
642 if (elem->positive_face_orientation(4))
643 {
644 // Case 3
645 eta = -zeta_saved;
646 zeta = eta_saved;
647 }
648 else
649 {
650 // Case 4
651 eta = eta_saved;
652 zeta = -zeta_saved;
653 }
654
655 else if (elem->point(7) == min_point)
656 if (elem->positive_face_orientation(4))
657 {
658 // Case 5
659 eta = -eta_saved;
660 zeta = -zeta_saved;
661 }
662 else
663 {
664 // Case 6
665 eta = -zeta_saved;
666 zeta = -eta_saved;
667 }
668
669 else if (elem->point(3) == min_point)
670 {
671 if (elem->positive_face_orientation(4))
672 {
673 // Case 7
674 eta = zeta_saved;
675 zeta = -eta_saved;
676 }
677 else
678 {
679 // Case 8
680 eta = -eta_saved;
681 zeta = zeta_saved;
682 }
683 }
684 }
685 // Face 5
686 else if (i < 8 + 12*e + 6*e*e)
687 {
688 unsigned int basisnum = i - 8 - 12*e - 5*e*e;
689 i0 = square_number_row[basisnum] + 2;
690 i1 = square_number_column[basisnum] + 2;
691 i2 = 1;
692 const Point min_point = get_min_point(elem, 4, 5, 6, 7);
693
694 if (elem->point(4) == min_point)
695 if (!elem->positive_face_orientation(5))
696 {
697 // Case 1
698 xi = xi_saved;
699 eta = eta_saved;
700 }
701 else
702 {
703 // Case 2
704 xi = eta_saved;
705 eta = xi_saved;
706 }
707
708 else if (elem->point(5) == min_point)
709 if (!elem->positive_face_orientation(5))
710 {
711 // Case 3
712 xi = eta_saved;
713 eta = -xi_saved;
714 }
715 else
716 {
717 // Case 4
718 xi = -xi_saved;
719 eta = eta_saved;
720 }
721
722 else if (elem->point(6) == min_point)
723 if (!elem->positive_face_orientation(5))
724 {
725 // Case 5
726 xi = -xi_saved;
727 eta = -eta_saved;
728 }
729 else
730 {
731 // Case 6
732 xi = -eta_saved;
733 eta = -xi_saved;
734 }
735
736 else if (elem->point(7) == min_point)
737 {
738 if (!elem->positive_face_orientation(5))
739 {
740 // Case 7
741 xi = -eta_saved;
742 eta = xi_saved;
743 }
744 else
745 {
746 // Case 8
747 xi = xi_saved;
748 eta = eta_saved;
749 }
750 }
751 }
752
753 // Internal DoFs
754 else
755 {
756 unsigned int basisnum = i - 8 - 12*e - 6*e*e;
757 i0 = cube_number_column[basisnum] + 2;
758 i1 = cube_number_row[basisnum] + 2;
759 i2 = cube_number_page[basisnum] + 2;
760 }
761}
762
763
764void prism_indices(const Elem * elem,
765 const unsigned int totalorder,
766 const unsigned int i,
767 Point & xi_eta, Real & zeta,
768 unsigned int & i01,
769 unsigned int & i2)
770{
771 // the number of DoFs per edge appears everywhere:
772 const unsigned int e = totalorder - 1u;
773
774 libmesh_assert_less (i, (totalorder+1u)*(totalorder+1u)*(totalorder+2u)/2u);
775
776 Point xi_eta_saved = xi_eta;
777 Real zeta_saved = zeta;
778
779 // Vertices:
780 if (i == 0)
781 {
782 i01 = 0;
783 i2 = 0;
784 }
785 else if (i == 1)
786 {
787 i01 = 1;
788 i2 = 0;
789 }
790 else if (i == 2)
791 {
792 i01 = 2;
793 i2 = 0;
794 }
795 else if (i == 3)
796 {
797 i01 = 0;
798 i2 = 1;
799 }
800 else if (i == 4)
801 {
802 i01 = 1;
803 i2 = 1;
804 }
805 else if (i == 5)
806 {
807 i01 = 2;
808 i2 = 1;
809 }
810 // Edges 0,1,2 (vertices 6,7,8)
811 else if (i < 6 + 3*e)
812 {
813 // The TRI code will handle any flips here
814 i01 = i - 3;
815 i2 = 0;
816 }
817 // Edge 3,4,5 (vertices 9,10,11)
818 else if (i < 6 + 6*e)
819 {
820 i01 = (i - 6 - 3*e)/e; // which tri DoF are we?
821 i2 = (i - 6 - 3*e)%e+2; // edge DoF? +2 to skip endpoints
822 // EDGE evaluations don't flip, so handle that here
823 if (elem->positive_edge_orientation(i01+3))
824 zeta = -zeta;
825 }
826 // Edge 6,7,8 (vertices 12,13,14)
827 else if (i < 6 + 9*e)
828 {
829 // The TRI code will handle any flips here
830 i01 = i - 3 - 6*e;
831 i2 = 1;
832 }
833 // Face 1, node 15 (*before* 0, via node 18 on prism20)
834 else if (i < 6 + 9*e + e*e)
835 {
836 unsigned int basisnum = i - 6 - 9*e;
837
838 // How wide is the stretch from one side to the other of the
839 // line in the xi-eta plane parallel to this face?
840 const Real xe_scale = 1 - xi_eta_saved(1);
841
842 // What percentage of the way along that stretch are we?
843 const Real xe_fraction = (xe_scale==0) ?
844 0 : xi_eta_saved(0)/xe_scale;
845
846 // indexes in edge numbering
847 unsigned int s0 = square_number_row[basisnum] + 2;
848 unsigned int s1 = square_number_column[basisnum] + 2;
849 const Point min_point = get_min_point(elem, 0, 1, 3, 4);
850
851 if (elem->point(0) == min_point)
852 {
853 if (!elem->positive_face_orientation(1))
854 {
855 // Case 1: no flips needed
856 i01 = s0+1; // edge to triangle side 0 numbering
857 i2 = s1;
858 }
859 else
860 {
861 // Case 2: flip about 0-4 diagonal
862 i01 = s1+1;
863 i2 = s0;
864 zeta = 2*xe_fraction-1;
865 xi_eta(0) = (zeta_saved+1)*xe_scale/2;
866 }
867 }
868 else if (elem->point(3) == min_point)
869 {
870 if (!elem->positive_face_orientation(1))
871 {
872 // Case 3: 0->3->4->1->0 rotation
873 i01 = s1+1;
874 i2 = s0;
875 zeta = 1-2*xe_fraction;
876 xi_eta(0) = (zeta_saved+1)*xe_scale/2;
877 }
878 else
879 {
880 // Case 4: flip about 9-10 midline
881 i01 = s0+1;
882 i2 = s1;
883 zeta = -zeta_saved;
884 }
885 }
886 else if (elem->point(1) == min_point)
887 {
888 if (!elem->positive_face_orientation(1))
889 {
890 // Case 5: 0->1->4->3->0 rotation
891 i01 = s1+1;
892 i2 = s0;
893 zeta = 2*xe_fraction-1;
894 xi_eta(0) = (1-zeta_saved)*xe_scale/2;
895 }
896 else
897 {
898 // Case 6: flip about 6-12 midline
899 i01 = s0+1;
900 i2 = s1;
901 xi_eta(0) = (1-xe_fraction)*xe_scale;
902 }
903 }
904 else if (elem->point(4) == min_point)
905 {
906 if (!elem->positive_face_orientation(1))
907 {
908 // Case 7: 180 degree rotation
909 i01 = s0+1;
910 i2 = s1;
911 xi_eta(0) = (1-xe_fraction)*xe_scale;
912 zeta = -zeta_saved;
913 }
914 else
915 {
916 // Case 8: flip about 1-3 diagonal
917 i01 = s1+1;
918 i2 = s0;
919 zeta = 1-2*xe_fraction;
920 xi_eta(0) = (1-zeta_saved)*xe_scale/2;
921 }
922 }
923 }
924 // Face 2, node 16
925 else if (i < 6 + 9*e + 2*e*e)
926 {
927 unsigned int basisnum = i - 6 - 9*e - e*e;
928
929 // How wide is the stretch from one side to the other of the
930 // line in the xi-eta plane parallel to this face?
931 const Real xe_scale = xi_eta_saved(0) + xi_eta_saved(1);
932
933 // What percentage of the way along that stretch are we?
934 const Real xe_fraction = (xe_scale==0) ?
935 0 : xi_eta_saved(0)/xe_scale;
936
937 // indexes in edge numbering
938 unsigned int s0 = square_number_row[basisnum] + 2;
939 unsigned int s1 = square_number_column[basisnum] + 2;
940 const Point min_point = get_min_point(elem, 1, 2, 4, 5);
941
942 if (elem->point(1) == min_point)
943 {
944 if (!elem->positive_face_orientation(2))
945 {
946 // Case 1: no flips needed
947 i01 = s0+1+3; // edge to triangle side 1 numbering
948 i2 = s1;
949 }
950 else
951 {
952 // Case 2: flip about 1-5 diagonal
953 i01 = s1+1+e;
954 i2 = s0;
955 zeta = 2*xe_fraction-1;
956 const Real xe = (zeta_saved+1)/2;
957 xi_eta(1) = xe*xe_scale;
958 xi_eta(0) = xe_scale - xi_eta(1);
959 }
960 }
961 else if (elem->point(4) == min_point)
962 {
963 if (!elem->positive_face_orientation(2))
964 {
965 // Case 3: 1->4->5->2->1 rotation
966 i01 = s1+1+e;
967 i2 = s0;
968 zeta = 1-2*xe_fraction;
969 const Real xe = (zeta_saved+1)/2;
970 xi_eta(1) = xe*xe_scale;
971 xi_eta(0) = xe_scale - xi_eta(1);
972 }
973 else
974 {
975 // Case 4: flip about 10-11 midline
976 i01 = s0+1+e;
977 i2 = s1;
978 zeta = -zeta_saved;
979 }
980 }
981 else if (elem->point(2) == min_point)
982 {
983 if (!elem->positive_face_orientation(2))
984 {
985 // Case 5: 1->2->5->4->1 rotation
986 i01 = s1+1+e;
987 i2 = s0;
988 zeta = 2*xe_fraction-1;
989 const Real xe = (1-zeta_saved)/2;
990 xi_eta(1) = xe*xe_scale;
991 xi_eta(0) = xe_scale - xi_eta(1);
992 }
993 else
994 {
995 // Case 6: flip about 7-13 midline
996 i01 = s0+1+e;
997 i2 = s1;
998 const Real xe = (1-xe_fraction);
999 xi_eta(1) = xe*xe_scale;
1000 xi_eta(0) = xe_scale - xi_eta(1);
1001 }
1002 }
1003 else if (elem->point(5) == min_point)
1004 {
1005 if (!elem->positive_face_orientation(2))
1006 {
1007 // Case 7: 180 degree rotation
1008 i01 = s0+1+e;
1009 i2 = s1;
1010 zeta = -zeta_saved;
1011 const Real xe = (1-xe_fraction);
1012 xi_eta(1) = xe*xe_scale;
1013 xi_eta(0) = xe_scale - xi_eta(1);
1014 }
1015 else
1016 {
1017 // Case 8: flip about 1-3 diagonal
1018 i01 = s1+1+e;
1019 i2 = s0;
1020 zeta = 1-2*xe_fraction;
1021 const Real xe = (1-zeta_saved)/2;
1022 xi_eta(1) = xe*xe_scale;
1023 xi_eta(0) = xe_scale - xi_eta(1);
1024 }
1025 }
1026 }
1027 // Face 3, node 17
1028 else if (i < 6 + 9*e + 3*e*e)
1029 {
1030 unsigned int basisnum = i - 6 - 9*e - 2*e*e;
1031
1032 // How wide is the stretch from one side to the other of the
1033 // line in the xi-eta plane parallel to this face?
1034 const Real xe_scale = 1 - xi_eta_saved(0);
1035
1036 // What percentage of the way along that stretch are we?
1037 const Real xe_fraction = (xe_scale==0) ?
1038 0 : (xe_scale - xi_eta_saved(1))/xe_scale;
1039
1040 // indexes in edge numbering
1041 unsigned int s0 = square_number_row[basisnum] + 2;
1042 unsigned int s1 = square_number_column[basisnum] + 2;
1043 const Point min_point = get_min_point(elem, 0, 2, 3, 5);
1044
1045 if (elem->point(2) == min_point)
1046 {
1047 if (!elem->positive_face_orientation(3))
1048 {
1049 // Case 1: no flips needed
1050 i01 = s0+1+2*e; // edge to triangle side 2 numbering
1051 i2 = s1;
1052 }
1053 else
1054 {
1055 // Case 2: flip about 2-3 diagonal
1056 i01 = s1+1+2*e;
1057 i2 = s0;
1058 zeta = 2*xe_fraction-1;
1059 const Real xe = (zeta_saved+1)/2;
1060 xi_eta(1) = xe_scale - xe*xe_scale;
1061 }
1062 }
1063 else if (elem->point(5) == min_point)
1064 {
1065 if (!elem->positive_face_orientation(3))
1066 {
1067 // Case 3: 2->5->3->0->2 rotation
1068 i01 = s1+1+2*e;
1069 i2 = s0;
1070 zeta = 1-2*xe_fraction;
1071 const Real xe = (zeta_saved+1)/2;
1072 xi_eta(1) = xe_scale - xe*xe_scale;
1073 }
1074 else
1075 {
1076 // Case 4: flip about 11-9 midline
1077 i01 = s0+1+2*e;
1078 i2 = s1;
1079 zeta = -zeta_saved;
1080 }
1081 }
1082 else if (elem->point(0) == min_point)
1083 {
1084 if (!elem->positive_face_orientation(3))
1085 {
1086 // Case 5: 2->0->3->5->2 rotation
1087 i01 = s1+1+2*e;
1088 i2 = s0;
1089 zeta = 2*xe_fraction-1;
1090 const Real xe = (1-zeta_saved)/2;
1091 xi_eta(1) = xe_scale - xe*xe_scale;
1092 }
1093 else
1094 {
1095 // Case 6: flip about 8-14 midline
1096 i01 = s0+1+2*e;
1097 i2 = s1;
1098 const Real xe = (1-xe_fraction);
1099 xi_eta(1) = xe_scale - xe*xe_scale;
1100 }
1101 }
1102 else if (elem->point(3) == min_point)
1103 {
1104 if (!elem->positive_face_orientation(3))
1105 {
1106 // Case 7: 180 degree rotation
1107 i01 = s0+1+2*e;
1108 i2 = s1;
1109 zeta = -zeta_saved;
1110 const Real xe = (1-xe_fraction);
1111 xi_eta(1) = xe_scale - xe*xe_scale;
1112 }
1113 else
1114 {
1115 // Case 8: flip about 0-5 diagonal
1116 i01 = s1+1+2*e;
1117 i2 = s0;
1118 zeta = 1-2*xe_fraction;
1119 const Real xe = (1-zeta_saved)/2;
1120 xi_eta(1) = xe_scale - xe*xe_scale;
1121 }
1122 }
1123 }
1124 // Face 0, node 18 - node order due to hierarchic numbering
1125 else if (i < 6 + 9*e + 3*e*e + e*(e-1)/2)
1126 {
1127 // The TRI code will handle any flips here
1128 i01 = i - 3 - 6*e - 3*e*e;
1129 i2 = 0;
1130 }
1131 // Face 4
1132 else if (i < 6 + 9*e + 3*e*e + e*(e-1))
1133 {
1134 // The TRI code will handle any flips here
1135 i01 = i - 3 - 6*e - 3*e*e - e*(e-1)/2;
1136 i2 = 1;
1137 }
1138 // Internal DoFs
1139 else
1140 {
1141 // We won't bother with any internal DoF reordering / flipping;
1142 // that's fine unless we ever get to 4D.
1143 unsigned int basisnum = i - 6 - 9*e - 3*e*e - e*(e-1);
1144 i01 = prism_number_triangle[basisnum] + 3 + 3*e;
1145 i2 = prism_number_page[basisnum] + 2;
1146 }
1147}
1148
1149#endif // LIBMESH_DIM > 2
1150
1151} // end anonymous namespace
1152
1153
1154
1155namespace libMesh
1156{
1157
1158
1159LIBMESH_DEFAULT_VECTORIZED_FE(3,HIERARCHIC)
1160LIBMESH_DEFAULT_VECTORIZED_FE(3,L2_HIERARCHIC)
1161LIBMESH_DEFAULT_VECTORIZED_FE(3,SIDE_HIERARCHIC)
1162
1163
1164template <>
1165Real FE<3,HIERARCHIC>::shape(const ElemType,
1166 const Order,
1167 const unsigned int,
1168 const Point &)
1169{
1170 libmesh_error_msg("Hierarchic shape functions require an Elem for edge/face orientation.");
1171 return 0.;
1172}
1173
1174
1175
1176template <>
1177Real FE<3,L2_HIERARCHIC>::shape(const ElemType,
1178 const Order,
1179 const unsigned int,
1180 const Point &)
1181{
1182 libmesh_error_msg("Hierarchic shape functions require an Elem for edge/face orientation.");
1183 return 0.;
1184}
1185
1186
1187
1188template <>
1189Real FE<3,SIDE_HIERARCHIC>::shape(const ElemType,
1190 const Order,
1191 const unsigned int,
1192 const Point &)
1193{
1194 libmesh_error_msg("Hierarchic shape functions require an Elem for edge/face orientation.");
1195 return 0.;
1196}
1197
1198
1199
1200template <>
1201Real FE<3,HIERARCHIC>::shape(const Elem * elem,
1202 const Order order,
1203 const unsigned int i,
1204 const Point & p,
1205 const bool add_p_level)
1206{
1207 return fe_hierarchic_3D_shape<HIERARCHIC>(elem, order, i, p, add_p_level);
1208}
1209
1210
1211template <>
1212Real FE<3,HIERARCHIC>::shape(const FEType fet,
1213 const Elem * elem,
1214 const unsigned int i,
1215 const Point & p,
1216 const bool add_p_level)
1217{
1218 return fe_hierarchic_3D_shape<HIERARCHIC>(elem, fet.order, i, p, add_p_level);
1219}
1220
1221
1222
1223
1224template <>
1225Real FE<3,L2_HIERARCHIC>::shape(const Elem * elem,
1226 const Order order,
1227 const unsigned int i,
1228 const Point & p,
1229 const bool add_p_level)
1230{
1231 return fe_hierarchic_3D_shape<L2_HIERARCHIC>(elem, order, i, p, add_p_level);
1232}
1233
1234
1235template <>
1236Real FE<3,L2_HIERARCHIC>::shape(const FEType fet,
1237 const Elem * elem,
1238 const unsigned int i,
1239 const Point & p,
1240 const bool add_p_level)
1241{
1242 return fe_hierarchic_3D_shape<L2_HIERARCHIC>(elem, fet.order, i, p, add_p_level);
1243}
1244
1245
1246
1247template <>
1248Real FE<3,SIDE_HIERARCHIC>::shape(const Elem * elem,
1249 const Order order,
1250 const unsigned int i,
1251 const Point & p,
1252 const bool add_p_level)
1253{
1254#if LIBMESH_DIM == 3
1255 libmesh_assert(elem);
1256 const ElemType type = elem->type();
1257
1258 const Order totalorder = order + add_p_level*elem->p_level();
1259
1260 switch (type)
1261 {
1262 case HEX27:
1263 {
1264 const unsigned int dofs_per_side = (totalorder+1u)*(totalorder+1u);
1265 libmesh_assert_less(i, 6*dofs_per_side);
1266
1267 const unsigned int sidenum = cube_side(p);
1268 if (sidenum > 5)
1269 return std::numeric_limits<Real>::quiet_NaN();
1270
1271 const unsigned int dof_offset = sidenum * dofs_per_side;
1272
1273 if (i < dof_offset) // i is on a previous side
1274 return 0;
1275
1276 if (i >= dof_offset + dofs_per_side) // i is on a later side
1277 return 0;
1278
1279 if (totalorder == 0) // special case since raw HIERARCHIC lacks CONSTANTs
1280 return 1;
1281
1282 unsigned int side_i = i - dof_offset;
1283
1284 std::unique_ptr<const Elem> side = elem->build_side_ptr(sidenum);
1285
1286 Point sidep = cube_side_point(sidenum, p);
1287
1288 cube_remap(side_i, *side, totalorder, sidep);
1289
1290 return FE<2,HIERARCHIC>::shape(side.get(), order, side_i, sidep, add_p_level);
1291 }
1292
1293 case TET14:
1294 {
1295 const unsigned int dofs_per_side = (totalorder+1u)*(totalorder+2u)/2u;
1296 libmesh_assert_less(i, 4*dofs_per_side);
1297
1298 const Real zeta[4] = { Real(1.) - p(0) - p(1) - p(2), p(0), p(1), p(2) };
1299
1300 unsigned int face_num = 0;
1301 if (zeta[0] > zeta[3] &&
1302 zeta[1] > zeta[3] &&
1303 zeta[2] > zeta[3])
1304 {
1305 face_num = 0;
1306 }
1307 else if (zeta[0] > zeta[2] &&
1308 zeta[1] > zeta[2] &&
1309 zeta[3] > zeta[2])
1310 {
1311 face_num = 1;
1312 }
1313 else if (zeta[1] > zeta[0] &&
1314 zeta[2] > zeta[0] &&
1315 zeta[3] > zeta[0])
1316 {
1317 face_num = 2;
1318 }
1319 else
1320 {
1321 // We'd better not be right between two faces
1322 libmesh_assert (zeta[0] > zeta[1] &&
1323 zeta[2] > zeta[1] &&
1324 zeta[3] > zeta[1]);
1325 face_num = 3;
1326 }
1327
1328 if (i < face_num * dofs_per_side ||
1329 i >= (face_num+1) * dofs_per_side)
1330 return 0;
1331
1332 if (totalorder == 0)
1333 return 1;
1334
1335 const std::array<unsigned int, 3> face_vertex =
1336 oriented_tet_nodes(*elem, face_num);
1337
1338 // We only need a Tri3 to evaluate L2_HIERARCHIC on the affine
1339 // master element
1340 Tri3 side;
1341
1342 // We pinky swear not to modify these nodes
1343 Elem & e = const_cast<Elem &>(*elem);
1344 side.set_node(0, e.node_ptr(face_vertex[0]));
1345 side.set_node(1, e.node_ptr(face_vertex[1]));
1346 side.set_node(2, e.node_ptr(face_vertex[2]));
1347
1348 const unsigned int basisnum = i - face_num*dofs_per_side;
1349
1350 Point sidep {zeta[face_vertex[1]], zeta[face_vertex[2]]};
1351
1352 return FE<2,L2_HIERARCHIC>::shape(&side, totalorder,
1353 basisnum, sidep, false);
1354 }
1355
1356 case PRISM20:
1357 case PRISM21:
1358 {
1359 const unsigned int dofs_per_quad = (totalorder+1u)*(totalorder+1u);
1360 const unsigned int dofs_per_tri = (totalorder+1u)*(totalorder+2u)/2u;
1361 libmesh_assert_less(i, 3*dofs_per_quad + 2*dofs_per_tri);
1362
1363 // We only need a Tri3 or Quad4 to evaluate L2_HIERARCHIC on
1364 // the affine master element
1365 Tri3 tri;
1366 Quad4 quad;
1367 Elem * side = &quad;
1368 unsigned int dofs_on_side = dofs_per_quad;
1369
1370 // We pinky swear not to modify the nodes we'll point to
1371 Elem & e = const_cast<Elem &>(*elem);
1372
1373 Point sidep;
1374
1375 // Face number calculation is tricky - the ordering of side
1376 // nodes on Prisms does *not* match the ordering of sides!
1377 // (the mid-triangle side nodes were added "later")
1378 // Here face_num will be the numbering that matches the side
1379 // number, but i_offset will have to consider the nodal
1380 // ordering.
1381 unsigned int face_num = 0;
1382 unsigned int i_offset = 0;
1383
1384 // Triangular coordinates
1385 const Real zeta[3] = { Real(1.) - p(0) - p(1), p(0), p(1) };
1386
1387 // Closeness to midplane
1388 const Real zmid = 1 - std::abs(p(2));
1389
1390 if (zeta[1] > zeta[2] && zeta[0] > zeta[2] &&
1391 zmid > 3*zeta[2]) // face 1, quad
1392 {
1393 face_num = 1;
1394 i_offset = 0;
1395 }
1396 else if (zeta[1] > zeta[0] && zeta[2] > zeta[0] &&
1397 zmid > 3*zeta[0]) // face 2, quad
1398 {
1399 face_num = 2;
1400 i_offset = dofs_per_quad;
1401 }
1402 else if (zeta[0] > zeta[1] && zeta[2] > zeta[1] &&
1403 zmid > 3*zeta[1]) // face 3, quad
1404 {
1405 face_num = 3;
1406 i_offset = 2*dofs_per_quad;
1407 }
1408 else if (p(2) + 1 < 3*zeta[0] &&
1409 p(2) + 1 < 3*zeta[1] &&
1410 p(2) + 1 < 3*zeta[2]) // face 0, tri
1411 {
1412 face_num = 0;
1413 i_offset = 3*dofs_per_quad;
1414 dofs_on_side = dofs_per_tri;
1415 side = &tri;
1416 }
1417 else if (1 - p(2) < 3*zeta[0] &&
1418 1 - p(2) < 3*zeta[1] &&
1419 1 - p(2) < 3*zeta[2]) // face 4, tri
1420 {
1421 face_num = 4;
1422 i_offset = dofs_per_tri + 3*dofs_per_quad;
1423 dofs_on_side = dofs_per_tri;
1424 side = &tri;
1425 }
1426 else
1427 {
1428 libmesh_error_msg("Evaluating SIDE_HIERARCHIC right between two Prism faces?");
1429 }
1430
1431 if (i < i_offset ||
1432 i >= i_offset + dofs_on_side)
1433 return 0;
1434
1435 if (totalorder == 0)
1436 return 1;
1437
1438 const std::array<unsigned int, 4> face_vertex =
1439 oriented_prism_nodes(*elem, face_num);
1440
1441 side->set_node(0, e.node_ptr(face_vertex[0]));
1442 side->set_node(1, e.node_ptr(face_vertex[1]));
1443 side->set_node(2, e.node_ptr(face_vertex[2]));
1444 if (face_vertex[3] < 21)
1445 side->set_node(3, e.node_ptr(face_vertex[3]));
1446
1447 if (face_num == 0 || face_num == 4)
1448 sidep = {zeta[face_vertex[1]%3], zeta[face_vertex[2]%3]};
1449 else
1450 {
1451 // Transform a coordinate from the master prism to the
1452 // master quad, based on two vertex indices defining the
1453 // coordinate's direction
1454 auto coord_val = [p](int v1, int v2){
1455 if (v2-v1 == 3)
1456 return p(2);
1457 else if (v2-v1 == -3)
1458 return -p(2);
1459 else if (v1%3 == 0 && v2%3 == 1)
1460 return 2*p(0)-1;
1461 else if (v2%3 == 0 && v1%3 == 1)
1462 return 1-2*p(0);
1463 else if (v1%3 == 1 && v2%3 == 2)
1464 return p(1)-p(0);
1465 else if (v2%3 == 1 && v1%3 == 2)
1466 return p(0)-p(1);
1467 else if (v1%3 == 2 && v2%3 == 0)
1468 return 1-2*p(1);
1469 else if (v2%3 == 2 && v1%3 == 0)
1470 return 2*p(1)-1;
1471 else
1472 libmesh_error();
1473 };
1474
1475 sidep = {coord_val(face_vertex[0], face_vertex[1]),
1476 coord_val(face_vertex[0], face_vertex[3])};
1477 }
1478
1479 const unsigned int basisnum = i - i_offset;
1480
1481 return FE<2,L2_HIERARCHIC>::shape(side, totalorder,
1482 basisnum, sidep, false);
1483 }
1484
1485
1486 default:
1487 libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
1488 }
1489
1490#else // LIBMESH_DIM != 3
1491 libmesh_ignore(elem, order, i, p, add_p_level);
1492 libmesh_not_implemented();
1493#endif
1494}
1495
1496
1497template <>
1498Real FE<3,SIDE_HIERARCHIC>::shape(const FEType fet,
1499 const Elem * elem,
1500 const unsigned int i,
1501 const Point & p,
1502 const bool add_p_level)
1503{
1504 return FE<3,SIDE_HIERARCHIC>::shape(elem,fet.order, i, p, add_p_level);
1505}
1506
1507
1508template <>
1509Real FE<3,HIERARCHIC>::shape_deriv(const ElemType,
1510 const Order,
1511 const unsigned int,
1512 const unsigned int,
1513 const Point & )
1514{
1515 libmesh_error_msg("Hierarchic shape functions require an Elem for edge/face orientation.");
1516 return 0.;
1517}
1518
1519
1520
1521template <>
1522Real FE<3,L2_HIERARCHIC>::shape_deriv(const ElemType,
1523 const Order,
1524 const unsigned int,
1525 const unsigned int,
1526 const Point & )
1527{
1528 libmesh_error_msg("Hierarchic shape functions require an Elem for edge/face orientation.");
1529 return 0.;
1530}
1531
1532
1533
1534template <>
1535Real FE<3,SIDE_HIERARCHIC>::shape_deriv(const ElemType,
1536 const Order,
1537 const unsigned int,
1538 const unsigned int,
1539 const Point & )
1540{
1541 libmesh_error_msg("Hierarchic shape functions require an Elem for edge/face orientation.");
1542 return 0.;
1543}
1544
1545
1546
1547template <>
1548Real FE<3,HIERARCHIC>::shape_deriv(const Elem * elem,
1549 const Order order,
1550 const unsigned int i,
1551 const unsigned int j,
1552 const Point & p,
1553 const bool add_p_level)
1554{
1555 return fe_hierarchic_3D_shape_deriv<HIERARCHIC>(elem, order, i, j, p, add_p_level);
1556}
1557
1558
1559template <>
1560Real FE<3,HIERARCHIC>::shape_deriv(const FEType fet,
1561 const Elem * elem,
1562 const unsigned int i,
1563 const unsigned int j,
1564 const Point & p,
1565 const bool add_p_level)
1566{
1567 return fe_hierarchic_3D_shape_deriv<HIERARCHIC>(elem, fet.order, i, j, p, add_p_level);
1568}
1569
1570
1571
1572template <>
1573Real FE<3,L2_HIERARCHIC>::shape_deriv(const Elem * elem,
1574 const Order order,
1575 const unsigned int i,
1576 const unsigned int j,
1577 const Point & p,
1578 const bool add_p_level)
1579{
1580 return fe_hierarchic_3D_shape_deriv<L2_HIERARCHIC>(elem, order, i, j, p, add_p_level);
1581}
1582
1583
1584template <>
1585Real FE<3,L2_HIERARCHIC>::shape_deriv(const FEType fet,
1586 const Elem * elem,
1587 const unsigned int i,
1588 const unsigned int j,
1589 const Point & p,
1590 const bool add_p_level)
1591{
1592 return fe_hierarchic_3D_shape_deriv<L2_HIERARCHIC>(elem, fet.order, i, j, p, add_p_level);
1593}
1594
1595
1596
1597template <>
1598Real FE<3,SIDE_HIERARCHIC>::shape_deriv(const Elem * elem,
1599 const Order order,
1600 const unsigned int i,
1601 const unsigned int j,
1602 const Point & p,
1603 const bool add_p_level)
1604{
1605#if LIBMESH_DIM == 3
1606 libmesh_assert(elem);
1607 const ElemType type = elem->type();
1608
1609 const Order totalorder = order + add_p_level*elem->p_level();
1610
1611 if (totalorder == 0) // special case since raw HIERARCHIC lacks CONSTANTs
1612 return 0; // constants have zero derivative
1613
1614 switch (type)
1615 {
1616 case HEX27:
1617 {
1618 // I need to debug the p>2 case here...
1619 if (totalorder > 2)
1620 return fe_fdm_deriv(elem, order, i, j, p, add_p_level, FE<3,SIDE_HIERARCHIC>::shape);
1621
1622 const unsigned int dofs_per_side = (totalorder+1u)*(totalorder+1u);
1623 libmesh_assert_less(i, 6*dofs_per_side);
1624
1625 const unsigned int sidenum = cube_side(p);
1626 if (sidenum > 5)
1627 return std::numeric_limits<Real>::quiet_NaN();
1628
1629 const unsigned int dof_offset = sidenum * dofs_per_side;
1630
1631 if (i < dof_offset) // i is on a previous side
1632 return 0;
1633
1634 if (i >= dof_offset + dofs_per_side) // i is on a later side
1635 return 0;
1636
1637 unsigned int side_i = i - dof_offset;
1638
1639 std::unique_ptr<const Elem> side = elem->build_side_ptr(sidenum);
1640
1641 Point sidep = cube_side_point(sidenum, p);
1642
1643 cube_remap(side_i, *side, totalorder, sidep);
1644
1645 // What direction on the side corresponds to the derivative
1646 // direction we want?
1647 unsigned int sidej = 100;
1648
1649 // Do we need a -1 here to flip that direction?
1650 Real f = 1.;
1651
1652 switch (j)
1653 {
1654 case 0: // d()/dxi
1655 {
1656 switch (sidenum)
1657 {
1658 case 0:
1659 sidej = 1;
1660 break;
1661 case 1:
1662 sidej = 0;
1663 break;
1664 case 2:
1665 return 0;
1666 case 3:
1667 sidej = 0;
1668 f = -1;
1669 break;
1670 case 4:
1671 return 0;
1672 case 5:
1673 sidej = 0;
1674 break;
1675 default:
1676 libmesh_error();
1677 }
1678 break;
1679 }
1680 case 1: // d()/deta
1681 {
1682 switch (sidenum)
1683 {
1684 case 0:
1685 sidej = 0;
1686 break;
1687 case 1:
1688 return 0;
1689 case 2:
1690 sidej = 0;
1691 break;
1692 case 3:
1693 return 0;
1694 case 4:
1695 sidej = 0;
1696 f = -1;
1697 break;
1698 case 5:
1699 sidej = 1;
1700 break;
1701 default:
1702 libmesh_error();
1703 }
1704 break;
1705 }
1706 case 2: // d()/dzeta
1707 {
1708 switch (sidenum)
1709 {
1710 case 0:
1711 return 0;
1712 case 1:
1713 case 2:
1714 case 3:
1715 case 4:
1716 sidej = 1;
1717 break;
1718 case 5:
1719 return 0;
1720 default:
1721 libmesh_error();
1722 }
1723 break;
1724 }
1725
1726 default:
1727 libmesh_error_msg("Invalid derivative index j = " << j);
1728 }
1729
1730 return f * FE<2,HIERARCHIC>::shape_deriv(side.get(), order,
1731 side_i, sidej, sidep,
1732 add_p_level);
1733 }
1734
1735 case TET14:
1736 case PRISM20:
1737 case PRISM21:
1738 {
1739 return fe_fdm_deriv(elem, order, i, j, p, add_p_level, FE<3,SIDE_HIERARCHIC>::shape);
1740 }
1741
1742 default:
1743 libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
1744 }
1745
1746#else // LIBMESH_DIM != 3
1747 libmesh_ignore(elem, order, i, j, p, add_p_level);
1748 libmesh_not_implemented();
1749#endif
1750}
1751
1752
1753template <>
1754Real FE<3,SIDE_HIERARCHIC>::shape_deriv(const FEType fet,
1755 const Elem * elem,
1756 const unsigned int i,
1757 const unsigned int j,
1758 const Point & p,
1759 const bool add_p_level)
1760{
1761 return FE<3,SIDE_HIERARCHIC>::shape_deriv(elem, fet.order, i, j, p, add_p_level);
1762}
1763
1764
1765#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1766
1767template <>
1768Real FE<3,HIERARCHIC>::shape_second_deriv(const ElemType,
1769 const Order,
1770 const unsigned int,
1771 const unsigned int,
1772 const Point & )
1773{
1774 libmesh_error_msg("Hierarchic shape functions require an Elem for edge/face orientation.");
1775 return 0.;
1776}
1777
1778
1779
1780template <>
1781Real FE<3,L2_HIERARCHIC>::shape_second_deriv(const ElemType,
1782 const Order,
1783 const unsigned int,
1784 const unsigned int,
1785 const Point & )
1786{
1787 libmesh_error_msg("Hierarchic shape functions require an Elem for edge/face orientation.");
1788 return 0.;
1789}
1790
1791
1792
1793template <>
1794Real FE<3,SIDE_HIERARCHIC>::shape_second_deriv(const ElemType,
1795 const Order,
1796 const unsigned int,
1797 const unsigned int,
1798 const Point & )
1799{
1800 libmesh_error_msg("Hierarchic shape functions require an Elem for edge/face orientation.");
1801 return 0.;
1802}
1803
1804
1805
1806template <>
1807Real FE<3,HIERARCHIC>::shape_second_deriv(const Elem * elem,
1808 const Order order,
1809 const unsigned int i,
1810 const unsigned int j,
1811 const Point & p,
1812 const bool add_p_level)
1813{
1814 return fe_hierarchic_3D_shape_second_deriv<HIERARCHIC>(elem, order, i, j, p, add_p_level);
1815}
1816
1817
1818
1819template <>
1820Real FE<3,HIERARCHIC>::shape_second_deriv(const FEType fet,
1821 const Elem * elem,
1822 const unsigned int i,
1823 const unsigned int j,
1824 const Point & p,
1825 const bool add_p_level)
1826{
1827 return fe_hierarchic_3D_shape_second_deriv<HIERARCHIC>(elem, fet.order, i, j, p, add_p_level);
1828}
1829
1830
1831
1832template <>
1833Real FE<3,L2_HIERARCHIC>::shape_second_deriv(const Elem * elem,
1834 const Order order,
1835 const unsigned int i,
1836 const unsigned int j,
1837 const Point & p,
1838 const bool add_p_level)
1839{
1840 return fe_hierarchic_3D_shape_second_deriv<L2_HIERARCHIC>(elem, order, i, j, p, add_p_level);
1841}
1842
1843
1844template <>
1845Real FE<3,L2_HIERARCHIC>::shape_second_deriv(const FEType fet,
1846 const Elem * elem,
1847 const unsigned int i,
1848 const unsigned int j,
1849 const Point & p,
1850 const bool add_p_level)
1851{
1852 return fe_hierarchic_3D_shape_second_deriv<L2_HIERARCHIC>(elem, fet.order, i, j, p, add_p_level);
1853}
1854
1855
1856template <>
1857Real FE<3,SIDE_HIERARCHIC>::shape_second_deriv(const Elem * elem,
1858 const Order order,
1859 const unsigned int i,
1860 const unsigned int j,
1861 const Point & p,
1862 const bool add_p_level)
1863{
1864#if LIBMESH_DIM == 3
1865 libmesh_assert(elem);
1866 const ElemType type = elem->type();
1867
1868 const Order totalorder = order + add_p_level*elem->p_level();
1869
1870 if (totalorder == 0) // special case since raw HIERARCHIC lacks CONSTANTs
1871 return 0; // constants have zero derivative
1872
1873 switch (type)
1874 {
1875 case HEX27:
1876 {
1877 // I need to debug the p>2 case here...
1878 if (totalorder > 2)
1879 return fe_fdm_second_deriv(elem, order, i, j, p, add_p_level,
1880 FE<3,SIDE_HIERARCHIC>::shape_deriv);
1881
1882 const unsigned int dofs_per_side = (totalorder+1u)*(totalorder+1u);
1883 libmesh_assert_less(i, 6*dofs_per_side);
1884
1885 const unsigned int sidenum = cube_side(p);
1886 if (sidenum > 5)
1887 return std::numeric_limits<Real>::quiet_NaN();
1888
1889 const unsigned int dof_offset = sidenum * dofs_per_side;
1890
1891 if (i < dof_offset) // i is on a previous side
1892 return 0;
1893
1894 if (i >= dof_offset + dofs_per_side) // i is on a later side
1895 return 0;
1896
1897 unsigned int side_i = i - dof_offset;
1898
1899 std::unique_ptr<const Elem> side = elem->build_side_ptr(sidenum);
1900
1901 Point sidep = cube_side_point(sidenum, p);
1902
1903 cube_remap(side_i, *side, totalorder, sidep);
1904
1905 // What second derivative or mixed derivative on the side
1906 // corresponds to the xi/eta/zeta mix we were asked for?
1907 unsigned int sidej = 100;
1908
1909 // Do we need a -1 here to flip the final derivative value?
1910 Real f = 1.;
1911
1912 switch (j)
1913 {
1914 case 0: // d^2()/dxi^2
1915 {
1916 switch (sidenum)
1917 {
1918 case 0:
1919 sidej = 2;
1920 break;
1921 case 1:
1922 sidej = 0;
1923 break;
1924 case 2:
1925 return 0;
1926 case 3:
1927 sidej = 0;
1928 break;
1929 case 4:
1930 return 0;
1931 case 5:
1932 sidej = 0;
1933 break;
1934 default:
1935 libmesh_error();
1936 }
1937 break;
1938 }
1939 case 1: // d^2()/dxideta
1940 {
1941 switch (sidenum)
1942 {
1943 case 0:
1944 sidej = 1;
1945 break;
1946 case 1:
1947 case 2:
1948 case 3:
1949 case 4:
1950 return 0;
1951 case 5:
1952 sidej = 1;
1953 break;
1954 default:
1955 libmesh_error();
1956 }
1957 break;
1958 }
1959 case 2: // d^2()/deta^2
1960 {
1961 switch (sidenum)
1962 {
1963 case 0:
1964 sidej = 0;
1965 break;
1966 case 1:
1967 return 0;
1968 case 2:
1969 sidej = 0;
1970 break;
1971 case 3:
1972 return 0;
1973 case 4:
1974 sidej = 0;
1975 break;
1976 case 5:
1977 sidej = 2;
1978 break;
1979 default:
1980 libmesh_error();
1981 }
1982 break;
1983 }
1984 case 3: // d^2()/dxidzeta
1985 {
1986 switch (sidenum)
1987 {
1988 case 0:
1989 return 0;
1990 case 1:
1991 sidej = 1;
1992 break;
1993 case 2:
1994 return 0;
1995 case 3:
1996 sidej = 1;
1997 f = -1;
1998 break;
1999 case 4:
2000 case 5:
2001 return 0;
2002 default:
2003 libmesh_error();
2004 }
2005 break;
2006 }
2007 case 4: // d^2()/detadzeta
2008 {
2009 switch (sidenum)
2010 {
2011 case 0:
2012 case 1:
2013 return 0;
2014 case 2:
2015 sidej = 1;
2016 break;
2017 case 3:
2018 return 0;
2019 case 4:
2020 sidej = 1;
2021 f = -1;
2022 break;
2023 case 5:
2024 return 0;
2025 default:
2026 libmesh_error();
2027 }
2028 break;
2029 }
2030 case 5: // d^2()/dzeta^2
2031 {
2032 switch (sidenum)
2033 {
2034 case 0:
2035 return 0;
2036 case 1:
2037 case 2:
2038 case 3:
2039 case 4:
2040 sidej = 2;
2041 break;
2042 case 5:
2043 return 0;
2044 default:
2045 libmesh_error();
2046 }
2047 break;
2048 }
2049
2050 default:
2051 libmesh_error_msg("Invalid derivative index j = " << j);
2052 }
2053
2054 return f * FE<2,HIERARCHIC>::shape_second_deriv(side.get(),
2055 order, side_i,
2056 sidej, sidep,
2057 add_p_level);
2058 }
2059
2060 case TET14:
2061 case PRISM20:
2062 case PRISM21:
2063 {
2064 return fe_fdm_second_deriv(elem, order, i, j, p, add_p_level,
2065 FE<3,SIDE_HIERARCHIC>::shape_deriv);
2066 }
2067
2068 default:
2069 libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
2070 }
2071
2072#else // LIBMESH_DIM != 3
2073 libmesh_ignore(elem, order, i, j, p, add_p_level);
2074 libmesh_not_implemented();
2075#endif
2076}
2077
2078
2079template <>
2080Real FE<3,SIDE_HIERARCHIC>::shape_second_deriv(const FEType fet,
2081 const Elem * elem,
2082 const unsigned int i,
2083 const unsigned int j,
2084 const Point & p,
2085 const bool add_p_level)
2086{
2087 return FE<3,SIDE_HIERARCHIC>::shape_second_deriv(elem, fet.order, i, j, p, add_p_level);
2088}
2089
2090#endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
2091
2092} // namespace libMesh
2093
2094
2095
2096namespace
2097{
2098using namespace libMesh;
2099
2100
2101unsigned int cube_side (const Point & p)
2102{
2103 const Real xi = p(0), eta = p(1), zeta = p(2);
2104 const Real absxi = std::abs(xi),
2105 abseta = std::abs(eta),
2106 abszeta = std::abs(zeta);
2107 const Real maxabs_xi_eta = std::max(absxi, abseta),
2108 maxabs_xi_zeta = std::max(absxi, abszeta),
2109 maxabs_eta_zeta = std::max(abseta, abszeta);
2110
2111 if (zeta < -maxabs_xi_eta)
2112 return 0;
2113 else if (eta < -maxabs_xi_zeta)
2114 return 1;
2115 else if (xi > maxabs_eta_zeta)
2116 return 2;
2117 else if (eta > maxabs_xi_zeta)
2118 return 3;
2119 else if (xi < -maxabs_eta_zeta)
2120 return 4;
2121 else if (zeta > maxabs_xi_eta)
2122 return 5;
2123
2124 // We need to be able to return invalid values for cases where
2125 // mixed FE are being evaluated together on edges and vertices
2126 return 65535;
2127}
2128
2129
2130
2131Point cube_side_point(unsigned int sidenum, const Point & p)
2132{
2133 Point sidep;
2134
2135 switch (sidenum)
2136 {
2137 case 0:
2138 sidep(0) = p(1);
2139 sidep(1) = p(0);
2140 break;
2141 case 1:
2142 sidep(0) = p(0);
2143 sidep(1) = p(2);
2144 break;
2145 case 2:
2146 sidep(0) = p(1);
2147 sidep(1) = p(2);
2148 break;
2149 case 3:
2150 sidep(0) = -p(0);
2151 sidep(1) = p(2);
2152 break;
2153 case 4:
2154 sidep(0) = -p(1);
2155 sidep(1) = p(2);
2156 break;
2157 case 5:
2158 sidep(0) = p(0);
2159 sidep(1) = p(1);
2160 break;
2161 default:
2162 libmesh_error();
2163 }
2164
2165 return sidep;
2166}
2167
2168
2169void orient_quad(const Elem & elem,
2170 std::array<unsigned int, 4> & face_vertex)
2171{
2172 // Sort the minimum point into face_vertex[0], the minimum of its
2173 // neighbors into face_vertex[1]. Keep the other two consistent; we
2174 // want to rotate or flip the quad but not to twist it.
2175
2176 const unsigned int min_pt =
2177 std::min_element(face_vertex.begin(), face_vertex.end(),
2178 [&elem](auto v1, auto v2)
2179 {return elem.point(v1)<elem.point(v2);}) -
2180 face_vertex.begin();
2181
2182 // Do we flip the quad?
2183 if (elem.point(face_vertex[(min_pt+3)%4]) <
2184 elem.point(face_vertex[(min_pt+1)%4]))
2185 face_vertex = { face_vertex[min_pt], face_vertex[(min_pt+3)%4],
2186 face_vertex[(min_pt+2)%4], face_vertex[(min_pt+1)%4] };
2187 else
2188 face_vertex = { face_vertex[min_pt], face_vertex[(min_pt+1)%4],
2189 face_vertex[(min_pt+2)%4], face_vertex[(min_pt+3)%4] };
2190}
2191
2192
2193void orient_triangle(const Elem & elem,
2194 unsigned int * face_vertex)
2195{
2196 // Reorient nodes to account for flipping and rotation.
2197 // We could try to identify indices with symmetric shape
2198 // functions, to skip this in those cases, if we really
2199 // need to optimize later.
2200 //
2201 // With only 3 items, we should bubble sort!
2202 // Programming-for-MechE's class pays off!
2203 bool lastcheck = true;
2204 if (elem.point(face_vertex[0]) > elem.point(face_vertex[1]))
2205 {
2206 std::swap(face_vertex[0], face_vertex[1]);
2207 lastcheck = true;
2208 }
2209 if (elem.point(face_vertex[1]) > elem.point(face_vertex[2]))
2210 std::swap(face_vertex[1], face_vertex[2]);
2211 if (lastcheck && elem.point(face_vertex[0]) > elem.point(face_vertex[1]))
2212 std::swap(face_vertex[0], face_vertex[1]);
2213}
2214
2215
2216std::array<unsigned int, 4> oriented_prism_nodes(const Elem & elem,
2217 unsigned int face_num)
2218{
2219 std::array<unsigned int, 4> face_vertex
2220 { Prism6::side_nodes_map[face_num][0],
2221 Prism6::side_nodes_map[face_num][1],
2222 Prism6::side_nodes_map[face_num][2],
2223 Prism6::side_nodes_map[face_num][3] };
2224
2225 if (face_num > 0 && face_num < 4)
2226 orient_quad(elem, face_vertex);
2227 else
2228 orient_triangle(elem, face_vertex.data());
2229
2230 return face_vertex;
2231}
2232
2233
2234std::array<unsigned int, 3> oriented_tet_nodes(const Elem & elem,
2235 unsigned int face_num)
2236{
2237 std::array<unsigned int, 3> face_vertex
2238 { Tet4::side_nodes_map[face_num][0],
2239 Tet4::side_nodes_map[face_num][1],
2240 Tet4::side_nodes_map[face_num][2] };
2241
2242 orient_triangle(elem, face_vertex.data());
2243
2244 return face_vertex;
2245}
2246
2247
2248template <FEFamily T>
2249Real fe_hierarchic_3D_shape(const Elem * elem,
2250 const Order order,
2251 const unsigned int i,
2252 const Point & p,
2253 const bool add_p_level)
2254{
2255#if LIBMESH_DIM == 3
2256
2257 libmesh_assert(elem);
2258 const ElemType type = elem->type();
2259
2260 const Order totalorder = order + add_p_level*elem->p_level();
2261
2262 switch (type)
2263 {
2264 case HEX8:
2265 case HEX20:
2266 libmesh_assert (T == L2_HIERARCHIC || totalorder < 2);
2267 libmesh_fallthrough();
2268 case HEX27:
2269 {
2270 libmesh_assert_less (i, (totalorder+1u)*(totalorder+1u)*(totalorder+1u));
2271
2272 // Compute hex shape functions as a tensor-product
2273 Real xi = p(0);
2274 Real eta = p(1);
2275 Real zeta = p(2);
2276
2277 unsigned int i0, i1, i2;
2278
2279 cube_indices(elem, totalorder, i, xi, eta, zeta, i0, i1, i2);
2280
2281 return (FE<1,T>::shape(EDGE3, totalorder, i0, xi)*
2282 FE<1,T>::shape(EDGE3, totalorder, i1, eta)*
2283 FE<1,T>::shape(EDGE3, totalorder, i2, zeta));
2284 }
2285
2286 case PRISM6:
2287 case PRISM15:
2288 libmesh_assert (T == L2_HIERARCHIC || totalorder < 2);
2289 libmesh_fallthrough();
2290 case PRISM18:
2291 libmesh_assert (T == L2_HIERARCHIC || totalorder < 3);
2292 libmesh_fallthrough();
2293 case PRISM20:
2294 case PRISM21:
2295 {
2296 libmesh_assert_less (i, (totalorder+1u)*(totalorder+1u)*(totalorder+2u)/2u);
2297
2298 // Compute prism shape functions as a tensor-product.
2299 // Non-const here, because prism_indices might need to do some
2300 // flips before evaluating edge or face DoFs
2301 Point xi_eta {p(0),p(1)};
2302 Real zeta = p(2);
2303
2304 unsigned int i01, i2;
2305
2306 prism_indices(elem, totalorder, i, xi_eta, zeta, i01, i2);
2307
2308 // We'll use the 2D Tri to handle any basis function flipping
2309 // needed in xi/eta for triangle face+edge bases.
2310 Tri3 tri;
2311
2312 // We pinky swear not to modify these nodes
2313 Elem & e = const_cast<Elem &>(*elem);
2314 if (i2 == 0)
2315 {
2316 tri.set_node(0, e.node_ptr(0));
2317 tri.set_node(1, e.node_ptr(1));
2318 tri.set_node(2, e.node_ptr(2));
2319 }
2320 else if (i2 == 1)
2321 {
2322 tri.set_node(0, e.node_ptr(3));
2323 tri.set_node(1, e.node_ptr(4));
2324 tri.set_node(2, e.node_ptr(5));
2325 }
2326 else
2327 {
2328 // For interior DoFs, no flipping is necessary or done; we
2329 // can just evaluate on any triangle ... but *not* the
2330 // obvious 9,10,11 triangle, because that might not exist
2331 // if we have L2_HIERARCHIC on Prism6.
2332 tri.set_node(0, e.node_ptr(0));
2333 tri.set_node(1, e.node_ptr(1));
2334 tri.set_node(2, e.node_ptr(2));
2335
2336 // For square face DoFs, prism_indices handles flipping,
2337 // and we *can't* override that in the tri shape call.
2338 if (i01 > 2 && i01 < 3u*totalorder)
2339 {
2340 // %(p-1) to find the edge number, %2 for even vs odd
2341 const bool odd_basis = ((i01-1)%(totalorder-1))%2;
2342 if (odd_basis)
2343 {
2344 const int tri_edge = (i01-3)/(totalorder-1);
2345 // Flip nodes now to avoid triggering a shape
2346 // function flip later
2347 if (tri.point(tri_edge) > tri.point((tri_edge+1)%3))
2348 {
2349 Node * n = tri.node_ptr(tri_edge);
2350 tri.set_node(tri_edge, tri.node_ptr((tri_edge+1)%3));
2351 tri.set_node((tri_edge+1)%3, n);
2352 }
2353 }
2354 }
2355 }
2356
2357 return (FE<2,L2_HIERARCHIC>::shape(&tri, totalorder, i01, xi_eta)*
2358 FE<1,L2_HIERARCHIC>::shape(EDGE2, totalorder, i2, zeta));
2359 }
2360
2361 case TET4:
2362 libmesh_assert (T == L2_HIERARCHIC || totalorder < 2);
2363 libmesh_fallthrough();
2364 case TET10:
2365 libmesh_assert (T == L2_HIERARCHIC || totalorder < 3);
2366 libmesh_fallthrough();
2367 case TET14:
2368 {
2369 const Real zeta[4] = { 1 - p(0) - p(1) - p(2), p(0), p(1), p(2) };
2370
2371 // Nodal DoFs
2372 if (i < 4)
2373 return zeta[i];
2374
2375 // Edge DoFs
2376 else if (i < 6u*totalorder - 2u)
2377 {
2378 const unsigned int edge_num = (i - 4) / (totalorder - 1u);
2379 // const int edge_node = edge_num + 4;
2380 const unsigned int basisorder = i - 2 - ((totalorder - 1u) * edge_num);
2381
2382 const unsigned int edgevertex0 = Tet4::edge_nodes_map[edge_num][0],
2383 edgevertex1 = Tet4::edge_nodes_map[edge_num][1];
2384
2385 // Get factors to account for edge-flipping
2386 Real flip = 1;
2387 if (basisorder%2 &&
2388 elem->positive_edge_orientation(edge_num))
2389 flip = -1;
2390
2391 const Real crossval = zeta[edgevertex0] + zeta[edgevertex1];
2392 const Real edgenumerator = zeta[edgevertex1] - zeta[edgevertex0];
2393
2394 if (crossval == 0.) // Yes, exact comparison; we seem numerically stable otherwise
2395 {
2396 unsigned int basisfactorial = 1.;
2397 for (unsigned int n=2; n <= basisorder; ++n)
2398 basisfactorial *= n;
2399
2400 return std::pow(edgenumerator, basisorder) / basisfactorial;
2401 }
2402
2403 const Real edgeval = edgenumerator / crossval;
2404 const Real crossfunc = std::pow(crossval, basisorder);
2405
2406 return flip * crossfunc *
2407 FE<1,HIERARCHIC>::shape(EDGE3, totalorder,
2408 basisorder, edgeval);
2409 }
2410
2411 // Face DoFs
2412 else if (i < 2u*totalorder*totalorder + 2u)
2413 {
2414 const int dofs_per_face = (totalorder - 1u) * (totalorder - 2u) / 2;
2415 const int face_num = (i - (6u*totalorder - 2u)) / dofs_per_face;
2416
2417 const std::array<unsigned int, 3> face_vertex =
2418 oriented_tet_nodes(*elem, face_num);
2419 const Real zeta0 = zeta[face_vertex[0]],
2420 zeta1 = zeta[face_vertex[1]],
2421 zeta2 = zeta[face_vertex[2]];
2422
2423 const unsigned int basisnum =
2424 i - 4 -
2425 (totalorder - 1u) * /*n_edges*/6 -
2426 (dofs_per_face * face_num);
2427
2428 const unsigned int exp0 = triangular_number_column[basisnum] + 1;
2429 const unsigned int exp1 = triangular_number_row[basisnum] + 1 -
2430 triangular_number_column[basisnum];
2431
2432 Real returnval = 1;
2433 for (unsigned int n = 0; n != exp0; ++n)
2434 returnval *= zeta0;
2435 for (unsigned int n = 0; n != exp1; ++n)
2436 returnval *= zeta1;
2437 returnval *= zeta2;
2438 return returnval;
2439 }
2440
2441 // Interior DoFs
2442 else
2443 {
2444 const unsigned int basisnum = i - 2u*totalorder*totalorder - 2u;
2445 const unsigned int exp0 = tetrahedral_number_column[basisnum] + 1;
2446 const unsigned int exp1 = tetrahedral_number_row[basisnum] + 1 -
2447 tetrahedral_number_column[basisnum] -
2448 tetrahedral_number_page[basisnum];
2449 const unsigned int exp2 = tetrahedral_number_page[basisnum] + 1;
2450
2451 Real returnval = 1;
2452 for (unsigned int n = 0; n != exp0; ++n)
2453 returnval *= zeta[0];
2454 for (unsigned int n = 0; n != exp1; ++n)
2455 returnval *= zeta[1];
2456 for (unsigned int n = 0; n != exp2; ++n)
2457 returnval *= zeta[2];
2458 returnval *= zeta[3];
2459 return returnval;
2460 }
2461 }
2462
2463 default:
2464 libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
2465 }
2466
2467#else // LIBMESH_DIM != 3
2468 libmesh_ignore(elem, order, i, p, add_p_level);
2469 libmesh_not_implemented();
2470#endif
2471}
2472
2473
2474
2475template <FEFamily T>
2476Real fe_hierarchic_3D_shape_deriv(const Elem * elem,
2477 const Order order,
2478 const unsigned int i,
2479 const unsigned int j,
2480 const Point & p,
2481 const bool add_p_level)
2482{
2483 return fe_fdm_deriv(elem, order, i, j, p, add_p_level, FE<3,T>::shape);
2484}
2485
2486
2487
2488#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
2489
2490template <FEFamily T>
2491Real fe_hierarchic_3D_shape_second_deriv(const Elem * elem,
2492 const Order order,
2493 const unsigned int i,
2494 const unsigned int j,
2495 const Point & p,
2496 const bool add_p_level)
2497{
2498 return fe_fdm_second_deriv(elem, order, i, j, p, add_p_level,
2499 FE<3,T>::shape_deriv);
2500}
2501
2502#endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
2503
2504
2505} // anonymous namespace
int
void ErrorVector unsigned int
Definition adjoints_ex3.C:360
libMesh::Elem
This is the base class from which all geometric element types are derived.
Definition elem.h:96
libMesh::Elem::set_node
virtual Node *& set_node(const unsigned int i)
Definition elem.h:2567
libMesh::Elem::point
const Point & point(const unsigned int i) const
Definition elem.h:2462
libMesh::Elem::positive_edge_orientation
bool positive_edge_orientation(const unsigned int i) const
Definition elem.C:3634
libMesh::Elem::positive_face_orientation
bool positive_face_orientation(const unsigned int i) const
Definition elem.C:3643
libMesh::Elem::node_ptr
const Node * node_ptr(const unsigned int i) const
Definition elem.h:2516
libMesh::Elem::type
virtual ElemType type() const =0
libMesh::Elem::p_level
unsigned int p_level() const
Definition elem.h:3125
libMesh::Elem::build_side_ptr
virtual std::unique_ptr< Elem > build_side_ptr(const unsigned int i)=0
libMesh::FEType
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition fe_type.h:197
libMesh::FEType::order
OrderWrapper order
The approximation order of the element (at 0 p-refinement level).
Definition fe_type.h:203
libMesh::FE
A specific instantiation of the FEBase class.
Definition fe.h:128
libMesh::FE::shape
static OutputShape shape(const ElemType t, const Order o, const unsigned int i, const Point &p)
libMesh::FE::shape_second_deriv
static OutputShape shape_second_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)
libMesh::FE::shape_deriv
static OutputShape shape_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)
libMesh::IntRange
The IntRange templated class is intended to make it easy to loop over integers which are indices of a...
Definition int_range.h:54
libMesh::Node
A Node is like a Point, but with more information.
Definition node.h:55
libMesh::Point
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition point.h:40
libMesh::Prism6::side_nodes_map
static const unsigned int side_nodes_map[num_sides][nodes_per_side]
This maps the node of the side to element node numbers.
Definition cell_prism6.h:172
libMesh::Quad4
The QUAD4 is an element in 2D composed of 4 nodes.
Definition face_quad4.h:54
libMesh::Tet4::side_nodes_map
static const unsigned int side_nodes_map[num_sides][nodes_per_side]
This maps the node of the side to element node numbers.
Definition cell_tet4.h:195
libMesh::Tet4::edge_nodes_map
static const unsigned int edge_nodes_map[num_edges][nodes_per_edge]
This maps the node of the edge to element node numbers.
Definition cell_tet4.h:201
libMesh::Tri3
The Tri3 is an element in 2D composed of 3 nodes.
Definition face_tri3.h:62
b
static const Real b
Definition disjoint_neighbor_test.C:24
libMesh::Utility::enum_to_string
std::string enum_to_string(const T e)
libMesh
The libMesh namespace provides an interface to certain functionality in the library.
libMesh::tetrahedral_number_column
const unsigned char tetrahedral_number_column[]
Definition number_lookups.C:421
libMesh::tetrahedral_number_row
const unsigned char tetrahedral_number_row[]
Definition number_lookups.C:488
libMesh::prism_number_triangle
const unsigned char prism_number_triangle[]
Definition number_lookups.C:621
libMesh::LIBMESH_DEFAULT_VECTORIZED_FE
LIBMESH_DEFAULT_VECTORIZED_FE(template<>Real FE< 0, BERNSTEIN)
Definition fe_bernstein_shape_0D.C:30
libMesh::prism_number_page
const unsigned char prism_number_page[]
Definition number_lookups.C:726
libMesh::fe_fdm_deriv
OutputShape fe_fdm_deriv(const Elem *elem, const Order order, const unsigned int i, const unsigned int j, const Point &p, const bool add_p_level, OutputShape(*shape_func)(const Elem *, const Order, const unsigned int, const Point &, const bool))
Helper functions for finite differenced derivatives in cases where analytical calculations haven't be...
Definition fe.C:911
libMesh::ElemType
ElemType
Defines an enum for geometric element types.
Definition enum_elem_type.h:33
libMesh::libmesh_ignore
void libmesh_ignore(const Args &...)
Definition libmesh_common.h:569
libMesh::fe_fdm_second_deriv
OutputShape fe_fdm_second_deriv(const Elem *elem, const Order order, const unsigned int i, const unsigned int j, const Point &p, const bool add_p_level, OutputShape(*deriv_func)(const Elem *, const Order, const unsigned int, const unsigned int, const Point &, const bool))
Definition fe.C:969
libMesh::libmesh_assert
libmesh_assert(ctx)
libMesh::invalid_uint
const unsigned int invalid_uint
A number which is used quite often to represent an invalid or uninitialized value for an unsigned int...
Definition libmesh.h:303
libMesh::triangular_number_row
const unsigned char triangular_number_row[]
Definition number_lookups.C:28
libMesh::cube_number_row
const unsigned char cube_number_row[]
Definition number_lookups.C:196
libMesh::cube_number_page
const unsigned char cube_number_page[]
Definition number_lookups.C:308
libMesh::cube_number_column
const unsigned char cube_number_column[]
Definition number_lookups.C:84
libMesh::square_number_row
const unsigned char square_number_row[]
Definition number_lookups.C:69
libMesh::square_number_column
const unsigned char square_number_column[]
Definition number_lookups.C:56
libMesh::triangular_number_column
const unsigned char triangular_number_column[]
Definition number_lookups.C:41
libMesh::tetrahedral_number_page
const unsigned char tetrahedral_number_page[]
Definition number_lookups.C:550
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition libmesh_common.h:144
Generated on Fri Aug 21 2026 13:32:28 for libMesh by doxygen 1.9.8