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