Line data Source code
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 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 45702 : unsigned int szabab_n_dofs(const ElemType t, const Order o)
107 : {
108 45702 : switch (o)
109 : {
110 : // Szabo-Babuska 1st-order polynomials.
111 6726 : case FIRST:
112 : {
113 6726 : switch (t)
114 : {
115 0 : case NODEELEM:
116 0 : return 1;
117 :
118 128 : case EDGE2:
119 : case EDGE3:
120 128 : return 2;
121 :
122 706 : case TRI3:
123 : case TRI6:
124 : case TRI7:
125 706 : return 3;
126 :
127 596 : case QUAD4:
128 : case QUAD8:
129 : case QUAD9:
130 596 : 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 13812 : case SECOND:
143 : {
144 13812 : switch (t)
145 : {
146 0 : case NODEELEM:
147 0 : return 1;
148 :
149 164 : case EDGE2:
150 : case EDGE3:
151 164 : return 3;
152 :
153 2280 : case TRI6:
154 : case TRI7:
155 2280 : return 6;
156 :
157 0 : case QUAD8:
158 0 : return 8;
159 776 : case QUAD9:
160 776 : 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 12582 : case THIRD:
173 : {
174 12582 : switch (t)
175 : {
176 0 : case NODEELEM:
177 0 : return 1;
178 :
179 164 : case EDGE2:
180 : case EDGE3:
181 164 : return 4;
182 :
183 2006 : case TRI6:
184 : case TRI7:
185 2006 : return 10;
186 :
187 3756 : case QUAD8:
188 : case QUAD9:
189 3756 : 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 12582 : case FOURTH:
202 : {
203 12582 : switch (t)
204 : {
205 0 : case NODEELEM:
206 0 : return 1;
207 :
208 164 : case EDGE2:
209 : case EDGE3:
210 164 : return 5;
211 :
212 8202 : case TRI6:
213 : case TRI7:
214 8202 : return 15;
215 :
216 3756 : case QUAD8:
217 : case QUAD9:
218 3756 : 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 42023 : unsigned int szabab_n_dofs(const Elem * e, const Order o)
323 : {
324 6863 : libmesh_assert(e);
325 45702 : return szabab_n_dofs(e->type(), o);
326 : }
327 :
328 :
329 :
330 321466 : unsigned int szabab_n_dofs_at_node(const ElemType t,
331 : const Order o,
332 : const unsigned int n)
333 : {
334 321466 : switch (o)
335 : {
336 : // The first-order Szabo-Babuska shape functions
337 16002 : case FIRST:
338 : {
339 16002 : switch (t)
340 : {
341 0 : case NODEELEM:
342 0 : return 1;
343 :
344 : // The 1D Szabo-Babuska defined on a three-noded edge
345 1908 : case EDGE2:
346 : case EDGE3:
347 : {
348 1908 : switch (n)
349 : {
350 189 : case 0:
351 : case 1:
352 189 : 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 8658 : case TRI3:
365 : case TRI6:
366 : case TRI7:
367 : {
368 7866 : switch (n)
369 : {
370 792 : case 0:
371 : case 1:
372 : case 2:
373 792 : 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 5436 : case QUAD4:
393 : case QUAD8:
394 : case QUAD9:
395 : {
396 4923 : switch (n)
397 : {
398 513 : case 0:
399 : case 1:
400 : case 2:
401 : case 3:
402 513 : 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 115888 : case SECOND:
430 : {
431 115888 : switch (t)
432 : {
433 0 : case NODEELEM:
434 0 : return 1;
435 :
436 : // The 1D Szabo-Babuska defined on a three-noded edge
437 3537 : case EDGE2:
438 : case EDGE3:
439 : {
440 3537 : switch (n)
441 : {
442 225 : case 0:
443 : case 1:
444 225 : return 1;
445 :
446 108 : case 2:
447 108 : 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 87106 : case TRI6:
457 : case TRI7:
458 : {
459 87106 : switch (n)
460 : {
461 3466 : case 0:
462 : case 1:
463 : case 2:
464 3466 : return 1;
465 :
466 3134 : case 3:
467 : case 4:
468 : case 5:
469 3134 : return 1;
470 :
471 810 : case 6:
472 405 : libmesh_assert(t == TRI7);
473 810 : 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 25245 : case QUAD8:
484 : case QUAD9:
485 : {
486 25245 : switch (n)
487 : {
488 1017 : case 0:
489 : case 1:
490 : case 2:
491 : case 3:
492 1017 : return 1;
493 :
494 945 : case 4:
495 : case 5:
496 : case 6:
497 : case 7:
498 945 : return 1;
499 :
500 243 : case 8:
501 243 : 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 94788 : case THIRD:
520 : {
521 94788 : switch (t)
522 : {
523 0 : case NODEELEM:
524 0 : return 1;
525 :
526 : // The 1D Szabo-Babuska defined on a three-noded edge
527 3537 : case EDGE2:
528 : case EDGE3:
529 : {
530 3537 : switch (n)
531 : {
532 225 : case 0:
533 : case 1:
534 225 : return 1;
535 :
536 216 : case 2:
537 216 : 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 35496 : case TRI7:
547 35496 : if (n == 6)
548 441 : return 1;
549 : libmesh_fallthrough();
550 : case TRI6:
551 : {
552 61020 : switch (n)
553 : {
554 2826 : case 0:
555 : case 1:
556 : case 2:
557 2826 : return 1;
558 :
559 5184 : case 3:
560 : case 4:
561 : case 5:
562 5184 : 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 25245 : case QUAD8:
573 : case QUAD9:
574 : {
575 23040 : switch (n)
576 : {
577 1017 : case 0:
578 : case 1:
579 : case 2:
580 : case 3:
581 1017 : return 1;
582 :
583 945 : case 4:
584 : case 5:
585 : case 6:
586 : case 7:
587 945 : return 2;
588 :
589 243 : case 8:
590 243 : 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 94788 : case FOURTH:
609 : {
610 94788 : switch (t)
611 : {
612 0 : case NODEELEM:
613 0 : return 1;
614 :
615 : // The 1D Szabo-Babuska defined on a three-noded edge
616 3537 : case EDGE2:
617 : case EDGE3:
618 : {
619 3537 : switch (n)
620 : {
621 225 : case 0:
622 : case 1:
623 225 : return 1;
624 :
625 216 : case 2:
626 216 : 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 35496 : case TRI7:
636 35496 : if (n == 6)
637 441 : return 3;
638 : libmesh_fallthrough();
639 : case TRI6:
640 : {
641 61020 : switch (n)
642 : {
643 2826 : case 0:
644 : case 1:
645 : case 2:
646 2826 : return 1;
647 :
648 5184 : case 3:
649 : case 4:
650 : case 5:
651 5184 : 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 25245 : case QUAD8:
662 : case QUAD9:
663 : {
664 23040 : switch (n)
665 : {
666 1017 : case 0:
667 : case 1:
668 : case 2:
669 : case 3:
670 1017 : return 1;
671 :
672 945 : case 4:
673 : case 5:
674 : case 6:
675 : case 7:
676 945 : return 3;
677 :
678 243 : case 8:
679 243 : 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 79972 : unsigned int szabab_n_dofs_at_node(const Elem & e,
969 : const Order o,
970 : const unsigned int n)
971 : {
972 87700 : return szabab_n_dofs_at_node(e.type(), o, n);
973 : }
974 :
975 :
976 :
977 44122 : unsigned int szabab_n_dofs_per_elem(const ElemType t, const Order o)
978 : {
979 44122 : switch (o)
980 : {
981 : // The first-order Szabo-Babuska shape functions
982 4752 : case FIRST:
983 : {
984 4320 : switch (t)
985 : {
986 0 : case NODEELEM:
987 0 : return 0;
988 :
989 : // The 1D Szabo-Babuska defined on an edge
990 81 : case EDGE2:
991 : case EDGE3:
992 81 : return 0;
993 :
994 : // The 2D Szabo-Babuska defined on a triangle
995 234 : case TRI3:
996 : case TRI6:
997 : case TRI7:
998 234 : return 0;
999 :
1000 : // The 2D tensor-product Szabo-Babuska defined on a
1001 : // quadrilateral.
1002 117 : case QUAD4:
1003 : case QUAD8:
1004 : case QUAD9:
1005 117 : 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 14998 : case SECOND:
1019 : {
1020 13764 : 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 99 : case EDGE3:
1031 99 : return 0;
1032 :
1033 : // The 2D Szabo-Babuska defined on a 6-noded triangle
1034 928 : case TRI6:
1035 : case TRI7:
1036 928 : 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 207 : case QUAD9:
1046 207 : 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 12186 : case THIRD:
1060 : {
1061 12186 : 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 99 : case EDGE3:
1072 99 : return 0;
1073 :
1074 : // The 2D Szabo-Babuska defined on a 6-noded triangle
1075 774 : case TRI6:
1076 774 : return 1;
1077 :
1078 387 : case TRI7:
1079 387 : 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 207 : case QUAD9:
1089 207 : 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 12186 : case FOURTH:
1103 : {
1104 12186 : 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 99 : case EDGE3:
1115 99 : return 0;
1116 :
1117 : // The 2D Szabo-Babuska defined on a 6-noded triangle
1118 774 : case TRI6:
1119 774 : return 3;
1120 :
1121 387 : case TRI7:
1122 387 : 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 207 : case QUAD9:
1132 207 : 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 40296 : unsigned int szabab_n_dofs_per_elem(const Elem & e, const Order o)
1279 : {
1280 44122 : 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 2370 : template <> unsigned int FE<1,SZABAB>::n_dofs(const Elem * e, const Order o) { return szabab_n_dofs(e, o); }
1299 43332 : 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 8874 : 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 224892 : 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 3645 : 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 84055 : 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 4140 : 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 39982 : 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 13260 : template <> FEContinuity FE<1,SZABAB>::get_continuity() const { return C_ZERO; }
1327 45874 : 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 2136 : void FE<2,SZABAB>::compute_constraints (DofConstraints & constraints,
1340 : DofMap & dof_map,
1341 : const unsigned int variable_number,
1342 : const Elem * elem)
1343 2136 : { 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 468 : template <> bool FE<1,SZABAB>::shapes_need_reinit() const { return true; }
1357 5883 : 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
|
