libMesh
Loading...
Searching...
No Matches
fe_hierarchic_shape_2D.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
25// Anonymous namespace for functions shared by HIERARCHIC and
26// L2_HIERARCHIC implementations. Implementations appear at the bottom
27// of this file.
28namespace
29{
30using namespace libMesh;
31
32Real fe_triangle_helper (const Elem & elem,
33 const Real edgenumerator,
34 const Real crossval,
35 const unsigned int basisorder,
36 const Order totalorder,
37 const unsigned int noden);
38
39template <FEFamily T>
40Real fe_hierarchic_2D_shape(const Elem * elem,
41 const Order order,
42 const unsigned int i,
43 const Point & p,
44 const bool add_p_level);
45
46template <FEFamily T>
47Real fe_hierarchic_2D_shape_deriv(const Elem * elem,
48 const Order order,
49 const unsigned int i,
50 const unsigned int j,
51 const Point & p,
52 const bool add_p_level);
53
54#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
55
56template <FEFamily T>
57Real fe_hierarchic_2D_shape_second_deriv(const Elem * elem,
58 const Order order,
59 const unsigned int i,
60 const unsigned int j,
61 const Point & p,
62 const bool add_p_level);
63
64#endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
65
66
67std::tuple<unsigned int, unsigned int, Real>
68quad_indices(const Elem * elem,
69 const unsigned int totalorder,
70 const unsigned int i)
71{
72 libmesh_assert_less (i, (totalorder+1u)*(totalorder+1u));
73
74 // Example i, i0, i1 values for totalorder = 5:
75 // 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 28 29 30 31 32 33 34 35
76 // static const unsigned int i0[] = {0, 1, 1, 0, 2, 3, 4, 5, 1, 1, 1, 1, 2, 3, 4, 5, 0, 0, 0, 0, 2, 3, 3, 2, 4, 4, 4, 3, 2, 5, 5, 5, 5, 4, 3, 2};
77 // static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 0, 0, 2, 3, 4, 5, 1, 1, 1, 1, 2, 3, 4, 5, 2, 2, 3, 3, 2, 3, 4, 4, 4, 2, 3, 4, 5, 5, 5, 5};
78
79 unsigned int i0, i1;
80
81 // Vertex DoFs
82 if (i == 0)
83 { i0 = 0; i1 = 0; }
84 else if (i == 1)
85 { i0 = 1; i1 = 0; }
86 else if (i == 2)
87 { i0 = 1; i1 = 1; }
88 else if (i == 3)
89 { i0 = 0; i1 = 1; }
90 // Edge DoFs
91 else if (i < totalorder + 3u)
92 { i0 = i - 2; i1 = 0; }
93 else if (i < 2u*totalorder + 2)
94 { i0 = 1; i1 = i - totalorder - 1; }
95 else if (i < 3u*totalorder + 1)
96 { i0 = i - 2u*totalorder; i1 = 1; }
97 else if (i < 4u*totalorder)
98 { i0 = 0; i1 = i - 3u*totalorder + 1; }
99 // Interior DoFs
100 else
101 {
102 unsigned int basisnum = i - 4*totalorder;
103 i0 = square_number_column[basisnum] + 2;
104 i1 = square_number_row[basisnum] + 2;
105 }
106
107 // Flip odd degree of freedom values if necessary
108 // to keep continuity on sides
109 Real f = 1.;
110
111 if ((i0%2) && (i0 > 2) && (i1 == 0))
112 f = elem->positive_edge_orientation(0)?-1.:1.;
113 else if ((i0%2) && (i0>2) && (i1 == 1))
114 f = !elem->positive_edge_orientation(2)?-1.:1.;
115 else if ((i0 == 0) && (i1%2) && (i1>2))
116 f = !elem->positive_edge_orientation(3)?-1.:1.;
117 else if ((i0 == 1) && (i1%2) && (i1>2))
118 f = elem->positive_edge_orientation(1)?-1.:1.;
119
120 return {i0, i1, f};
121}
122
123} // anonymous namespace
124
125
126
127namespace libMesh
128{
129
130
131LIBMESH_DEFAULT_VECTORIZED_FE(2,HIERARCHIC)
132LIBMESH_DEFAULT_VECTORIZED_FE(2,L2_HIERARCHIC)
133LIBMESH_DEFAULT_VECTORIZED_FE(2,SIDE_HIERARCHIC)
134
135
136template <>
137Real FE<2,HIERARCHIC>::shape(const ElemType,
138 const Order,
139 const unsigned int,
140 const Point &)
141{
142 libmesh_error_msg("Hierarchic shape functions require an Elem for edge orientation.");
143 return 0.;
144}
145
146
147
148template <>
149Real FE<2,L2_HIERARCHIC>::shape(const ElemType,
150 const Order,
151 const unsigned int,
152 const Point &)
153{
154 libmesh_error_msg("Hierarchic shape functions require an Elem for edge orientation.");
155 return 0.;
156}
157
158
159
160template <>
161Real FE<2,SIDE_HIERARCHIC>::shape(const ElemType,
162 const Order,
163 const unsigned int,
164 const Point &)
165{
166 libmesh_error_msg("Hierarchic shape functions require an Elem for edge orientation.");
167 return 0.;
168}
169
170
171
172template <>
173Real FE<2,HIERARCHIC>::shape(const Elem * elem,
174 const Order order,
175 const unsigned int i,
176 const Point & p,
177 const bool add_p_level)
178{
179 return fe_hierarchic_2D_shape<HIERARCHIC>(elem, order, i, p, add_p_level);
180}
181
182
183
184template <>
185Real FE<2,HIERARCHIC>::shape(const FEType fet,
186 const Elem * elem,
187 const unsigned int i,
188 const Point & p,
189 const bool add_p_level)
190{
191 return fe_hierarchic_2D_shape<HIERARCHIC>(elem, fet.order, i, p, add_p_level);
192}
193
194
195template <>
196Real FE<2,L2_HIERARCHIC>::shape(const Elem * elem,
197 const Order order,
198 const unsigned int i,
199 const Point & p,
200 const bool add_p_level)
201{
202 return fe_hierarchic_2D_shape<L2_HIERARCHIC>(elem, order, i, p, add_p_level);
203}
204
205
206template <>
207Real FE<2,L2_HIERARCHIC>::shape(const FEType fet,
208 const Elem * elem,
209 const unsigned int i,
210 const Point & p,
211 const bool add_p_level)
212{
213 return fe_hierarchic_2D_shape<L2_HIERARCHIC>(elem, fet.order, i, p, add_p_level);
214}
215
216
217template <>
218Real FE<2,SIDE_HIERARCHIC>::shape(const Elem * elem,
219 const Order order,
220 const unsigned int i,
221 const Point & p,
222 const bool add_p_level)
223{
224 libmesh_assert(elem);
225 const ElemType type = elem->type();
226
227 const Order totalorder = order + add_p_level*elem->p_level();
228
229 const unsigned int dofs_per_side = totalorder+1u;
230
231 switch (type)
232 {
233 case TRI6:
234 case TRI7:
235 {
236 libmesh_assert_less(i, 3*dofs_per_side);
237
238 // Flip odd degree of freedom values if necessary
239 // to keep continuity on sides. We'll flip xi/eta rather than
240 // flipping phi, so that we can use this to handle the "nodal"
241 // degrees of freedom too.
242 Real f = 1.;
243
244 const Real zeta1 = p(0);
245 const Real zeta2 = p(1);
246 const Real zeta0 = 1. - zeta1 - zeta2;
247
248 if (zeta1 > zeta2 && zeta0 > zeta2) // side 0
249 {
250 if (i >= dofs_per_side)
251 return 0;
252
253 if (totalorder == 0) // special case since raw HIERARCHIC lacks CONSTANTs
254 return 1;
255
256 if ((i < 2 || i % 2) &&
257 elem->positive_edge_orientation(0))
258 f = -1;
259
260 return FE<1,HIERARCHIC>::shape(EDGE3, totalorder, i, f*(zeta1-zeta0));
261 }
262 else if (zeta1 > zeta0 && zeta2 > zeta0) // side 1
263 {
264 if (i < dofs_per_side ||
265 i >= 2*dofs_per_side)
266 return 0;
267
268 if (totalorder == 0) // special case since raw HIERARCHIC lacks CONSTANTs
269 return 1;
270
271 const unsigned int side_i = i - dofs_per_side;
272
273 if ((side_i < 2 || side_i % 2) &&
274 elem->positive_edge_orientation(1))
275 f = -1;
276
277 return FE<1,HIERARCHIC>::shape(EDGE3, totalorder, side_i, f*(zeta2-zeta1));
278 }
279 else // side 2
280 {
281 libmesh_assert (zeta2 >= zeta1 && zeta0 >= zeta1); // On a corner???
282
283 if (i < 2*dofs_per_side)
284 return 0;
285
286 if (totalorder == 0) // special case since raw HIERARCHIC lacks CONSTANTs
287 return 1;
288
289 const unsigned int side_i = i - 2*dofs_per_side;
290
291 if ((side_i < 2 || side_i % 2) &&
292 elem->positive_edge_orientation(2))
293 f = -1;
294
295 return FE<1,HIERARCHIC>::shape(EDGE3, totalorder, side_i, f*(zeta0-zeta2));
296 }
297 }
298 case QUAD8:
299 case QUADSHELL8:
300 case QUAD9:
301 case QUADSHELL9:
302 {
303 libmesh_assert_less(i, 4*dofs_per_side);
304
305 // Flip odd degree of freedom values if necessary
306 // to keep continuity on sides. We'll flip xi/eta rather than
307 // flipping phi, so that we can use this to handle the "nodal"
308 // degrees of freedom too.
309 Real f = 1.;
310
311 const Real xi = p(0), eta = p(1);
312 if (eta < xi)
313 {
314 if (eta < -xi) // side 0
315 {
316 if (i >= dofs_per_side)
317 return 0;
318
319 if (totalorder == 0) // special case since raw HIERARCHIC lacks CONSTANTs
320 return 1;
321
322 if ((i < 2 || i % 2) &&
323 elem->positive_edge_orientation(0))
324 f = -1;
325
326 return FE<1,HIERARCHIC>::shape(EDGE3, totalorder, i, f*xi);
327 }
328 else // side 1
329 {
330 if (i < dofs_per_side ||
331 i >= 2*dofs_per_side)
332 return 0;
333
334 if (totalorder == 0) // special case since raw HIERARCHIC lacks CONSTANTs
335 return 1;
336
337 const unsigned int side_i = i - dofs_per_side;
338
339 if ((side_i < 2 || side_i % 2) &&
340 elem->positive_edge_orientation(1))
341 f = -1;
342
343 return FE<1,HIERARCHIC>::shape(EDGE3, totalorder, side_i, f*eta);
344 }
345 }
346 else // xi < eta
347 {
348 if (eta > -xi) // side 2
349 {
350 if (i < 2*dofs_per_side ||
351 i >= 3*dofs_per_side)
352 return 0;
353
354 if (totalorder == 0) // special case since raw HIERARCHIC lacks CONSTANTs
355 return 1;
356
357 const unsigned int side_i = i - 2*dofs_per_side;
358
359 if ((side_i < 2 || side_i % 2) &&
360 !elem->positive_edge_orientation(2))
361 f = -1;
362
363 return FE<1,HIERARCHIC>::shape(EDGE3, totalorder, side_i, f*xi);
364 }
365 else // side 3
366 {
367 if (i < 3*dofs_per_side)
368 return 0;
369
370 if (totalorder == 0) // special case since raw HIERARCHIC lacks CONSTANTs
371 return 1;
372
373 const unsigned int side_i = i - 3*dofs_per_side;
374
375 if ((side_i < 2 || side_i % 2) &&
376 !elem->positive_edge_orientation(3))
377 f = -1;
378
379 return FE<1,HIERARCHIC>::shape(EDGE3, totalorder, side_i, f*eta);
380 }
381 }
382 }
383 default:
384 libmesh_error_msg("ERROR: Unsupported element type = " << Utility::enum_to_string(elem->type()));
385 }
386 return 0;
387}
388
389
390template <>
391Real FE<2,SIDE_HIERARCHIC>::shape(const FEType fet,
392 const Elem * elem,
393 const unsigned int i,
394 const Point & p,
395 const bool add_p_level)
396{
397 return FE<2,SIDE_HIERARCHIC>::shape(elem, fet.order, i, p, add_p_level);
398}
399
400
401template <>
402Real FE<2,HIERARCHIC>::shape_deriv(const ElemType,
403 const Order,
404 const unsigned int,
405 const unsigned int,
406 const Point &)
407{
408 libmesh_error_msg("Hierarchic shape functions require an Elem for edge orientation.");
409 return 0.;
410}
411
412
413
414template <>
415Real FE<2,L2_HIERARCHIC>::shape_deriv(const ElemType,
416 const Order,
417 const unsigned int,
418 const unsigned int,
419 const Point &)
420{
421 libmesh_error_msg("Hierarchic shape functions require an Elem for edge orientation.");
422 return 0.;
423}
424
425
426
427template <>
428Real FE<2,SIDE_HIERARCHIC>::shape_deriv(const ElemType,
429 const Order,
430 const unsigned int,
431 const unsigned int,
432 const Point &)
433{
434 libmesh_error_msg("Hierarchic shape functions require an Elem for edge orientation.");
435 return 0.;
436}
437
438
439
440template <>
441Real FE<2,HIERARCHIC>::shape_deriv(const Elem * elem,
442 const Order order,
443 const unsigned int i,
444 const unsigned int j,
445 const Point & p,
446 const bool add_p_level)
447{
448 return fe_hierarchic_2D_shape_deriv<HIERARCHIC>(elem, order, i, j, p, add_p_level);
449}
450
451
452template <>
453Real FE<2,HIERARCHIC>::shape_deriv(const FEType fet,
454 const Elem * elem,
455 const unsigned int i,
456 const unsigned int j,
457 const Point & p,
458 const bool add_p_level)
459{
460 return fe_hierarchic_2D_shape_deriv<HIERARCHIC>(elem, fet.order, i, j, p, add_p_level);
461}
462
463
464
465
466template <>
467Real FE<2,L2_HIERARCHIC>::shape_deriv(const Elem * elem,
468 const Order order,
469 const unsigned int i,
470 const unsigned int j,
471 const Point & p,
472 const bool add_p_level)
473{
474 return fe_hierarchic_2D_shape_deriv<L2_HIERARCHIC>(elem, order, i, j, p, add_p_level);
475}
476
477
478template <>
479Real FE<2,L2_HIERARCHIC>::shape_deriv(const FEType fet,
480 const Elem * elem,
481 const unsigned int i,
482 const unsigned int j,
483 const Point & p,
484 const bool add_p_level)
485{
486 return fe_hierarchic_2D_shape_deriv<L2_HIERARCHIC>(elem, fet.order, i, j, p, add_p_level);
487}
488
489
490
491template <>
492Real FE<2,SIDE_HIERARCHIC>::shape_deriv(const Elem * elem,
493 const Order order,
494 const unsigned int i,
495 const unsigned int j,
496 const Point & p,
497 const bool add_p_level)
498{
499 libmesh_assert(elem);
500
501 const ElemType type = elem->type();
502
503 const Order totalorder = order + add_p_level*elem->p_level();
504
505 if (totalorder == 0) // special case since raw HIERARCHIC lacks CONSTANTs
506 return 0;
507
508 const unsigned int dofs_per_side = totalorder+1u;
509
510 switch (type)
511 {
512 case TRI6:
513 case TRI7:
514 {
515 return fe_fdm_deriv(elem, order, i, j, p, add_p_level, FE<2,SIDE_HIERARCHIC>::shape);
516 }
517#if 0
518 {
519 libmesh_assert_less(i, 3*dofs_per_side);
520 libmesh_assert_less (j, 2);
521
522 // Flip odd degree of freedom values if necessary
523 // to keep continuity on sides. We'll flip xi/eta rather than
524 // flipping phi, so that we can use this to handle the "nodal"
525 // degrees of freedom too.
526 Real f = 1.;
527
528 const Real zeta1 = p(0);
529 const Real zeta2 = p(1);
530 const Real zeta0 = 1. - zeta1 - zeta2;
531
532 if (zeta1 > zeta2 && zeta0 > zeta2) // side 0
533 {
534 if (j == 1) // d/deta is perpendicular here
535 return 0;
536
537 if (i >= dofs_per_side)
538 return 0;
539
540 if (totalorder == 0) // special case since raw HIERARCHIC lacks CONSTANTs
541 return 0;
542
543 if ((i < 2 || i % 2) &&
544 elem->positive_edge_orientation(0))
545 f = -1;
546
547 return f*FE<1,HIERARCHIC>::shape_deriv(EDGE3, totalorder, i, 0, f*(zeta1-zeta0));
548 }
549 else if (zeta1 > zeta0 && zeta2 > zeta0) // side 1
550 {
551 if (i < dofs_per_side ||
552 i >= 2*dofs_per_side)
553 return 0;
554
555 if (totalorder == 0) // special case since raw HIERARCHIC lacks CONSTANTs
556 return 0;
557
558 const unsigned int side_i = i - dofs_per_side;
559
560 if ((side_i < 2 || side_i % 2) &&
561 elem->positive_edge_orientation(1))
562 f = -1;
563
564 Real g = 1;
565 if (j == 0) // 2D d/dxi is in the opposite direction on this edge
566 g = -1;
567
568 return f*g*FE<1,HIERARCHIC>::shape_deriv(EDGE3, totalorder, side_i, 0, f*(zeta2-zeta1));
569 }
570 else // side 2
571 {
572 libmesh_assert (zeta2 >= zeta1 && zeta0 >= zeta1); // On a corner???
573
574 if (j == 0) // d/dxi is perpendicular here
575 return 0;
576
577 if (i < 2*dofs_per_side)
578 return 0;
579
580 if (totalorder == 0) // special case since raw HIERARCHIC lacks CONSTANTs
581 return 0;
582
583 const unsigned int side_i = i - 2*dofs_per_side;
584
585 if ((side_i < 2 || side_i % 2) &&
586 elem->positive_edge_orientation(2))
587 f = -1;
588
589 return -f*FE<1,HIERARCHIC>::shape_deriv(EDGE3, totalorder, side_i, 0, f*(zeta0-zeta2));
590 }
591 }
592#endif
593 case QUAD8:
594 case QUADSHELL8:
595 case QUAD9:
596 case QUADSHELL9:
597 {
598 libmesh_assert_less(i, 4*dofs_per_side);
599
600 // Flip odd degree of freedom values if necessary
601 // to keep continuity on sides. We'll flip xi/eta rather than
602 // flipping phi, so that we can use this to handle the "nodal"
603 // degrees of freedom too.
604 Real f = 1.;
605
606 const Real xi = p(0), eta = p(1);
607 if (eta < xi)
608 {
609 if (eta < -xi) // side 0
610 {
611 if (i >= dofs_per_side)
612 return 0;
613 if (j != 0)
614 return 0;
615 if ((i < 2 || i % 2) &&
616 elem->positive_edge_orientation(0))
617 f = -1;
618
619 return f*FE<1,HIERARCHIC>::shape_deriv(EDGE3, totalorder, i, 0, f*xi);
620 }
621 else // side 1
622 {
623 if (i < dofs_per_side ||
624 i >= 2*dofs_per_side)
625 return 0;
626 if (j != 1)
627 return 0;
628
629 const unsigned int side_i = i - dofs_per_side;
630
631 if ((side_i < 2 || side_i % 2) &&
632 elem->positive_edge_orientation(1))
633 f = -1;
634
635 return f*FE<1,HIERARCHIC>::shape_deriv(EDGE3, totalorder, side_i, 0, f*eta);
636 }
637 }
638 else // xi < eta
639 {
640 if (eta > -xi) // side 2
641 {
642 if (i < 2*dofs_per_side ||
643 i >= 3*dofs_per_side)
644 return 0;
645 if (j != 0)
646 return 0;
647
648 const unsigned int side_i = i - 2*dofs_per_side;
649
650 if ((side_i < 2 || side_i % 2) &&
651 !elem->positive_edge_orientation(2))
652 f = -1;
653
654 return f*FE<1,HIERARCHIC>::shape_deriv(EDGE3, totalorder, side_i, 0, f*xi);
655 }
656 else // side 3
657 {
658 if (i < 3*dofs_per_side)
659 return 0;
660 if (j != 1)
661 return 0;
662
663 const unsigned int side_i = i - 3*dofs_per_side;
664
665 if ((side_i < 2 || side_i % 2) &&
666 !elem->positive_edge_orientation(3))
667 f = -1;
668
669 return f*FE<1,HIERARCHIC>::shape_deriv(EDGE3, totalorder, side_i, 0, f*eta);
670 }
671 }
672 }
673 default:
674 libmesh_error_msg("ERROR: Unsupported element type = " << Utility::enum_to_string(elem->type()));
675 }
676 return 0;
677}
678
679
680template <>
681Real FE<2,SIDE_HIERARCHIC>::shape_deriv(const FEType fet,
682 const Elem * elem,
683 const unsigned int i,
684 const unsigned int j,
685 const Point & p,
686 const bool add_p_level)
687{
688 return FE<2,SIDE_HIERARCHIC>::shape_deriv(elem, fet.order, i, j, p, add_p_level);
689}
690
691
692#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
693
694template <>
695Real FE<2,HIERARCHIC>::shape_second_deriv(const ElemType,
696 const Order,
697 const unsigned int,
698 const unsigned int,
699 const Point &)
700{
701 libmesh_error_msg("Hierarchic shape functions require an Elem for edge orientation.");
702 return 0.;
703}
704
705
706
707template <>
708Real FE<2,L2_HIERARCHIC>::shape_second_deriv(const ElemType,
709 const Order,
710 const unsigned int,
711 const unsigned int,
712 const Point &)
713{
714 libmesh_error_msg("Hierarchic shape functions require an Elem for edge orientation.");
715 return 0.;
716}
717
718
719
720template <>
721Real FE<2,SIDE_HIERARCHIC>::shape_second_deriv(const ElemType,
722 const Order,
723 const unsigned int,
724 const unsigned int,
725 const Point &)
726{
727 libmesh_error_msg("Hierarchic shape functions require an Elem for edge orientation.");
728 return 0.;
729}
730
731
732
733template <>
734Real FE<2,HIERARCHIC>::shape_second_deriv(const Elem * elem,
735 const Order order,
736 const unsigned int i,
737 const unsigned int j,
738 const Point & p,
739 const bool add_p_level)
740{
741 return fe_hierarchic_2D_shape_second_deriv<HIERARCHIC>(elem, order, i, j, p, add_p_level);
742}
743
744
745template <>
746Real FE<2,HIERARCHIC>::shape_second_deriv(const FEType fet,
747 const Elem * elem,
748 const unsigned int i,
749 const unsigned int j,
750 const Point & p,
751 const bool add_p_level)
752{
753 return fe_hierarchic_2D_shape_second_deriv<HIERARCHIC>(elem, fet.order, i, j, p, add_p_level);
754}
755
756
757template <>
758Real FE<2,L2_HIERARCHIC>::shape_second_deriv(const Elem * elem,
759 const Order order,
760 const unsigned int i,
761 const unsigned int j,
762 const Point & p,
763 const bool add_p_level)
764{
765 return fe_hierarchic_2D_shape_second_deriv<L2_HIERARCHIC>(elem, order, i, j, p, add_p_level);
766}
767
768
769template <>
770Real FE<2,L2_HIERARCHIC>::shape_second_deriv(const FEType fet,
771 const Elem * elem,
772 const unsigned int i,
773 const unsigned int j,
774 const Point & p,
775 const bool add_p_level)
776{
777 return fe_hierarchic_2D_shape_second_deriv<L2_HIERARCHIC>(elem, fet.order, i, j, p, add_p_level);
778}
779
780
781template <>
782Real FE<2,SIDE_HIERARCHIC>::shape_second_deriv(const Elem * elem,
783 const Order order,
784 const unsigned int i,
785 const unsigned int j,
786 const Point & p,
787 const bool add_p_level)
788{
789 libmesh_assert(elem);
790 const ElemType type = elem->type();
791
792 const Order totalorder = order + add_p_level*elem->p_level();
793
794 if (totalorder == 0) // special case since raw HIERARCHIC lacks CONSTANTs
795 return 0;
796
797 const unsigned int dofs_per_side = totalorder+1u;
798
799 switch (type)
800 {
801 case TRI6:
802 case TRI7:
803 {
804 return fe_fdm_second_deriv(elem, order, i, j, p, add_p_level,
805 FE<2,SIDE_HIERARCHIC>::shape_deriv);
806 }
807 case QUAD8:
808 case QUADSHELL8:
809 case QUAD9:
810 case QUADSHELL9:
811 {
812 libmesh_assert_less(i, 4*dofs_per_side);
813
814 // Flip odd degree of freedom values if necessary
815 // to keep continuity on sides. We'll flip xi/eta rather than
816 // flipping phi, so that we can use this to handle the "nodal"
817 // degrees of freedom too.
818 Real f = 1.;
819
820 const Real xi = p(0), eta = p(1);
821 if (eta < xi)
822 {
823 if (eta < -xi) // side 0
824 {
825 if (i >= dofs_per_side)
826 return 0;
827 if (j != 0)
828 return 0;
829 if ((i < 2 || i % 2) &&
830 elem->positive_edge_orientation(0))
831 f = -1;
832
833 return FE<1,HIERARCHIC>::shape_second_deriv(EDGE3, totalorder, i, 0, f*xi);
834 }
835 else // side 1
836 {
837 if (i < dofs_per_side ||
838 i >= 2*dofs_per_side)
839 return 0;
840 if (j != 2)
841 return 0;
842
843 const unsigned int side_i = i - dofs_per_side;
844
845 if ((side_i < 2 || side_i % 2) &&
846 elem->positive_edge_orientation(1))
847 f = -1;
848
849 return FE<1,HIERARCHIC>::shape_second_deriv(EDGE3, totalorder, side_i, 0, f*eta);
850 }
851 }
852 else // xi < eta
853 {
854 if (eta > -xi) // side 2
855 {
856 if (i < 2*dofs_per_side ||
857 i >= 3*dofs_per_side)
858 return 0;
859 if (j != 0)
860 return 0;
861
862 const unsigned int side_i = i - 2*dofs_per_side;
863
864 if ((side_i < 2 || side_i % 2) &&
865 !elem->positive_edge_orientation(2))
866 f = -1;
867
868 return FE<1,HIERARCHIC>::shape_second_deriv(EDGE3, totalorder, side_i, 0, f*xi);
869 }
870 else // side 3
871 {
872 if (i < 3*dofs_per_side)
873 return 0;
874 if (j != 2)
875 return 0;
876
877 const unsigned int side_i = i - 3*dofs_per_side;
878
879 if ((side_i < 2 || side_i % 2) &&
880 !elem->positive_edge_orientation(3))
881 f = -1;
882
883 return FE<1,HIERARCHIC>::shape_second_deriv(EDGE3, totalorder, side_i, 0, f*eta);
884 }
885 }
886 }
887 default:
888 libmesh_error_msg("ERROR: Unsupported element type = " << Utility::enum_to_string(elem->type()));
889 }
890 return 0;
891}
892
893
894template <>
895Real FE<2,SIDE_HIERARCHIC>::shape_second_deriv(const FEType fet,
896 const Elem * elem,
897 const unsigned int i,
898 const unsigned int j,
899 const Point & p,
900 const bool add_p_level)
901{
902 return FE<2,SIDE_HIERARCHIC>::shape_second_deriv(elem, fet.order, i, j, p, add_p_level);
903}
904
905#endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
906
907} // namespace libMesh
908
909
910
911namespace
912{
913using namespace libMesh;
914
915Real fe_triangle_helper (const Elem & elem,
916 const Real edgenumerator,
917 const Real crossval,
918 const unsigned int basisorder,
919 const Order totalorder,
920 const unsigned int noden)
921{
922 // Get factors to account for edge-flipping
923 Real flip = 1;
924 if (basisorder%2 && (elem.point(noden) > elem.point((noden+1)%3)))
925 flip = -1.;
926
927 // Avoid NaN around vertices ... but we still have to match the true
928 // function, even when we're *outside* the triangle (crossval==0 on
929 // a line, not just at a point!), to handle imprecise queries and
930 // FDM derivatives correctly!
931 if (crossval == 0.)
932 {
933 unsigned int basisfactorial = 1.;
934 for (unsigned int n=2; n <= basisorder; ++n)
935 basisfactorial *= n;
936
937 return std::pow(edgenumerator, basisorder) / basisfactorial;
938 }
939 // Experimentally, as c -> 0, n propto c, I'm still seeing good
940 // behavior from the default implementation below:
941
942 const Real edgeval = edgenumerator / crossval;
943 const Real crossfunc = std::pow(crossval, basisorder);
944
945 return flip * crossfunc *
946 FE<1,HIERARCHIC>::shape(EDGE3, totalorder,
947 basisorder, edgeval);
948}
949
950template <FEFamily T>
951Real fe_hierarchic_2D_shape(const Elem * elem,
952 const Order order,
953 const unsigned int i,
954 const Point & p,
955 const bool add_p_level)
956{
957 libmesh_assert(elem);
958
959 const Order totalorder = order + add_p_level*elem->p_level();
960 libmesh_assert_greater (totalorder, 0);
961
962 switch (elem->type())
963 {
964 case TRI3:
965 case TRISHELL3:
966 case TRI6:
967 case TRI7:
968 {
969 const Real zeta1 = p(0);
970 const Real zeta2 = p(1);
971 const Real zeta0 = 1. - zeta1 - zeta2;
972
973 libmesh_assert_less (i, (totalorder+1u)*(totalorder+2u)/2);
974 libmesh_assert (T == L2_HIERARCHIC || elem->type() == TRI6 ||
975 elem->type() == TRI7 || totalorder < 2);
976
977 // Vertex DoFs
978 if (i == 0)
979 return zeta0;
980 else if (i == 1)
981 return zeta1;
982 else if (i == 2)
983 return zeta2;
984 // Edge DoFs
985 else if (i < totalorder + 2u)
986 {
987 const unsigned int basisorder = i - 1;
988
989 const Real crossval = zeta0 + zeta1;
990 const Real edgenumerator = zeta1 - zeta0;
991
992 return fe_triangle_helper(*elem, edgenumerator, crossval,
993 basisorder, totalorder, 0);
994 }
995 else if (i < 2u*totalorder + 1)
996 {
997 const unsigned int basisorder = i - totalorder;
998
999 const Real crossval = zeta2 + zeta1;
1000 const Real edgenumerator = zeta2 - zeta1;
1001
1002 return fe_triangle_helper(*elem, edgenumerator, crossval,
1003 basisorder, totalorder, 1);
1004 }
1005 else if (i < 3u*totalorder)
1006 {
1007 const unsigned int basisorder = i - (2u*totalorder) + 1;
1008
1009 const Real crossval = zeta0 + zeta2;
1010 const Real edgenumerator = zeta0 - zeta2;
1011
1012 return fe_triangle_helper(*elem, edgenumerator, crossval,
1013 basisorder, totalorder, 2);
1014 }
1015 // Interior DoFs
1016 else
1017 {
1018 const unsigned int basisnum = i - (3u*totalorder);
1019 unsigned int exp0 = triangular_number_column[basisnum] + 1;
1020 unsigned int exp1 = triangular_number_row[basisnum] + 1 -
1021 triangular_number_column[basisnum];
1022
1023 Real returnval = 1;
1024 for (unsigned int n = 0; n != exp0; ++n)
1025 returnval *= zeta0;
1026 for (unsigned int n = 0; n != exp1; ++n)
1027 returnval *= zeta1;
1028 returnval *= zeta2;
1029 return returnval;
1030 }
1031 }
1032
1033 // Hierarchic shape functions on the quadrilateral.
1034 case QUAD4:
1035 case QUADSHELL4:
1036 libmesh_assert (T == L2_HIERARCHIC || totalorder < 2);
1037 libmesh_fallthrough();
1038 case QUAD8:
1039 case QUADSHELL8:
1040 case QUAD9:
1041 case QUADSHELL9:
1042 {
1043 // Compute quad shape functions as a tensor-product
1044 auto [i0, i1, f] = quad_indices(elem, totalorder, i);
1045
1046 return f*(FE<1,T>::shape(EDGE3, totalorder, i0, p(0))*
1047 FE<1,T>::shape(EDGE3, totalorder, i1, p(1)));
1048 }
1049
1050 default:
1051 libmesh_error_msg("ERROR: Unsupported element type = " << Utility::enum_to_string(elem->type()));
1052 }
1053
1054 return 0.;
1055}
1056
1057
1058template <FEFamily T>
1059Real fe_hierarchic_2D_shape_deriv(const Elem * elem,
1060 const Order order,
1061 const unsigned int i,
1062 const unsigned int j,
1063 const Point & p,
1064 const bool add_p_level)
1065{
1066 libmesh_assert(elem);
1067
1068 const ElemType type = elem->type();
1069
1070 const Order totalorder = order + add_p_level*elem->p_level();
1071
1072 libmesh_assert_greater (totalorder, 0);
1073
1074 switch (type)
1075 {
1076 // 1st & 2nd-order Hierarchics.
1077 case TRI3:
1078 case TRISHELL3:
1079 case TRI6:
1080 case TRI7:
1081 {
1082 return fe_fdm_deriv(elem, order, i, j, p, add_p_level, FE<2,T>::shape);
1083 }
1084
1085 case QUAD4:
1086 case QUADSHELL4:
1087 libmesh_assert (T == L2_HIERARCHIC || totalorder < 2);
1088 libmesh_fallthrough();
1089 case QUAD8:
1090 case QUADSHELL8:
1091 case QUAD9:
1092 case QUADSHELL9:
1093 {
1094 // Compute quad shape functions as a tensor-product
1095 auto [i0, i1, f] = quad_indices(elem, totalorder, i);
1096
1097 switch (j)
1098 {
1099 // d()/dxi
1100 case 0:
1101 return f*(FE<1,T>::shape_deriv(EDGE3, totalorder, i0, 0, p(0))*
1102 FE<1,T>::shape (EDGE3, totalorder, i1, p(1)));
1103
1104 // d()/deta
1105 case 1:
1106 return f*(FE<1,T>::shape (EDGE3, totalorder, i0, p(0))*
1107 FE<1,T>::shape_deriv(EDGE3, totalorder, i1, 0, p(1)));
1108
1109 default:
1110 libmesh_error_msg("Invalid derivative index j = " << j);
1111 }
1112 }
1113
1114 default:
1115 libmesh_error_msg("ERROR: Unsupported element type = " << Utility::enum_to_string(type));
1116 }
1117
1118 return 0.;
1119}
1120
1121
1122
1123#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1124
1125template <FEFamily T>
1126Real fe_hierarchic_2D_shape_second_deriv(const Elem * elem,
1127 const Order order,
1128 const unsigned int i,
1129 const unsigned int j,
1130 const Point & p,
1131 const bool add_p_level)
1132{
1133 return fe_fdm_second_deriv(elem, order, i, j, p, add_p_level,
1134 FE<2,T>::shape_deriv);
1135}
1136
1137#endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
1138
1139} // anonymous namespace
libMesh::Elem
This is the base class from which all geometric element types are derived.
Definition elem.h:96
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::type
virtual ElemType type() const =0
libMesh::Elem::p_level
unsigned int p_level() const
Definition elem.h:3125
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::Point
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition point.h:40
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::LIBMESH_DEFAULT_VECTORIZED_FE
LIBMESH_DEFAULT_VECTORIZED_FE(template<>Real FE< 0, BERNSTEIN)
Definition fe_bernstein_shape_0D.C:30
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::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::triangular_number_row
const unsigned char triangular_number_row[]
Definition number_lookups.C:28
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::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