libMesh
Loading...
Searching...
No Matches
fe_szabab.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
20// Local includes
21#include "libmesh/libmesh_config.h"
22#ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
23
24#include "libmesh/elem.h"
25#include "libmesh/enum_to_string.h"
26#include "libmesh/fe.h"
27#include "libmesh/fe_interface.h"
28#include "libmesh/fe_macro.h"
29
30namespace libMesh
31{
32
33// ------------------------------------------------------------
34// Szabo-Babuska-specific implementations, Steffen Petersen 2004
35
36// Anonymous namespace for local helper functions
37namespace {
38
39void szabab_nodal_soln(const Elem * elem,
40 const Order order,
41 const std::vector<Number> & elem_soln,
42 std::vector<Number> & nodal_soln,
43 const bool add_p_level)
44{
45 const unsigned int n_nodes = elem->n_nodes();
46
47 const ElemType elem_type = elem->type();
48
49 nodal_soln.resize(n_nodes);
50
51 const Order totalorder = order + add_p_level*elem->p_level();
52
53 // FEType object to be passed to various FEInterface functions below.
54 FEType fe_type(order, SZABAB);
55
56 switch (totalorder)
57 {
58 // Constant shape functions
59 case CONSTANT:
60 {
61 libmesh_assert_equal_to (elem_soln.size(), 1);
62
63 std::fill(nodal_soln.begin(), nodal_soln.end(), elem_soln[0]);
64
65 return;
66 }
67
68
69 // For other bases do interpolation at the nodes
70 // explicitly.
71 case FIRST:
72 case SECOND:
73 case THIRD:
74 case FOURTH:
75 case FIFTH:
76 case SIXTH:
77 case SEVENTH:
78 {
79 const unsigned int n_sf =
80 FEInterface::n_shape_functions(fe_type, elem);
81
82 std::vector<Point> refspace_nodes;
83 FEBase::get_refspace_nodes(elem_type,refspace_nodes);
84 libmesh_assert_equal_to (refspace_nodes.size(), n_nodes);
85 libmesh_assert_equal_to (elem_soln.size(), n_sf);
86
87 // Zero before summation
88 std::fill(nodal_soln.begin(), nodal_soln.end(), 0);
89
90 for (unsigned int n=0; n<n_nodes; n++)
91 // u_i = Sum (alpha_i phi_i)
92 for (unsigned int i=0; i<n_sf; i++)
93 nodal_soln[n] += elem_soln[i] *
94 FEInterface::shape(fe_type, elem, i, refspace_nodes[n]);
95
96 return;
97 }
98
99 default:
100 libmesh_error_msg("ERROR: Invalid total order " << totalorder);
101 }
102} // szabab_nodal_soln()
103
104
105
106unsigned int szabab_n_dofs(const ElemType t, const Order o)
107{
108 switch (o)
109 {
110 // Szabo-Babuska 1st-order polynomials.
111 case FIRST:
112 {
113 switch (t)
114 {
115 case NODEELEM:
116 return 1;
117
118 case EDGE2:
119 case EDGE3:
120 return 2;
121
122 case TRI3:
123 case TRI6:
124 case TRI7:
125 return 3;
126
127 case QUAD4:
128 case QUAD8:
129 case QUAD9:
130 return 4;
131
132 case INVALID_ELEM:
133 return 0;
134
135 default:
136 libmesh_error_msg("ERROR: Invalid ElemType " << Utility::enum_to_string(t) << " selected for SZABAB FE family!");
137 }
138 }
139
140
141 // Szabo-Babuska 2nd-order polynomials.
142 case SECOND:
143 {
144 switch (t)
145 {
146 case NODEELEM:
147 return 1;
148
149 case EDGE2:
150 case EDGE3:
151 return 3;
152
153 case TRI6:
154 case TRI7:
155 return 6;
156
157 case QUAD8:
158 return 8;
159 case QUAD9:
160 return 9;
161
162 case INVALID_ELEM:
163 return 0;
164
165 default:
166 libmesh_error_msg("ERROR: Invalid ElemType " << Utility::enum_to_string(t) << " selected for SZABAB FE family!");
167 }
168 }
169
170
171 // Szabo-Babuska 3rd-order polynomials.
172 case THIRD:
173 {
174 switch (t)
175 {
176 case NODEELEM:
177 return 1;
178
179 case EDGE2:
180 case EDGE3:
181 return 4;
182
183 case TRI6:
184 case TRI7:
185 return 10;
186
187 case QUAD8:
188 case QUAD9:
189 return 16;
190
191 case INVALID_ELEM:
192 return 0;
193
194 default:
195 libmesh_error_msg("ERROR: Invalid ElemType " << Utility::enum_to_string(t) << " selected for SZABAB FE family!");
196 }
197 }
198
199
200 // Szabo-Babuska 4th-order polynomials.
201 case FOURTH:
202 {
203 switch (t)
204 {
205 case NODEELEM:
206 return 1;
207
208 case EDGE2:
209 case EDGE3:
210 return 5;
211
212 case TRI6:
213 case TRI7:
214 return 15;
215
216 case QUAD8:
217 case QUAD9:
218 return 25;
219
220 case INVALID_ELEM:
221 return 0;
222
223 default:
224 libmesh_error_msg("ERROR: Invalid ElemType " << Utility::enum_to_string(t) << " selected for SZABAB FE family!");
225 }
226 }
227
228
229 // Szabo-Babuska 5th-order polynomials.
230 case FIFTH:
231 {
232 switch (t)
233 {
234 case NODEELEM:
235 return 1;
236
237 case EDGE2:
238 case EDGE3:
239 return 6;
240
241 case TRI6:
242 case TRI7:
243 return 21;
244
245 case QUAD8:
246 case QUAD9:
247 return 36;
248
249 case INVALID_ELEM:
250 return 0;
251
252 default:
253 libmesh_error_msg("ERROR: Invalid ElemType " << Utility::enum_to_string(t) << " selected for SZABAB FE family!");
254 }
255 }
256
257
258 // Szabo-Babuska 6th-order polynomials.
259 case SIXTH:
260 {
261 switch (t)
262 {
263 case NODEELEM:
264 return 1;
265
266 case EDGE2:
267 case EDGE3:
268 return 7;
269
270 case TRI6:
271 case TRI7:
272 return 28;
273
274 case QUAD8:
275 case QUAD9:
276 return 49;
277
278 case INVALID_ELEM:
279 return 0;
280
281 default:
282 libmesh_error_msg("ERROR: Invalid ElemType " << Utility::enum_to_string(t) << " selected for SZABAB FE family!");
283 }
284 }
285
286 // Szabo-Babuska 7th-order polynomials.
287 case SEVENTH:
288 {
289 switch (t)
290 {
291 case NODEELEM:
292 return 1;
293
294 case EDGE2:
295 case EDGE3:
296 return 8;
297
298 case TRI6:
299 case TRI7:
300 return 36;
301
302 case QUAD8:
303 case QUAD9:
304 return 64;
305
306 case INVALID_ELEM:
307 return 0;
308
309 default:
310 libmesh_error_msg("ERROR: Invalid ElemType " << Utility::enum_to_string(t) << " selected for SZABAB FE family!");
311 }
312 }
313
314
315 default:
316 libmesh_error_msg("ERROR: Invalid Order " << Utility::enum_to_string(o) << " selected for SZABAB FE family!");
317 }
318} // szabab_n_dofs()
319
320
321
322unsigned int szabab_n_dofs(const Elem * e, const Order o)
323{
324 libmesh_assert(e);
325 return szabab_n_dofs(e->type(), o);
326}
327
328
329
330unsigned int szabab_n_dofs_at_node(const ElemType t,
331 const Order o,
332 const unsigned int n)
333{
334 switch (o)
335 {
336 // The first-order Szabo-Babuska shape functions
337 case FIRST:
338 {
339 switch (t)
340 {
341 case NODEELEM:
342 return 1;
343
344 // The 1D Szabo-Babuska defined on a three-noded edge
345 case EDGE2:
346 case EDGE3:
347 {
348 switch (n)
349 {
350 case 0:
351 case 1:
352 return 1;
353
354 case 2:
355 return 0;
356
357 default:
358 libmesh_error_msg("ERROR: Invalid node ID " << n << " selected for EDGE2/3!");
359 }
360 }
361
362
363 // The 2D Szabo-Babuska defined on a 6-noded triangle
364 case TRI3:
365 case TRI6:
366 case TRI7:
367 {
368 switch (n)
369 {
370 case 0:
371 case 1:
372 case 2:
373 return 1;
374
375 case 3:
376 case 4:
377 case 5:
378 return 0;
379
380 case 6:
381 libmesh_assert(t == TRI7);
382 return 0;
383
384 default:
385 libmesh_error_msg("ERROR: Invalid node ID " << n << " selected for TRI3/6/7!");
386 }
387 }
388
389
390 // The 2D tensor-product Szabo-Babuska defined on a
391 // nine-noded quadrilateral.
392 case QUAD4:
393 case QUAD8:
394 case QUAD9:
395 {
396 switch (n)
397 {
398 case 0:
399 case 1:
400 case 2:
401 case 3:
402 return 1;
403
404 case 4:
405 case 5:
406 case 6:
407 case 7:
408 return 0;
409
410 case 8:
411 return 0;
412
413 default:
414 libmesh_error_msg("ERROR: Invalid node ID " << n << " selected for QUAD4/8/9!");
415 }
416 }
417
418 case INVALID_ELEM:
419 return 0;
420
421 default:
422 libmesh_error_msg("ERROR: Invalid ElemType " << Utility::enum_to_string(t) << " selected for SZABAB FE family!");
423 }
424 }
425
426
427
428 // The second-order Szabo-Babuska shape functions
429 case SECOND:
430 {
431 switch (t)
432 {
433 case NODEELEM:
434 return 1;
435
436 // The 1D Szabo-Babuska defined on a three-noded edge
437 case EDGE2:
438 case EDGE3:
439 {
440 switch (n)
441 {
442 case 0:
443 case 1:
444 return 1;
445
446 case 2:
447 return 1;
448
449 default:
450 libmesh_error_msg("ERROR: Invalid node ID " << n << " selected for EDGE2/3!");
451 }
452 }
453
454
455 // The 2D Szabo-Babuska defined on a 6-noded triangle
456 case TRI6:
457 case TRI7:
458 {
459 switch (n)
460 {
461 case 0:
462 case 1:
463 case 2:
464 return 1;
465
466 case 3:
467 case 4:
468 case 5:
469 return 1;
470
471 case 6:
472 libmesh_assert(t == TRI7);
473 return 0;
474
475 default:
476 libmesh_error_msg("ERROR: Invalid node ID " << n << " selected for TRI6/7!");
477 }
478 }
479
480
481 // The 2D tensor-product Szabo-Babuska defined on a
482 // nine-noded quadrilateral.
483 case QUAD8:
484 case QUAD9:
485 {
486 switch (n)
487 {
488 case 0:
489 case 1:
490 case 2:
491 case 3:
492 return 1;
493
494 case 4:
495 case 5:
496 case 6:
497 case 7:
498 return 1;
499
500 case 8:
501 return 1;
502
503 default:
504 libmesh_error_msg("ERROR: Invalid node ID " << n << " selected for QUAD8/9!");
505 }
506 }
507
508 case INVALID_ELEM:
509 return 0;
510
511 default:
512 libmesh_error_msg("ERROR: Invalid ElemType " << Utility::enum_to_string(t) << " selected for SZABAB FE family!");
513 }
514 }
515
516
517
518 // The third-order Szabo-Babuska shape functions
519 case THIRD:
520 {
521 switch (t)
522 {
523 case NODEELEM:
524 return 1;
525
526 // The 1D Szabo-Babuska defined on a three-noded edge
527 case EDGE2:
528 case EDGE3:
529 {
530 switch (n)
531 {
532 case 0:
533 case 1:
534 return 1;
535
536 case 2:
537 return 2;
538
539 default:
540 libmesh_error_msg("ERROR: Invalid node ID " << n << " selected for EDGE2/3!");
541 }
542 }
543
544
545 // The 2D Szabo-Babuska defined on a 6-noded triangle
546 case TRI7:
547 if (n == 6)
548 return 1;
549 libmesh_fallthrough();
550 case TRI6:
551 {
552 switch (n)
553 {
554 case 0:
555 case 1:
556 case 2:
557 return 1;
558
559 case 3:
560 case 4:
561 case 5:
562 return 2;
563
564 default:
565 libmesh_error_msg("ERROR: Invalid node ID " << n << " selected for TRI6/7!");
566 }
567 }
568
569
570 // The 2D tensor-product Szabo-Babuska defined on a
571 // nine-noded quadrilateral.
572 case QUAD8:
573 case QUAD9:
574 {
575 switch (n)
576 {
577 case 0:
578 case 1:
579 case 2:
580 case 3:
581 return 1;
582
583 case 4:
584 case 5:
585 case 6:
586 case 7:
587 return 2;
588
589 case 8:
590 return 4;
591
592 default:
593 libmesh_error_msg("ERROR: Invalid node ID " << n << " selected for QUAD8/9!");
594 }
595 }
596
597 case INVALID_ELEM:
598 return 0;
599
600 default:
601 libmesh_error_msg("ERROR: Invalid ElemType " << Utility::enum_to_string(t) << " selected for SZABAB FE family!");
602 }
603 }
604
605
606
607 // The fourth-order Szabo-Babuska shape functions
608 case FOURTH:
609 {
610 switch (t)
611 {
612 case NODEELEM:
613 return 1;
614
615 // The 1D Szabo-Babuska defined on a three-noded edge
616 case EDGE2:
617 case EDGE3:
618 {
619 switch (n)
620 {
621 case 0:
622 case 1:
623 return 1;
624
625 case 2:
626 return 3;
627
628 default:
629 libmesh_error_msg("ERROR: Invalid node ID " << n << " selected for EDGE2/3!");
630 }
631 }
632
633
634 // The 2D Szabo-Babuska defined on a 6/7-noded triangle
635 case TRI7:
636 if (n == 6)
637 return 3;
638 libmesh_fallthrough();
639 case TRI6:
640 {
641 switch (n)
642 {
643 case 0:
644 case 1:
645 case 2:
646 return 1;
647
648 case 3:
649 case 4:
650 case 5:
651 return 3;
652
653 default:
654 libmesh_error_msg("ERROR: Invalid node ID " << n << " selected for TRI6!");
655 }
656 }
657
658
659 // The 2D tensor-product Szabo-Babuska defined on a
660 // nine-noded quadrilateral.
661 case QUAD8:
662 case QUAD9:
663 {
664 switch (n)
665 {
666 case 0:
667 case 1:
668 case 2:
669 case 3:
670 return 1;
671
672 case 4:
673 case 5:
674 case 6:
675 case 7:
676 return 3;
677
678 case 8:
679 return 9;
680
681 default:
682 libmesh_error_msg("ERROR: Invalid node ID " << n << " selected for QUAD8/9!");
683 }
684 }
685
686 case INVALID_ELEM:
687 return 0;
688
689 default:
690 libmesh_error_msg("ERROR: Invalid ElemType " << Utility::enum_to_string(t) << " selected for SZABAB FE family!");
691 }
692 }
693
694
695
696 // The fifth-order Szabo-Babuska shape functions
697 case FIFTH:
698 {
699 switch (t)
700 {
701 case NODEELEM:
702 return 1;
703
704 // The 1D Szabo-Babuska defined on a three-noded edge
705 case EDGE2:
706 case EDGE3:
707 {
708 switch (n)
709 {
710 case 0:
711 case 1:
712 return 1;
713
714 case 2:
715 return 4;
716
717 default:
718 libmesh_error_msg("ERROR: Invalid node ID " << n << " selected for EDGE2/3!");
719 }
720 }
721
722
723 // The 2D Szabo-Babuska defined on a 6/7-noded triangle
724 case TRI7:
725 if (n == 6)
726 return 6;
727 libmesh_fallthrough();
728 case TRI6:
729 {
730 switch (n)
731 {
732 case 0:
733 case 1:
734 case 2:
735 return 1;
736
737 case 3:
738 case 4:
739 case 5:
740 return 4;
741
742 default:
743 libmesh_error_msg("ERROR: Invalid node ID " << n << " selected for TRI6!");
744 }
745 }
746
747
748 // The 2D tensor-product Szabo-Babuska defined on a
749 // nine-noded quadrilateral.
750 case QUAD8:
751 case QUAD9:
752 {
753 switch (n)
754 {
755 case 0:
756 case 1:
757 case 2:
758 case 3:
759 return 1;
760
761 case 4:
762 case 5:
763 case 6:
764 case 7:
765 return 4;
766
767 case 8:
768 return 16;
769
770 default:
771 libmesh_error_msg("ERROR: Invalid node ID " << n << " selected for QUAD8/9!");
772 }
773 }
774
775 case INVALID_ELEM:
776 return 0;
777
778 default:
779 libmesh_error_msg("ERROR: Invalid ElemType " << Utility::enum_to_string(t) << " selected for SZABAB FE family!");
780 }
781 }
782
783
784
785 // The sixth-order Szabo-Babuska shape functions
786 case SIXTH:
787 {
788 switch (t)
789 {
790 case NODEELEM:
791 return 1;
792
793 // The 1D Szabo-Babuska defined on a three-noded edge
794 case EDGE2:
795 case EDGE3:
796 {
797 switch (n)
798 {
799 case 0:
800 case 1:
801 return 1;
802
803 case 2:
804 return 5;
805
806 default:
807 libmesh_error_msg("ERROR: Invalid node ID " << n << " selected for EDGE2/3!");
808 }
809 }
810
811
812 // The 2D Szabo-Babuska defined on a 6/7-noded triangle
813 case TRI7:
814 if (n == 6)
815 return 10;
816 libmesh_fallthrough();
817 case TRI6:
818 {
819 switch (n)
820 {
821 case 0:
822 case 1:
823 case 2:
824 return 1;
825
826 case 3:
827 case 4:
828 case 5:
829 return 5;
830
831 default:
832 libmesh_error_msg("ERROR: Invalid node ID " << n << " selected for TRI6!");
833 }
834 }
835
836
837 // The 2D tensor-product Szabo-Babuska defined on a
838 // nine-noded quadrilateral.
839 case QUAD8:
840 case QUAD9:
841 {
842 switch (n)
843 {
844 case 0:
845 case 1:
846 case 2:
847 case 3:
848 return 1;
849
850 case 4:
851 case 5:
852 case 6:
853 case 7:
854 return 5;
855
856 case 8:
857 return 25;
858
859 default:
860 libmesh_error_msg("ERROR: Invalid node ID " << n << " selected for QUAD8/9!");
861 }
862 }
863
864 case INVALID_ELEM:
865 return 0;
866
867 default:
868 libmesh_error_msg("ERROR: Invalid ElemType " << Utility::enum_to_string(t) << " selected for SZABAB FE family!");
869 }
870 }
871
872
873 // The seventh-order Szabo-Babuska shape functions
874 case SEVENTH:
875 {
876 switch (t)
877 {
878 case NODEELEM:
879 return 1;
880
881 // The 1D Szabo-Babuska defined on a three-noded edge
882 case EDGE2:
883 case EDGE3:
884 {
885 switch (n)
886 {
887 case 0:
888 case 1:
889 return 1;
890
891 case 2:
892 return 6;
893
894 default:
895 libmesh_error_msg("ERROR: Invalid node ID " << n << " selected for EDGE2/3!");
896 }
897 }
898
899
900 // The 2D Szabo-Babuska defined on a 6/7-noded triangle
901 case TRI7:
902 if (n == 6)
903 return 15;
904 libmesh_fallthrough();
905 case TRI6:
906 {
907 switch (n)
908 {
909 case 0:
910 case 1:
911 case 2:
912 return 1;
913
914 case 3:
915 case 4:
916 case 5:
917 return 6;
918
919 default:
920 libmesh_error_msg("ERROR: Invalid node ID " << n << " selected for TRI6!");
921 }
922 }
923
924
925 // The 2D tensor-product Szabo-Babuska defined on a
926 // nine-noded quadrilateral.
927 case QUAD8:
928 case QUAD9:
929 {
930 switch (n)
931 {
932 case 0:
933 case 1:
934 case 2:
935 case 3:
936 return 1;
937
938 case 4:
939 case 5:
940 case 6:
941 case 7:
942 return 6;
943
944 case 8:
945 return 36;
946
947 default:
948 libmesh_error_msg("ERROR: Invalid node ID " << n << " selected for QUAD8/9!");
949 }
950 }
951
952 case INVALID_ELEM:
953 return 0;
954
955 default:
956 libmesh_error_msg("ERROR: Invalid ElemType " << Utility::enum_to_string(t) << " selected for SZABAB FE family!");
957 }
958 }
959
960
961 default:
962 libmesh_error_msg("ERROR: Invalid Order " << Utility::enum_to_string(o) << " selected for SZABAB FE family!");
963 }
964} // szabab_n_dofs_at_node()
965
966
967
968unsigned int szabab_n_dofs_at_node(const Elem & e,
969 const Order o,
970 const unsigned int n)
971{
972 return szabab_n_dofs_at_node(e.type(), o, n);
973}
974
975
976
977unsigned int szabab_n_dofs_per_elem(const ElemType t, const Order o)
978{
979 switch (o)
980 {
981 // The first-order Szabo-Babuska shape functions
982 case FIRST:
983 {
984 switch (t)
985 {
986 case NODEELEM:
987 return 0;
988
989 // The 1D Szabo-Babuska defined on an edge
990 case EDGE2:
991 case EDGE3:
992 return 0;
993
994 // The 2D Szabo-Babuska defined on a triangle
995 case TRI3:
996 case TRI6:
997 case TRI7:
998 return 0;
999
1000 // The 2D tensor-product Szabo-Babuska defined on a
1001 // quadrilateral.
1002 case QUAD4:
1003 case QUAD8:
1004 case QUAD9:
1005 return 0;
1006
1007 case INVALID_ELEM:
1008 return 0;
1009
1010 default:
1011 libmesh_error_msg("ERROR: Invalid ElemType " << Utility::enum_to_string(t) << " selected for SZABAB FE family!");
1012 }
1013 }
1014
1015
1016
1017 // The second-order Szabo-Babuska shape functions
1018 case SECOND:
1019 {
1020 switch (t)
1021 {
1022 case NODEELEM:
1023 return 0;
1024
1025 // The 1D Szabo-Babuska defined on a two-noded edge
1026 case EDGE2:
1027 return 1;
1028
1029 // The 1D Szabo-Babuska defined on a three-noded edge
1030 case EDGE3:
1031 return 0;
1032
1033 // The 2D Szabo-Babuska defined on a 6-noded triangle
1034 case TRI6:
1035 case TRI7:
1036 return 0;
1037
1038 // The 2D tensor-product Szabo-Babuska defined on a
1039 // eight-noded quadrilateral.
1040 case QUAD8:
1041 return 0;
1042
1043 // The 2D tensor-product Szabo-Babuska defined on a
1044 // nine-noded quadrilateral.
1045 case QUAD9:
1046 return 0;
1047
1048 case INVALID_ELEM:
1049 return 0;
1050
1051 default:
1052 libmesh_error_msg("ERROR: Invalid ElemType " << Utility::enum_to_string(t) << " selected for SZABAB FE family!");
1053 }
1054 }
1055
1056
1057
1058 // The third-order Szabo-Babuska shape functions
1059 case THIRD:
1060 {
1061 switch (t)
1062 {
1063 case NODEELEM:
1064 return 0;
1065
1066 // The 1D Szabo-Babuska defined on a two-noded edge
1067 case EDGE2:
1068 return 2;
1069
1070 // The 1D Szabo-Babuska defined on a three-noded edge
1071 case EDGE3:
1072 return 0;
1073
1074 // The 2D Szabo-Babuska defined on a 6-noded triangle
1075 case TRI6:
1076 return 1;
1077
1078 case TRI7:
1079 return 0;
1080
1081 // The 2D tensor-product Szabo-Babuska defined on a
1082 // eight-noded quadrilateral.
1083 case QUAD8:
1084 return 4;
1085
1086 // The 2D tensor-product Szabo-Babuska defined on a
1087 // nine-noded quadrilateral.
1088 case QUAD9:
1089 return 0;
1090
1091 case INVALID_ELEM:
1092 return 0;
1093
1094 default:
1095 libmesh_error_msg("ERROR: Invalid ElemType " << Utility::enum_to_string(t) << " selected for SZABAB FE family!");
1096 }
1097 }
1098
1099
1100
1101 // The fourth-order Szabo-Babuska shape functions
1102 case FOURTH:
1103 {
1104 switch (t)
1105 {
1106 case NODEELEM:
1107 return 0;
1108
1109 // The 1D Szabo-Babuska defined on a two-noded edge
1110 case EDGE2:
1111 return 3;
1112
1113 // The 1D Szabo-Babuska defined on a three-noded edge
1114 case EDGE3:
1115 return 0;
1116
1117 // The 2D Szabo-Babuska defined on a 6-noded triangle
1118 case TRI6:
1119 return 3;
1120
1121 case TRI7:
1122 return 0;
1123
1124 // The 2D tensor-product Szabo-Babuska defined on a
1125 // eight-noded quadrilateral.
1126 case QUAD8:
1127 return 9;
1128
1129 // The 2D tensor-product Szabo-Babuska defined on a
1130 // nine-noded quadrilateral.
1131 case QUAD9:
1132 return 0;
1133
1134 case INVALID_ELEM:
1135 return 0;
1136
1137 default:
1138 libmesh_error_msg("ERROR: Invalid ElemType " << Utility::enum_to_string(t) << " selected for SZABAB FE family!");
1139 }
1140 }
1141
1142
1143
1144 // The fifth-order Szabo-Babuska shape functions
1145 case FIFTH:
1146 {
1147 switch (t)
1148 {
1149 case NODEELEM:
1150 return 0;
1151
1152 // The 1D Szabo-Babuska defined on a two-noded edge
1153 case EDGE2:
1154 return 4;
1155
1156 // The 1D Szabo-Babuska defined on a three-noded edge
1157 case EDGE3:
1158 return 0;
1159
1160 // The 2D Szabo-Babuska defined on a 6-noded triangle
1161 case TRI6:
1162 return 6;
1163
1164 case TRI7:
1165 return 0;
1166
1167 // The 2D tensor-product Szabo-Babuska defined on a
1168 // eight-noded quadrilateral.
1169 case QUAD8:
1170 return 16;
1171
1172 // The 2D tensor-product Szabo-Babuska defined on a
1173 // nine-noded quadrilateral.
1174 case QUAD9:
1175 return 0;
1176
1177 case INVALID_ELEM:
1178 return 0;
1179
1180 default:
1181 libmesh_error_msg("ERROR: Invalid ElemType " << Utility::enum_to_string(t) << " selected for SZABAB FE family!");
1182 }
1183 }
1184
1185
1186 // The sixth-order Szabo-Babuska shape functions
1187 case SIXTH:
1188 {
1189 switch (t)
1190 {
1191 case NODEELEM:
1192 return 0;
1193
1194 // The 1D Szabo-Babuska defined on a two-noded edge
1195 case EDGE2:
1196 return 5;
1197
1198 // The 1D Szabo-Babuska defined on a three-noded edge
1199 case EDGE3:
1200 return 0;
1201
1202 // The 2D Szabo-Babuska defined on a 6-noded triangle
1203 case TRI6:
1204 return 10;
1205
1206 case TRI7:
1207 return 0;
1208
1209 // The 2D tensor-product Szabo-Babuska defined on a
1210 // eight-noded quadrilateral.
1211 case QUAD8:
1212 return 25;
1213
1214 // The 2D tensor-product Szabo-Babuska defined on a
1215 // nine-noded quadrilateral.
1216 case QUAD9:
1217 return 0;
1218
1219 case INVALID_ELEM:
1220 return 0;
1221
1222 default:
1223 libmesh_error_msg("ERROR: Invalid ElemType " << Utility::enum_to_string(t) << " selected for SZABAB FE family!");
1224 }
1225 }
1226
1227
1228 // The seventh-order Szabo-Babuska shape functions
1229 case SEVENTH:
1230 {
1231 switch (t)
1232 {
1233 case NODEELEM:
1234 return 0;
1235
1236 // The 1D Szabo-Babuska defined on a two-noded edge
1237 case EDGE2:
1238 return 6;
1239
1240 // The 1D Szabo-Babuska defined on a three-noded edge
1241 case EDGE3:
1242 return 0;
1243
1244 // The 2D Szabo-Babuska defined on a 6-noded triangle
1245 case TRI6:
1246 return 15;
1247
1248 case TRI7:
1249 return 0;
1250
1251 // The 2D tensor-product Szabo-Babuska defined on a
1252 // eight-noded quadrilateral.
1253 case QUAD8:
1254 return 36;
1255
1256 // The 2D tensor-product Szabo-Babuska defined on a
1257 // nine-noded quadrilateral.
1258 case QUAD9:
1259 return 0;
1260
1261 case INVALID_ELEM:
1262 return 0;
1263
1264 default:
1265 libmesh_error_msg("ERROR: Invalid ElemType " << Utility::enum_to_string(t) << " selected for SZABAB FE family!");
1266 }
1267 }
1268
1269
1270 // Otherwise no DOFS per element
1271 default:
1272 return 0;
1273 }
1274} // szabab_n_dofs_per_elem
1275
1276
1277
1278unsigned int szabab_n_dofs_per_elem(const Elem & e, const Order o)
1279{
1280 return szabab_n_dofs_per_elem(e.type(), o);
1281}
1282
1283} // anonymous namespace
1284
1285
1286// Instantiate (side_) nodal_soln() function for every dimension
1287LIBMESH_FE_NODAL_SOLN(SZABAB, szabab_nodal_soln)
1288LIBMESH_FE_SIDE_NODAL_SOLN(SZABAB)
1289
1290
1291// Full specialization of n_dofs() function for every dimension
1292template <> unsigned int FE<0,SZABAB>::n_dofs(const ElemType t, const Order o) { return szabab_n_dofs(t, o); }
1293template <> unsigned int FE<1,SZABAB>::n_dofs(const ElemType t, const Order o) { return szabab_n_dofs(t, o); }
1294template <> unsigned int FE<2,SZABAB>::n_dofs(const ElemType t, const Order o) { return szabab_n_dofs(t, o); }
1295template <> unsigned int FE<3,SZABAB>::n_dofs(const ElemType t, const Order o) { return szabab_n_dofs(t, o); }
1296
1297template <> unsigned int FE<0,SZABAB>::n_dofs(const Elem * e, const Order o) { return szabab_n_dofs(e, o); }
1298template <> unsigned int FE<1,SZABAB>::n_dofs(const Elem * e, const Order o) { return szabab_n_dofs(e, o); }
1299template <> unsigned int FE<2,SZABAB>::n_dofs(const Elem * e, const Order o) { return szabab_n_dofs(e, o); }
1300template <> unsigned int FE<3,SZABAB>::n_dofs(const Elem * e, const Order o) { return szabab_n_dofs(e, o); }
1301
1302// Full specialization of n_dofs_at_node() function for every dimension.
1303template <> unsigned int FE<0,SZABAB>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return szabab_n_dofs_at_node(t, o, n); }
1304template <> unsigned int FE<1,SZABAB>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return szabab_n_dofs_at_node(t, o, n); }
1305template <> unsigned int FE<2,SZABAB>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return szabab_n_dofs_at_node(t, o, n); }
1306template <> unsigned int FE<3,SZABAB>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return szabab_n_dofs_at_node(t, o, n); }
1307
1308template <> unsigned int FE<0,SZABAB>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return szabab_n_dofs_at_node(e, o, n); }
1309template <> unsigned int FE<1,SZABAB>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return szabab_n_dofs_at_node(e, o, n); }
1310template <> unsigned int FE<2,SZABAB>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return szabab_n_dofs_at_node(e, o, n); }
1311template <> unsigned int FE<3,SZABAB>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return szabab_n_dofs_at_node(e, o, n); }
1312
1313// Full specialization of n_dofs_per_elem() function for every dimension.
1314template <> unsigned int FE<0,SZABAB>::n_dofs_per_elem(const ElemType t, const Order o) { return szabab_n_dofs_per_elem(t, o); }
1315template <> unsigned int FE<1,SZABAB>::n_dofs_per_elem(const ElemType t, const Order o) { return szabab_n_dofs_per_elem(t, o); }
1316template <> unsigned int FE<2,SZABAB>::n_dofs_per_elem(const ElemType t, const Order o) { return szabab_n_dofs_per_elem(t, o); }
1317template <> unsigned int FE<3,SZABAB>::n_dofs_per_elem(const ElemType t, const Order o) { return szabab_n_dofs_per_elem(t, o); }
1318
1319template <> unsigned int FE<0,SZABAB>::n_dofs_per_elem(const Elem & e, const Order o) { return szabab_n_dofs_per_elem(e, o); }
1320template <> unsigned int FE<1,SZABAB>::n_dofs_per_elem(const Elem & e, const Order o) { return szabab_n_dofs_per_elem(e, o); }
1321template <> unsigned int FE<2,SZABAB>::n_dofs_per_elem(const Elem & e, const Order o) { return szabab_n_dofs_per_elem(e, o); }
1322template <> unsigned int FE<3,SZABAB>::n_dofs_per_elem(const Elem & e, const Order o) { return szabab_n_dofs_per_elem(e, o); }
1323
1324// Szabab FEMs are C^0 continuous
1325template <> FEContinuity FE<0,SZABAB>::get_continuity() const { return C_ZERO; }
1326template <> FEContinuity FE<1,SZABAB>::get_continuity() const { return C_ZERO; }
1327template <> FEContinuity FE<2,SZABAB>::get_continuity() const { return C_ZERO; }
1328template <> FEContinuity FE<3,SZABAB>::get_continuity() const { return C_ZERO; }
1329
1330// Szabab FEMs are hierarchic
1331template <> bool FE<0,SZABAB>::is_hierarchic() const { return true; }
1332template <> bool FE<1,SZABAB>::is_hierarchic() const { return true; }
1333template <> bool FE<2,SZABAB>::is_hierarchic() const { return true; }
1334template <> bool FE<3,SZABAB>::is_hierarchic() const { return true; }
1335
1336#ifdef LIBMESH_ENABLE_AMR
1337// compute_constraints() specializations are only needed for 2 and 3D
1338template <>
1339void FE<2,SZABAB>::compute_constraints (DofConstraints & constraints,
1340 DofMap & dof_map,
1341 const unsigned int variable_number,
1342 const Elem * elem)
1343{ compute_proj_constraints(constraints, dof_map, variable_number, elem); }
1344
1345template <>
1346void FE<3,SZABAB>::compute_constraints (DofConstraints & constraints,
1347 DofMap & dof_map,
1348 const unsigned int variable_number,
1349 const Elem * elem)
1350{ compute_proj_constraints(constraints, dof_map, variable_number, elem); }
1351#endif // #ifdef LIBMESH_ENABLE_AMR
1352
1353// Szabab shapes need reinit only for approximation orders >= 3,
1354// but we might reach that with p refinement
1355template <> bool FE<0,SZABAB>::shapes_need_reinit() const { return true; }
1356template <> bool FE<1,SZABAB>::shapes_need_reinit() const { return true; }
1357template <> bool FE<2,SZABAB>::shapes_need_reinit() const { return true; }
1358template <> bool FE<3,SZABAB>::shapes_need_reinit() const { return true; }
1359
1360} // namespace libMesh
1361
1362#endif //LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
libMesh::DofConstraints
The constraint matrix storage format.
Definition dof_map.h:112
libMesh::DofMap
This class handles the numbering of degrees of freedom on a mesh.
Definition dof_map.h:181
libMesh::Elem
This is the base class from which all geometric element types are derived.
Definition elem.h:96
libMesh::FEAbstract::get_refspace_nodes
static void get_refspace_nodes(const ElemType t, std::vector< Point > &nodes)
Definition fe_abstract.C:400
libMesh::FEInterface::shape
static Real shape(const unsigned int dim, const FEType &fe_t, const ElemType t, const unsigned int i, const Point &p)
Definition fe_interface.C:761
libMesh::FEInterface::n_shape_functions
static unsigned int n_shape_functions(const unsigned int dim, const FEType &fe_t, const ElemType t)
Definition fe_interface.C:272
libMesh::FE::is_hierarchic
virtual bool is_hierarchic() const override
libMesh::FE::n_dofs_per_elem
static unsigned int n_dofs_per_elem(const ElemType t, const Order o)
libMesh::FE::get_continuity
virtual FEContinuity get_continuity() const override
libMesh::FE::compute_constraints
static void compute_constraints(DofConstraints &constraints, DofMap &dof_map, const unsigned int variable_number, const Elem *elem)
Computes the constraint matrix contributions (for non-conforming adapted meshes) corresponding to var...
libMesh::FE::n_dofs
static unsigned int n_dofs(const ElemType t, const Order o)
libMesh::FE::n_dofs_at_node
static unsigned int n_dofs_at_node(const ElemType t, const Order o, const unsigned int n)
libMesh::FE::shapes_need_reinit
virtual bool shapes_need_reinit() const override
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_FE_NODAL_SOLN
LIBMESH_FE_NODAL_SOLN(BERNSTEIN, bernstein_nodal_soln)
Definition fe_bernstein.C:415
libMesh::ElemType
ElemType
Defines an enum for geometric element types.
Definition enum_elem_type.h:33
libMesh::libmesh_assert
libmesh_assert(ctx)
libMesh::LIBMESH_FE_SIDE_NODAL_SOLN
LIBMESH_FE_SIDE_NODAL_SOLN(HIERARCHIC_VEC)
Definition fe_hierarchic_vec.C:117
n_nodes
const dof_id_type n_nodes
Definition tecplot_io.C:67
Generated on Fri Aug 21 2026 13:32:29 for libMesh by doxygen 1.9.8