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