libMesh
fe_monomial_shape_3D.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 // C++ includes
20 
21 // Local includes
22 #include "libmesh/fe.h"
23 #include "libmesh/elem.h"
24 
25 
26 namespace libMesh
27 {
28 
29 
30 LIBMESH_DEFAULT_VECTORIZED_FE(3,MONOMIAL)
31 
32 
33 template <>
34 Real 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 
189 template <>
190 Real 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 
204 template <>
205 Real 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 
218 template <>
219 Real 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 
658 template <>
659 Real 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 
673 template <>
674 Real 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 
689 template <>
690 Real 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 
1339 template <>
1340 Real 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 
1354 template <>
1355 Real 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::FEType
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:196
libMesh::ElemType
ElemType
Defines an enum for geometric element types.
Definition: enum_elem_type.h:33
libMesh::Order
Order
defines an enum for polynomial orders.
Definition: enum_order.h:40
libMesh::FE::shape
static OutputShape shape(const ElemType t, const Order o, const unsigned int i, const Point &p)
libMesh::Elem
This is the base class from which all geometric element types are derived.
Definition: elem.h:94
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::Elem::p_level
unsigned int p_level() const
Definition: elem.h:3108
libMesh::FEType::order
OrderWrapper order
The approximation order of the element.
Definition: fe_type.h:215
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::libmesh_ignore
void libmesh_ignore(const Args &...)
Definition: libmesh_common.h:515
libMesh::libmesh_assert
libmesh_assert(ctx)
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:144
libMesh::Elem::type
virtual ElemType type() const =0
libMesh::Point
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
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)
Generated on Thu Jul 17 2025 01:29:07 for libMesh by   doxygen 1.8.14