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