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