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