fe_monomial_shape_3D.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2014 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++ inlcludes
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();
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 
628 #endif
629 
630  libmesh_error();
631  return 0.;
632 }
633 
634 
635 
636 template <>
638  const Order order,
639  const unsigned int i,
640  const unsigned int j,
641  const Point& p)
642 {
643  libmesh_assert(elem);
644 
645  // call the orientation-independent shape function derivatives
646  return FE<3,MONOMIAL>::shape_deriv(elem->type(), static_cast<Order>(order + elem->p_level()), i, j, p);
647 }
648 
649 
650 
651 template <>
653  const Order libmesh_dbg_var(order),
654  const unsigned int i,
655  const unsigned int j,
656  const Point& p)
657 {
658 #if LIBMESH_DIM == 3
659 
660  libmesh_assert_less (j, 6);
661 
662  libmesh_assert_less (i, (static_cast<unsigned int>(order)+1)*
663  (static_cast<unsigned int>(order)+2)*
664  (static_cast<unsigned int>(order)+3)/6);
665 
666  const Real xi = p(0);
667  const Real eta = p(1);
668  const Real zeta = p(2);
669 
670  // monomials. since they are hierarchic we only need one case block.
671  switch (j)
672  {
673  // d^2()/dxi^2
674  case 0:
675  {
676  switch (i)
677  {
678  // constant
679  case 0:
680 
681  // linear
682  case 1:
683  case 2:
684  case 3:
685  return 0.;
686 
687  // quadratic
688  case 4:
689  return 2.;
690 
691  case 5:
692  case 6:
693  case 7:
694  case 8:
695  case 9:
696  return 0.;
697 
698  // cubic
699  case 10:
700  return 6.*xi;
701 
702  case 11:
703  return 2.*eta;
704 
705  case 12:
706  case 13:
707  return 0.;
708 
709  case 14:
710  return 2.*zeta;
711 
712  case 15:
713  case 16:
714  case 17:
715  case 18:
716  case 19:
717  return 0.;
718 
719  // quartics
720  case 20:
721  return 12.*xi*xi;
722 
723  case 21:
724  return 6.*xi*eta;
725 
726  case 22:
727  return 2.*eta*eta;
728 
729  case 23:
730  case 24:
731  return 0.;
732 
733  case 25:
734  return 6.*xi*zeta;
735 
736  case 26:
737  return 2.*eta*zeta;
738 
739  case 27:
740  case 28:
741  return 0.;
742 
743  case 29:
744  return 2.*zeta*zeta;
745 
746  case 30:
747  case 31:
748  case 32:
749  case 33:
750  case 34:
751  return 0.;
752 
753  default:
754  unsigned int o = 0;
755  for (; i >= (o+1)*(o+2)*(o+3)/6; o++) { }
756  unsigned int i2 = i - (o*(o+1)*(o+2)/6);
757  unsigned int block=o, nz = 0;
758  for (; block < i2; block += (o-nz+1)) { nz++; }
759  const unsigned int nx = block - i2;
760  const unsigned int ny = o - nx - nz;
761  Real val = nx * (nx - 1);
762  for (unsigned int index=2; index < nx; index++)
763  val *= xi;
764  for (unsigned int index=0; index != ny; index++)
765  val *= eta;
766  for (unsigned int index=0; index != nz; index++)
767  val *= zeta;
768  return val;
769  }
770  }
771 
772 
773  // d^2()/dxideta
774  case 1:
775  {
776  switch (i)
777  {
778  // constant
779  case 0:
780 
781  // linear
782  case 1:
783  case 2:
784  case 3:
785  return 0.;
786 
787  // quadratic
788  case 4:
789  return 0.;
790 
791  case 5:
792  return 1.;
793 
794  case 6:
795  case 7:
796  case 8:
797  case 9:
798  return 0.;
799 
800  // cubic
801  case 10:
802  return 0.;
803 
804  case 11:
805  return 2.*xi;
806 
807  case 12:
808  return 2.*eta;
809 
810  case 13:
811  case 14:
812  return 0.;
813 
814  case 15:
815  return zeta;
816 
817  case 16:
818  case 17:
819  case 18:
820  case 19:
821  return 0.;
822 
823  // quartics
824  case 20:
825  return 0.;
826 
827  case 21:
828  return 3.*xi*xi;
829 
830  case 22:
831  return 4.*xi*eta;
832 
833  case 23:
834  return 3.*eta*eta;
835 
836  case 24:
837  case 25:
838  return 0.;
839 
840  case 26:
841  return 2.*xi*zeta;
842 
843  case 27:
844  return 2.*eta*zeta;
845 
846  case 28:
847  case 29:
848  return 0.;
849 
850  case 30:
851  return zeta*zeta;
852 
853  case 31:
854  case 32:
855  case 33:
856  case 34:
857  return 0.;
858 
859  default:
860  unsigned int o = 0;
861  for (; i >= (o+1)*(o+2)*(o+3)/6; o++) { }
862  unsigned int i2 = i - (o*(o+1)*(o+2)/6);
863  unsigned int block=o, nz = 0;
864  for (; block < i2; block += (o-nz+1)) { nz++; }
865  const unsigned int nx = block - i2;
866  const unsigned int ny = o - nx - nz;
867  Real val = nx * ny;
868  for (unsigned int index=1; index < nx; index++)
869  val *= xi;
870  for (unsigned int index=1; index < ny; index++)
871  val *= eta;
872  for (unsigned int index=0; index != nz; index++)
873  val *= zeta;
874  return val;
875  }
876  }
877 
878 
879  // d^2()/deta^2
880  case 2:
881  {
882  switch (i)
883  {
884  // constant
885  case 0:
886 
887  // linear
888  case 1:
889  case 2:
890  case 3:
891  return 0.;
892 
893  // quadratic
894  case 4:
895  case 5:
896  return 0.;
897 
898  case 6:
899  return 2.;
900 
901  case 7:
902  case 8:
903  case 9:
904  return 0.;
905 
906  // cubic
907  case 10:
908  case 11:
909  return 0.;
910 
911  case 12:
912  return 2.*xi;
913  case 13:
914  return 6.*eta;
915 
916  case 14:
917  case 15:
918  return 0.;
919 
920  case 16:
921  return 2.*zeta;
922 
923  case 17:
924  case 18:
925  case 19:
926  return 0.;
927 
928  // quartics
929  case 20:
930  case 21:
931  return 0.;
932 
933  case 22:
934  return 2.*xi*xi;
935 
936  case 23:
937  return 6.*xi*eta;
938 
939  case 24:
940  return 12.*eta*eta;
941 
942  case 25:
943  case 26:
944  return 0.;
945 
946  case 27:
947  return 2.*xi*zeta;
948 
949  case 28:
950  return 6.*eta*zeta;
951 
952  case 29:
953  case 30:
954  return 0.;
955 
956  case 31:
957  return 2.*zeta*zeta;
958 
959  case 32:
960  case 33:
961  case 34:
962  return 0.;
963 
964  default:
965  unsigned int o = 0;
966  for (; i >= (o+1)*(o+2)*(o+3)/6; o++) { }
967  unsigned int i2 = i - (o*(o+1)*(o+2)/6);
968  unsigned int block=o, nz = 0;
969  for (; block < i2; block += (o-nz+1)) { nz++; }
970  const unsigned int nx = block - i2;
971  const unsigned int ny = o - nx - nz;
972  Real val = ny * (ny - 1);
973  for (unsigned int index=0; index != nx; index++)
974  val *= xi;
975  for (unsigned int index=2; index < ny; index++)
976  val *= eta;
977  for (unsigned int index=0; index != nz; index++)
978  val *= zeta;
979  return val;
980  }
981  }
982 
983 
984  // d^2()/dxidzeta
985  case 3:
986  {
987  switch (i)
988  {
989  // constant
990  case 0:
991 
992  // linear
993  case 1:
994  case 2:
995  case 3:
996  return 0.;
997 
998  // quadratic
999  case 4:
1000  case 5:
1001  case 6:
1002  return 0.;
1003 
1004  case 7:
1005  return 1.;
1006 
1007  case 8:
1008  case 9:
1009  return 0.;
1010 
1011  // cubic
1012  case 10:
1013  case 11:
1014  case 12:
1015  case 13:
1016  return 0.;
1017 
1018  case 14:
1019  return 2.*xi;
1020 
1021  case 15:
1022  return eta;
1023 
1024  case 16:
1025  return 0.;
1026 
1027  case 17:
1028  return 2.*zeta;
1029 
1030  case 18:
1031  case 19:
1032  return 0.;
1033 
1034  // quartics
1035  case 20:
1036  case 21:
1037  case 22:
1038  case 23:
1039  case 24:
1040  return 0.;
1041 
1042  case 25:
1043  return 3.*xi*xi;
1044 
1045  case 26:
1046  return 2.*xi*eta;
1047 
1048  case 27:
1049  return eta*eta;
1050 
1051  case 28:
1052  return 0.;
1053 
1054  case 29:
1055  return 4.*xi*zeta;
1056 
1057  case 30:
1058  return 2.*eta*zeta;
1059 
1060  case 31:
1061  return 0.;
1062 
1063  case 32:
1064  return 3.*zeta*zeta;
1065 
1066  case 33:
1067  case 34:
1068  return 0.;
1069 
1070  default:
1071  unsigned int o = 0;
1072  for (; i >= (o+1)*(o+2)*(o+3)/6; o++) { }
1073  unsigned int i2 = i - (o*(o+1)*(o+2)/6);
1074  unsigned int block=o, nz = 0;
1075  for (; block < i2; block += (o-nz+1)) { nz++; }
1076  const unsigned int nx = block - i2;
1077  const unsigned int ny = o - nx - nz;
1078  Real val = nx * nz;
1079  for (unsigned int index=1; index < nx; index++)
1080  val *= xi;
1081  for (unsigned int index=0; index != ny; index++)
1082  val *= eta;
1083  for (unsigned int index=1; index < nz; index++)
1084  val *= zeta;
1085  return val;
1086  }
1087  }
1088 
1089  // d^2()/detadzeta
1090  case 4:
1091  {
1092  switch (i)
1093  {
1094  // constant
1095  case 0:
1096 
1097  // linear
1098  case 1:
1099  case 2:
1100  case 3:
1101  return 0.;
1102 
1103  // quadratic
1104  case 4:
1105  case 5:
1106  case 6:
1107  case 7:
1108  return 0.;
1109 
1110  case 8:
1111  return 1.;
1112 
1113  case 9:
1114  return 0.;
1115 
1116  // cubic
1117  case 10:
1118  case 11:
1119  case 12:
1120  case 13:
1121  case 14:
1122  return 0.;
1123 
1124  case 15:
1125  return xi;
1126 
1127  case 16:
1128  return 2.*eta;
1129 
1130  case 17:
1131  return 0.;
1132 
1133  case 18:
1134  return 2.*zeta;
1135 
1136  case 19:
1137  return 0.;
1138 
1139  // quartics
1140  case 20:
1141  case 21:
1142  case 22:
1143  case 23:
1144  case 24:
1145  case 25:
1146  return 0.;
1147 
1148  case 26:
1149  return xi*xi;
1150 
1151  case 27:
1152  return 2.*xi*eta;
1153 
1154  case 28:
1155  return 3.*eta*eta;
1156 
1157  case 29:
1158  return 0.;
1159 
1160  case 30:
1161  return 2.*xi*zeta;
1162 
1163  case 31:
1164  return 4.*eta*zeta;
1165 
1166  case 32:
1167  return 0.;
1168 
1169  case 33:
1170  return 3.*zeta*zeta;
1171 
1172  case 34:
1173  return 0.;
1174 
1175  default:
1176  unsigned int o = 0;
1177  for (; i >= (o+1)*(o+2)*(o+3)/6; o++) { }
1178  unsigned int i2 = i - (o*(o+1)*(o+2)/6);
1179  unsigned int block=o, nz = 0;
1180  for (; block < i2; block += (o-nz+1)) { nz++; }
1181  const unsigned int nx = block - i2;
1182  const unsigned int ny = o - nx - nz;
1183  Real val = ny * nz;
1184  for (unsigned int index=0; index != nx; index++)
1185  val *= xi;
1186  for (unsigned int index=1; index < ny; index++)
1187  val *= eta;
1188  for (unsigned int index=1; index < nz; index++)
1189  val *= zeta;
1190  return val;
1191  }
1192  }
1193 
1194 
1195  // d^2()/dzeta^2
1196  case 5:
1197  {
1198  switch (i)
1199  {
1200  // constant
1201  case 0:
1202 
1203  // linear
1204  case 1:
1205  case 2:
1206  case 3:
1207  return 0.;
1208 
1209  // quadratic
1210  case 4:
1211  case 5:
1212  case 6:
1213  case 7:
1214  case 8:
1215  return 0.;
1216 
1217  case 9:
1218  return 2.;
1219 
1220  // cubic
1221  case 10:
1222  case 11:
1223  case 12:
1224  case 13:
1225  case 14:
1226  case 15:
1227  case 16:
1228  return 0.;
1229 
1230  case 17:
1231  return 2.*xi;
1232 
1233  case 18:
1234  return 2.*eta;
1235 
1236  case 19:
1237  return 6.*zeta;
1238 
1239  // quartics
1240  case 20:
1241  case 21:
1242  case 22:
1243  case 23:
1244  case 24:
1245  case 25:
1246  case 26:
1247  case 27:
1248  case 28:
1249  return 0.;
1250 
1251  case 29:
1252  return 2.*xi*xi;
1253 
1254  case 30:
1255  return 2.*xi*eta;
1256 
1257  case 31:
1258  return 2.*eta*eta;
1259 
1260  case 32:
1261  return 6.*xi*zeta;
1262 
1263  case 33:
1264  return 6.*eta*zeta;
1265 
1266  case 34:
1267  return 12.*zeta*zeta;
1268 
1269  default:
1270  unsigned int o = 0;
1271  for (; i >= (o+1)*(o+2)*(o+3)/6; o++) { }
1272  unsigned int i2 = i - (o*(o+1)*(o+2)/6);
1273  unsigned int block=o, nz = 0;
1274  for (; block < i2; block += (o-nz+1)) { nz++; }
1275  const unsigned int nx = block - i2;
1276  const unsigned int ny = o - nx - nz;
1277  Real val = nz * (nz - 1);
1278  for (unsigned int index=0; index != nx; index++)
1279  val *= xi;
1280  for (unsigned int index=0; index != ny; index++)
1281  val *= eta;
1282  for (unsigned int index=2; index < nz; index++)
1283  val *= zeta;
1284  return val;
1285  }
1286  }
1287 
1288 
1289  default:
1290  libmesh_error();
1291  }
1292 
1293 #endif
1294 
1295  libmesh_error();
1296  return 0.;
1297 }
1298 
1299 
1300 
1301 template <>
1303  const Order order,
1304  const unsigned int i,
1305  const unsigned int j,
1306  const Point& p)
1307 {
1308  libmesh_assert(elem);
1309 
1310  // call the orientation-independent shape function derivatives
1311  return FE<3,MONOMIAL>::shape_second_deriv(elem->type(), static_cast<Order>(order + elem->p_level()), i, j, p);
1312 }
1313 
1314 } // namespace libMesh

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:05 UTC

Hosted By:
SourceForge.net Logo