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

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

Hosted By:
SourceForge.net Logo