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