Line data Source code
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 :
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 576 : 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 576 : const unsigned int n_nodes = elem->n_nodes();
46 :
47 576 : const ElemType elem_type = elem->type();
48 :
49 576 : nodal_soln.resize(n_nodes);
50 :
51 624 : const Order totalorder = order + add_p_level*elem->p_level();
52 :
53 : // FEType object to be passed to various FEInterface functions below.
54 48 : FEType fe_type(order, SZABAB);
55 :
56 576 : switch (totalorder)
57 : {
58 : // Constant shape functions
59 0 : case CONSTANT:
60 : {
61 0 : libmesh_assert_equal_to (elem_soln.size(), 1);
62 :
63 0 : std::fill(nodal_soln.begin(), nodal_soln.end(), elem_soln[0]);
64 :
65 0 : return;
66 : }
67 :
68 :
69 : // For other bases do interpolation at the nodes
70 : // explicitly.
71 576 : 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 576 : FEInterface::n_shape_functions(fe_type, elem);
81 :
82 96 : std::vector<Point> refspace_nodes;
83 576 : FEBase::get_refspace_nodes(elem_type,refspace_nodes);
84 48 : libmesh_assert_equal_to (refspace_nodes.size(), n_nodes);
85 48 : libmesh_assert_equal_to (elem_soln.size(), n_sf);
86 :
87 : // Zero before summation
88 48 : std::fill(nodal_soln.begin(), nodal_soln.end(), 0);
89 :
90 4032 : for (unsigned int n=0; n<n_nodes; n++)
91 : // u_i = Sum (alpha_i phi_i)
92 24192 : for (unsigned int i=0; i<n_sf; i++)
93 22464 : nodal_soln[n] += elem_soln[i] *
94 22464 : FEInterface::shape(fe_type, elem, i, refspace_nodes[n]);
95 :
96 48 : return;
97 : }
98 :
99 0 : default:
100 0 : libmesh_error_msg("ERROR: Invalid total order " << totalorder);
101 : }
102 : } // szabab_nodal_soln()
103 :
104 :
105 :
106 46768 : unsigned int szabab_n_dofs(const ElemType t, const Order o)
107 : {
108 46768 : switch (o)
109 : {
110 : // Szabo-Babuska 1st-order polynomials.
111 6762 : case FIRST:
112 : {
113 6762 : switch (t)
114 : {
115 0 : case NODEELEM:
116 0 : return 1;
117 :
118 134 : case EDGE2:
119 : case EDGE3:
120 134 : return 2;
121 :
122 726 : case TRI3:
123 : case TRI6:
124 : case TRI7:
125 726 : return 3;
126 :
127 606 : case QUAD4:
128 : case QUAD8:
129 : case QUAD9:
130 606 : return 4;
131 :
132 0 : case INVALID_ELEM:
133 0 : return 0;
134 :
135 0 : default:
136 0 : 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 14118 : case SECOND:
143 : {
144 14118 : switch (t)
145 : {
146 0 : case NODEELEM:
147 0 : return 1;
148 :
149 174 : case EDGE2:
150 : case EDGE3:
151 174 : return 3;
152 :
153 2370 : case TRI6:
154 : case TRI7:
155 2370 : return 6;
156 :
157 0 : case QUAD8:
158 0 : return 8;
159 806 : case QUAD9:
160 806 : return 9;
161 :
162 0 : case INVALID_ELEM:
163 0 : return 0;
164 :
165 0 : default:
166 0 : 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 12944 : case THIRD:
173 : {
174 12944 : switch (t)
175 : {
176 0 : case NODEELEM:
177 0 : return 1;
178 :
179 174 : case EDGE2:
180 : case EDGE3:
181 174 : return 4;
182 :
183 2112 : case TRI6:
184 : case TRI7:
185 2112 : return 10;
186 :
187 3836 : case QUAD8:
188 : case QUAD9:
189 3836 : return 16;
190 :
191 0 : case INVALID_ELEM:
192 0 : return 0;
193 :
194 0 : default:
195 0 : 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 12944 : case FOURTH:
202 : {
203 12944 : switch (t)
204 : {
205 0 : case NODEELEM:
206 0 : return 1;
207 :
208 174 : case EDGE2:
209 : case EDGE3:
210 174 : return 5;
211 :
212 8464 : case TRI6:
213 : case TRI7:
214 8464 : return 15;
215 :
216 3836 : case QUAD8:
217 : case QUAD9:
218 3836 : return 25;
219 :
220 0 : case INVALID_ELEM:
221 0 : return 0;
222 :
223 0 : default:
224 0 : 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 0 : case FIFTH:
231 : {
232 0 : switch (t)
233 : {
234 0 : case NODEELEM:
235 0 : return 1;
236 :
237 0 : case EDGE2:
238 : case EDGE3:
239 0 : return 6;
240 :
241 0 : case TRI6:
242 : case TRI7:
243 0 : return 21;
244 :
245 0 : case QUAD8:
246 : case QUAD9:
247 0 : return 36;
248 :
249 0 : case INVALID_ELEM:
250 0 : return 0;
251 :
252 0 : default:
253 0 : 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 0 : case SIXTH:
260 : {
261 0 : switch (t)
262 : {
263 0 : case NODEELEM:
264 0 : return 1;
265 :
266 0 : case EDGE2:
267 : case EDGE3:
268 0 : return 7;
269 :
270 0 : case TRI6:
271 : case TRI7:
272 0 : return 28;
273 :
274 0 : case QUAD8:
275 : case QUAD9:
276 0 : return 49;
277 :
278 0 : case INVALID_ELEM:
279 0 : return 0;
280 :
281 0 : default:
282 0 : 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 0 : case SEVENTH:
288 : {
289 0 : switch (t)
290 : {
291 0 : case NODEELEM:
292 0 : return 1;
293 :
294 0 : case EDGE2:
295 : case EDGE3:
296 0 : return 8;
297 :
298 0 : case TRI6:
299 : case TRI7:
300 0 : return 36;
301 :
302 0 : case QUAD8:
303 : case QUAD9:
304 0 : return 64;
305 :
306 0 : case INVALID_ELEM:
307 0 : return 0;
308 :
309 0 : default:
310 0 : libmesh_error_msg("ERROR: Invalid ElemType " << Utility::enum_to_string(t) << " selected for SZABAB FE family!");
311 : }
312 : }
313 :
314 :
315 0 : default:
316 0 : 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 43027 : unsigned int szabab_n_dofs(const Elem * e, const Order o)
323 : {
324 7259 : libmesh_assert(e);
325 46768 : return szabab_n_dofs(e->type(), o);
326 : }
327 :
328 :
329 :
330 354216 : unsigned int szabab_n_dofs_at_node(const ElemType t,
331 : const Order o,
332 : const unsigned int n)
333 : {
334 354216 : switch (o)
335 : {
336 : // The first-order Szabo-Babuska shape functions
337 17780 : case FIRST:
338 : {
339 17780 : switch (t)
340 : {
341 0 : case NODEELEM:
342 0 : return 1;
343 :
344 : // The 1D Szabo-Babuska defined on a three-noded edge
345 2120 : case EDGE2:
346 : case EDGE3:
347 : {
348 2120 : switch (n)
349 : {
350 210 : case 0:
351 : case 1:
352 210 : return 1;
353 :
354 0 : case 2:
355 0 : return 0;
356 :
357 0 : default:
358 0 : 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 9620 : case TRI3:
365 : case TRI6:
366 : case TRI7:
367 : {
368 8740 : switch (n)
369 : {
370 880 : case 0:
371 : case 1:
372 : case 2:
373 880 : return 1;
374 :
375 0 : case 3:
376 : case 4:
377 : case 5:
378 0 : return 0;
379 :
380 0 : case 6:
381 0 : libmesh_assert(t == TRI7);
382 0 : return 0;
383 :
384 0 : default:
385 0 : 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 6040 : case QUAD4:
393 : case QUAD8:
394 : case QUAD9:
395 : {
396 5470 : switch (n)
397 : {
398 570 : case 0:
399 : case 1:
400 : case 2:
401 : case 3:
402 570 : return 1;
403 :
404 0 : case 4:
405 : case 5:
406 : case 6:
407 : case 7:
408 0 : return 0;
409 :
410 0 : case 8:
411 0 : return 0;
412 :
413 0 : default:
414 0 : libmesh_error_msg("ERROR: Invalid node ID " << n << " selected for QUAD4/8/9!");
415 : }
416 : }
417 :
418 0 : case INVALID_ELEM:
419 0 : return 0;
420 :
421 0 : default:
422 0 : 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 125796 : case SECOND:
430 : {
431 125796 : switch (t)
432 : {
433 0 : case NODEELEM:
434 0 : return 1;
435 :
436 : // The 1D Szabo-Babuska defined on a three-noded edge
437 3930 : case EDGE2:
438 : case EDGE3:
439 : {
440 3930 : switch (n)
441 : {
442 250 : case 0:
443 : case 1:
444 250 : return 1;
445 :
446 120 : case 2:
447 120 : return 1;
448 :
449 0 : default:
450 0 : 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 93816 : case TRI6:
457 : case TRI7:
458 : {
459 93816 : switch (n)
460 : {
461 3756 : case 0:
462 : case 1:
463 : case 2:
464 3756 : return 1;
465 :
466 3398 : case 3:
467 : case 4:
468 : case 5:
469 3398 : return 1;
470 :
471 900 : case 6:
472 450 : libmesh_assert(t == TRI7);
473 900 : return 0;
474 :
475 0 : default:
476 0 : 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 28050 : case QUAD8:
484 : case QUAD9:
485 : {
486 28050 : switch (n)
487 : {
488 1130 : case 0:
489 : case 1:
490 : case 2:
491 : case 3:
492 1130 : return 1;
493 :
494 1050 : case 4:
495 : case 5:
496 : case 6:
497 : case 7:
498 1050 : return 1;
499 :
500 270 : case 8:
501 270 : return 1;
502 :
503 0 : default:
504 0 : libmesh_error_msg("ERROR: Invalid node ID " << n << " selected for QUAD8/9!");
505 : }
506 : }
507 :
508 0 : case INVALID_ELEM:
509 0 : return 0;
510 :
511 0 : default:
512 0 : 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 105320 : case THIRD:
520 : {
521 105320 : switch (t)
522 : {
523 0 : case NODEELEM:
524 0 : return 1;
525 :
526 : // The 1D Szabo-Babuska defined on a three-noded edge
527 3930 : case EDGE2:
528 : case EDGE3:
529 : {
530 3930 : switch (n)
531 : {
532 250 : case 0:
533 : case 1:
534 250 : return 1;
535 :
536 240 : case 2:
537 240 : return 2;
538 :
539 0 : default:
540 0 : 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 39440 : case TRI7:
547 39440 : if (n == 6)
548 490 : return 1;
549 : libmesh_fallthrough();
550 : case TRI6:
551 : {
552 67800 : switch (n)
553 : {
554 3140 : case 0:
555 : case 1:
556 : case 2:
557 3140 : return 1;
558 :
559 5760 : case 3:
560 : case 4:
561 : case 5:
562 5760 : return 2;
563 :
564 0 : default:
565 0 : 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 28050 : case QUAD8:
573 : case QUAD9:
574 : {
575 25600 : switch (n)
576 : {
577 1130 : case 0:
578 : case 1:
579 : case 2:
580 : case 3:
581 1130 : return 1;
582 :
583 1050 : case 4:
584 : case 5:
585 : case 6:
586 : case 7:
587 1050 : return 2;
588 :
589 270 : case 8:
590 270 : return 4;
591 :
592 0 : default:
593 0 : libmesh_error_msg("ERROR: Invalid node ID " << n << " selected for QUAD8/9!");
594 : }
595 : }
596 :
597 0 : case INVALID_ELEM:
598 0 : return 0;
599 :
600 0 : default:
601 0 : 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 105320 : case FOURTH:
609 : {
610 105320 : switch (t)
611 : {
612 0 : case NODEELEM:
613 0 : return 1;
614 :
615 : // The 1D Szabo-Babuska defined on a three-noded edge
616 3930 : case EDGE2:
617 : case EDGE3:
618 : {
619 3930 : switch (n)
620 : {
621 250 : case 0:
622 : case 1:
623 250 : return 1;
624 :
625 240 : case 2:
626 240 : return 3;
627 :
628 0 : default:
629 0 : 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 39440 : case TRI7:
636 39440 : if (n == 6)
637 490 : return 3;
638 : libmesh_fallthrough();
639 : case TRI6:
640 : {
641 67800 : switch (n)
642 : {
643 3140 : case 0:
644 : case 1:
645 : case 2:
646 3140 : return 1;
647 :
648 5760 : case 3:
649 : case 4:
650 : case 5:
651 5760 : return 3;
652 :
653 0 : default:
654 0 : 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 28050 : case QUAD8:
662 : case QUAD9:
663 : {
664 25600 : switch (n)
665 : {
666 1130 : case 0:
667 : case 1:
668 : case 2:
669 : case 3:
670 1130 : return 1;
671 :
672 1050 : case 4:
673 : case 5:
674 : case 6:
675 : case 7:
676 1050 : return 3;
677 :
678 270 : case 8:
679 270 : return 9;
680 :
681 0 : default:
682 0 : libmesh_error_msg("ERROR: Invalid node ID " << n << " selected for QUAD8/9!");
683 : }
684 : }
685 :
686 0 : case INVALID_ELEM:
687 0 : return 0;
688 :
689 0 : default:
690 0 : 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 0 : case FIFTH:
698 : {
699 0 : switch (t)
700 : {
701 0 : case NODEELEM:
702 0 : return 1;
703 :
704 : // The 1D Szabo-Babuska defined on a three-noded edge
705 0 : case EDGE2:
706 : case EDGE3:
707 : {
708 0 : switch (n)
709 : {
710 0 : case 0:
711 : case 1:
712 0 : return 1;
713 :
714 0 : case 2:
715 0 : return 4;
716 :
717 0 : default:
718 0 : 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 0 : case TRI7:
725 0 : if (n == 6)
726 0 : return 6;
727 : libmesh_fallthrough();
728 : case TRI6:
729 : {
730 0 : switch (n)
731 : {
732 0 : case 0:
733 : case 1:
734 : case 2:
735 0 : return 1;
736 :
737 0 : case 3:
738 : case 4:
739 : case 5:
740 0 : return 4;
741 :
742 0 : default:
743 0 : 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 0 : case QUAD8:
751 : case QUAD9:
752 : {
753 0 : switch (n)
754 : {
755 0 : case 0:
756 : case 1:
757 : case 2:
758 : case 3:
759 0 : return 1;
760 :
761 0 : case 4:
762 : case 5:
763 : case 6:
764 : case 7:
765 0 : return 4;
766 :
767 0 : case 8:
768 0 : return 16;
769 :
770 0 : default:
771 0 : libmesh_error_msg("ERROR: Invalid node ID " << n << " selected for QUAD8/9!");
772 : }
773 : }
774 :
775 0 : case INVALID_ELEM:
776 0 : return 0;
777 :
778 0 : default:
779 0 : 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 0 : case SIXTH:
787 : {
788 0 : switch (t)
789 : {
790 0 : case NODEELEM:
791 0 : return 1;
792 :
793 : // The 1D Szabo-Babuska defined on a three-noded edge
794 0 : case EDGE2:
795 : case EDGE3:
796 : {
797 0 : switch (n)
798 : {
799 0 : case 0:
800 : case 1:
801 0 : return 1;
802 :
803 0 : case 2:
804 0 : return 5;
805 :
806 0 : default:
807 0 : 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 0 : case TRI7:
814 0 : if (n == 6)
815 0 : return 10;
816 : libmesh_fallthrough();
817 : case TRI6:
818 : {
819 0 : switch (n)
820 : {
821 0 : case 0:
822 : case 1:
823 : case 2:
824 0 : return 1;
825 :
826 0 : case 3:
827 : case 4:
828 : case 5:
829 0 : return 5;
830 :
831 0 : default:
832 0 : 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 0 : case QUAD8:
840 : case QUAD9:
841 : {
842 0 : switch (n)
843 : {
844 0 : case 0:
845 : case 1:
846 : case 2:
847 : case 3:
848 0 : return 1;
849 :
850 0 : case 4:
851 : case 5:
852 : case 6:
853 : case 7:
854 0 : return 5;
855 :
856 0 : case 8:
857 0 : return 25;
858 :
859 0 : default:
860 0 : libmesh_error_msg("ERROR: Invalid node ID " << n << " selected for QUAD8/9!");
861 : }
862 : }
863 :
864 0 : case INVALID_ELEM:
865 0 : return 0;
866 :
867 0 : default:
868 0 : 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 0 : case SEVENTH:
875 : {
876 0 : switch (t)
877 : {
878 0 : case NODEELEM:
879 0 : return 1;
880 :
881 : // The 1D Szabo-Babuska defined on a three-noded edge
882 0 : case EDGE2:
883 : case EDGE3:
884 : {
885 0 : switch (n)
886 : {
887 0 : case 0:
888 : case 1:
889 0 : return 1;
890 :
891 0 : case 2:
892 0 : return 6;
893 :
894 0 : default:
895 0 : 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 0 : case TRI7:
902 0 : if (n == 6)
903 0 : return 15;
904 : libmesh_fallthrough();
905 : case TRI6:
906 : {
907 0 : switch (n)
908 : {
909 0 : case 0:
910 : case 1:
911 : case 2:
912 0 : return 1;
913 :
914 0 : case 3:
915 : case 4:
916 : case 5:
917 0 : return 6;
918 :
919 0 : default:
920 0 : 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 0 : case QUAD8:
928 : case QUAD9:
929 : {
930 0 : switch (n)
931 : {
932 0 : case 0:
933 : case 1:
934 : case 2:
935 : case 3:
936 0 : return 1;
937 :
938 0 : case 4:
939 : case 5:
940 : case 6:
941 : case 7:
942 0 : return 6;
943 :
944 0 : case 8:
945 0 : return 36;
946 :
947 0 : default:
948 0 : libmesh_error_msg("ERROR: Invalid node ID " << n << " selected for QUAD8/9!");
949 : }
950 : }
951 :
952 0 : case INVALID_ELEM:
953 0 : return 0;
954 :
955 0 : default:
956 0 : libmesh_error_msg("ERROR: Invalid ElemType " << Utility::enum_to_string(t) << " selected for SZABAB FE family!");
957 : }
958 : }
959 :
960 :
961 0 : default:
962 0 : 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 87786 : unsigned int szabab_n_dofs_at_node(const Elem & e,
969 : const Order o,
970 : const unsigned int n)
971 : {
972 96312 : return szabab_n_dofs_at_node(e.type(), o, n);
973 : }
974 :
975 :
976 :
977 48616 : unsigned int szabab_n_dofs_per_elem(const ElemType t, const Order o)
978 : {
979 48616 : switch (o)
980 : {
981 : // The first-order Szabo-Babuska shape functions
982 5280 : case FIRST:
983 : {
984 4800 : switch (t)
985 : {
986 0 : case NODEELEM:
987 0 : return 0;
988 :
989 : // The 1D Szabo-Babuska defined on an edge
990 90 : case EDGE2:
991 : case EDGE3:
992 90 : return 0;
993 :
994 : // The 2D Szabo-Babuska defined on a triangle
995 260 : case TRI3:
996 : case TRI6:
997 : case TRI7:
998 260 : return 0;
999 :
1000 : // The 2D tensor-product Szabo-Babuska defined on a
1001 : // quadrilateral.
1002 130 : case QUAD4:
1003 : case QUAD8:
1004 : case QUAD9:
1005 130 : return 0;
1006 :
1007 0 : case INVALID_ELEM:
1008 0 : return 0;
1009 :
1010 0 : default:
1011 0 : 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 16256 : case SECOND:
1019 : {
1020 14910 : switch (t)
1021 : {
1022 0 : case NODEELEM:
1023 0 : return 0;
1024 :
1025 : // The 1D Szabo-Babuska defined on a two-noded edge
1026 0 : case EDGE2:
1027 0 : return 1;
1028 :
1029 : // The 1D Szabo-Babuska defined on a three-noded edge
1030 110 : case EDGE3:
1031 110 : return 0;
1032 :
1033 : // The 2D Szabo-Babuska defined on a 6-noded triangle
1034 1006 : case TRI6:
1035 : case TRI7:
1036 1006 : return 0;
1037 :
1038 : // The 2D tensor-product Szabo-Babuska defined on a
1039 : // eight-noded quadrilateral.
1040 0 : case QUAD8:
1041 0 : return 0;
1042 :
1043 : // The 2D tensor-product Szabo-Babuska defined on a
1044 : // nine-noded quadrilateral.
1045 230 : case QUAD9:
1046 230 : return 0;
1047 :
1048 0 : case INVALID_ELEM:
1049 0 : return 0;
1050 :
1051 0 : default:
1052 0 : 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 13540 : case THIRD:
1060 : {
1061 13540 : switch (t)
1062 : {
1063 0 : case NODEELEM:
1064 0 : return 0;
1065 :
1066 : // The 1D Szabo-Babuska defined on a two-noded edge
1067 0 : case EDGE2:
1068 0 : return 2;
1069 :
1070 : // The 1D Szabo-Babuska defined on a three-noded edge
1071 110 : case EDGE3:
1072 110 : return 0;
1073 :
1074 : // The 2D Szabo-Babuska defined on a 6-noded triangle
1075 860 : case TRI6:
1076 860 : return 1;
1077 :
1078 430 : case TRI7:
1079 430 : return 0;
1080 :
1081 : // The 2D tensor-product Szabo-Babuska defined on a
1082 : // eight-noded quadrilateral.
1083 0 : case QUAD8:
1084 0 : return 4;
1085 :
1086 : // The 2D tensor-product Szabo-Babuska defined on a
1087 : // nine-noded quadrilateral.
1088 230 : case QUAD9:
1089 230 : return 0;
1090 :
1091 0 : case INVALID_ELEM:
1092 0 : return 0;
1093 :
1094 0 : default:
1095 0 : 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 13540 : case FOURTH:
1103 : {
1104 13540 : switch (t)
1105 : {
1106 0 : case NODEELEM:
1107 0 : return 0;
1108 :
1109 : // The 1D Szabo-Babuska defined on a two-noded edge
1110 0 : case EDGE2:
1111 0 : return 3;
1112 :
1113 : // The 1D Szabo-Babuska defined on a three-noded edge
1114 110 : case EDGE3:
1115 110 : return 0;
1116 :
1117 : // The 2D Szabo-Babuska defined on a 6-noded triangle
1118 860 : case TRI6:
1119 860 : return 3;
1120 :
1121 430 : case TRI7:
1122 430 : return 0;
1123 :
1124 : // The 2D tensor-product Szabo-Babuska defined on a
1125 : // eight-noded quadrilateral.
1126 0 : case QUAD8:
1127 0 : return 9;
1128 :
1129 : // The 2D tensor-product Szabo-Babuska defined on a
1130 : // nine-noded quadrilateral.
1131 230 : case QUAD9:
1132 230 : return 0;
1133 :
1134 0 : case INVALID_ELEM:
1135 0 : return 0;
1136 :
1137 0 : default:
1138 0 : 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 0 : case FIFTH:
1146 : {
1147 0 : switch (t)
1148 : {
1149 0 : case NODEELEM:
1150 0 : return 0;
1151 :
1152 : // The 1D Szabo-Babuska defined on a two-noded edge
1153 0 : case EDGE2:
1154 0 : return 4;
1155 :
1156 : // The 1D Szabo-Babuska defined on a three-noded edge
1157 0 : case EDGE3:
1158 0 : return 0;
1159 :
1160 : // The 2D Szabo-Babuska defined on a 6-noded triangle
1161 0 : case TRI6:
1162 0 : return 6;
1163 :
1164 0 : case TRI7:
1165 0 : return 0;
1166 :
1167 : // The 2D tensor-product Szabo-Babuska defined on a
1168 : // eight-noded quadrilateral.
1169 0 : case QUAD8:
1170 0 : return 16;
1171 :
1172 : // The 2D tensor-product Szabo-Babuska defined on a
1173 : // nine-noded quadrilateral.
1174 0 : case QUAD9:
1175 0 : return 0;
1176 :
1177 0 : case INVALID_ELEM:
1178 0 : return 0;
1179 :
1180 0 : default:
1181 0 : 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 0 : case SIXTH:
1188 : {
1189 0 : switch (t)
1190 : {
1191 0 : case NODEELEM:
1192 0 : return 0;
1193 :
1194 : // The 1D Szabo-Babuska defined on a two-noded edge
1195 0 : case EDGE2:
1196 0 : return 5;
1197 :
1198 : // The 1D Szabo-Babuska defined on a three-noded edge
1199 0 : case EDGE3:
1200 0 : return 0;
1201 :
1202 : // The 2D Szabo-Babuska defined on a 6-noded triangle
1203 0 : case TRI6:
1204 0 : return 10;
1205 :
1206 0 : case TRI7:
1207 0 : return 0;
1208 :
1209 : // The 2D tensor-product Szabo-Babuska defined on a
1210 : // eight-noded quadrilateral.
1211 0 : case QUAD8:
1212 0 : return 25;
1213 :
1214 : // The 2D tensor-product Szabo-Babuska defined on a
1215 : // nine-noded quadrilateral.
1216 0 : case QUAD9:
1217 0 : return 0;
1218 :
1219 0 : case INVALID_ELEM:
1220 0 : return 0;
1221 :
1222 0 : default:
1223 0 : 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 0 : case SEVENTH:
1230 : {
1231 0 : switch (t)
1232 : {
1233 0 : case NODEELEM:
1234 0 : return 0;
1235 :
1236 : // The 1D Szabo-Babuska defined on a two-noded edge
1237 0 : case EDGE2:
1238 0 : return 6;
1239 :
1240 : // The 1D Szabo-Babuska defined on a three-noded edge
1241 0 : case EDGE3:
1242 0 : return 0;
1243 :
1244 : // The 2D Szabo-Babuska defined on a 6-noded triangle
1245 0 : case TRI6:
1246 0 : return 15;
1247 :
1248 0 : case TRI7:
1249 0 : return 0;
1250 :
1251 : // The 2D tensor-product Szabo-Babuska defined on a
1252 : // eight-noded quadrilateral.
1253 0 : case QUAD8:
1254 0 : return 36;
1255 :
1256 : // The 2D tensor-product Szabo-Babuska defined on a
1257 : // nine-noded quadrilateral.
1258 0 : case QUAD9:
1259 0 : return 0;
1260 :
1261 0 : case INVALID_ELEM:
1262 0 : return 0;
1263 :
1264 0 : default:
1265 0 : 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 0 : default:
1272 0 : return 0;
1273 : }
1274 : } // szabab_n_dofs_per_elem
1275 :
1276 :
1277 :
1278 44390 : unsigned int szabab_n_dofs_per_elem(const Elem & e, const Order o)
1279 : {
1280 48616 : 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 576 : LIBMESH_FE_NODAL_SOLN(SZABAB, szabab_nodal_soln)
1288 384 : LIBMESH_FE_SIDE_NODAL_SOLN(SZABAB)
1289 :
1290 :
1291 : // Full specialization of n_dofs() function for every dimension
1292 0 : template <> unsigned int FE<0,SZABAB>::n_dofs(const ElemType t, const Order o) { return szabab_n_dofs(t, o); }
1293 0 : template <> unsigned int FE<1,SZABAB>::n_dofs(const ElemType t, const Order o) { return szabab_n_dofs(t, o); }
1294 0 : template <> unsigned int FE<2,SZABAB>::n_dofs(const ElemType t, const Order o) { return szabab_n_dofs(t, o); }
1295 0 : template <> unsigned int FE<3,SZABAB>::n_dofs(const ElemType t, const Order o) { return szabab_n_dofs(t, o); }
1296 :
1297 0 : template <> unsigned int FE<0,SZABAB>::n_dofs(const Elem * e, const Order o) { return szabab_n_dofs(e, o); }
1298 2436 : template <> unsigned int FE<1,SZABAB>::n_dofs(const Elem * e, const Order o) { return szabab_n_dofs(e, o); }
1299 44332 : template <> unsigned int FE<2,SZABAB>::n_dofs(const Elem * e, const Order o) { return szabab_n_dofs(e, o); }
1300 0 : 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 0 : 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 9860 : 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 248044 : 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 0 : 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 0 : 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 4050 : 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 92262 : 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 0 : 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 0 : 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 0 : 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 0 : 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 0 : 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 0 : 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 4600 : 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 44016 : 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 0 : 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 0 : template <> FEContinuity FE<0,SZABAB>::get_continuity() const { return C_ZERO; }
1326 14608 : template <> FEContinuity FE<1,SZABAB>::get_continuity() const { return C_ZERO; }
1327 49764 : template <> FEContinuity FE<2,SZABAB>::get_continuity() const { return C_ZERO; }
1328 0 : template <> FEContinuity FE<3,SZABAB>::get_continuity() const { return C_ZERO; }
1329 :
1330 : // Szabab FEMs are hierarchic
1331 0 : template <> bool FE<0,SZABAB>::is_hierarchic() const { return true; }
1332 48 : template <> bool FE<1,SZABAB>::is_hierarchic() const { return true; }
1333 202 : template <> bool FE<2,SZABAB>::is_hierarchic() const { return true; }
1334 0 : 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 2352 : void FE<2,SZABAB>::compute_constraints (DofConstraints & constraints,
1340 : DofMap & dof_map,
1341 : const unsigned int variable_number,
1342 : const Elem * elem)
1343 2352 : { compute_proj_constraints(constraints, dof_map, variable_number, elem); }
1344 :
1345 : template <>
1346 0 : void FE<3,SZABAB>::compute_constraints (DofConstraints & constraints,
1347 : DofMap & dof_map,
1348 : const unsigned int variable_number,
1349 : const Elem * elem)
1350 0 : { 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 0 : template <> bool FE<0,SZABAB>::shapes_need_reinit() const { return true; }
1356 504 : template <> bool FE<1,SZABAB>::shapes_need_reinit() const { return true; }
1357 6836 : template <> bool FE<2,SZABAB>::shapes_need_reinit() const { return true; }
1358 0 : template <> bool FE<3,SZABAB>::shapes_need_reinit() const { return true; }
1359 :
1360 : } // namespace libMesh
1361 :
1362 : #endif //LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
|
