libMesh
fe_lagrange_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 // Local includes
20 #include "libmesh/fe.h"
21 #include "libmesh/elem.h"
22 #include "libmesh/fe_lagrange_shape_1D.h"
23 
24 // Anonymous namespace for functions shared by LAGRANGE and
25 // L2_LAGRANGE implementations. Implementations appear at the bottom
26 // of this file.
27 namespace
28 {
29 using namespace libMesh;
30 
31 Real fe_lagrange_3D_shape(const ElemType,
32  const Order order,
33  const unsigned int i,
34  const Point & p);
35 
36 Real fe_lagrange_3D_shape_deriv(const ElemType type,
37  const Order order,
38  const unsigned int i,
39  const unsigned int j,
40  const Point & p);
41 
42 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
43 
44 Real fe_lagrange_3D_shape_second_deriv(const ElemType type,
45  const Order order,
46  const unsigned int i,
47  const unsigned int j,
48  const Point & p);
49 
50 #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
51 
52 } // anonymous namespace
53 
54 namespace libMesh
55 {
56 
57 template <>
58 Real FE<3,LAGRANGE>::shape(const ElemType type,
59  const Order order,
60  const unsigned int i,
61  const Point & p)
62 {
63  return fe_lagrange_3D_shape(type, order, i, p);
64 }
65 
66 
67 
68 template <>
69 Real FE<3,L2_LAGRANGE>::shape(const ElemType type,
70  const Order order,
71  const unsigned int i,
72  const Point & p)
73 {
74  return fe_lagrange_3D_shape(type, order, i, p);
75 }
76 
77 
78 
79 template <>
80 Real FE<3,LAGRANGE>::shape(const Elem * elem,
81  const Order order,
82  const unsigned int i,
83  const Point & p,
84  const bool add_p_level)
85 {
86  libmesh_assert(elem);
87 
88  // call the orientation-independent shape functions
89  return fe_lagrange_3D_shape(elem->type(), static_cast<Order>(order + add_p_level * elem->p_level()), i, p);
90 }
91 
92 
93 
94 template <>
95 Real FE<3,L2_LAGRANGE>::shape(const Elem * elem,
96  const Order order,
97  const unsigned int i,
98  const Point & p,
99  const bool add_p_level)
100 {
101  libmesh_assert(elem);
102 
103  // call the orientation-independent shape functions
104  return fe_lagrange_3D_shape(elem->type(), static_cast<Order>(order + add_p_level * elem->p_level()), i, p);
105 }
106 
107 
108 
109 template <>
110 Real FE<3,LAGRANGE>::shape_deriv(const ElemType type,
111  const Order order,
112  const unsigned int i,
113  const unsigned int j,
114  const Point & p)
115 {
116  return fe_lagrange_3D_shape_deriv(type, order, i, j, p);
117 }
118 
119 
120 
121 template <>
122 Real FE<3,L2_LAGRANGE>::shape_deriv(const ElemType type,
123  const Order order,
124  const unsigned int i,
125  const unsigned int j,
126  const Point & p)
127 {
128  return fe_lagrange_3D_shape_deriv(type, order, i, j, p);
129 }
130 
131 
132 
133 template <>
134 Real FE<3,LAGRANGE>::shape_deriv(const Elem * elem,
135  const Order order,
136  const unsigned int i,
137  const unsigned int j,
138  const Point & p,
139  const bool add_p_level)
140 {
141  libmesh_assert(elem);
142 
143  // call the orientation-independent shape function derivatives
144  return fe_lagrange_3D_shape_deriv(elem->type(), static_cast<Order>(order + add_p_level * elem->p_level()), i, j, p);
145 }
146 
147 
148 template <>
149 Real FE<3,L2_LAGRANGE>::shape_deriv(const Elem * elem,
150  const Order order,
151  const unsigned int i,
152  const unsigned int j,
153  const Point & p,
154  const bool add_p_level)
155 {
156  libmesh_assert(elem);
157 
158  // call the orientation-independent shape function derivatives
159  return fe_lagrange_3D_shape_deriv(elem->type(), static_cast<Order>(order + add_p_level * elem->p_level()), i, j, p);
160 }
161 
162 
163 
164 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
165 
166 template <>
167 Real FE<3,LAGRANGE>::shape_second_deriv(const ElemType type,
168  const Order order,
169  const unsigned int i,
170  const unsigned int j,
171  const Point & p)
172 {
173  return fe_lagrange_3D_shape_second_deriv(type, order, i, j, p);
174 }
175 
176 
177 
178 template <>
179 Real FE<3,L2_LAGRANGE>::shape_second_deriv(const ElemType type,
180  const Order order,
181  const unsigned int i,
182  const unsigned int j,
183  const Point & p)
184 {
185  return fe_lagrange_3D_shape_second_deriv(type, order, i, j, p);
186 }
187 
188 
189 
190 template <>
191 Real FE<3,LAGRANGE>::shape_second_deriv(const Elem * elem,
192  const Order order,
193  const unsigned int i,
194  const unsigned int j,
195  const Point & p,
196  const bool add_p_level)
197 {
198  libmesh_assert(elem);
199 
200  // call the orientation-independent shape function derivatives
201  return fe_lagrange_3D_shape_second_deriv
202  (elem->type(), static_cast<Order>(order + add_p_level * elem->p_level()), i, j, p);
203 }
204 
205 
206 
207 template <>
208 Real FE<3,L2_LAGRANGE>::shape_second_deriv(const Elem * elem,
209  const Order order,
210  const unsigned int i,
211  const unsigned int j,
212  const Point & p,
213  const bool add_p_level)
214 {
215  libmesh_assert(elem);
216 
217  // call the orientation-independent shape function derivatives
218  return fe_lagrange_3D_shape_second_deriv
219  (elem->type(), static_cast<Order>(order + add_p_level * elem->p_level()), i, j, p);
220 }
221 
222 #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
223 
224 } // namespace libMesh
225 
226 
227 
228 namespace
229 {
230 using namespace libMesh;
231 
232 Real fe_lagrange_3D_shape(const ElemType type,
233  const Order order,
234  const unsigned int i,
235  const Point & p)
236 {
237 #if LIBMESH_DIM == 3
238 
239  switch (order)
240  {
241  // linear Lagrange shape functions
242  case FIRST:
243  {
244  switch (type)
245  {
246  // trilinear hexahedral shape functions
247  case HEX8:
248  case HEX20:
249  case HEX27:
250  {
251  libmesh_assert_less (i, 8);
252 
253  // Compute hex shape functions as a tensor-product
254  const Real xi = p(0);
255  const Real eta = p(1);
256  const Real zeta = p(2);
257 
258  // 0 1 2 3 4 5 6 7
259  static const unsigned int i0[] = {0, 1, 1, 0, 0, 1, 1, 0};
260  static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 1, 1};
261  static const unsigned int i2[] = {0, 0, 0, 0, 1, 1, 1, 1};
262 
263  return (fe_lagrange_1D_linear_shape(i0[i], xi)*
264  fe_lagrange_1D_linear_shape(i1[i], eta)*
265  fe_lagrange_1D_linear_shape(i2[i], zeta));
266  }
267 
268  // linear tetrahedral shape functions
269  case TET4:
270  case TET10:
271  {
272  libmesh_assert_less (i, 4);
273 
274  // Area coordinates, pg. 205, Vol. I, Carey, Oden, Becker FEM
275  const Real zeta1 = p(0);
276  const Real zeta2 = p(1);
277  const Real zeta3 = p(2);
278  const Real zeta0 = 1. - zeta1 - zeta2 - zeta3;
279 
280  switch(i)
281  {
282  case 0:
283  return zeta0;
284 
285  case 1:
286  return zeta1;
287 
288  case 2:
289  return zeta2;
290 
291  case 3:
292  return zeta3;
293 
294  default:
295  libmesh_error_msg("Invalid i = " << i);
296  }
297  }
298 
299  // linear prism shape functions
300  case PRISM6:
301  case PRISM15:
302  case PRISM18:
303  {
304  libmesh_assert_less (i, 6);
305 
306  // Compute prism shape functions as a tensor-product
307  // of a triangle and an edge
308 
309  Point p2d(p(0),p(1));
310  Real p1d = p(2);
311 
312  // 0 1 2 3 4 5
313  static const unsigned int i0[] = {0, 0, 0, 1, 1, 1};
314  static const unsigned int i1[] = {0, 1, 2, 0, 1, 2};
315 
316  return (FE<2,LAGRANGE>::shape(TRI3, FIRST, i1[i], p2d)*
317  fe_lagrange_1D_linear_shape(i0[i], p1d));
318  }
319 
320  // linear pyramid shape functions
321  case PYRAMID5:
322  case PYRAMID13:
323  case PYRAMID14:
324  {
325  libmesh_assert_less (i, 5);
326 
327  const Real xi = p(0);
328  const Real eta = p(1);
329  const Real zeta = p(2);
330  const Real eps = 1.e-35;
331 
332  switch(i)
333  {
334  case 0:
335  return .25*(zeta + xi - 1.)*(zeta + eta - 1.)/((1. - zeta) + eps);
336 
337  case 1:
338  return .25*(zeta - xi - 1.)*(zeta + eta - 1.)/((1. - zeta) + eps);
339 
340  case 2:
341  return .25*(zeta - xi - 1.)*(zeta - eta - 1.)/((1. - zeta) + eps);
342 
343  case 3:
344  return .25*(zeta + xi - 1.)*(zeta - eta - 1.)/((1. - zeta) + eps);
345 
346  case 4:
347  return zeta;
348 
349  default:
350  libmesh_error_msg("Invalid i = " << i);
351  }
352  }
353 
354 
355  default:
356  libmesh_error_msg("ERROR: Unsupported 3D element type!: " << type);
357  }
358  }
359 
360 
361  // quadratic Lagrange shape functions
362  case SECOND:
363  {
364  switch (type)
365  {
366 
367  // serendipity hexahedral quadratic shape functions
368  case HEX20:
369  {
370  libmesh_assert_less (i, 20);
371 
372  const Real xi = p(0);
373  const Real eta = p(1);
374  const Real zeta = p(2);
375 
376  // these functions are defined for (x,y,z) in [0,1]^3
377  // so transform the locations
378  const Real x = .5*(xi + 1.);
379  const Real y = .5*(eta + 1.);
380  const Real z = .5*(zeta + 1.);
381 
382  switch (i)
383  {
384  case 0:
385  return (1. - x)*(1. - y)*(1. - z)*(1. - 2.*x - 2.*y - 2.*z);
386 
387  case 1:
388  return x*(1. - y)*(1. - z)*(2.*x - 2.*y - 2.*z - 1.);
389 
390  case 2:
391  return x*y*(1. - z)*(2.*x + 2.*y - 2.*z - 3.);
392 
393  case 3:
394  return (1. - x)*y*(1. - z)*(2.*y - 2.*x - 2.*z - 1.);
395 
396  case 4:
397  return (1. - x)*(1. - y)*z*(2.*z - 2.*x - 2.*y - 1.);
398 
399  case 5:
400  return x*(1. - y)*z*(2.*x - 2.*y + 2.*z - 3.);
401 
402  case 6:
403  return x*y*z*(2.*x + 2.*y + 2.*z - 5.);
404 
405  case 7:
406  return (1. - x)*y*z*(2.*y - 2.*x + 2.*z - 3.);
407 
408  case 8:
409  return 4.*x*(1. - x)*(1. - y)*(1. - z);
410 
411  case 9:
412  return 4.*x*y*(1. - y)*(1. - z);
413 
414  case 10:
415  return 4.*x*(1. - x)*y*(1. - z);
416 
417  case 11:
418  return 4.*(1. - x)*y*(1. - y)*(1. - z);
419 
420  case 12:
421  return 4.*(1. - x)*(1. - y)*z*(1. - z);
422 
423  case 13:
424  return 4.*x*(1. - y)*z*(1. - z);
425 
426  case 14:
427  return 4.*x*y*z*(1. - z);
428 
429  case 15:
430  return 4.*(1. - x)*y*z*(1. - z);
431 
432  case 16:
433  return 4.*x*(1. - x)*(1. - y)*z;
434 
435  case 17:
436  return 4.*x*y*(1. - y)*z;
437 
438  case 18:
439  return 4.*x*(1. - x)*y*z;
440 
441  case 19:
442  return 4.*(1. - x)*y*(1. - y)*z;
443 
444  default:
445  libmesh_error_msg("Invalid i = " << i);
446  }
447  }
448 
449  // triquadratic hexahedral shape functions
450  case HEX27:
451  {
452  libmesh_assert_less (i, 27);
453 
454  // Compute hex shape functions as a tensor-product
455  const Real xi = p(0);
456  const Real eta = p(1);
457  const Real zeta = p(2);
458 
459  // The only way to make any sense of this
460  // is to look at the mgflo/mg2/mgf documentation
461  // and make the cut-out cube!
462  // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26
463  static const unsigned int i0[] = {0, 1, 1, 0, 0, 1, 1, 0, 2, 1, 2, 0, 0, 1, 1, 0, 2, 1, 2, 0, 2, 2, 1, 2, 0, 2, 2};
464  static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 1, 1, 0, 2, 1, 2, 0, 0, 1, 1, 0, 2, 1, 2, 2, 0, 2, 1, 2, 2, 2};
465  static const unsigned int i2[] = {0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 2, 2, 2, 2, 1, 1, 1, 1, 0, 2, 2, 2, 2, 1, 2};
466 
467  return (fe_lagrange_1D_quadratic_shape(i0[i], xi)*
468  fe_lagrange_1D_quadratic_shape(i1[i], eta)*
469  fe_lagrange_1D_quadratic_shape(i2[i], zeta));
470  }
471 
472  // quadratic tetrahedral shape functions
473  case TET10:
474  {
475  libmesh_assert_less (i, 10);
476 
477  // Area coordinates, pg. 205, Vol. I, Carey, Oden, Becker FEM
478  const Real zeta1 = p(0);
479  const Real zeta2 = p(1);
480  const Real zeta3 = p(2);
481  const Real zeta0 = 1. - zeta1 - zeta2 - zeta3;
482 
483  switch(i)
484  {
485  case 0:
486  return zeta0*(2.*zeta0 - 1.);
487 
488  case 1:
489  return zeta1*(2.*zeta1 - 1.);
490 
491  case 2:
492  return zeta2*(2.*zeta2 - 1.);
493 
494  case 3:
495  return zeta3*(2.*zeta3 - 1.);
496 
497  case 4:
498  return 4.*zeta0*zeta1;
499 
500  case 5:
501  return 4.*zeta1*zeta2;
502 
503  case 6:
504  return 4.*zeta2*zeta0;
505 
506  case 7:
507  return 4.*zeta0*zeta3;
508 
509  case 8:
510  return 4.*zeta1*zeta3;
511 
512  case 9:
513  return 4.*zeta2*zeta3;
514 
515  default:
516  libmesh_error_msg("Invalid i = " << i);
517  }
518  }
519 
520  // "serendipity" prism
521  case PRISM15:
522  {
523  libmesh_assert_less (i, 15);
524 
525  const Real xi = p(0);
526  const Real eta = p(1);
527  const Real zeta = p(2);
528 
529  switch(i)
530  {
531  case 0:
532  return (1. - zeta)*(xi + eta - 1.)*(xi + eta + 0.5*zeta);
533 
534  case 1:
535  return (1. - zeta)*xi*(xi - 1. - 0.5*zeta);
536 
537  case 2: // phi1 with xi <- eta
538  return (1. - zeta)*eta*(eta - 1. - 0.5*zeta);
539 
540  case 3: // phi0 with zeta <- (-zeta)
541  return (1. + zeta)*(xi + eta - 1.)*(xi + eta - 0.5*zeta);
542 
543  case 4: // phi1 with zeta <- (-zeta)
544  return (1. + zeta)*xi*(xi - 1. + 0.5*zeta);
545 
546  case 5: // phi4 with xi <- eta
547  return (1. + zeta)*eta*(eta - 1. + 0.5*zeta);
548 
549  case 6:
550  return 2.*(1. - zeta)*xi*(1. - xi - eta);
551 
552  case 7:
553  return 2.*(1. - zeta)*xi*eta;
554 
555  case 8:
556  return 2.*(1. - zeta)*eta*(1. - xi - eta);
557 
558  case 9:
559  return (1. - zeta)*(1. + zeta)*(1. - xi - eta);
560 
561  case 10:
562  return (1. - zeta)*(1. + zeta)*xi;
563 
564  case 11: // phi10 with xi <-> eta
565  return (1. - zeta)*(1. + zeta)*eta;
566 
567  case 12: // phi6 with zeta <- (-zeta)
568  return 2.*(1. + zeta)*xi*(1. - xi - eta);
569 
570  case 13: // phi7 with zeta <- (-zeta)
571  return 2.*(1. + zeta)*xi*eta;
572 
573  case 14: // phi8 with zeta <- (-zeta)
574  return 2.*(1. + zeta)*eta*(1. - xi - eta);
575 
576  default:
577  libmesh_error_msg("Invalid i = " << i);
578  }
579  }
580 
581  // quadratic prism shape functions
582  case PRISM18:
583  {
584  libmesh_assert_less (i, 18);
585 
586  // Compute prism shape functions as a tensor-product
587  // of a triangle and an edge
588 
589  Point p2d(p(0),p(1));
590  Real p1d = p(2);
591 
592  // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
593  static const unsigned int i0[] = {0, 0, 0, 1, 1, 1, 0, 0, 0, 2, 2, 2, 1, 1, 1, 2, 2, 2};
594  static const unsigned int i1[] = {0, 1, 2, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 3, 4, 5};
595 
596  return (FE<2,LAGRANGE>::shape(TRI6, SECOND, i1[i], p2d)*
597  fe_lagrange_1D_quadratic_shape(i0[i], p1d));
598  }
599 
600  // G. Bedrosian, "Shape functions and integration formulas for
601  // three-dimensional finite element analysis", Int. J. Numerical
602  // Methods Engineering, vol 35, p. 95-108, 1992.
603  case PYRAMID13:
604  {
605  libmesh_assert_less (i, 13);
606 
607  const Real xi = p(0);
608  const Real eta = p(1);
609  const Real zeta = p(2);
610  const Real eps = 1.e-35;
611 
612  // Denominators are perturbed by epsilon to avoid
613  // divide-by-zero issues.
614  Real den = (1. - zeta + eps);
615 
616  switch(i)
617  {
618  case 0:
619  return 0.25*(-xi - eta - 1.)*((1. - xi)*(1. - eta) - zeta + xi*eta*zeta/den);
620 
621  case 1:
622  return 0.25*(-eta + xi - 1.)*((1. + xi)*(1. - eta) - zeta - xi*eta*zeta/den);
623 
624  case 2:
625  return 0.25*(xi + eta - 1.)*((1. + xi)*(1. + eta) - zeta + xi*eta*zeta/den);
626 
627  case 3:
628  return 0.25*(eta - xi - 1.)*((1. - xi)*(1. + eta) - zeta - xi*eta*zeta/den);
629 
630  case 4:
631  return zeta*(2.*zeta - 1.);
632 
633  case 5:
634  return 0.5*(1. + xi - zeta)*(1. - xi - zeta)*(1. - eta - zeta)/den;
635 
636  case 6:
637  return 0.5*(1. + eta - zeta)*(1. - eta - zeta)*(1. + xi - zeta)/den;
638 
639  case 7:
640  return 0.5*(1. + xi - zeta)*(1. - xi - zeta)*(1. + eta - zeta)/den;
641 
642  case 8:
643  return 0.5*(1. + eta - zeta)*(1. - eta - zeta)*(1. - xi - zeta)/den;
644 
645  case 9:
646  return zeta*(1. - xi - zeta)*(1. - eta - zeta)/den;
647 
648  case 10:
649  return zeta*(1. + xi - zeta)*(1. - eta - zeta)/den;
650 
651  case 11:
652  return zeta*(1. + eta - zeta)*(1. + xi - zeta)/den;
653 
654  case 12:
655  return zeta*(1. - xi - zeta)*(1. + eta - zeta)/den;
656 
657  default:
658  libmesh_error_msg("Invalid i = " << i);
659  }
660  }
661 
662  // Quadratic shape functions, as defined in R. Graglia, "Higher order
663  // bases on pyramidal elements", IEEE Trans Antennas and Propagation,
664  // vol 47, no 5, May 1999.
665  case PYRAMID14:
666  {
667  libmesh_assert_less (i, 14);
668 
669  const Real xi = p(0);
670  const Real eta = p(1);
671  const Real zeta = p(2);
672  const Real eps = 1.e-35;
673 
674  // The "normalized coordinates" defined by Graglia. These are
675  // the planes which define the faces of the pyramid.
676  Real
677  p1 = 0.5*(1. - eta - zeta), // back
678  p2 = 0.5*(1. + xi - zeta), // left
679  p3 = 0.5*(1. + eta - zeta), // front
680  p4 = 0.5*(1. - xi - zeta); // right
681 
682  // Denominators are perturbed by epsilon to avoid
683  // divide-by-zero issues.
684  Real
685  den = (-1. + zeta + eps),
686  den2 = den*den;
687 
688  switch(i)
689  {
690  case 0:
691  return p4*p1*(xi*eta - zeta + zeta*zeta)/den2;
692 
693  case 1:
694  return -p1*p2*(xi*eta + zeta - zeta*zeta)/den2;
695 
696  case 2:
697  return p2*p3*(xi*eta - zeta + zeta*zeta)/den2;
698 
699  case 3:
700  return -p3*p4*(xi*eta + zeta - zeta*zeta)/den2;
701 
702  case 4:
703  return zeta*(2.*zeta - 1.);
704 
705  case 5:
706  return -4.*p2*p1*p4*eta/den2;
707 
708  case 6:
709  return 4.*p1*p2*p3*xi/den2;
710 
711  case 7:
712  return 4.*p2*p3*p4*eta/den2;
713 
714  case 8:
715  return -4.*p3*p4*p1*xi/den2;
716 
717  case 9:
718  return -4.*p1*p4*zeta/den;
719 
720  case 10:
721  return -4.*p2*p1*zeta/den;
722 
723  case 11:
724  return -4.*p3*p2*zeta/den;
725 
726  case 12:
727  return -4.*p4*p3*zeta/den;
728 
729  case 13:
730  return 16.*p1*p2*p3*p4/den2;
731 
732  default:
733  libmesh_error_msg("Invalid i = " << i);
734  }
735  }
736 
737 
738  default:
739  libmesh_error_msg("ERROR: Unsupported 3D element type!: " << type);
740  }
741  }
742 
743 
744  // unsupported order
745  default:
746  libmesh_error_msg("ERROR: Unsupported 3D FE order!: " << order);
747  }
748 
749 #else // LIBMESH_DIM != 3
750  libmesh_ignore(type, order, i, p);
751  libmesh_not_implemented();
752 #endif
753 }
754 
755 
756 
757 Real fe_lagrange_3D_shape_deriv(const ElemType type,
758  const Order order,
759  const unsigned int i,
760  const unsigned int j,
761  const Point & p)
762 {
763 #if LIBMESH_DIM == 3
764 
765  libmesh_assert_less (j, 3);
766 
767  switch (order)
768  {
769  // linear Lagrange shape functions
770  case FIRST:
771  {
772  switch (type)
773  {
774  // trilinear hexahedral shape functions
775  case HEX8:
776  case HEX20:
777  case HEX27:
778  {
779  libmesh_assert_less (i, 8);
780 
781  // Compute hex shape functions as a tensor-product
782  const Real xi = p(0);
783  const Real eta = p(1);
784  const Real zeta = p(2);
785 
786  static const unsigned int i0[] = {0, 1, 1, 0, 0, 1, 1, 0};
787  static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 1, 1};
788  static const unsigned int i2[] = {0, 0, 0, 0, 1, 1, 1, 1};
789 
790  switch(j)
791  {
792  case 0:
793  return (fe_lagrange_1D_linear_shape_deriv(i0[i], 0, xi)*
794  fe_lagrange_1D_linear_shape (i1[i], eta)*
795  fe_lagrange_1D_linear_shape (i2[i], zeta));
796 
797  case 1:
798  return (fe_lagrange_1D_linear_shape (i0[i], xi)*
799  fe_lagrange_1D_linear_shape_deriv(i1[i], 0, eta)*
800  fe_lagrange_1D_linear_shape (i2[i], zeta));
801 
802  case 2:
803  return (fe_lagrange_1D_linear_shape (i0[i], xi)*
804  fe_lagrange_1D_linear_shape (i1[i], eta)*
805  fe_lagrange_1D_linear_shape_deriv(i2[i], 0, zeta));
806 
807  default:
808  libmesh_error_msg("Invalid j = " << j);
809  }
810  }
811 
812  // linear tetrahedral shape functions
813  case TET4:
814  case TET10:
815  {
816  libmesh_assert_less (i, 4);
817 
818  // Area coordinates, pg. 205, Vol. I, Carey, Oden, Becker FEM
819  const Real dzeta0dxi = -1.;
820  const Real dzeta1dxi = 1.;
821  const Real dzeta2dxi = 0.;
822  const Real dzeta3dxi = 0.;
823 
824  const Real dzeta0deta = -1.;
825  const Real dzeta1deta = 0.;
826  const Real dzeta2deta = 1.;
827  const Real dzeta3deta = 0.;
828 
829  const Real dzeta0dzeta = -1.;
830  const Real dzeta1dzeta = 0.;
831  const Real dzeta2dzeta = 0.;
832  const Real dzeta3dzeta = 1.;
833 
834  switch (j)
835  {
836  // d()/dxi
837  case 0:
838  {
839  switch(i)
840  {
841  case 0:
842  return dzeta0dxi;
843 
844  case 1:
845  return dzeta1dxi;
846 
847  case 2:
848  return dzeta2dxi;
849 
850  case 3:
851  return dzeta3dxi;
852 
853  default:
854  libmesh_error_msg("Invalid i = " << i);
855  }
856  }
857 
858  // d()/deta
859  case 1:
860  {
861  switch(i)
862  {
863  case 0:
864  return dzeta0deta;
865 
866  case 1:
867  return dzeta1deta;
868 
869  case 2:
870  return dzeta2deta;
871 
872  case 3:
873  return dzeta3deta;
874 
875  default:
876  libmesh_error_msg("Invalid i = " << i);
877  }
878  }
879 
880  // d()/dzeta
881  case 2:
882  {
883  switch(i)
884  {
885  case 0:
886  return dzeta0dzeta;
887 
888  case 1:
889  return dzeta1dzeta;
890 
891  case 2:
892  return dzeta2dzeta;
893 
894  case 3:
895  return dzeta3dzeta;
896 
897  default:
898  libmesh_error_msg("Invalid i = " << i);
899  }
900  }
901 
902  default:
903  libmesh_error_msg("Invalid shape function derivative j = " << j);
904  }
905  }
906 
907  // linear prism shape functions
908  case PRISM6:
909  case PRISM15:
910  case PRISM18:
911  {
912  libmesh_assert_less (i, 6);
913 
914  // Compute prism shape functions as a tensor-product
915  // of a triangle and an edge
916 
917  Point p2d(p(0),p(1));
918  Real p1d = p(2);
919 
920  // 0 1 2 3 4 5
921  static const unsigned int i0[] = {0, 0, 0, 1, 1, 1};
922  static const unsigned int i1[] = {0, 1, 2, 0, 1, 2};
923 
924  switch (j)
925  {
926  // d()/dxi
927  case 0:
928  return (FE<2,LAGRANGE>::shape_deriv(TRI3, FIRST, i1[i], 0, p2d)*
929  fe_lagrange_1D_linear_shape(i0[i], p1d));
930 
931  // d()/deta
932  case 1:
933  return (FE<2,LAGRANGE>::shape_deriv(TRI3, FIRST, i1[i], 1, p2d)*
934  fe_lagrange_1D_linear_shape(i0[i], p1d));
935 
936  // d()/dzeta
937  case 2:
938  return (FE<2,LAGRANGE>::shape(TRI3, FIRST, i1[i], p2d)*
939  fe_lagrange_1D_linear_shape_deriv(i0[i], 0, p1d));
940 
941  default:
942  libmesh_error_msg("Invalid shape function derivative j = " << j);
943  }
944  }
945 
946  // linear pyramid shape functions
947  case PYRAMID5:
948  case PYRAMID13:
949  case PYRAMID14:
950  {
951  libmesh_assert_less (i, 5);
952 
953  const Real xi = p(0);
954  const Real eta = p(1);
955  const Real zeta = p(2);
956  const Real eps = 1.e-35;
957 
958  switch (j)
959  {
960  // d/dxi
961  case 0:
962  switch(i)
963  {
964  case 0:
965  return .25*(zeta + eta - 1.)/((1. - zeta) + eps);
966 
967  case 1:
968  return -.25*(zeta + eta - 1.)/((1. - zeta) + eps);
969 
970  case 2:
971  return -.25*(zeta - eta - 1.)/((1. - zeta) + eps);
972 
973  case 3:
974  return .25*(zeta - eta - 1.)/((1. - zeta) + eps);
975 
976  case 4:
977  return 0;
978 
979  default:
980  libmesh_error_msg("Invalid i = " << i);
981  }
982 
983 
984  // d/deta
985  case 1:
986  switch(i)
987  {
988  case 0:
989  return .25*(zeta + xi - 1.)/((1. - zeta) + eps);
990 
991  case 1:
992  return .25*(zeta - xi - 1.)/((1. - zeta) + eps);
993 
994  case 2:
995  return -.25*(zeta - xi - 1.)/((1. - zeta) + eps);
996 
997  case 3:
998  return -.25*(zeta + xi - 1.)/((1. - zeta) + eps);
999 
1000  case 4:
1001  return 0;
1002 
1003  default:
1004  libmesh_error_msg("Invalid i = " << i);
1005  }
1006 
1007 
1008  // d/dzeta
1009  case 2:
1010  {
1011  // We computed the derivatives with general eps and
1012  // then let eps tend to zero in the numerators...
1013  Real
1014  num = zeta*(2. - zeta) - 1.,
1015  den = (1. - zeta + eps)*(1. - zeta + eps);
1016 
1017  switch(i)
1018  {
1019  case 0:
1020  case 2:
1021  return .25*(num + xi*eta)/den;
1022 
1023  case 1:
1024  case 3:
1025  return .25*(num - xi*eta)/den;
1026 
1027  case 4:
1028  return 1.;
1029 
1030  default:
1031  libmesh_error_msg("Invalid i = " << i);
1032  }
1033  }
1034 
1035  default:
1036  libmesh_error_msg("Invalid j = " << j);
1037  }
1038  }
1039 
1040 
1041  default:
1042  libmesh_error_msg("ERROR: Unsupported 3D element type!: " << type);
1043  }
1044  }
1045 
1046 
1047  // quadratic Lagrange shape functions
1048  case SECOND:
1049  {
1050  switch (type)
1051  {
1052 
1053  // serendipity hexahedral quadratic shape functions
1054  case HEX20:
1055  {
1056  libmesh_assert_less (i, 20);
1057 
1058  const Real xi = p(0);
1059  const Real eta = p(1);
1060  const Real zeta = p(2);
1061 
1062  // these functions are defined for (x,y,z) in [0,1]^3
1063  // so transform the locations
1064  const Real x = .5*(xi + 1.);
1065  const Real y = .5*(eta + 1.);
1066  const Real z = .5*(zeta + 1.);
1067 
1068  // and don't forget the chain rule!
1069 
1070  switch (j)
1071  {
1072 
1073  // d/dx*dx/dxi
1074  case 0:
1075  switch (i)
1076  {
1077  case 0:
1078  return .5*(1. - y)*(1. - z)*((1. - x)*(-2.) +
1079  (-1.)*(1. - 2.*x - 2.*y - 2.*z));
1080 
1081  case 1:
1082  return .5*(1. - y)*(1. - z)*(x*(2.) +
1083  (1.)*(2.*x - 2.*y - 2.*z - 1.));
1084 
1085  case 2:
1086  return .5*y*(1. - z)*(x*(2.) +
1087  (1.)*(2.*x + 2.*y - 2.*z - 3.));
1088 
1089  case 3:
1090  return .5*y*(1. - z)*((1. - x)*(-2.) +
1091  (-1.)*(2.*y - 2.*x - 2.*z - 1.));
1092 
1093  case 4:
1094  return .5*(1. - y)*z*((1. - x)*(-2.) +
1095  (-1.)*(2.*z - 2.*x - 2.*y - 1.));
1096 
1097  case 5:
1098  return .5*(1. - y)*z*(x*(2.) +
1099  (1.)*(2.*x - 2.*y + 2.*z - 3.));
1100 
1101  case 6:
1102  return .5*y*z*(x*(2.) +
1103  (1.)*(2.*x + 2.*y + 2.*z - 5.));
1104 
1105  case 7:
1106  return .5*y*z*((1. - x)*(-2.) +
1107  (-1.)*(2.*y - 2.*x + 2.*z - 3.));
1108 
1109  case 8:
1110  return 2.*(1. - y)*(1. - z)*(1. - 2.*x);
1111 
1112  case 9:
1113  return 2.*y*(1. - y)*(1. - z);
1114 
1115  case 10:
1116  return 2.*y*(1. - z)*(1. - 2.*x);
1117 
1118  case 11:
1119  return 2.*y*(1. - y)*(1. - z)*(-1.);
1120 
1121  case 12:
1122  return 2.*(1. - y)*z*(1. - z)*(-1.);
1123 
1124  case 13:
1125  return 2.*(1. - y)*z*(1. - z);
1126 
1127  case 14:
1128  return 2.*y*z*(1. - z);
1129 
1130  case 15:
1131  return 2.*y*z*(1. - z)*(-1.);
1132 
1133  case 16:
1134  return 2.*(1. - y)*z*(1. - 2.*x);
1135 
1136  case 17:
1137  return 2.*y*(1. - y)*z;
1138 
1139  case 18:
1140  return 2.*y*z*(1. - 2.*x);
1141 
1142  case 19:
1143  return 2.*y*(1. - y)*z*(-1.);
1144 
1145  default:
1146  libmesh_error_msg("Invalid i = " << i);
1147  }
1148 
1149 
1150  // d/dy*dy/deta
1151  case 1:
1152  switch (i)
1153  {
1154  case 0:
1155  return .5*(1. - x)*(1. - z)*((1. - y)*(-2.) +
1156  (-1.)*(1. - 2.*x - 2.*y - 2.*z));
1157 
1158  case 1:
1159  return .5*x*(1. - z)*((1. - y)*(-2.) +
1160  (-1.)*(2.*x - 2.*y - 2.*z - 1.));
1161 
1162  case 2:
1163  return .5*x*(1. - z)*(y*(2.) +
1164  (1.)*(2.*x + 2.*y - 2.*z - 3.));
1165 
1166  case 3:
1167  return .5*(1. - x)*(1. - z)*(y*(2.) +
1168  (1.)*(2.*y - 2.*x - 2.*z - 1.));
1169 
1170  case 4:
1171  return .5*(1. - x)*z*((1. - y)*(-2.) +
1172  (-1.)*(2.*z - 2.*x - 2.*y - 1.));
1173 
1174  case 5:
1175  return .5*x*z*((1. - y)*(-2.) +
1176  (-1.)*(2.*x - 2.*y + 2.*z - 3.));
1177 
1178  case 6:
1179  return .5*x*z*(y*(2.) +
1180  (1.)*(2.*x + 2.*y + 2.*z - 5.));
1181 
1182  case 7:
1183  return .5*(1. - x)*z*(y*(2.) +
1184  (1.)*(2.*y - 2.*x + 2.*z - 3.));
1185 
1186  case 8:
1187  return 2.*x*(1. - x)*(1. - z)*(-1.);
1188 
1189  case 9:
1190  return 2.*x*(1. - z)*(1. - 2.*y);
1191 
1192  case 10:
1193  return 2.*x*(1. - x)*(1. - z);
1194 
1195  case 11:
1196  return 2.*(1. - x)*(1. - z)*(1. - 2.*y);
1197 
1198  case 12:
1199  return 2.*(1. - x)*z*(1. - z)*(-1.);
1200 
1201  case 13:
1202  return 2.*x*z*(1. - z)*(-1.);
1203 
1204  case 14:
1205  return 2.*x*z*(1. - z);
1206 
1207  case 15:
1208  return 2.*(1. - x)*z*(1. - z);
1209 
1210  case 16:
1211  return 2.*x*(1. - x)*z*(-1.);
1212 
1213  case 17:
1214  return 2.*x*z*(1. - 2.*y);
1215 
1216  case 18:
1217  return 2.*x*(1. - x)*z;
1218 
1219  case 19:
1220  return 2.*(1. - x)*z*(1. - 2.*y);
1221 
1222  default:
1223  libmesh_error_msg("Invalid i = " << i);
1224  }
1225 
1226 
1227  // d/dz*dz/dzeta
1228  case 2:
1229  switch (i)
1230  {
1231  case 0:
1232  return .5*(1. - x)*(1. - y)*((1. - z)*(-2.) +
1233  (-1.)*(1. - 2.*x - 2.*y - 2.*z));
1234 
1235  case 1:
1236  return .5*x*(1. - y)*((1. - z)*(-2.) +
1237  (-1.)*(2.*x - 2.*y - 2.*z - 1.));
1238 
1239  case 2:
1240  return .5*x*y*((1. - z)*(-2.) +
1241  (-1.)*(2.*x + 2.*y - 2.*z - 3.));
1242 
1243  case 3:
1244  return .5*(1. - x)*y*((1. - z)*(-2.) +
1245  (-1.)*(2.*y - 2.*x - 2.*z - 1.));
1246 
1247  case 4:
1248  return .5*(1. - x)*(1. - y)*(z*(2.) +
1249  (1.)*(2.*z - 2.*x - 2.*y - 1.));
1250 
1251  case 5:
1252  return .5*x*(1. - y)*(z*(2.) +
1253  (1.)*(2.*x - 2.*y + 2.*z - 3.));
1254 
1255  case 6:
1256  return .5*x*y*(z*(2.) +
1257  (1.)*(2.*x + 2.*y + 2.*z - 5.));
1258 
1259  case 7:
1260  return .5*(1. - x)*y*(z*(2.) +
1261  (1.)*(2.*y - 2.*x + 2.*z - 3.));
1262 
1263  case 8:
1264  return 2.*x*(1. - x)*(1. - y)*(-1.);
1265 
1266  case 9:
1267  return 2.*x*y*(1. - y)*(-1.);
1268 
1269  case 10:
1270  return 2.*x*(1. - x)*y*(-1.);
1271 
1272  case 11:
1273  return 2.*(1. - x)*y*(1. - y)*(-1.);
1274 
1275  case 12:
1276  return 2.*(1. - x)*(1. - y)*(1. - 2.*z);
1277 
1278  case 13:
1279  return 2.*x*(1. - y)*(1. - 2.*z);
1280 
1281  case 14:
1282  return 2.*x*y*(1. - 2.*z);
1283 
1284  case 15:
1285  return 2.*(1. - x)*y*(1. - 2.*z);
1286 
1287  case 16:
1288  return 2.*x*(1. - x)*(1. - y);
1289 
1290  case 17:
1291  return 2.*x*y*(1. - y);
1292 
1293  case 18:
1294  return 2.*x*(1. - x)*y;
1295 
1296  case 19:
1297  return 2.*(1. - x)*y*(1. - y);
1298 
1299  default:
1300  libmesh_error_msg("Invalid i = " << i);
1301  }
1302 
1303  default:
1304  libmesh_error_msg("Invalid shape function derivative j = " << j);
1305  }
1306  }
1307 
1308  // triquadratic hexahedral shape functions
1309  case HEX27:
1310  {
1311  libmesh_assert_less (i, 27);
1312 
1313  // Compute hex shape functions as a tensor-product
1314  const Real xi = p(0);
1315  const Real eta = p(1);
1316  const Real zeta = p(2);
1317 
1318  // The only way to make any sense of this
1319  // is to look at the mgflo/mg2/mgf documentation
1320  // and make the cut-out cube!
1321  // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26
1322  static const unsigned int i0[] = {0, 1, 1, 0, 0, 1, 1, 0, 2, 1, 2, 0, 0, 1, 1, 0, 2, 1, 2, 0, 2, 2, 1, 2, 0, 2, 2};
1323  static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 1, 1, 0, 2, 1, 2, 0, 0, 1, 1, 0, 2, 1, 2, 2, 0, 2, 1, 2, 2, 2};
1324  static const unsigned int i2[] = {0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 2, 2, 2, 2, 1, 1, 1, 1, 0, 2, 2, 2, 2, 1, 2};
1325 
1326  switch(j)
1327  {
1328  case 0:
1329  return (fe_lagrange_1D_quadratic_shape_deriv(i0[i], 0, xi)*
1330  fe_lagrange_1D_quadratic_shape (i1[i], eta)*
1331  fe_lagrange_1D_quadratic_shape (i2[i], zeta));
1332 
1333  case 1:
1334  return (fe_lagrange_1D_quadratic_shape (i0[i], xi)*
1335  fe_lagrange_1D_quadratic_shape_deriv(i1[i], 0, eta)*
1336  fe_lagrange_1D_quadratic_shape (i2[i], zeta));
1337 
1338  case 2:
1339  return (fe_lagrange_1D_quadratic_shape (i0[i], xi)*
1340  fe_lagrange_1D_quadratic_shape (i1[i], eta)*
1341  fe_lagrange_1D_quadratic_shape_deriv(i2[i], 0, zeta));
1342 
1343  default:
1344  libmesh_error_msg("Invalid j = " << j);
1345  }
1346  }
1347 
1348  // quadratic tetrahedral shape functions
1349  case TET10:
1350  {
1351  libmesh_assert_less (i, 10);
1352 
1353  // Area coordinates, pg. 205, Vol. I, Carey, Oden, Becker FEM
1354  const Real zeta1 = p(0);
1355  const Real zeta2 = p(1);
1356  const Real zeta3 = p(2);
1357  const Real zeta0 = 1. - zeta1 - zeta2 - zeta3;
1358 
1359  const Real dzeta0dxi = -1.;
1360  const Real dzeta1dxi = 1.;
1361  const Real dzeta2dxi = 0.;
1362  const Real dzeta3dxi = 0.;
1363 
1364  const Real dzeta0deta = -1.;
1365  const Real dzeta1deta = 0.;
1366  const Real dzeta2deta = 1.;
1367  const Real dzeta3deta = 0.;
1368 
1369  const Real dzeta0dzeta = -1.;
1370  const Real dzeta1dzeta = 0.;
1371  const Real dzeta2dzeta = 0.;
1372  const Real dzeta3dzeta = 1.;
1373 
1374  switch (j)
1375  {
1376  // d()/dxi
1377  case 0:
1378  {
1379  switch(i)
1380  {
1381  case 0:
1382  return (4.*zeta0 - 1.)*dzeta0dxi;
1383 
1384  case 1:
1385  return (4.*zeta1 - 1.)*dzeta1dxi;
1386 
1387  case 2:
1388  return (4.*zeta2 - 1.)*dzeta2dxi;
1389 
1390  case 3:
1391  return (4.*zeta3 - 1.)*dzeta3dxi;
1392 
1393  case 4:
1394  return 4.*(zeta0*dzeta1dxi + dzeta0dxi*zeta1);
1395 
1396  case 5:
1397  return 4.*(zeta1*dzeta2dxi + dzeta1dxi*zeta2);
1398 
1399  case 6:
1400  return 4.*(zeta0*dzeta2dxi + dzeta0dxi*zeta2);
1401 
1402  case 7:
1403  return 4.*(zeta0*dzeta3dxi + dzeta0dxi*zeta3);
1404 
1405  case 8:
1406  return 4.*(zeta1*dzeta3dxi + dzeta1dxi*zeta3);
1407 
1408  case 9:
1409  return 4.*(zeta2*dzeta3dxi + dzeta2dxi*zeta3);
1410 
1411  default:
1412  libmesh_error_msg("Invalid i = " << i);
1413  }
1414  }
1415 
1416  // d()/deta
1417  case 1:
1418  {
1419  switch(i)
1420  {
1421  case 0:
1422  return (4.*zeta0 - 1.)*dzeta0deta;
1423 
1424  case 1:
1425  return (4.*zeta1 - 1.)*dzeta1deta;
1426 
1427  case 2:
1428  return (4.*zeta2 - 1.)*dzeta2deta;
1429 
1430  case 3:
1431  return (4.*zeta3 - 1.)*dzeta3deta;
1432 
1433  case 4:
1434  return 4.*(zeta0*dzeta1deta + dzeta0deta*zeta1);
1435 
1436  case 5:
1437  return 4.*(zeta1*dzeta2deta + dzeta1deta*zeta2);
1438 
1439  case 6:
1440  return 4.*(zeta0*dzeta2deta + dzeta0deta*zeta2);
1441 
1442  case 7:
1443  return 4.*(zeta0*dzeta3deta + dzeta0deta*zeta3);
1444 
1445  case 8:
1446  return 4.*(zeta1*dzeta3deta + dzeta1deta*zeta3);
1447 
1448  case 9:
1449  return 4.*(zeta2*dzeta3deta + dzeta2deta*zeta3);
1450 
1451  default:
1452  libmesh_error_msg("Invalid i = " << i);
1453  }
1454  }
1455 
1456  // d()/dzeta
1457  case 2:
1458  {
1459  switch(i)
1460  {
1461  case 0:
1462  return (4.*zeta0 - 1.)*dzeta0dzeta;
1463 
1464  case 1:
1465  return (4.*zeta1 - 1.)*dzeta1dzeta;
1466 
1467  case 2:
1468  return (4.*zeta2 - 1.)*dzeta2dzeta;
1469 
1470  case 3:
1471  return (4.*zeta3 - 1.)*dzeta3dzeta;
1472 
1473  case 4:
1474  return 4.*(zeta0*dzeta1dzeta + dzeta0dzeta*zeta1);
1475 
1476  case 5:
1477  return 4.*(zeta1*dzeta2dzeta + dzeta1dzeta*zeta2);
1478 
1479  case 6:
1480  return 4.*(zeta0*dzeta2dzeta + dzeta0dzeta*zeta2);
1481 
1482  case 7:
1483  return 4.*(zeta0*dzeta3dzeta + dzeta0dzeta*zeta3);
1484 
1485  case 8:
1486  return 4.*(zeta1*dzeta3dzeta + dzeta1dzeta*zeta3);
1487 
1488  case 9:
1489  return 4.*(zeta2*dzeta3dzeta + dzeta2dzeta*zeta3);
1490 
1491  default:
1492  libmesh_error_msg("Invalid i = " << i);
1493  }
1494  }
1495 
1496  default:
1497  libmesh_error_msg("Invalid j = " << j);
1498  }
1499  }
1500 
1501 
1502  // "serendipity" prism
1503  case PRISM15:
1504  {
1505  libmesh_assert_less (i, 15);
1506 
1507  const Real xi = p(0);
1508  const Real eta = p(1);
1509  const Real zeta = p(2);
1510 
1511  switch (j)
1512  {
1513  // d()/dxi
1514  case 0:
1515  {
1516  switch(i)
1517  {
1518  case 0:
1519  return (2.*xi + 2.*eta + 0.5*zeta - 1.)*(1. - zeta);
1520  case 1:
1521  return (2.*xi - 1. - 0.5*zeta)*(1. - zeta);
1522  case 2:
1523  return 0.;
1524  case 3:
1525  return (2.*xi + 2.*eta - 0.5*zeta - 1.)*(1. + zeta);
1526  case 4:
1527  return (2.*xi - 1. + 0.5*zeta)*(1. + zeta);
1528  case 5:
1529  return 0.;
1530  case 6:
1531  return (4.*xi + 2.*eta - 2.)*(zeta - 1.);
1532  case 7:
1533  return -2.*(zeta - 1.)*eta;
1534  case 8:
1535  return 2.*(zeta - 1.)*eta;
1536  case 9:
1537  return (zeta - 1.)*(1. + zeta);
1538  case 10:
1539  return (1. - zeta)*(1. + zeta);
1540  case 11:
1541  return 0.;
1542  case 12:
1543  return (-4.*xi - 2.*eta + 2.)*(1. + zeta);
1544  case 13:
1545  return 2.*(1. + zeta)*eta;
1546  case 14:
1547  return -2.*(1. + zeta)*eta;
1548  default:
1549  libmesh_error_msg("Invalid i = " << i);
1550  }
1551  }
1552 
1553  // d()/deta
1554  case 1:
1555  {
1556  switch(i)
1557  {
1558  case 0:
1559  return (2.*xi + 2.*eta + 0.5*zeta - 1.)*(1. - zeta);
1560  case 1:
1561  return 0.;
1562  case 2:
1563  return (2.*eta - 1. - 0.5*zeta)*(1. - zeta);
1564  case 3:
1565  return (2.*xi + 2.*eta - 0.5*zeta - 1.)*(1. + zeta);
1566  case 4:
1567  return 0.;
1568  case 5:
1569  return (2.*eta - 1. + 0.5*zeta)*(1. + zeta);
1570  case 6:
1571  return 2.*(zeta - 1.)*xi;
1572  case 7:
1573  return 2.*(1. - zeta)*xi;
1574  case 8:
1575  return (2.*xi + 4.*eta - 2.)*(zeta - 1.);
1576  case 9:
1577  return (zeta - 1.)*(1. + zeta);
1578  case 10:
1579  return 0.;
1580  case 11:
1581  return (1. - zeta)*(1. + zeta);
1582  case 12:
1583  return -2.*(1. + zeta)*xi;
1584  case 13:
1585  return 2.*(1. + zeta)*xi;
1586  case 14:
1587  return (-2.*xi - 4.*eta + 2.)*(1. + zeta);
1588 
1589  default:
1590  libmesh_error_msg("Invalid i = " << i);
1591  }
1592  }
1593 
1594  // d()/dzeta
1595  case 2:
1596  {
1597  switch(i)
1598  {
1599  case 0:
1600  return (-xi - eta - zeta + 0.5)*(xi + eta - 1.);
1601  case 1:
1602  return -0.5*xi*(2.*xi - 1. - 2.*zeta);
1603  case 2:
1604  return -0.5*eta*(2.*eta - 1. - 2.*zeta);
1605  case 3:
1606  return (xi + eta - zeta - 0.5)*(xi + eta - 1.);
1607  case 4:
1608  return 0.5*xi*(2.*xi - 1. + 2.*zeta);
1609  case 5:
1610  return 0.5*eta*(2.*eta - 1. + 2.*zeta);
1611  case 6:
1612  return 2.*xi*(xi + eta - 1.);
1613  case 7:
1614  return -2.*xi*eta;
1615  case 8:
1616  return 2.*eta*(xi + eta - 1.);
1617  case 9:
1618  return 2.*zeta*(xi + eta - 1.);
1619  case 10:
1620  return -2.*xi*zeta;
1621  case 11:
1622  return -2.*eta*zeta;
1623  case 12:
1624  return 2.*xi*(1. - xi - eta);
1625  case 13:
1626  return 2.*xi*eta;
1627  case 14:
1628  return 2.*eta*(1. - xi - eta);
1629 
1630  default:
1631  libmesh_error_msg("Invalid i = " << i);
1632  }
1633  }
1634 
1635  default:
1636  libmesh_error_msg("Invalid j = " << j);
1637  }
1638  }
1639 
1640 
1641 
1642  // quadratic prism shape functions
1643  case PRISM18:
1644  {
1645  libmesh_assert_less (i, 18);
1646 
1647  // Compute prism shape functions as a tensor-product
1648  // of a triangle and an edge
1649 
1650  Point p2d(p(0),p(1));
1651  Real p1d = p(2);
1652 
1653  // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
1654  static const unsigned int i0[] = {0, 0, 0, 1, 1, 1, 0, 0, 0, 2, 2, 2, 1, 1, 1, 2, 2, 2};
1655  static const unsigned int i1[] = {0, 1, 2, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 3, 4, 5};
1656 
1657  switch (j)
1658  {
1659  // d()/dxi
1660  case 0:
1661  return (FE<2,LAGRANGE>::shape_deriv(TRI6, SECOND, i1[i], 0, p2d)*
1662  fe_lagrange_1D_quadratic_shape(i0[i], p1d));
1663 
1664  // d()/deta
1665  case 1:
1666  return (FE<2,LAGRANGE>::shape_deriv(TRI6, SECOND, i1[i], 1, p2d)*
1667  fe_lagrange_1D_quadratic_shape(i0[i], p1d));
1668 
1669  // d()/dzeta
1670  case 2:
1671  return (FE<2,LAGRANGE>::shape(TRI6, SECOND, i1[i], p2d)*
1672  fe_lagrange_1D_quadratic_shape_deriv(i0[i], 0, p1d));
1673 
1674  default:
1675  libmesh_error_msg("Invalid shape function derivative j = " << j);
1676  }
1677  }
1678 
1679  // G. Bedrosian, "Shape functions and integration formulas for
1680  // three-dimensional finite element analysis", Int. J. Numerical
1681  // Methods Engineering, vol 35, p. 95-108, 1992.
1682  case PYRAMID13:
1683  {
1684  libmesh_assert_less (i, 13);
1685 
1686  const Real xi = p(0);
1687  const Real eta = p(1);
1688  const Real zeta = p(2);
1689  const Real eps = 1.e-35;
1690 
1691  // Denominators are perturbed by epsilon to avoid
1692  // divide-by-zero issues.
1693  Real
1694  den = (-1. + zeta + eps),
1695  den2 = den*den,
1696  xi2 = xi*xi,
1697  eta2 = eta*eta,
1698  zeta2 = zeta*zeta,
1699  zeta3 = zeta2*zeta;
1700 
1701  switch (j)
1702  {
1703  // d/dxi
1704  case 0:
1705  switch(i)
1706  {
1707  case 0:
1708  return 0.25*(-zeta - eta + 2.*eta*zeta - 2.*xi + 2.*zeta*xi + 2.*eta*xi + zeta2 + eta2)/den;
1709 
1710  case 1:
1711  return -0.25*(-zeta - eta + 2.*eta*zeta + 2.*xi - 2.*zeta*xi - 2.*eta*xi + zeta2 + eta2)/den;
1712 
1713  case 2:
1714  return -0.25*(-zeta + eta - 2.*eta*zeta + 2.*xi - 2.*zeta*xi + 2.*eta*xi + zeta2 + eta2)/den;
1715 
1716  case 3:
1717  return 0.25*(-zeta + eta - 2.*eta*zeta - 2.*xi + 2.*zeta*xi - 2.*eta*xi + zeta2 + eta2)/den;
1718 
1719  case 4:
1720  return 0.;
1721 
1722  case 5:
1723  return -(-1. + eta + zeta)*xi/den;
1724 
1725  case 6:
1726  return 0.5*(-1. + eta + zeta)*(1. + eta - zeta)/den;
1727 
1728  case 7:
1729  return (1. + eta - zeta)*xi/den;
1730 
1731  case 8:
1732  return -0.5*(-1. + eta + zeta)*(1. + eta - zeta)/den;
1733 
1734  case 9:
1735  return -(-1. + eta + zeta)*zeta/den;
1736 
1737  case 10:
1738  return (-1. + eta + zeta)*zeta/den;
1739 
1740  case 11:
1741  return -(1. + eta - zeta)*zeta/den;
1742 
1743  case 12:
1744  return (1. + eta - zeta)*zeta/den;
1745 
1746  default:
1747  libmesh_error_msg("Invalid i = " << i);
1748  }
1749 
1750  // d/deta
1751  case 1:
1752  switch(i)
1753  {
1754  case 0:
1755  return 0.25*(-zeta - 2.*eta + 2.*eta*zeta - xi + 2.*zeta*xi + 2.*eta*xi + zeta2 + xi2)/den;
1756 
1757  case 1:
1758  return -0.25*(zeta + 2.*eta - 2.*eta*zeta - xi + 2.*zeta*xi + 2.*eta*xi - zeta2 - xi2)/den;
1759 
1760  case 2:
1761  return -0.25*(-zeta + 2.*eta - 2.*eta*zeta + xi - 2.*zeta*xi + 2.*eta*xi + zeta2 + xi2)/den;
1762 
1763  case 3:
1764  return 0.25*(zeta - 2.*eta + 2.*eta*zeta + xi - 2.*zeta*xi + 2.*eta*xi - zeta2 - xi2)/den;
1765 
1766  case 4:
1767  return 0.;
1768 
1769  case 5:
1770  return -0.5*(-1. + xi + zeta)*(1. + xi - zeta)/den;
1771 
1772  case 6:
1773  return (1. + xi - zeta)*eta/den;
1774 
1775  case 7:
1776  return 0.5*(-1. + xi + zeta)*(1. + xi - zeta)/den;
1777 
1778  case 8:
1779  return -(-1. + xi + zeta)*eta/den;
1780 
1781  case 9:
1782  return -(-1. + xi + zeta)*zeta/den;
1783 
1784  case 10:
1785  return (1. + xi - zeta)*zeta/den;
1786 
1787  case 11:
1788  return -(1. + xi - zeta)*zeta/den;
1789 
1790  case 12:
1791  return (-1. + xi + zeta)*zeta/den;
1792 
1793  default:
1794  libmesh_error_msg("Invalid i = " << i);
1795  }
1796 
1797  // d/dzeta
1798  case 2:
1799  {
1800  switch(i)
1801  {
1802  case 0:
1803  return -0.25*(xi + eta + 1.)*(-1. + 2.*zeta - zeta2 + eta*xi)/den2;
1804 
1805  case 1:
1806  return 0.25*(eta - xi + 1.)*(1. - 2.*zeta + zeta2 + eta*xi)/den2;
1807 
1808  case 2:
1809  return 0.25*(xi + eta - 1.)*(-1. + 2.*zeta - zeta2 + eta*xi)/den2;
1810 
1811  case 3:
1812  return -0.25*(eta - xi - 1.)*(1. - 2.*zeta + zeta2 + eta*xi)/den2;
1813 
1814  case 4:
1815  return 4.*zeta - 1.;
1816 
1817  case 5:
1818  return 0.5*(-2 + eta + 6.*zeta + eta*xi2 + eta*zeta2 - 6.*zeta2 + 2.*zeta3 - 2.*eta*zeta)/den2;
1819 
1820  case 6:
1821  return -0.5*(2 - 6.*zeta + xi + xi*zeta2 + eta2*xi + 6.*zeta2 - 2.*zeta3 - 2.*zeta*xi)/den2;
1822 
1823  case 7:
1824  return -0.5*(2 + eta - 6.*zeta + eta*xi2 + eta*zeta2 + 6.*zeta2 - 2.*zeta3 - 2.*eta*zeta)/den2;
1825 
1826  case 8:
1827  return 0.5*(-2 + 6.*zeta + xi + xi*zeta2 + eta2*xi - 6.*zeta2 + 2.*zeta3 - 2.*zeta*xi)/den2;
1828 
1829  case 9:
1830  return (1. - eta - 4.*zeta - xi - xi*zeta2 - eta*zeta2 + eta*xi + 5.*zeta2 - 2.*zeta3 + 2.*eta*zeta + 2.*zeta*xi)/den2;
1831 
1832  case 10:
1833  return -(-1. + eta + 4.*zeta - xi - xi*zeta2 + eta*zeta2 + eta*xi - 5.*zeta2 + 2.*zeta3 - 2.*eta*zeta + 2.*zeta*xi)/den2;
1834 
1835  case 11:
1836  return (1. + eta - 4.*zeta + xi + xi*zeta2 + eta*zeta2 + eta*xi + 5.*zeta2 - 2.*zeta3 - 2.*eta*zeta - 2.*zeta*xi)/den2;
1837 
1838  case 12:
1839  return -(-1. - eta + 4.*zeta + xi + xi*zeta2 - eta*zeta2 + eta*xi - 5.*zeta2 + 2.*zeta3 + 2.*eta*zeta - 2.*zeta*xi)/den2;
1840 
1841  default:
1842  libmesh_error_msg("Invalid i = " << i);
1843  }
1844  }
1845 
1846  default:
1847  libmesh_error_msg("Invalid j = " << j);
1848  }
1849  }
1850 
1851  // Quadratic shape functions, as defined in R. Graglia, "Higher order
1852  // bases on pyramidal elements", IEEE Trans Antennas and Propagation,
1853  // vol 47, no 5, May 1999.
1854  case PYRAMID14:
1855  {
1856  libmesh_assert_less (i, 14);
1857 
1858  const Real xi = p(0);
1859  const Real eta = p(1);
1860  const Real zeta = p(2);
1861  const Real eps = 1.e-35;
1862 
1863  // The "normalized coordinates" defined by Graglia. These are
1864  // the planes which define the faces of the pyramid.
1865  Real
1866  p1 = 0.5*(1. - eta - zeta), // back
1867  p2 = 0.5*(1. + xi - zeta), // left
1868  p3 = 0.5*(1. + eta - zeta), // front
1869  p4 = 0.5*(1. - xi - zeta); // right
1870 
1871  // Denominators are perturbed by epsilon to avoid
1872  // divide-by-zero issues.
1873  Real
1874  den = (-1. + zeta + eps),
1875  den2 = den*den,
1876  den3 = den2*den;
1877 
1878  switch (j)
1879  {
1880  // d/dxi
1881  case 0:
1882  switch(i)
1883  {
1884  case 0:
1885  return 0.5*p1*(-xi*eta + zeta - zeta*zeta + 2.*p4*eta)/den2;
1886 
1887  case 1:
1888  return -0.5*p1*(xi*eta + zeta - zeta*zeta + 2.*p2*eta)/den2;
1889 
1890  case 2:
1891  return 0.5*p3*(xi*eta - zeta + zeta*zeta + 2.*p2*eta)/den2;
1892 
1893  case 3:
1894  return -0.5*p3*(-xi*eta - zeta + zeta*zeta + 2.*p4*eta)/den2;
1895 
1896  case 4:
1897  return 0.;
1898 
1899  case 5:
1900  return 2.*p1*eta*xi/den2;
1901 
1902  case 6:
1903  return 2.*p1*p3*(xi + 2.*p2)/den2;
1904 
1905  case 7:
1906  return -2.*p3*eta*xi/den2;
1907 
1908  case 8:
1909  return -2.*p1*p3*(-xi + 2.*p4)/den2;
1910 
1911  case 9:
1912  return 2.*p1*zeta/den;
1913 
1914  case 10:
1915  return -2.*p1*zeta/den;
1916 
1917  case 11:
1918  return -2.*p3*zeta/den;
1919 
1920  case 12:
1921  return 2.*p3*zeta/den;
1922 
1923  case 13:
1924  return -8.*p1*p3*xi/den2;
1925 
1926  default:
1927  libmesh_error_msg("Invalid i = " << i);
1928  }
1929 
1930  // d/deta
1931  case 1:
1932  switch(i)
1933  {
1934  case 0:
1935  return -0.5*p4*(xi*eta - zeta + zeta*zeta - 2.*p1*xi)/den2;
1936 
1937  case 1:
1938  return 0.5*p2*(xi*eta + zeta - zeta*zeta - 2.*p1*xi)/den2;
1939 
1940  case 2:
1941  return 0.5*p2*(xi*eta - zeta + zeta*zeta + 2.*p3*xi)/den2;
1942 
1943  case 3:
1944  return -0.5*p4*(xi*eta + zeta - zeta*zeta + 2.*p3*xi)/den2;
1945 
1946  case 4:
1947  return 0.;
1948 
1949  case 5:
1950  return 2.*p2*p4*(eta - 2.*p1)/den2;
1951 
1952  case 6:
1953  return -2.*p2*xi*eta/den2;
1954 
1955  case 7:
1956  return 2.*p2*p4*(eta + 2.*p3)/den2;
1957 
1958  case 8:
1959  return 2.*p4*xi*eta/den2;
1960 
1961  case 9:
1962  return 2.*p4*zeta/den;
1963 
1964  case 10:
1965  return 2.*p2*zeta/den;
1966 
1967  case 11:
1968  return -2.*p2*zeta/den;
1969 
1970  case 12:
1971  return -2.*p4*zeta/den;
1972 
1973  case 13:
1974  return -8.*p2*p4*eta/den2;
1975 
1976  default:
1977  libmesh_error_msg("Invalid i = " << i);
1978  }
1979 
1980 
1981  // d/dzeta
1982  case 2:
1983  {
1984  switch(i)
1985  {
1986  case 0:
1987  return -0.5*p1*(xi*eta - zeta + zeta*zeta)/den2
1988  - 0.5*p4*(xi*eta - zeta + zeta*zeta)/den2
1989  + p4*p1*(2.*zeta - 1)/den2
1990  - 2.*p4*p1*(xi*eta - zeta + zeta*zeta)/den3;
1991 
1992  case 1:
1993  return 0.5*p2*(xi*eta + zeta - zeta*zeta)/den2
1994  + 0.5*p1*(xi*eta + zeta - zeta*zeta)/den2
1995  - p1*p2*(1 - 2.*zeta)/den2
1996  + 2.*p1*p2*(xi*eta + zeta - zeta*zeta)/den3;
1997 
1998  case 2:
1999  return -0.5*p3*(xi*eta - zeta + zeta*zeta)/den2
2000  - 0.5*p2*(xi*eta - zeta + zeta*zeta)/den2
2001  + p2*p3*(2.*zeta - 1)/den2
2002  - 2.*p2*p3*(xi*eta - zeta + zeta*zeta)/den3;
2003 
2004  case 3:
2005  return 0.5*p4*(xi*eta + zeta - zeta*zeta)/den2
2006  + 0.5*p3*(xi*eta + zeta - zeta*zeta)/den2
2007  - p3*p4*(1 - 2.*zeta)/den2
2008  + 2.*p3*p4*(xi*eta + zeta - zeta*zeta)/den3;
2009 
2010  case 4:
2011  return 4.*zeta - 1.;
2012 
2013  case 5:
2014  return 2.*p4*p1*eta/den2
2015  + 2.*p2*p4*eta/den2
2016  + 2.*p1*p2*eta/den2
2017  + 8.*p2*p1*p4*eta/den3;
2018 
2019  case 6:
2020  return -2.*p2*p3*xi/den2
2021  - 2.*p1*p3*xi/den2
2022  - 2.*p1*p2*xi/den2
2023  - 8.*p1*p2*p3*xi/den3;
2024 
2025  case 7:
2026  return -2.*p3*p4*eta/den2
2027  - 2.*p2*p4*eta/den2
2028  - 2.*p2*p3*eta/den2
2029  - 8.*p2*p3*p4*eta/den3;
2030 
2031  case 8:
2032  return 2.*p4*p1*xi/den2
2033  + 2.*p1*p3*xi/den2
2034  + 2.*p3*p4*xi/den2
2035  + 8.*p3*p4*p1*xi/den3;
2036 
2037  case 9:
2038  return 2.*p4*zeta/den
2039  + 2.*p1*zeta/den
2040  - 4.*p1*p4/den
2041  + 4.*p1*p4*zeta/den2;
2042 
2043  case 10:
2044  return 2.*p1*zeta/den
2045  + 2.*p2*zeta/den
2046  - 4.*p2*p1/den
2047  + 4.*p2*p1*zeta/den2;
2048 
2049  case 11:
2050  return 2.*p2*zeta/den
2051  + 2.*p3*zeta/den
2052  - 4.*p3*p2/den
2053  + 4.*p3*p2*zeta/den2;
2054 
2055  case 12:
2056  return 2.*p3*zeta/den
2057  + 2.*p4*zeta/den
2058  - 4.*p4*p3/den
2059  + 4.*p4*p3*zeta/den2;
2060 
2061  case 13:
2062  return -8.*p2*p3*p4/den2
2063  - 8.*p3*p4*p1/den2
2064  - 8.*p2*p1*p4/den2
2065  - 8.*p1*p2*p3/den2
2066  - 32.*p1*p2*p3*p4/den3;
2067 
2068  default:
2069  libmesh_error_msg("Invalid i = " << i);
2070  }
2071  }
2072 
2073  default:
2074  libmesh_error_msg("Invalid j = " << j);
2075  }
2076  }
2077 
2078 
2079  default:
2080  libmesh_error_msg("ERROR: Unsupported 3D element type!: " << type);
2081  }
2082  }
2083 
2084 
2085  // unsupported order
2086  default:
2087  libmesh_error_msg("ERROR: Unsupported 3D FE order!: " << order);
2088  }
2089 
2090 #else // LIBMESH_DIM != 3
2091  libmesh_ignore(type, order, i, j, p);
2092  libmesh_not_implemented();
2093 #endif
2094 }
2095 
2096 
2097 
2098 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
2099 
2100 Real fe_lagrange_3D_shape_second_deriv(const ElemType type,
2101  const Order order,
2102  const unsigned int i,
2103  const unsigned int j,
2104  const Point & p)
2105 {
2106 #if LIBMESH_DIM == 3
2107 
2108  libmesh_assert_less (j, 6);
2109 
2110  switch (order)
2111  {
2112  // linear Lagrange shape functions
2113  case FIRST:
2114  {
2115  switch (type)
2116  {
2117  // Linear tets have all second derivatives = 0
2118  case TET4:
2119  case TET10:
2120  {
2121  return 0.;
2122  }
2123 
2124  // The following elements use either tensor product or
2125  // rational basis functions, and therefore probably have
2126  // second derivatives, but we have not implemented them
2127  // yet...
2128  case PRISM6:
2129  case PRISM15:
2130  case PRISM18:
2131  {
2132  libmesh_assert_less (i, 6);
2133 
2134  // Compute prism shape functions as a tensor-product
2135  // of a triangle and an edge
2136 
2137  Point p2d(p(0),p(1));
2138  Real p1d = p(2);
2139 
2140  // 0 1 2 3 4 5
2141  static const unsigned int i0[] = {0, 0, 0, 1, 1, 1};
2142  static const unsigned int i1[] = {0, 1, 2, 0, 1, 2};
2143 
2144  switch (j)
2145  {
2146  // All repeated second derivatives and the xi-eta derivative are zero on PRISMs
2147  case 0: // d^2()/dxi^2
2148  case 1: // d^2()/dxideta
2149  case 2: // d^2()/deta^2
2150  case 5: // d^2()/dzeta^2
2151  {
2152  return 0.;
2153  }
2154 
2155  case 3: // d^2()/dxidzeta
2156  return (FE<2,LAGRANGE>::shape_deriv(TRI3, FIRST, i1[i], 0, p2d)*
2157  fe_lagrange_1D_linear_shape_deriv(i0[i], 0, p1d));
2158 
2159  case 4: // d^2()/detadzeta
2160  return (FE<2,LAGRANGE>::shape_deriv(TRI3, FIRST, i1[i], 1, p2d)*
2161  fe_lagrange_1D_linear_shape_deriv(i0[i], 0, p1d));
2162 
2163  default:
2164  libmesh_error_msg("Invalid j = " << j);
2165  }
2166  }
2167 
2168  case PYRAMID5:
2169  case PYRAMID13:
2170  case PYRAMID14:
2171  {
2172  libmesh_assert_less (i, 5);
2173 
2174  const Real xi = p(0);
2175  const Real eta = p(1);
2176  const Real zeta = p(2);
2177  const Real eps = 1.e-35;
2178 
2179  switch (j)
2180  {
2181  // xi-xi and eta-eta derivatives are all zero for PYRAMID5.
2182  case 0: // d^2()/dxi^2
2183  case 2: // d^2()/deta^2
2184  return 0.;
2185 
2186  case 1: // d^2()/dxideta
2187  {
2188  switch (i)
2189  {
2190  case 0:
2191  case 2:
2192  return 0.25/(1. - zeta + eps);
2193  case 1:
2194  case 3:
2195  return -0.25/(1. - zeta + eps);
2196  case 4:
2197  return 0.;
2198  default:
2199  libmesh_error_msg("Invalid i = " << i);
2200  }
2201  }
2202 
2203  case 3: // d^2()/dxidzeta
2204  {
2205  Real den = (1. - zeta + eps)*(1. - zeta + eps);
2206 
2207  switch (i)
2208  {
2209  case 0:
2210  case 2:
2211  return 0.25*eta/den;
2212  case 1:
2213  case 3:
2214  return -0.25*eta/den;
2215  case 4:
2216  return 0.;
2217  default:
2218  libmesh_error_msg("Invalid i = " << i);
2219  }
2220  }
2221 
2222  case 4: // d^2()/detadzeta
2223  {
2224  Real den = (1. - zeta + eps)*(1. - zeta + eps);
2225 
2226  switch (i)
2227  {
2228  case 0:
2229  case 2:
2230  return 0.25*xi/den;
2231  case 1:
2232  case 3:
2233  return -0.25*xi/den;
2234  case 4:
2235  return 0.;
2236  default:
2237  libmesh_error_msg("Invalid i = " << i);
2238  }
2239  }
2240 
2241  case 5: // d^2()/dzeta^2
2242  {
2243  Real den = (1. - zeta + eps)*(1. - zeta + eps)*(1. - zeta + eps);
2244 
2245  switch (i)
2246  {
2247  case 0:
2248  case 2:
2249  return 0.5*xi*eta/den;
2250  case 1:
2251  case 3:
2252  return -0.5*xi*eta/den;
2253  case 4:
2254  return 0.;
2255  default:
2256  libmesh_error_msg("Invalid i = " << i);
2257  }
2258  }
2259 
2260  default:
2261  libmesh_error_msg("Invalid j = " << j);
2262  }
2263  }
2264 
2265  // Trilinear shape functions on HEX8s have nonzero mixed second derivatives
2266  case HEX8:
2267  case HEX20:
2268  case HEX27:
2269  {
2270  libmesh_assert_less (i, 8);
2271 
2272  // Compute hex shape functions as a tensor-product
2273  const Real xi = p(0);
2274  const Real eta = p(1);
2275  const Real zeta = p(2);
2276 
2277  static const unsigned int i0[] = {0, 1, 1, 0, 0, 1, 1, 0};
2278  static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 1, 1};
2279  static const unsigned int i2[] = {0, 0, 0, 0, 1, 1, 1, 1};
2280 
2281  switch (j)
2282  {
2283  // All repeated second derivatives are zero on HEX8
2284  case 0: // d^2()/dxi^2
2285  case 2: // d^2()/deta^2
2286  case 5: // d^2()/dzeta^2
2287  {
2288  return 0.;
2289  }
2290 
2291  case 1: // d^2()/dxideta
2292  return (fe_lagrange_1D_linear_shape_deriv(i0[i], 0, xi)*
2293  fe_lagrange_1D_linear_shape_deriv(i1[i], 0, eta)*
2294  fe_lagrange_1D_linear_shape (i2[i], zeta));
2295 
2296  case 3: // d^2()/dxidzeta
2297  return (fe_lagrange_1D_linear_shape_deriv(i0[i], 0, xi)*
2298  fe_lagrange_1D_linear_shape (i1[i], eta)*
2299  fe_lagrange_1D_linear_shape_deriv(i2[i], 0, zeta));
2300 
2301  case 4: // d^2()/detadzeta
2302  return (fe_lagrange_1D_linear_shape (i0[i], xi)*
2303  fe_lagrange_1D_linear_shape_deriv(i1[i], 0, eta)*
2304  fe_lagrange_1D_linear_shape_deriv(i2[i], 0, zeta));
2305 
2306  default:
2307  libmesh_error_msg("Invalid j = " << j);
2308  }
2309  }
2310 
2311  default:
2312  libmesh_error_msg("ERROR: Unsupported 3D element type!: " << type);
2313  }
2314 
2315  }
2316 
2317  // quadratic Lagrange shape functions
2318  case SECOND:
2319  {
2320  switch (type)
2321  {
2322 
2323  // serendipity hexahedral quadratic shape functions
2324  case HEX20:
2325  {
2326  libmesh_assert_less (i, 20);
2327 
2328  const Real xi = p(0);
2329  const Real eta = p(1);
2330  const Real zeta = p(2);
2331 
2332  // these functions are defined for (x,y,z) in [0,1]^3
2333  // so transform the locations
2334  const Real x = .5*(xi + 1.);
2335  const Real y = .5*(eta + 1.);
2336  const Real z = .5*(zeta + 1.);
2337 
2338  switch(j)
2339  {
2340  case 0: // d^2()/dxi^2
2341  {
2342  switch(i)
2343  {
2344  case 0:
2345  case 1:
2346  return (1. - y) * (1. - z);
2347  case 2:
2348  case 3:
2349  return y * (1. - z);
2350  case 4:
2351  case 5:
2352  return (1. - y) * z;
2353  case 6:
2354  case 7:
2355  return y * z;
2356  case 8:
2357  return -2. * (1. - y) * (1. - z);
2358  case 10:
2359  return -2. * y * (1. - z);
2360  case 16:
2361  return -2. * (1. - y) * z;
2362  case 18:
2363  return -2. * y * z;
2364  case 9:
2365  case 11:
2366  case 12:
2367  case 13:
2368  case 14:
2369  case 15:
2370  case 17:
2371  case 19:
2372  return 0;
2373  default:
2374  libmesh_error_msg("Invalid i = " << i);
2375  }
2376  }
2377  case 1: // d^2()/dxideta
2378  {
2379  switch(i)
2380  {
2381  case 0:
2382  return (1.25 - x - y - .5*z) * (1. - z);
2383  case 1:
2384  return (-x + y + .5*z - .25) * (1. - z);
2385  case 2:
2386  return (x + y - .5*z - .75) * (1. - z);
2387  case 3:
2388  return (-y + x + .5*z - .25) * (1. - z);
2389  case 4:
2390  return -.25*z * (4.*x + 4.*y - 2.*z - 3);
2391  case 5:
2392  return -.25*z * (-4.*y + 4.*x + 2.*z - 1.);
2393  case 6:
2394  return .25*z * (-5 + 4.*x + 4.*y + 2.*z);
2395  case 7:
2396  return .25*z * (4.*x - 4.*y - 2.*z + 1.);
2397  case 8:
2398  return (-1. + 2.*x) * (1. - z);
2399  case 9:
2400  return (1. - 2.*y) * (1. - z);
2401  case 10:
2402  return (1. - 2.*x) * (1. - z);
2403  case 11:
2404  return (-1. + 2.*y) * (1. - z);
2405  case 12:
2406  return z * (1. - z);
2407  case 13:
2408  return -z * (1. - z);
2409  case 14:
2410  return z * (1. - z);
2411  case 15:
2412  return -z * (1. - z);
2413  case 16:
2414  return (-1. + 2.*x) * z;
2415  case 17:
2416  return (1. - 2.*y) * z;
2417  case 18:
2418  return (1. - 2.*x) * z;
2419  case 19:
2420  return (-1. + 2.*y) * z;
2421  default:
2422  libmesh_error_msg("Invalid i = " << i);
2423  }
2424  }
2425  case 2: // d^2()/deta^2
2426  switch(i)
2427  {
2428  case 0:
2429  case 3:
2430  return (1. - x) * (1. - z);
2431  case 1:
2432  case 2:
2433  return x * (1. - z);
2434  case 4:
2435  case 7:
2436  return (1. - x) * z;
2437  case 5:
2438  case 6:
2439  return x * z;
2440  case 9:
2441  return -2. * x * (1. - z);
2442  case 11:
2443  return -2. * (1. - x) * (1. - z);
2444  case 17:
2445  return -2. * x * z;
2446  case 19:
2447  return -2. * (1. - x) * z;
2448  case 8:
2449  case 10:
2450  case 12:
2451  case 13:
2452  case 14:
2453  case 15:
2454  case 16:
2455  case 18:
2456  return 0.;
2457  default:
2458  libmesh_error_msg("Invalid i = " << i);
2459  }
2460  case 3: // d^2()/dxidzeta
2461  switch(i)
2462  {
2463  case 0:
2464  return (1.25 - x - .5*y - z) * (1. - y);
2465  case 1:
2466  return (-x + .5*y + z - .25) * (1. - y);
2467  case 2:
2468  return -.25*y * (2.*y + 4.*x - 4.*z - 1.);
2469  case 3:
2470  return -.25*y * (-2.*y + 4.*x + 4.*z - 3);
2471  case 4:
2472  return (-z + x + .5*y - .25) * (1. - y);
2473  case 5:
2474  return (x - .5*y + z - .75) * (1. - y);
2475  case 6:
2476  return .25*y * (2.*y + 4.*x + 4.*z - 5);
2477  case 7:
2478  return .25*y * (-2.*y + 4.*x - 4.*z + 1.);
2479  case 8:
2480  return (-1. + 2.*x) * (1. - y);
2481  case 9:
2482  return -y * (1. - y);
2483  case 10:
2484  return (-1. + 2.*x) * y;
2485  case 11:
2486  return y * (1. - y);
2487  case 12:
2488  return (-1. + 2.*z) * (1. - y);
2489  case 13:
2490  return (1. - 2.*z) * (1. - y);
2491  case 14:
2492  return (1. - 2.*z) * y;
2493  case 15:
2494  return (-1. + 2.*z) * y;
2495  case 16:
2496  return (1. - 2.*x) * (1. - y);
2497  case 17:
2498  return y * (1. - y);
2499  case 18:
2500  return (1. - 2.*x) * y;
2501  case 19:
2502  return -y * (1. - y);
2503  default:
2504  libmesh_error_msg("Invalid i = " << i);
2505  }
2506  case 4: // d^2()/detadzeta
2507  switch(i)
2508  {
2509  case 0:
2510  return (1.25 - .5*x - y - z) * (1. - x);
2511  case 1:
2512  return .25*x * (2.*x - 4.*y - 4.*z + 3.);
2513  case 2:
2514  return -.25*x * (2.*x + 4.*y - 4.*z - 1.);
2515  case 3:
2516  return (-y + .5*x + z - .25) * (1. - x);
2517  case 4:
2518  return (-z + .5*x + y - .25) * (1. - x);
2519  case 5:
2520  return -.25*x * (2.*x - 4.*y + 4.*z - 1.);
2521  case 6:
2522  return .25*x * (2.*x + 4.*y + 4.*z - 5);
2523  case 7:
2524  return (y - .5*x + z - .75) * (1. - x);
2525  case 8:
2526  return x * (1. - x);
2527  case 9:
2528  return (-1. + 2.*y) * x;
2529  case 10:
2530  return -x * (1. - x);
2531  case 11:
2532  return (-1. + 2.*y) * (1. - x);
2533  case 12:
2534  return (-1. + 2.*z) * (1. - x);
2535  case 13:
2536  return (-1. + 2.*z) * x;
2537  case 14:
2538  return (1. - 2.*z) * x;
2539  case 15:
2540  return (1. - 2.*z) * (1. - x);
2541  case 16:
2542  return -x * (1. - x);
2543  case 17:
2544  return (1. - 2.*y) * x;
2545  case 18:
2546  return x * (1. - x);
2547  case 19:
2548  return (1. - 2.*y) * (1. - x);
2549  default:
2550  libmesh_error_msg("Invalid i = " << i);
2551  }
2552  case 5: // d^2()/dzeta^2
2553  switch(i)
2554  {
2555  case 0:
2556  case 4:
2557  return (1. - x) * (1. - y);
2558  case 1:
2559  case 5:
2560  return x * (1. - y);
2561  case 2:
2562  case 6:
2563  return x * y;
2564  case 3:
2565  case 7:
2566  return (1. - x) * y;
2567  case 12:
2568  return -2. * (1. - x) * (1. - y);
2569  case 13:
2570  return -2. * x * (1. - y);
2571  case 14:
2572  return -2. * x * y;
2573  case 15:
2574  return -2. * (1. - x) * y;
2575  case 8:
2576  case 9:
2577  case 10:
2578  case 11:
2579  case 16:
2580  case 17:
2581  case 18:
2582  case 19:
2583  return 0.;
2584  default:
2585  libmesh_error_msg("Invalid i = " << i);
2586  }
2587  default:
2588  libmesh_error_msg("Invalid j = " << j);
2589  }
2590  }
2591 
2592  // triquadratic hexahedral shape functions
2593  case HEX27:
2594  {
2595  libmesh_assert_less (i, 27);
2596 
2597  // Compute hex shape functions as a tensor-product
2598  const Real xi = p(0);
2599  const Real eta = p(1);
2600  const Real zeta = p(2);
2601 
2602  // The only way to make any sense of this
2603  // is to look at the mgflo/mg2/mgf documentation
2604  // and make the cut-out cube!
2605  // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26
2606  static const unsigned int i0[] = {0, 1, 1, 0, 0, 1, 1, 0, 2, 1, 2, 0, 0, 1, 1, 0, 2, 1, 2, 0, 2, 2, 1, 2, 0, 2, 2};
2607  static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 1, 1, 0, 2, 1, 2, 0, 0, 1, 1, 0, 2, 1, 2, 2, 0, 2, 1, 2, 2, 2};
2608  static const unsigned int i2[] = {0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 2, 2, 2, 2, 1, 1, 1, 1, 0, 2, 2, 2, 2, 1, 2};
2609 
2610  switch(j)
2611  {
2612  // d^2()/dxi^2
2613  case 0:
2614  return (fe_lagrange_1D_quadratic_shape_second_deriv(i0[i], 0, xi)*
2615  fe_lagrange_1D_quadratic_shape (i1[i], eta)*
2616  fe_lagrange_1D_quadratic_shape (i2[i], zeta));
2617 
2618  // d^2()/dxideta
2619  case 1:
2620  return (fe_lagrange_1D_quadratic_shape_deriv(i0[i], 0, xi)*
2621  fe_lagrange_1D_quadratic_shape_deriv(i1[i], 0, eta)*
2622  fe_lagrange_1D_quadratic_shape (i2[i], zeta));
2623 
2624  // d^2()/deta^2
2625  case 2:
2626  return (fe_lagrange_1D_quadratic_shape (i0[i], xi)*
2627  fe_lagrange_1D_quadratic_shape_second_deriv(i1[i], 0, eta)*
2628  fe_lagrange_1D_quadratic_shape (i2[i], zeta));
2629 
2630  // d^2()/dxidzeta
2631  case 3:
2632  return (fe_lagrange_1D_quadratic_shape_deriv(i0[i], 0, xi)*
2633  fe_lagrange_1D_quadratic_shape (i1[i], eta)*
2634  fe_lagrange_1D_quadratic_shape_deriv(i2[i], 0, zeta));
2635 
2636  // d^2()/detadzeta
2637  case 4:
2638  return (fe_lagrange_1D_quadratic_shape (i0[i], xi)*
2639  fe_lagrange_1D_quadratic_shape_deriv(i1[i], 0, eta)*
2640  fe_lagrange_1D_quadratic_shape_deriv(i2[i], 0, zeta));
2641 
2642  // d^2()/dzeta^2
2643  case 5:
2644  return (fe_lagrange_1D_quadratic_shape (i0[i], xi)*
2645  fe_lagrange_1D_quadratic_shape (i1[i], eta)*
2646  fe_lagrange_1D_quadratic_shape_second_deriv(i2[i], 0, zeta));
2647 
2648  default:
2649  libmesh_error_msg("Invalid j = " << j);
2650  }
2651  }
2652 
2653  // quadratic tetrahedral shape functions
2654  case TET10:
2655  {
2656  // The area coordinates are the same as used for the
2657  // shape() and shape_deriv() functions.
2658  // const Real zeta0 = 1. - zeta1 - zeta2 - zeta3;
2659  // const Real zeta1 = p(0);
2660  // const Real zeta2 = p(1);
2661  // const Real zeta3 = p(2);
2662  static const Real dzetadxi[4][3] =
2663  {
2664  {-1., -1., -1.},
2665  {1., 0., 0.},
2666  {0., 1., 0.},
2667  {0., 0., 1.}
2668  };
2669 
2670  // Convert from j -> (j,k) indices for independent variable
2671  // (0=xi, 1=eta, 2=zeta)
2672  static const unsigned short int independent_var_indices[6][2] =
2673  {
2674  {0, 0}, // d^2 phi / dxi^2
2675  {0, 1}, // d^2 phi / dxi deta
2676  {1, 1}, // d^2 phi / deta^2
2677  {0, 2}, // d^2 phi / dxi dzeta
2678  {1, 2}, // d^2 phi / deta dzeta
2679  {2, 2} // d^2 phi / dzeta^2
2680  };
2681 
2682  // Convert from i -> zeta indices. Each quadratic shape
2683  // function for the Tet10 depends on up to two of the zeta
2684  // area coordinate functions (see the shape() function above).
2685  // This table just tells which two area coords it uses.
2686  static const unsigned short int zeta_indices[10][2] =
2687  {
2688  {0, 0},
2689  {1, 1},
2690  {2, 2},
2691  {3, 3},
2692  {0, 1},
2693  {1, 2},
2694  {2, 0},
2695  {0, 3},
2696  {1, 3},
2697  {2, 3},
2698  };
2699 
2700  // Look up the independent variable indices for this value of j.
2701  const unsigned int my_j = independent_var_indices[j][0];
2702  const unsigned int my_k = independent_var_indices[j][1];
2703 
2704  if (i<4)
2705  {
2706  return 4.*dzetadxi[i][my_j]*dzetadxi[i][my_k];
2707  }
2708 
2709  else if (i<10)
2710  {
2711  const unsigned short int my_m = zeta_indices[i][0];
2712  const unsigned short int my_n = zeta_indices[i][1];
2713 
2714  return 4.*(dzetadxi[my_n][my_j]*dzetadxi[my_m][my_k] +
2715  dzetadxi[my_m][my_j]*dzetadxi[my_n][my_k] );
2716  }
2717  else
2718  libmesh_error_msg("Invalid shape function index " << i);
2719  }
2720 
2721 
2722 
2723  // "serendipity" prism
2724  case PRISM15:
2725  {
2726  libmesh_assert_less (i, 15);
2727 
2728  const Real xi = p(0);
2729  const Real eta = p(1);
2730  const Real zeta = p(2);
2731 
2732  switch (j)
2733  {
2734  // d^2()/dxi^2
2735  case 0:
2736  {
2737  switch(i)
2738  {
2739  case 0:
2740  case 1:
2741  return 2.*(1. - zeta);
2742  case 2:
2743  case 5:
2744  case 7:
2745  case 8:
2746  case 9:
2747  case 10:
2748  case 11:
2749  case 13:
2750  case 14:
2751  return 0.;
2752  case 3:
2753  case 4:
2754  return 2.*(1. + zeta);
2755  case 6:
2756  return 4.*(zeta - 1);
2757  case 12:
2758  return -4.*(1. + zeta);
2759  default:
2760  libmesh_error_msg("Invalid i = " << i);
2761  }
2762  }
2763 
2764  // d^2()/dxideta
2765  case 1:
2766  {
2767  switch(i)
2768  {
2769  case 0:
2770  case 7:
2771  return 2.*(1. - zeta);
2772  case 1:
2773  case 2:
2774  case 4:
2775  case 5:
2776  case 9:
2777  case 10:
2778  case 11:
2779  return 0.;
2780  case 3:
2781  case 13:
2782  return 2.*(1. + zeta);
2783  case 6:
2784  case 8:
2785  return 2.*(zeta - 1.);
2786  case 12:
2787  case 14:
2788  return -2.*(1. + zeta);
2789  default:
2790  libmesh_error_msg("Invalid i = " << i);
2791  }
2792  }
2793 
2794  // d^2()/deta^2
2795  case 2:
2796  {
2797  switch(i)
2798  {
2799  case 0:
2800  case 2:
2801  return 2.*(1. - zeta);
2802  case 1:
2803  case 4:
2804  case 6:
2805  case 7:
2806  case 9:
2807  case 10:
2808  case 11:
2809  case 12:
2810  case 13:
2811  return 0.;
2812  case 3:
2813  case 5:
2814  return 2.*(1. + zeta);
2815  case 8:
2816  return 4.*(zeta - 1.);
2817  case 14:
2818  return -4.*(1. + zeta);
2819  default:
2820  libmesh_error_msg("Invalid i = " << i);
2821  }
2822  }
2823 
2824  // d^2()/dxidzeta
2825  case 3:
2826  {
2827  switch(i)
2828  {
2829  case 0:
2830  return 1.5 - zeta - 2.*xi - 2.*eta;
2831  case 1:
2832  return 0.5 + zeta - 2.*xi;
2833  case 2:
2834  case 5:
2835  case 11:
2836  return 0.;
2837  case 3:
2838  return -1.5 - zeta + 2.*xi + 2.*eta;
2839  case 4:
2840  return -0.5 + zeta + 2.*xi;
2841  case 6:
2842  return 4.*xi + 2.*eta - 2.;
2843  case 7:
2844  return -2.*eta;
2845  case 8:
2846  return 2.*eta;
2847  case 9:
2848  return 2.*zeta;
2849  case 10:
2850  return -2.*zeta;
2851  case 12:
2852  return -4.*xi - 2.*eta + 2.;
2853  case 13:
2854  return 2.*eta;
2855  case 14:
2856  return -2.*eta;
2857  default:
2858  libmesh_error_msg("Invalid i = " << i);
2859  }
2860  }
2861 
2862  // d^2()/detadzeta
2863  case 4:
2864  {
2865  switch(i)
2866  {
2867  case 0:
2868  return 1.5 - zeta - 2.*xi - 2.*eta;
2869  case 1:
2870  case 4:
2871  case 10:
2872  return 0.;
2873  case 2:
2874  return .5 + zeta - 2.*eta;
2875  case 3:
2876  return -1.5 - zeta + 2.*xi + 2.*eta;
2877  case 5:
2878  return -.5 + zeta + 2.*eta;
2879  case 6:
2880  return 2.*xi;
2881  case 7:
2882  return -2.*xi;
2883  case 8:
2884  return 2.*xi + 4.*eta - 2.;
2885  case 9:
2886  return 2.*zeta;
2887  case 11:
2888  return -2.*zeta;
2889  case 12:
2890  return -2.*xi;
2891  case 13:
2892  return 2.*xi;
2893  case 14:
2894  return -2.*xi - 4.*eta + 2.;
2895  default:
2896  libmesh_error_msg("Invalid i = " << i);
2897  }
2898  }
2899 
2900  // d^2()/dzeta^2
2901  case 5:
2902  {
2903  switch(i)
2904  {
2905  case 0:
2906  case 3:
2907  return 1. - xi - eta;
2908  case 1:
2909  case 4:
2910  return xi;
2911  case 2:
2912  case 5:
2913  return eta;
2914  case 6:
2915  case 7:
2916  case 8:
2917  case 12:
2918  case 13:
2919  case 14:
2920  return 0.;
2921  case 9:
2922  return 2.*xi + 2.*eta - 2.;
2923  case 10:
2924  return -2.*xi;
2925  case 11:
2926  return -2.*eta;
2927  default:
2928  libmesh_error_msg("Invalid i = " << i);
2929  }
2930  }
2931 
2932  default:
2933  libmesh_error_msg("Invalid j = " << j);
2934  }
2935  }
2936 
2937 
2938 
2939  // quadratic prism shape functions
2940  case PRISM18:
2941  {
2942  libmesh_assert_less (i, 18);
2943 
2944  // Compute prism shape functions as a tensor-product
2945  // of a triangle and an edge
2946 
2947  Point p2d(p(0),p(1));
2948  Real p1d = p(2);
2949 
2950  // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
2951  static const unsigned int i0[] = {0, 0, 0, 1, 1, 1, 0, 0, 0, 2, 2, 2, 1, 1, 1, 2, 2, 2};
2952  static const unsigned int i1[] = {0, 1, 2, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 3, 4, 5};
2953 
2954  switch (j)
2955  {
2956  // d^2()/dxi^2
2957  case 0:
2958  return (FE<2,LAGRANGE>::shape_second_deriv(TRI6, SECOND, i1[i], 0, p2d)*
2959  fe_lagrange_1D_quadratic_shape(i0[i], p1d));
2960 
2961  // d^2()/dxideta
2962  case 1:
2963  return (FE<2,LAGRANGE>::shape_second_deriv(TRI6, SECOND, i1[i], 1, p2d)*
2964  fe_lagrange_1D_quadratic_shape(i0[i], p1d));
2965 
2966  // d^2()/deta^2
2967  case 2:
2968  return (FE<2,LAGRANGE>::shape_second_deriv(TRI6, SECOND, i1[i], 2, p2d)*
2969  fe_lagrange_1D_quadratic_shape(i0[i], p1d));
2970 
2971  // d^2()/dxidzeta
2972  case 3:
2973  return (FE<2,LAGRANGE>::shape_deriv(TRI6, SECOND, i1[i], 0, p2d)*
2974  fe_lagrange_1D_quadratic_shape_deriv(i0[i], 0, p1d));
2975 
2976  // d^2()/detadzeta
2977  case 4:
2978  return (FE<2,LAGRANGE>::shape_deriv(TRI6, SECOND, i1[i], 1, p2d)*
2979  fe_lagrange_1D_quadratic_shape_deriv(i0[i], 0, p1d));
2980 
2981  // d^2()/dzeta^2
2982  case 5:
2983  return (FE<2,LAGRANGE>::shape(TRI6, SECOND, i1[i], p2d)*
2984  fe_lagrange_1D_quadratic_shape_second_deriv(i0[i], 0, p1d));
2985 
2986  default:
2987  libmesh_error_msg("Invalid shape function derivative j = " << j);
2988  }
2989  }
2990 
2991 
2992  // Quadratic shape functions, as defined in R. Graglia, "Higher order
2993  // bases on pyramidal elements", IEEE Trans Antennas and Propagation,
2994  // vol 47, no 5, May 1999.
2995  case PYRAMID14:
2996  {
2997  libmesh_assert_less (i, 14);
2998 
2999  const Real xi = p(0);
3000  const Real eta = p(1);
3001  const Real zeta = p(2);
3002  const Real eps = 1.e-35;
3003 
3004  // The "normalized coordinates" defined by Graglia. These are
3005  // the planes which define the faces of the pyramid.
3006  Real
3007  p1 = 0.5*(1. - eta - zeta), // back
3008  p2 = 0.5*(1. + xi - zeta), // left
3009  p3 = 0.5*(1. + eta - zeta), // front
3010  p4 = 0.5*(1. - xi - zeta); // right
3011 
3012  // Denominators are perturbed by epsilon to avoid
3013  // divide-by-zero issues.
3014  Real
3015  den = (-1. + zeta + eps),
3016  den2 = den*den,
3017  den3 = den2*den,
3018  den4 = den2*den2;
3019 
3020  // These terms are used in several of the derivatives
3021  Real
3022  numer_mp = xi*eta - zeta + zeta*zeta,
3023  numer_pm = xi*eta + zeta - zeta*zeta;
3024 
3025  switch (j)
3026  {
3027  case 0: // d^2()/dxi^2
3028  {
3029  switch(i)
3030  {
3031  case 0:
3032  case 1:
3033  return -p1*eta/den2;
3034  case 2:
3035  case 3:
3036  return p3*eta/den2;
3037  case 4:
3038  case 9:
3039  case 10:
3040  case 11:
3041  case 12:
3042  return 0.;
3043  case 5:
3044  return 2.*p1*eta/den2;
3045  case 6:
3046  case 8:
3047  return 4.*p1*p3/den2;
3048  case 7:
3049  return -2.*p3*eta/den2;
3050  case 13:
3051  return -8.*p1*p3/den2;
3052  default:
3053  libmesh_error_msg("Invalid i = " << i);
3054  }
3055  }
3056 
3057  case 1: // d^2()/dxideta
3058  {
3059  switch(i)
3060  {
3061  case 0:
3062  return 0.25*numer_mp/den2
3063  - 0.5*p1*xi/den2
3064  - 0.5*p4*eta/den2
3065  + p4*p1/den2;
3066 
3067  case 1:
3068  return 0.25*numer_pm/den2
3069  - 0.5*p1*xi/den2
3070  + 0.5*p2*eta/den2
3071  - p1*p2/den2;
3072 
3073  case 2:
3074  return 0.25*numer_mp/den2
3075  + 0.5*p3*xi/den2
3076  + 0.5*p2*eta/den2
3077  + p2*p3/den2;
3078 
3079  case 3:
3080  return 0.25*numer_pm/den2
3081  + 0.5*p3*xi/den2
3082  - 0.5*p4*eta/den2
3083  - p3*p4/den2;
3084 
3085  case 4:
3086  return 0.;
3087 
3088  case 5:
3089  return p4*eta/den2
3090  - 2.*p4*p1/den2
3091  - p2*eta/den2
3092  + 2.*p1*p2/den2;
3093 
3094  case 6:
3095  return -p3*xi/den2
3096  + p1*xi/den2
3097  - 2.*p2*p3/den2
3098  + 2.*p1*p2/den2;
3099 
3100  case 7:
3101  return p4*eta/den2
3102  + 2.*p3*p4/den2
3103  - p2*eta/den2
3104  - 2.*p2*p3/den2;
3105 
3106  case 8:
3107  return -p3*xi/den2
3108  + p1*xi/den2
3109  - 2.*p4*p1/den2
3110  + 2.*p3*p4/den2;
3111 
3112  case 9:
3113  case 11:
3114  return -zeta/den;
3115 
3116  case 10:
3117  case 12:
3118  return zeta/den;
3119 
3120  case 13:
3121  return 4.*p4*p1/den2
3122  - 4.*p3*p4/den2
3123  + 4.*p2*p3/den2
3124  - 4.*p1*p2/den2;
3125 
3126  default:
3127  libmesh_error_msg("Invalid i = " << i);
3128  }
3129  }
3130 
3131 
3132  case 2: // d^2()/deta^2
3133  {
3134  switch(i)
3135  {
3136  case 0:
3137  case 3:
3138  return -p4*xi/den2;
3139  case 1:
3140  case 2:
3141  return p2*xi/den2;
3142  case 4:
3143  case 9:
3144  case 10:
3145  case 11:
3146  case 12:
3147  return 0.;
3148  case 5:
3149  case 7:
3150  return 4.*p2*p4/den2;
3151  case 6:
3152  return -2.*p2*xi/den2;
3153  case 8:
3154  return 2.*p4*xi/den2;
3155  case 13:
3156  return -8.*p2*p4/den2;
3157  default:
3158  libmesh_error_msg("Invalid i = " << i);
3159  }
3160  }
3161 
3162 
3163  case 3: // d^2()/dxidzeta
3164  {
3165  switch(i)
3166  {
3167  case 0:
3168  return 0.25*numer_mp/den2
3169  - 0.5*p1*(2.*zeta - 1.)/den2
3170  + p1*numer_mp/den3
3171  - 0.5*p1*eta/den2
3172  - 0.5*p4*eta/den2
3173  - 2.*p4*p1*eta/den3;
3174 
3175  case 1:
3176  return 0.25*numer_pm/den2
3177  - 0.5*p1*(1 - 2.*zeta)/den2
3178  + p1*numer_pm/den3
3179  + 0.5*p2*eta/den2
3180  + 0.5*p1*eta/den2
3181  + 2.*p1*p2*eta/den3;
3182 
3183  case 2:
3184  return -0.25*numer_mp/den2
3185  + 0.5*p3*(2.*zeta - 1.)/den2
3186  - p3*numer_mp/den3
3187  - 0.5*p3*eta/den2
3188  - 0.5*p2*eta/den2
3189  - 2.*p2*p3*eta/den3;
3190 
3191  case 3:
3192  return -0.25*numer_pm/den2
3193  + 0.5*p3*(1 - 2.*zeta)/den2
3194  - p3*numer_pm/den3
3195  + 0.5*p4*eta/den2
3196  + 0.5*p3*eta/den2
3197  + 2.*p3*p4*eta/den3;
3198 
3199  case 4:
3200  return 0.;
3201 
3202  case 5:
3203  return p4*eta/den2
3204  + 4.*p4*p1*eta/den3
3205  - p2*eta/den2
3206  - 4.*p1*p2*eta/den3;
3207 
3208  case 6:
3209  return -p3*xi/den2
3210  - p1*xi/den2
3211  - 4.*p1*p3*xi/den3
3212  - 2.*p2*p3/den2
3213  - 2.*p1*p3/den2
3214  - 2.*p1*p2/den2
3215  - 8.*p1*p2*p3/den3;
3216 
3217  case 7:
3218  return -p4*eta/den2
3219  - 4.*p3*p4*eta/den3
3220  + p2*eta/den2
3221  + 4.*p2*p3*eta/den3;
3222 
3223  case 8:
3224  return -p3*xi/den2
3225  - p1*xi/den2
3226  - 4.*p1*p3*xi/den3
3227  + 2.*p4*p1/den2
3228  + 2.*p1*p3/den2
3229  + 2.*p3*p4/den2
3230  + 8.*p3*p4*p1/den3;
3231 
3232  case 9:
3233  return -zeta/den
3234  + 2.*p1/den
3235  - 2.*p1*zeta/den2;
3236 
3237  case 10:
3238  return zeta/den
3239  - 2.*p1/den
3240  + 2.*p1*zeta/den2;
3241 
3242  case 11:
3243  return zeta/den
3244  - 2.*p3/den
3245  + 2.*p3*zeta/den2;
3246 
3247  case 12:
3248  return -zeta/den
3249  + 2.*p3/den
3250  - 2.*p3*zeta/den2;
3251 
3252  case 13:
3253  return -4.*p4*p1/den2
3254  - 4.*p3*p4/den2
3255  - 16.*p3*p4*p1/den3
3256  + 4.*p2*p3/den2
3257  + 4.*p1*p2/den2
3258  + 16.*p1*p2*p3/den3;
3259 
3260  default:
3261  libmesh_error_msg("Invalid i = " << i);
3262  }
3263  }
3264 
3265  case 4: // d^2()/detadzeta
3266  {
3267  switch(i)
3268  {
3269  case 0:
3270  return 0.25*numer_mp/den2
3271  - 0.5*p4*(2.*zeta - 1.)/den2
3272  + p4*numer_mp/den3
3273  - 0.5*p1*xi/den2
3274  - 0.5*p4*xi/den2
3275  - 2.*p4*p1*xi/den3;
3276 
3277  case 1:
3278  return -0.25*numer_pm/den2
3279  + 0.5*p2*(1. - 2.*zeta)/den2
3280  - p2*numer_pm/den3
3281  + 0.5*p2*xi/den2
3282  + 0.5*p1*xi/den2
3283  + 2.*p1*p2*xi/den3;
3284 
3285  case 2:
3286  return -0.25*numer_mp/den2
3287  + 0.5*p2*(2.*zeta - 1.)/den2
3288  - p2*numer_mp/den3
3289  - 0.5*p3*xi/den2
3290  - 0.5*p2*xi/den2
3291  - 2.*p2*p3*xi/den3;
3292 
3293  case 3:
3294  return 0.25*numer_pm/den2
3295  - 0.5*p4*(1. - 2.*zeta)/den2
3296  + p4*numer_pm/den3
3297  + 0.5*p4*xi/den2
3298  + 0.5*p3*xi/den2
3299  + 2.*p3*p4*xi/den3;
3300 
3301  case 4:
3302  return 0.;
3303 
3304  case 5:
3305  return -p4*eta/den2
3306  - p2*eta/den2
3307  - 4.*p2*p4*eta/den3
3308  + 2.*p4*p1/den2
3309  + 2.*p2*p4/den2
3310  + 2.*p1*p2/den2
3311  + 8.*p2*p1*p4/den3;
3312 
3313  case 6:
3314  return p3*xi/den2
3315  + 4.*p2*p3*xi/den3
3316  - p1*xi/den2
3317  - 4.*p1*p2*xi/den3;
3318 
3319  case 7:
3320  return -p4*eta/den2
3321  - p2*eta/den2
3322  - 4.*p2*p4*eta/den3
3323  - 2.*p3*p4/den2
3324  - 2.*p2*p4/den2
3325  - 2.*p2*p3/den2
3326  - 8.*p2*p3*p4/den3;
3327 
3328  case 8:
3329  return p1*xi/den2
3330  + 4.*p4*p1*xi/den3
3331  - p3*xi/den2
3332  - 4.*p3*p4*xi/den3;
3333 
3334  case 9:
3335  return -zeta/den
3336  + 2.*p4/den
3337  - 2.*p4*zeta/den2;
3338 
3339  case 10:
3340  return -zeta/den
3341  + 2.*p2/den
3342  - 2.*p2*zeta/den2;
3343 
3344  case 11:
3345  return zeta/den
3346  - 2.*p2/den
3347  + 2.*p2*zeta/den2;
3348 
3349  case 12:
3350  return zeta/den
3351  - 2.*p4/den
3352  + 2.*p4*zeta/den2;
3353 
3354  case 13:
3355  return 4.*p3*p4/den2
3356  + 4.*p2*p3/den2
3357  + 16.*p2*p3*p4/den3
3358  - 4.*p4*p1/den2
3359  - 4.*p1*p2/den2
3360  - 16.*p2*p1*p4/den3;
3361 
3362  default:
3363  libmesh_error_msg("Invalid i = " << i);
3364  }
3365  }
3366 
3367  case 5: // d^2()/dzeta^2
3368  {
3369  switch(i)
3370  {
3371  case 0:
3372  return 0.5*numer_mp/den2
3373  - p1*(2.*zeta - 1.)/den2
3374  + 2.*p1*numer_mp/den3
3375  - p4*(2.*zeta - 1.)/den2
3376  + 2.*p4*numer_mp/den3
3377  + 2.*p4*p1/den2
3378  - 4.*p4*p1*(2.*zeta - 1.)/den3
3379  + 6.*p4*p1*numer_mp/den4;
3380 
3381  case 1:
3382  return -0.5*numer_pm/den2
3383  + p2*(1 - 2.*zeta)/den2
3384  - 2.*p2*numer_pm/den3
3385  + p1*(1 - 2.*zeta)/den2
3386  - 2.*p1*numer_pm/den3
3387  + 2.*p1*p2/den2
3388  + 4.*p1*p2*(1 - 2.*zeta)/den3
3389  - 6.*p1*p2*numer_pm/den4;
3390 
3391  case 2:
3392  return 0.5*numer_mp/den2
3393  - p3*(2.*zeta - 1.)/den2
3394  + 2.*p3*numer_mp/den3
3395  - p2*(2.*zeta - 1.)/den2
3396  + 2.*p2*numer_mp/den3
3397  + 2.*p2*p3/den2
3398  - 4.*p2*p3*(2.*zeta - 1.)/den3
3399  + 6.*p2*p3*numer_mp/den4;
3400 
3401  case 3:
3402  return -0.5*numer_pm/den2
3403  + p4*(1 - 2.*zeta)/den2
3404  - 2.*p4*numer_pm/den3
3405  + p3*(1 - 2.*zeta)/den2
3406  - 2.*p3*numer_pm/den3
3407  + 2.*p3*p4/den2
3408  + 4.*p3*p4*(1 - 2.*zeta)/den3
3409  - 6.*p3*p4*numer_pm/den4;
3410 
3411  case 4:
3412  return 4.;
3413 
3414  case 5:
3415  return -2.*p1*eta/den2
3416  - 2.*p4*eta/den2
3417  - 8.*p4*p1*eta/den3
3418  - 2.*p2*eta/den2
3419  - 8.*p2*p4*eta/den3
3420  - 8.*p1*p2*eta/den3
3421  - 24.*p2*p1*p4*eta/den4;
3422 
3423  case 6:
3424  return 2.*p3*xi/den2
3425  + 2.*p2*xi/den2
3426  + 8.*p2*p3*xi/den3
3427  + 2.*p1*xi/den2
3428  + 8.*p1*p3*xi/den3
3429  + 8.*p1*p2*xi/den3
3430  + 24.*p1*p2*p3*xi/den4;
3431 
3432  case 7:
3433  return 2.*p4*eta/den2
3434  + 2.*p3*eta/den2
3435  + 8.*p3*p4*eta/den3
3436  + 2.*p2*eta/den2
3437  + 8.*p2*p4*eta/den3
3438  + 8.*p2*p3*eta/den3
3439  + 24.*p2*p3*p4*eta/den4;
3440 
3441  case 8:
3442  return -2.*p1*xi/den2
3443  - 2.*p4*xi/den2
3444  - 8.*p4*p1*xi/den3
3445  - 2.*p3*xi/den2
3446  - 8.*p1*p3*xi/den3
3447  - 8.*p3*p4*xi/den3
3448  - 24.*p3*p4*p1*xi/den4;
3449 
3450  case 9:
3451  return -2.*zeta/den
3452  + 4.*p4/den
3453  - 4.*p4*zeta/den2
3454  + 4.*p1/den
3455  - 4.*p1*zeta/den2
3456  + 8.*p4*p1/den2
3457  - 8.*p1*p4*zeta/den3;
3458 
3459  case 10:
3460  return -2.*zeta/den
3461  + 4.*p1/den
3462  - 4.*p1*zeta/den2
3463  + 4.*p2/den
3464  - 4.*p2*zeta/den2
3465  + 8.*p1*p2/den2
3466  - 8.*p2*p1*zeta/den3;
3467 
3468  case 11:
3469  return -2.*zeta/den
3470  + 4.*p2/den
3471  - 4.*p2*zeta/den2
3472  + 4.*p3/den
3473  - 4.*p3*zeta/den2
3474  + 8.*p2*p3/den2
3475  - 8.*p3*p2*zeta/den3;
3476 
3477  case 12:
3478  return -2.*zeta/den
3479  + 4.*p3/den
3480  - 4.*p3*zeta/den2
3481  + 4.*p4/den
3482  - 4.*p4*zeta/den2
3483  + 8.*p3*p4/den2
3484  - 8.*p4*p3*zeta/den3;
3485 
3486  case 13:
3487  return 8.*p3*p4/den2
3488  + 8.*p2*p4/den2
3489  + 8.*p2*p3/den2
3490  + 32.*p2*p3*p4/den3
3491  + 8.*p4*p1/den2
3492  + 8.*p1*p3/den2
3493  + 32.*p3*p4*p1/den3
3494  + 8.*p1*p2/den2
3495  + 32.*p2*p1*p4/den3
3496  + 32.*p1*p2*p3/den3
3497  + 96.*p1*p2*p3*p4/den4;
3498 
3499  default:
3500  libmesh_error_msg("Invalid i = " << i);
3501  }
3502  }
3503 
3504  default:
3505  libmesh_error_msg("Invalid j = " << j);
3506  }
3507  }
3508 
3509  // G. Bedrosian, "Shape functions and integration formulas for
3510  // three-dimensional finite element analysis", Int. J. Numerical
3511  // Methods Engineering, vol 35, p. 95-108, 1992.
3512  case PYRAMID13:
3513  {
3514  libmesh_assert_less (i, 13);
3515 
3516  const Real xi = p(0);
3517  const Real eta = p(1);
3518  const Real zeta = p(2);
3519  const Real eps = 1.e-35;
3520 
3521  // Denominators are perturbed by epsilon to avoid
3522  // divide-by-zero issues.
3523  Real
3524  den = (-1. + zeta + eps),
3525  den2 = den*den,
3526  den3 = den2*den,
3527  xi2 = xi*xi,
3528  eta2 = eta*eta,
3529  zeta2 = zeta*zeta,
3530  zeta3 = zeta2*zeta;
3531 
3532  switch (j)
3533  {
3534  case 0: // d^2()/dxi^2
3535  {
3536  switch(i)
3537  {
3538  case 0:
3539  case 1:
3540  return 0.5*(-1. + zeta + eta)/den;
3541 
3542  case 2:
3543  case 3:
3544  return 0.5*(-1. + zeta - eta)/den;
3545 
3546  case 4:
3547  case 6:
3548  case 8:
3549  case 9:
3550  case 10:
3551  case 11:
3552  case 12:
3553  return 0.;
3554 
3555  case 5:
3556  return (1. - eta - zeta)/den;
3557 
3558  case 7:
3559  return (1. + eta - zeta)/den;
3560 
3561  default:
3562  libmesh_error_msg("Invalid i = " << i);
3563  }
3564  }
3565 
3566  case 1: // d^2()/dxideta
3567  {
3568  switch(i)
3569  {
3570  case 0:
3571  return 0.25*(-1. + 2.*zeta + 2.*xi + 2.*eta)/den;
3572 
3573  case 1:
3574  return -0.25*(-1. + 2.*zeta - 2.*xi + 2.*eta)/den;
3575 
3576  case 2:
3577  return -0.25*(1. - 2.*zeta + 2.*xi + 2.*eta)/den;
3578 
3579  case 3:
3580  return 0.25*(1. - 2.*zeta - 2.*xi + 2.*eta)/den;
3581 
3582  case 4:
3583  return 0.;
3584 
3585  case 5:
3586  return -xi/den;
3587 
3588  case 6:
3589  return eta/den;
3590 
3591  case 7:
3592  return xi/den;
3593 
3594  case 8:
3595  return -eta/den;
3596 
3597  case 9:
3598  return -zeta/den;
3599 
3600  case 10:
3601  return zeta/den;
3602 
3603  case 11:
3604  return -zeta/den;
3605 
3606  case 12:
3607  return zeta/den;
3608 
3609  default:
3610  libmesh_error_msg("Invalid i = " << i);
3611  }
3612  }
3613 
3614 
3615  case 2: // d^2()/deta^2
3616  {
3617  switch(i)
3618  {
3619  case 0:
3620  case 3:
3621  return 0.5*(-1. + zeta + xi)/den;
3622 
3623  case 1:
3624  case 2:
3625  return 0.5*(-1. + zeta - xi)/den;
3626 
3627  case 4:
3628  case 5:
3629  case 7:
3630  case 9:
3631  case 10:
3632  case 11:
3633  case 12:
3634  return 0.;
3635 
3636  case 6:
3637  return (1. + xi - zeta)/den;
3638 
3639  case 8:
3640  return (1. - xi - zeta)/den;
3641 
3642  default:
3643  libmesh_error_msg("Invalid i = " << i);
3644  }
3645  }
3646 
3647 
3648  case 3: // d^2()/dxidzeta
3649  {
3650  switch(i)
3651  {
3652  case 0:
3653  return -0.25*(-1. + 2.*zeta - zeta2 + eta + 2.*eta*xi + eta2)/den2;
3654 
3655  case 1:
3656  return 0.25*(-1. + 2.*zeta - zeta2 + eta - 2.*eta*xi + eta2)/den2;
3657 
3658  case 2:
3659  return 0.25*(-1. + 2.*zeta - zeta2 - eta + 2.*eta*xi + eta2)/den2;
3660 
3661  case 3:
3662  return -0.25*(-1. + 2.*zeta - zeta2 - eta - 2.*eta*xi + eta2)/den2;
3663 
3664  case 4:
3665  return 0.;
3666 
3667  case 5:
3668  return eta*xi/den2;
3669 
3670  case 6:
3671  return -0.5*(1. + zeta2 + eta2 - 2.*zeta)/den2;
3672 
3673  case 7:
3674  return -eta*xi/den2;
3675 
3676  case 8:
3677  return 0.5*(1. + zeta2 + eta2 - 2.*zeta)/den2;
3678 
3679  case 9:
3680  return (-1. - zeta2 + eta + 2.*zeta)/den2;
3681 
3682  case 10:
3683  return -(-1. - zeta2 + eta + 2.*zeta)/den2;
3684 
3685  case 11:
3686  return (1. + zeta2 + eta - 2.*zeta)/den2;
3687 
3688  case 12:
3689  return -(1. + zeta2 + eta - 2.*zeta)/den2;
3690 
3691  default:
3692  libmesh_error_msg("Invalid i = " << i);
3693  }
3694  }
3695 
3696  case 4: // d^2()/detadzeta
3697  {
3698  switch(i)
3699  {
3700  case 0:
3701  return -0.25*(-1. + 2.*zeta - zeta2 + xi + 2.*eta*xi + xi2)/den2;
3702 
3703  case 1:
3704  return 0.25*(1. - 2.*zeta + zeta2 + xi + 2.*eta*xi - xi2)/den2;
3705 
3706  case 2:
3707  return 0.25*(-1. + 2.*zeta - zeta2 - xi + 2.*eta*xi + xi2)/den2;
3708 
3709  case 3:
3710  return -0.25*(1. - 2.*zeta + zeta2 - xi + 2.*eta*xi - xi2)/den2;
3711 
3712  case 4:
3713  return 0.;
3714 
3715  case 5:
3716  return 0.5*(1. + xi2 + zeta2 - 2.*zeta)/den2;
3717 
3718  case 6:
3719  return -eta*xi/den2;
3720 
3721  case 7:
3722  return -0.5*(1. + xi2 + zeta2 - 2.*zeta)/den2;
3723 
3724  case 8:
3725  return eta*xi/den2;
3726 
3727  case 9:
3728  return (-1. - zeta2 + xi + 2.*zeta)/den2;
3729 
3730  case 10:
3731  return -(1. + zeta2 + xi - 2.*zeta)/den2;
3732 
3733  case 11:
3734  return (1. + zeta2 + xi - 2.*zeta)/den2;
3735 
3736  case 12:
3737  return -(-1. - zeta2 + xi + 2.*zeta)/den2;
3738 
3739  default:
3740  libmesh_error_msg("Invalid i = " << i);
3741  }
3742  }
3743 
3744  case 5: // d^2()/dzeta^2
3745  {
3746  switch(i)
3747  {
3748  case 0:
3749  return 0.5*(xi + eta + 1.)*eta*xi/den3;
3750 
3751  case 1:
3752  return -0.5*(eta - xi + 1.)*eta*xi/den3;
3753 
3754  case 2:
3755  return -0.5*(xi + eta - 1.)*eta*xi/den3;
3756 
3757  case 3:
3758  return 0.5*(eta - xi - 1.)*eta*xi/den3;
3759 
3760  case 4:
3761  return 4.;
3762 
3763  case 5:
3764  return -(1. - 3.*zeta + 3.*zeta2 - zeta3 + eta*xi2)/den3;
3765 
3766  case 6:
3767  return (-1. + 3.*zeta - 3.*zeta2 + zeta3 + eta2*xi)/den3;
3768 
3769  case 7:
3770  return (-1. + 3.*zeta - 3.*zeta2 + zeta3 + eta*xi2)/den3;
3771 
3772  case 8:
3773  return -(1. - 3.*zeta + 3.*zeta2 - zeta3 + eta2*xi)/den3;
3774 
3775  case 9:
3776  return -2.*(-1. + 3.*zeta - 3.*zeta2 + zeta3 + eta*xi)/den3;
3777 
3778  case 10:
3779  return 2.*(1. - 3.*zeta + 3.*zeta2 - zeta3 + eta*xi)/den3;
3780 
3781  case 11:
3782  return -2.*(-1. + 3.*zeta - 3.*zeta2 + zeta3 + eta*xi)/den3;
3783 
3784  case 12:
3785  return 2.*(1. - 3.*zeta + 3.*zeta2 - zeta3 + eta*xi)/den3;
3786 
3787  default:
3788  libmesh_error_msg("Invalid i = " << i);
3789  }
3790  }
3791 
3792  default:
3793  libmesh_error_msg("Invalid j = " << j);
3794  }
3795  }
3796 
3797  default:
3798  libmesh_error_msg("ERROR: Unsupported 3D element type!: " << type);
3799  }
3800  }
3801 
3802 
3803  // unsupported order
3804  default:
3805  libmesh_error_msg("ERROR: Unsupported 3D FE order!: " << order);
3806  }
3807 
3808 #else // LIBMESH_DIM != 3
3809  libmesh_ignore(type, order, i, j, p);
3810  libmesh_not_implemented();
3811 #endif
3812 }
3813 
3814 #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
3815 
3816 
3817 } // anonymous namespace
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::HEX20
Definition: enum_elem_type.h:48
libMesh::fe_lagrange_1D_linear_shape
Real fe_lagrange_1D_linear_shape(const unsigned int i, const Real xi)
Definition: fe_lagrange_shape_1D.h:30
libMesh::PRISM6
Definition: enum_elem_type.h:50
libMesh::HEX8
Definition: enum_elem_type.h:47
libMesh
The libMesh namespace provides an interface to certain functionality in the library.
Definition: factoryfunction.C:55
libMesh::TET10
Definition: enum_elem_type.h:46
libMesh::Order
Order
Definition: enum_order.h:40
libMesh::Elem::p_level
unsigned int p_level() const
Definition: elem.h:2512
libMesh::fe_lagrange_1D_quadratic_shape
Real fe_lagrange_1D_quadratic_shape(const unsigned int i, const Real xi)
Definition: fe_lagrange_shape_1D.h:51
libMesh::SECOND
Definition: enum_order.h:43
libMesh::TET4
Definition: enum_elem_type.h:45
libMesh::PRISM15
Definition: enum_elem_type.h:51
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::HEX27
Definition: enum_elem_type.h:49
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::FE
A specific instantiation of the FEBase class.
Definition: fe.h:95
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::TRI3
Definition: enum_elem_type.h:39
libMesh::TRI6
Definition: enum_elem_type.h:40
libMesh::PYRAMID5
Definition: enum_elem_type.h:53
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::PRISM18
Definition: enum_elem_type.h:52
libMesh::Elem::type
virtual ElemType type() const =0
libMesh::FIRST
Definition: enum_order.h:42
libMesh::fe_lagrange_1D_quadratic_shape_deriv
Real fe_lagrange_1D_quadratic_shape_deriv(const unsigned int i, const unsigned int libmesh_dbg_var(j), const Real xi)
Definition: fe_lagrange_shape_1D.h:152
libMesh::fe_lagrange_1D_quadratic_shape_second_deriv
Real fe_lagrange_1D_quadratic_shape_second_deriv(const unsigned int i, const unsigned int libmesh_dbg_var(j), const Real)
Definition: fe_lagrange_shape_1D.h:238
libMesh::fe_lagrange_1D_linear_shape_deriv
Real fe_lagrange_1D_linear_shape_deriv(const unsigned int i, const unsigned int libmesh_dbg_var(j), const Real)
Definition: fe_lagrange_shape_1D.h:128
libMesh::PYRAMID13
Definition: enum_elem_type.h:54
libMesh::PYRAMID14
Definition: enum_elem_type.h:55
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