libMesh
fe_szabab.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2025 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 
30 namespace libMesh
31 {
32 
33 // ------------------------------------------------------------
34 // Szabo-Babuska-specific implementations, Steffen Petersen 2004
35 
36 // Anonymous namespace for local helper functions
37 namespace {
38 
39 void 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 
106 unsigned 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 
322 unsigned 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 
330 unsigned 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 
968 unsigned 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 
977 unsigned 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 
1278 unsigned 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
1287 LIBMESH_FE_NODAL_SOLN(SZABAB, szabab_nodal_soln)
1288 LIBMESH_FE_SIDE_NODAL_SOLN(SZABAB)
1289 
1290 
1291 // Full specialization of n_dofs() function for every dimension
1292 template <> unsigned int FE<0,SZABAB>::n_dofs(const ElemType t, const Order o) { return szabab_n_dofs(t, o); }
1293 template <> unsigned int FE<1,SZABAB>::n_dofs(const ElemType t, const Order o) { return szabab_n_dofs(t, o); }
1294 template <> unsigned int FE<2,SZABAB>::n_dofs(const ElemType t, const Order o) { return szabab_n_dofs(t, o); }
1295 template <> unsigned int FE<3,SZABAB>::n_dofs(const ElemType t, const Order o) { return szabab_n_dofs(t, o); }
1296 
1297 template <> unsigned int FE<0,SZABAB>::n_dofs(const Elem * e, const Order o) { return szabab_n_dofs(e, o); }
1298 template <> unsigned int FE<1,SZABAB>::n_dofs(const Elem * e, const Order o) { return szabab_n_dofs(e, o); }
1299 template <> unsigned int FE<2,SZABAB>::n_dofs(const Elem * e, const Order o) { return szabab_n_dofs(e, o); }
1300 template <> 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.
1303 template <> 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); }
1304 template <> 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); }
1305 template <> 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); }
1306 template <> 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 
1308 template <> 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); }
1309 template <> 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); }
1310 template <> 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); }
1311 template <> 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.
1314 template <> unsigned int FE<0,SZABAB>::n_dofs_per_elem(const ElemType t, const Order o) { return szabab_n_dofs_per_elem(t, o); }
1315 template <> unsigned int FE<1,SZABAB>::n_dofs_per_elem(const ElemType t, const Order o) { return szabab_n_dofs_per_elem(t, o); }
1316 template <> unsigned int FE<2,SZABAB>::n_dofs_per_elem(const ElemType t, const Order o) { return szabab_n_dofs_per_elem(t, o); }
1317 template <> unsigned int FE<3,SZABAB>::n_dofs_per_elem(const ElemType t, const Order o) { return szabab_n_dofs_per_elem(t, o); }
1318 
1319 template <> unsigned int FE<0,SZABAB>::n_dofs_per_elem(const Elem & e, const Order o) { return szabab_n_dofs_per_elem(e, o); }
1320 template <> unsigned int FE<1,SZABAB>::n_dofs_per_elem(const Elem & e, const Order o) { return szabab_n_dofs_per_elem(e, o); }
1321 template <> unsigned int FE<2,SZABAB>::n_dofs_per_elem(const Elem & e, const Order o) { return szabab_n_dofs_per_elem(e, o); }
1322 template <> 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
1325 template <> FEContinuity FE<0,SZABAB>::get_continuity() const { return C_ZERO; }
1326 template <> FEContinuity FE<1,SZABAB>::get_continuity() const { return C_ZERO; }
1327 template <> FEContinuity FE<2,SZABAB>::get_continuity() const { return C_ZERO; }
1328 template <> FEContinuity FE<3,SZABAB>::get_continuity() const { return C_ZERO; }
1329 
1330 // Szabab FEMs are hierarchic
1331 template <> bool FE<0,SZABAB>::is_hierarchic() const { return true; }
1332 template <> bool FE<1,SZABAB>::is_hierarchic() const { return true; }
1333 template <> bool FE<2,SZABAB>::is_hierarchic() const { return true; }
1334 template <> 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
1338 template <>
1339 void 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 
1345 template <>
1346 void 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
1355 template <> bool FE<0,SZABAB>::shapes_need_reinit() const { return true; }
1356 template <> bool FE<1,SZABAB>::shapes_need_reinit() const { return true; }
1357 template <> bool FE<2,SZABAB>::shapes_need_reinit() const { return true; }
1358 template <> bool FE<3,SZABAB>::shapes_need_reinit() const { return true; }
1359 
1360 } // namespace libMesh
1361 
1362 #endif //LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
libMesh::FE::n_dofs
static unsigned int n_dofs(const ElemType t, const Order o)
libMesh::ElemType
ElemType
Defines an enum for geometric element types.
Definition: enum_elem_type.h:33
libMesh::Order
Order
defines an enum for polynomial orders.
Definition: enum_order.h:40
libMesh::Elem
This is the base class from which all geometric element types are derived.
Definition: elem.h:94
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
The libMesh namespace provides an interface to certain functionality in the library.
libMesh::LIBMESH_FE_NODAL_SOLN
LIBMESH_FE_NODAL_SOLN(LIBMESH_FE_SIDE_NODAL_SOLN() LIBMESH_DEFAULT_NDOFS(BERNSTEIN) template<> FEContinuity FE< 0 BERNSTEIN, bernstein_nodal_soln)
Definition: fe_bernstein.C:415
libMesh::FE::is_hierarchic
virtual bool is_hierarchic() const override
libMesh::DofMap
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:179
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
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::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:760
libMesh::FEAbstract::get_refspace_nodes
static void get_refspace_nodes(const ElemType t, std::vector< Point > &nodes)
Definition: fe_abstract.C:400
libMesh::libmesh_assert
libmesh_assert(ctx)
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::Utility::enum_to_string
std::string enum_to_string(const T e)
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::FEContinuity
FEContinuity
defines an enum for finite element types to libmesh_assert a certain level (or type? Hcurl?) of continuity.
Definition: enum_fe_family.h:84
libMesh::DofConstraints
The constraint matrix storage format.
Definition: dof_map.h:108
Generated on Thu Jul 17 2025 01:29:08 for libMesh by   doxygen 1.8.14