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