libMesh
fe_l2_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,L2_LAGRANGE>::shape(EDGE2, FIRST, i0[i], xi)*
65  FE<1,L2_LAGRANGE>::shape(EDGE2, FIRST, i1[i], eta)*
66  FE<1,L2_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,L2_LAGRANGE>::shape(TRI3, FIRST, i1[i], p2d)*
118  FE<1,L2_LAGRANGE>::shape(EDGE2, FIRST, i0[i], p1d));
119  }
120 
121  // linear pyramid shape functions
122  case PYRAMID5:
123  {
124  libmesh_assert_less (i, 5);
125 
126  const Real xi = p(0);
127  const Real eta = p(1);
128  const Real zeta = p(2);
129  const Real eps = 1.e-35;
130 
131  switch(i)
132  {
133  case 0:
134  return .25*(zeta + xi - 1.)*(zeta + eta - 1.)/((1. - zeta) + eps);
135 
136  case 1:
137  return .25*(zeta - xi - 1.)*(zeta + eta - 1.)/((1. - zeta) + eps);
138 
139  case 2:
140  return .25*(zeta - xi - 1.)*(zeta - eta - 1.)/((1. - zeta) + eps);
141 
142  case 3:
143  return .25*(zeta + xi - 1.)*(zeta - eta - 1.)/((1. - zeta) + eps);
144 
145  case 4:
146  return zeta;
147 
148  default:
149  libmesh_error_msg("Invalid i = " << i);
150  }
151  }
152 
153 
154  default:
155  libmesh_error_msg("ERROR: Unsupported 3D element type!: " << type);
156  }
157  }
158 
159 
160  // quadratic Lagrange shape functions
161  case SECOND:
162  {
163  switch (type)
164  {
165 
166  // serendipity hexahedral quadratic shape functions
167  case HEX20:
168  {
169  libmesh_assert_less (i, 20);
170 
171  const Real xi = p(0);
172  const Real eta = p(1);
173  const Real zeta = p(2);
174 
175  // these functions are defined for (x,y,z) in [0,1]^3
176  // so transform the locations
177  const Real x = .5*(xi + 1.);
178  const Real y = .5*(eta + 1.);
179  const Real z = .5*(zeta + 1.);
180 
181  switch (i)
182  {
183  case 0:
184  return (1. - x)*(1. - y)*(1. - z)*(1. - 2.*x - 2.*y - 2.*z);
185 
186  case 1:
187  return x*(1. - y)*(1. - z)*(2.*x - 2.*y - 2.*z - 1.);
188 
189  case 2:
190  return x*y*(1. - z)*(2.*x + 2.*y - 2.*z - 3.);
191 
192  case 3:
193  return (1. - x)*y*(1. - z)*(2.*y - 2.*x - 2.*z - 1.);
194 
195  case 4:
196  return (1. - x)*(1. - y)*z*(2.*z - 2.*x - 2.*y - 1.);
197 
198  case 5:
199  return x*(1. - y)*z*(2.*x - 2.*y + 2.*z - 3.);
200 
201  case 6:
202  return x*y*z*(2.*x + 2.*y + 2.*z - 5.);
203 
204  case 7:
205  return (1. - x)*y*z*(2.*y - 2.*x + 2.*z - 3.);
206 
207  case 8:
208  return 4.*x*(1. - x)*(1. - y)*(1. - z);
209 
210  case 9:
211  return 4.*x*y*(1. - y)*(1. - z);
212 
213  case 10:
214  return 4.*x*(1. - x)*y*(1. - z);
215 
216  case 11:
217  return 4.*(1. - x)*y*(1. - y)*(1. - z);
218 
219  case 12:
220  return 4.*(1. - x)*(1. - y)*z*(1. - z);
221 
222  case 13:
223  return 4.*x*(1. - y)*z*(1. - z);
224 
225  case 14:
226  return 4.*x*y*z*(1. - z);
227 
228  case 15:
229  return 4.*(1. - x)*y*z*(1. - z);
230 
231  case 16:
232  return 4.*x*(1. - x)*(1. - y)*z;
233 
234  case 17:
235  return 4.*x*y*(1. - y)*z;
236 
237  case 18:
238  return 4.*x*(1. - x)*y*z;
239 
240  case 19:
241  return 4.*(1. - x)*y*(1. - y)*z;
242 
243  default:
244  libmesh_error_msg("Invalid i = " << i);
245  }
246  }
247 
248  // triquadratic hexahedral shape functions
249  case HEX27:
250  {
251  libmesh_assert_less (i, 27);
252 
253  // Compute hex shape functions as a tensor-product
254  const Real xi = p(0);
255  const Real eta = p(1);
256  const Real zeta = p(2);
257 
258  // The only way to make any sense of this
259  // is to look at the mgflo/mg2/mgf documentation
260  // and make the cut-out cube!
261  // 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
262  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};
263  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};
264  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};
265 
266  return (FE<1,L2_LAGRANGE>::shape(EDGE3, SECOND, i0[i], xi)*
267  FE<1,L2_LAGRANGE>::shape(EDGE3, SECOND, i1[i], eta)*
268  FE<1,L2_LAGRANGE>::shape(EDGE3, SECOND, i2[i], zeta));
269  }
270 
271  // quadratic tetrahedral shape functions
272  case TET10:
273  {
274  libmesh_assert_less (i, 10);
275 
276  // Area coordinates, pg. 205, Vol. I, Carey, Oden, Becker FEM
277  const Real zeta1 = p(0);
278  const Real zeta2 = p(1);
279  const Real zeta3 = p(2);
280  const Real zeta0 = 1. - zeta1 - zeta2 - zeta3;
281 
282  switch(i)
283  {
284  case 0:
285  return zeta0*(2.*zeta0 - 1.);
286 
287  case 1:
288  return zeta1*(2.*zeta1 - 1.);
289 
290  case 2:
291  return zeta2*(2.*zeta2 - 1.);
292 
293  case 3:
294  return zeta3*(2.*zeta3 - 1.);
295 
296  case 4:
297  return 4.*zeta0*zeta1;
298 
299  case 5:
300  return 4.*zeta1*zeta2;
301 
302  case 6:
303  return 4.*zeta2*zeta0;
304 
305  case 7:
306  return 4.*zeta0*zeta3;
307 
308  case 8:
309  return 4.*zeta1*zeta3;
310 
311  case 9:
312  return 4.*zeta2*zeta3;
313 
314  default:
315  libmesh_error_msg("Invalid i = " << i);
316  }
317  }
318 
319  // quadratic prism shape functions
320  case PRISM18:
321  {
322  libmesh_assert_less (i, 18);
323 
324  // Compute prism shape functions as a tensor-product
325  // of a triangle and an edge
326 
327  Point p2d(p(0),p(1));
328  Point p1d(p(2));
329 
330  // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
331  static const unsigned int i0[] = {0, 0, 0, 1, 1, 1, 0, 0, 0, 2, 2, 2, 1, 1, 1, 2, 2, 2};
332  static const unsigned int i1[] = {0, 1, 2, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 3, 4, 5};
333 
334  return (FE<2,L2_LAGRANGE>::shape(TRI6, SECOND, i1[i], p2d)*
335  FE<1,L2_LAGRANGE>::shape(EDGE3, SECOND, i0[i], p1d));
336  }
337 
338 
339  default:
340  libmesh_error_msg("ERROR: Unsupported 3D element type!: " << type);
341  }
342  }
343 
344 
345  // unsupported order
346  default:
347  libmesh_error_msg("ERROR: Unsupported 3D FE order!: " << order);
348  }
349 
350 #endif
351 
352  libmesh_error_msg("We'll never get here!");
353  return 0.;
354 }
355 
356 
357 
358 template <>
360  const Order order,
361  const unsigned int i,
362  const Point & p)
363 {
364  libmesh_assert(elem);
365 
366  // call the orientation-independent shape functions
367  return FE<3,L2_LAGRANGE>::shape(elem->type(), static_cast<Order>(order + elem->p_level()), i, p);
368 }
369 
370 
371 
372 
373 template <>
375  const Order order,
376  const unsigned int i,
377  const unsigned int j,
378  const Point & p)
379 {
380 #if LIBMESH_DIM == 3
381 
382  libmesh_assert_less (j, 3);
383 
384  switch (order)
385  {
386  // linear Lagrange shape functions
387  case FIRST:
388  {
389  switch (type)
390  {
391  // trilinear hexahedral shape functions
392  case HEX8:
393  case HEX20:
394  case HEX27:
395  {
396  libmesh_assert_less (i, 8);
397 
398  // Compute hex shape functions as a tensor-product
399  const Real xi = p(0);
400  const Real eta = p(1);
401  const Real zeta = p(2);
402 
403  static const unsigned int i0[] = {0, 1, 1, 0, 0, 1, 1, 0};
404  static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 1, 1};
405  static const unsigned int i2[] = {0, 0, 0, 0, 1, 1, 1, 1};
406 
407  switch(j)
408  {
409  case 0:
410  return (FE<1,L2_LAGRANGE>::shape_deriv(EDGE2, FIRST, i0[i], 0, xi)*
411  FE<1,L2_LAGRANGE>::shape (EDGE2, FIRST, i1[i], eta)*
412  FE<1,L2_LAGRANGE>::shape (EDGE2, FIRST, i2[i], zeta));
413 
414  case 1:
415  return (FE<1,L2_LAGRANGE>::shape (EDGE2, FIRST, i0[i], xi)*
416  FE<1,L2_LAGRANGE>::shape_deriv(EDGE2, FIRST, i1[i], 0, eta)*
417  FE<1,L2_LAGRANGE>::shape (EDGE2, FIRST, i2[i], zeta));
418 
419  case 2:
420  return (FE<1,L2_LAGRANGE>::shape (EDGE2, FIRST, i0[i], xi)*
421  FE<1,L2_LAGRANGE>::shape (EDGE2, FIRST, i1[i], eta)*
422  FE<1,L2_LAGRANGE>::shape_deriv(EDGE2, FIRST, i2[i], 0, zeta));
423 
424  default:
425  libmesh_error_msg("Invalid j = " << j);
426  }
427  }
428 
429  // linear tetrahedral shape functions
430  case TET4:
431  case TET10:
432  {
433  libmesh_assert_less (i, 4);
434 
435  // Area coordinates, pg. 205, Vol. I, Carey, Oden, Becker FEM
436  const Real dzeta0dxi = -1.;
437  const Real dzeta1dxi = 1.;
438  const Real dzeta2dxi = 0.;
439  const Real dzeta3dxi = 0.;
440 
441  const Real dzeta0deta = -1.;
442  const Real dzeta1deta = 0.;
443  const Real dzeta2deta = 1.;
444  const Real dzeta3deta = 0.;
445 
446  const Real dzeta0dzeta = -1.;
447  const Real dzeta1dzeta = 0.;
448  const Real dzeta2dzeta = 0.;
449  const Real dzeta3dzeta = 1.;
450 
451  switch (j)
452  {
453  // d()/dxi
454  case 0:
455  {
456  switch(i)
457  {
458  case 0:
459  return dzeta0dxi;
460 
461  case 1:
462  return dzeta1dxi;
463 
464  case 2:
465  return dzeta2dxi;
466 
467  case 3:
468  return dzeta3dxi;
469 
470  default:
471  libmesh_error_msg("Invalid i = " << i);
472  }
473  }
474 
475  // d()/deta
476  case 1:
477  {
478  switch(i)
479  {
480  case 0:
481  return dzeta0deta;
482 
483  case 1:
484  return dzeta1deta;
485 
486  case 2:
487  return dzeta2deta;
488 
489  case 3:
490  return dzeta3deta;
491 
492  default:
493  libmesh_error_msg("Invalid i = " << i);
494  }
495  }
496 
497  // d()/dzeta
498  case 2:
499  {
500  switch(i)
501  {
502  case 0:
503  return dzeta0dzeta;
504 
505  case 1:
506  return dzeta1dzeta;
507 
508  case 2:
509  return dzeta2dzeta;
510 
511  case 3:
512  return dzeta3dzeta;
513 
514  default:
515  libmesh_error_msg("Invalid i = " << i);
516  }
517  }
518 
519  default:
520  libmesh_error_msg("Invalid shape function derivative j = " << j);
521  }
522  }
523 
524  // linear prism shape functions
525  case PRISM6:
526  case PRISM15:
527  case PRISM18:
528  {
529  libmesh_assert_less (i, 6);
530 
531  // Compute prism shape functions as a tensor-product
532  // of a triangle and an edge
533 
534  Point p2d(p(0),p(1));
535  Point p1d(p(2));
536 
537  // 0 1 2 3 4 5
538  static const unsigned int i0[] = {0, 0, 0, 1, 1, 1};
539  static const unsigned int i1[] = {0, 1, 2, 0, 1, 2};
540 
541  switch (j)
542  {
543  // d()/dxi
544  case 0:
545  return (FE<2,L2_LAGRANGE>::shape_deriv(TRI3, FIRST, i1[i], 0, p2d)*
546  FE<1,L2_LAGRANGE>::shape(EDGE2, FIRST, i0[i], p1d));
547 
548  // d()/deta
549  case 1:
550  return (FE<2,L2_LAGRANGE>::shape_deriv(TRI3, FIRST, i1[i], 1, p2d)*
551  FE<1,L2_LAGRANGE>::shape(EDGE2, FIRST, i0[i], p1d));
552 
553  // d()/dzeta
554  case 2:
555  return (FE<2,L2_LAGRANGE>::shape(TRI3, FIRST, i1[i], p2d)*
556  FE<1,L2_LAGRANGE>::shape_deriv(EDGE2, FIRST, i0[i], 0, p1d));
557 
558  default:
559  libmesh_error_msg("Invalid shape function derivative j = " << j);
560  }
561  }
562 
563  // linear pyramid shape functions
564  case PYRAMID5:
565  {
566  libmesh_assert_less (i, 5);
567 
568  const Real xi = p(0);
569  const Real eta = p(1);
570  const Real zeta = p(2);
571  const Real eps = 1.e-35;
572 
573  switch (j)
574  {
575  // d/dxi
576  case 0:
577  switch(i)
578  {
579  case 0:
580  return .25*(zeta + eta - 1.)/((1. - zeta) + eps);
581 
582  case 1:
583  return -.25*(zeta + eta - 1.)/((1. - zeta) + eps);
584 
585  case 2:
586  return -.25*(zeta - eta - 1.)/((1. - zeta) + eps);
587 
588  case 3:
589  return .25*(zeta - eta - 1.)/((1. - zeta) + eps);
590 
591  case 4:
592  return 0;
593 
594  default:
595  libmesh_error_msg("Invalid i = " << i);
596  }
597 
598 
599  // d/deta
600  case 1:
601  switch(i)
602  {
603  case 0:
604  return .25*(zeta + xi - 1.)/((1. - zeta) + eps);
605 
606  case 1:
607  return .25*(zeta - xi - 1.)/((1. - zeta) + eps);
608 
609  case 2:
610  return -.25*(zeta - xi - 1.)/((1. - zeta) + eps);
611 
612  case 3:
613  return -.25*(zeta + xi - 1.)/((1. - zeta) + eps);
614 
615  case 4:
616  return 0;
617 
618  default:
619  libmesh_error_msg("Invalid i = " << i);
620  }
621 
622 
623  // d/dzeta
624  case 2:
625  switch(i)
626  {
627  case 0:
628  {
629  const Real a=1.;
630  const Real b=1.;
631 
632  return .25*(((zeta + a*xi - 1.)*(zeta + b*eta - 1.) +
633  (1. - zeta)*((zeta + a*xi -1.) + (zeta + b*eta - 1.)))/
634  ((1. - zeta)*(1. - zeta) + eps));
635  }
636 
637  case 1:
638  {
639  const Real a=-1.;
640  const Real b=1.;
641 
642  return .25*(((zeta + a*xi - 1.)*(zeta + b*eta - 1.) +
643  (1. - zeta)*((zeta + a*xi -1.) + (zeta + b*eta - 1.)))/
644  ((1. - zeta)*(1. - zeta) + eps));
645  }
646 
647  case 2:
648  {
649  const Real a=-1.;
650  const Real b=-1.;
651 
652  return .25*(((zeta + a*xi - 1.)*(zeta + b*eta - 1.) +
653  (1. - zeta)*((zeta + a*xi -1.) + (zeta + b*eta - 1.)))/
654  ((1. - zeta)*(1. - zeta) + eps));
655  }
656 
657  case 3:
658  {
659  const Real a=1.;
660  const Real b=-1.;
661 
662  return .25*(((zeta + a*xi - 1.)*(zeta + b*eta - 1.) +
663  (1. - zeta)*((zeta + a*xi -1.) + (zeta + b*eta - 1.)))/
664  ((1. - zeta)*(1. - zeta) + eps));
665  }
666 
667  case 4:
668  return 1.;
669 
670  default:
671  libmesh_error_msg("Invalid i = " << i);
672  }
673 
674 
675  default:
676  libmesh_error_msg("Invalid j = " << j);
677  }
678  }
679 
680 
681  default:
682  libmesh_error_msg("ERROR: Unsupported 3D element type!: " << type);
683  }
684  }
685 
686 
687  // quadratic Lagrange shape functions
688  case SECOND:
689  {
690  switch (type)
691  {
692 
693  // serendipity hexahedral quadratic shape functions
694  case HEX20:
695  {
696  libmesh_assert_less (i, 20);
697 
698  const Real xi = p(0);
699  const Real eta = p(1);
700  const Real zeta = p(2);
701 
702  // these functions are defined for (x,y,z) in [0,1]^3
703  // so transform the locations
704  const Real x = .5*(xi + 1.);
705  const Real y = .5*(eta + 1.);
706  const Real z = .5*(zeta + 1.);
707 
708  // and don't forget the chain rule!
709 
710  switch (j)
711  {
712 
713  // d/dx*dx/dxi
714  case 0:
715  switch (i)
716  {
717  case 0:
718  return .5*(1. - y)*(1. - z)*((1. - x)*(-2.) +
719  (-1.)*(1. - 2.*x - 2.*y - 2.*z));
720 
721  case 1:
722  return .5*(1. - y)*(1. - z)*(x*(2.) +
723  (1.)*(2.*x - 2.*y - 2.*z - 1.));
724 
725  case 2:
726  return .5*y*(1. - z)*(x*(2.) +
727  (1.)*(2.*x + 2.*y - 2.*z - 3.));
728 
729  case 3:
730  return .5*y*(1. - z)*((1. - x)*(-2.) +
731  (-1.)*(2.*y - 2.*x - 2.*z - 1.));
732 
733  case 4:
734  return .5*(1. - y)*z*((1. - x)*(-2.) +
735  (-1.)*(2.*z - 2.*x - 2.*y - 1.));
736 
737  case 5:
738  return .5*(1. - y)*z*(x*(2.) +
739  (1.)*(2.*x - 2.*y + 2.*z - 3.));
740 
741  case 6:
742  return .5*y*z*(x*(2.) +
743  (1.)*(2.*x + 2.*y + 2.*z - 5.));
744 
745  case 7:
746  return .5*y*z*((1. - x)*(-2.) +
747  (-1.)*(2.*y - 2.*x + 2.*z - 3.));
748 
749  case 8:
750  return 2.*(1. - y)*(1. - z)*(1. - 2.*x);
751 
752  case 9:
753  return 2.*y*(1. - y)*(1. - z);
754 
755  case 10:
756  return 2.*y*(1. - z)*(1. - 2.*x);
757 
758  case 11:
759  return 2.*y*(1. - y)*(1. - z)*(-1.);
760 
761  case 12:
762  return 2.*(1. - y)*z*(1. - z)*(-1.);
763 
764  case 13:
765  return 2.*(1. - y)*z*(1. - z);
766 
767  case 14:
768  return 2.*y*z*(1. - z);
769 
770  case 15:
771  return 2.*y*z*(1. - z)*(-1.);
772 
773  case 16:
774  return 2.*(1. - y)*z*(1. - 2.*x);
775 
776  case 17:
777  return 2.*y*(1. - y)*z;
778 
779  case 18:
780  return 2.*y*z*(1. - 2.*x);
781 
782  case 19:
783  return 2.*y*(1. - y)*z*(-1.);
784 
785  default:
786  libmesh_error_msg("Invalid i = " << i);
787  }
788 
789 
790  // d/dy*dy/deta
791  case 1:
792  switch (i)
793  {
794  case 0:
795  return .5*(1. - x)*(1. - z)*((1. - y)*(-2.) +
796  (-1.)*(1. - 2.*x - 2.*y - 2.*z));
797 
798  case 1:
799  return .5*x*(1. - z)*((1. - y)*(-2.) +
800  (-1.)*(2.*x - 2.*y - 2.*z - 1.));
801 
802  case 2:
803  return .5*x*(1. - z)*(y*(2.) +
804  (1.)*(2.*x + 2.*y - 2.*z - 3.));
805 
806  case 3:
807  return .5*(1. - x)*(1. - z)*(y*(2.) +
808  (1.)*(2.*y - 2.*x - 2.*z - 1.));
809 
810  case 4:
811  return .5*(1. - x)*z*((1. - y)*(-2.) +
812  (-1.)*(2.*z - 2.*x - 2.*y - 1.));
813 
814  case 5:
815  return .5*x*z*((1. - y)*(-2.) +
816  (-1.)*(2.*x - 2.*y + 2.*z - 3.));
817 
818  case 6:
819  return .5*x*z*(y*(2.) +
820  (1.)*(2.*x + 2.*y + 2.*z - 5.));
821 
822  case 7:
823  return .5*(1. - x)*z*(y*(2.) +
824  (1.)*(2.*y - 2.*x + 2.*z - 3.));
825 
826  case 8:
827  return 2.*x*(1. - x)*(1. - z)*(-1.);
828 
829  case 9:
830  return 2.*x*(1. - z)*(1. - 2.*y);
831 
832  case 10:
833  return 2.*x*(1. - x)*(1. - z);
834 
835  case 11:
836  return 2.*(1. - x)*(1. - z)*(1. - 2.*y);
837 
838  case 12:
839  return 2.*(1. - x)*z*(1. - z)*(-1.);
840 
841  case 13:
842  return 2.*x*z*(1. - z)*(-1.);
843 
844  case 14:
845  return 2.*x*z*(1. - z);
846 
847  case 15:
848  return 2.*(1. - x)*z*(1. - z);
849 
850  case 16:
851  return 2.*x*(1. - x)*z*(-1.);
852 
853  case 17:
854  return 2.*x*z*(1. - 2.*y);
855 
856  case 18:
857  return 2.*x*(1. - x)*z;
858 
859  case 19:
860  return 2.*(1. - x)*z*(1. - 2.*y);
861 
862  default:
863  libmesh_error_msg("Invalid i = " << i);
864  }
865 
866 
867  // d/dz*dz/dzeta
868  case 2:
869  switch (i)
870  {
871  case 0:
872  return .5*(1. - x)*(1. - y)*((1. - z)*(-2.) +
873  (-1.)*(1. - 2.*x - 2.*y - 2.*z));
874 
875  case 1:
876  return .5*x*(1. - y)*((1. - z)*(-2.) +
877  (-1.)*(2.*x - 2.*y - 2.*z - 1.));
878 
879  case 2:
880  return .5*x*y*((1. - z)*(-2.) +
881  (-1.)*(2.*x + 2.*y - 2.*z - 3.));
882 
883  case 3:
884  return .5*(1. - x)*y*((1. - z)*(-2.) +
885  (-1.)*(2.*y - 2.*x - 2.*z - 1.));
886 
887  case 4:
888  return .5*(1. - x)*(1. - y)*(z*(2.) +
889  (1.)*(2.*z - 2.*x - 2.*y - 1.));
890 
891  case 5:
892  return .5*x*(1. - y)*(z*(2.) +
893  (1.)*(2.*x - 2.*y + 2.*z - 3.));
894 
895  case 6:
896  return .5*x*y*(z*(2.) +
897  (1.)*(2.*x + 2.*y + 2.*z - 5.));
898 
899  case 7:
900  return .5*(1. - x)*y*(z*(2.) +
901  (1.)*(2.*y - 2.*x + 2.*z - 3.));
902 
903  case 8:
904  return 2.*x*(1. - x)*(1. - y)*(-1.);
905 
906  case 9:
907  return 2.*x*y*(1. - y)*(-1.);
908 
909  case 10:
910  return 2.*x*(1. - x)*y*(-1.);
911 
912  case 11:
913  return 2.*(1. - x)*y*(1. - y)*(-1.);
914 
915  case 12:
916  return 2.*(1. - x)*(1. - y)*(1. - 2.*z);
917 
918  case 13:
919  return 2.*x*(1. - y)*(1. - 2.*z);
920 
921  case 14:
922  return 2.*x*y*(1. - 2.*z);
923 
924  case 15:
925  return 2.*(1. - x)*y*(1. - 2.*z);
926 
927  case 16:
928  return 2.*x*(1. - x)*(1. - y);
929 
930  case 17:
931  return 2.*x*y*(1. - y);
932 
933  case 18:
934  return 2.*x*(1. - x)*y;
935 
936  case 19:
937  return 2.*(1. - x)*y*(1. - y);
938 
939  default:
940  libmesh_error_msg("Invalid i = " << i);
941  }
942 
943  default:
944  libmesh_error_msg("Invalid shape function derivative j = " << j);
945  }
946  }
947 
948  // triquadratic hexahedral shape functions
949  case HEX27:
950  {
951  libmesh_assert_less (i, 27);
952 
953  // Compute hex shape functions as a tensor-product
954  const Real xi = p(0);
955  const Real eta = p(1);
956  const Real zeta = p(2);
957 
958  // The only way to make any sense of this
959  // is to look at the mgflo/mg2/mgf documentation
960  // and make the cut-out cube!
961  // 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
962  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};
963  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};
964  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};
965 
966  switch(j)
967  {
968  case 0:
969  return (FE<1,L2_LAGRANGE>::shape_deriv(EDGE3, SECOND, i0[i], 0, xi)*
970  FE<1,L2_LAGRANGE>::shape (EDGE3, SECOND, i1[i], eta)*
971  FE<1,L2_LAGRANGE>::shape (EDGE3, SECOND, i2[i], zeta));
972 
973  case 1:
974  return (FE<1,L2_LAGRANGE>::shape (EDGE3, SECOND, i0[i], xi)*
976  FE<1,L2_LAGRANGE>::shape (EDGE3, SECOND, i2[i], zeta));
977 
978  case 2:
979  return (FE<1,L2_LAGRANGE>::shape (EDGE3, SECOND, i0[i], xi)*
980  FE<1,L2_LAGRANGE>::shape (EDGE3, SECOND, i1[i], eta)*
981  FE<1,L2_LAGRANGE>::shape_deriv(EDGE3, SECOND, i2[i], 0, zeta));
982 
983  default:
984  libmesh_error_msg("Invalid j = " << j);
985  }
986  }
987 
988  // quadratic tetrahedral shape functions
989  case TET10:
990  {
991  libmesh_assert_less (i, 10);
992 
993  // Area coordinates, pg. 205, Vol. I, Carey, Oden, Becker FEM
994  const Real zeta1 = p(0);
995  const Real zeta2 = p(1);
996  const Real zeta3 = p(2);
997  const Real zeta0 = 1. - zeta1 - zeta2 - zeta3;
998 
999  const Real dzeta0dxi = -1.;
1000  const Real dzeta1dxi = 1.;
1001  const Real dzeta2dxi = 0.;
1002  const Real dzeta3dxi = 0.;
1003 
1004  const Real dzeta0deta = -1.;
1005  const Real dzeta1deta = 0.;
1006  const Real dzeta2deta = 1.;
1007  const Real dzeta3deta = 0.;
1008 
1009  const Real dzeta0dzeta = -1.;
1010  const Real dzeta1dzeta = 0.;
1011  const Real dzeta2dzeta = 0.;
1012  const Real dzeta3dzeta = 1.;
1013 
1014  switch (j)
1015  {
1016  // d()/dxi
1017  case 0:
1018  {
1019  switch(i)
1020  {
1021  case 0:
1022  return (4.*zeta0 - 1.)*dzeta0dxi;
1023 
1024  case 1:
1025  return (4.*zeta1 - 1.)*dzeta1dxi;
1026 
1027  case 2:
1028  return (4.*zeta2 - 1.)*dzeta2dxi;
1029 
1030  case 3:
1031  return (4.*zeta3 - 1.)*dzeta3dxi;
1032 
1033  case 4:
1034  return 4.*(zeta0*dzeta1dxi + dzeta0dxi*zeta1);
1035 
1036  case 5:
1037  return 4.*(zeta1*dzeta2dxi + dzeta1dxi*zeta2);
1038 
1039  case 6:
1040  return 4.*(zeta0*dzeta2dxi + dzeta0dxi*zeta2);
1041 
1042  case 7:
1043  return 4.*(zeta0*dzeta3dxi + dzeta0dxi*zeta3);
1044 
1045  case 8:
1046  return 4.*(zeta1*dzeta3dxi + dzeta1dxi*zeta3);
1047 
1048  case 9:
1049  return 4.*(zeta2*dzeta3dxi + dzeta2dxi*zeta3);
1050 
1051  default:
1052  libmesh_error_msg("Invalid i = " << i);
1053  }
1054  }
1055 
1056  // d()/deta
1057  case 1:
1058  {
1059  switch(i)
1060  {
1061  case 0:
1062  return (4.*zeta0 - 1.)*dzeta0deta;
1063 
1064  case 1:
1065  return (4.*zeta1 - 1.)*dzeta1deta;
1066 
1067  case 2:
1068  return (4.*zeta2 - 1.)*dzeta2deta;
1069 
1070  case 3:
1071  return (4.*zeta3 - 1.)*dzeta3deta;
1072 
1073  case 4:
1074  return 4.*(zeta0*dzeta1deta + dzeta0deta*zeta1);
1075 
1076  case 5:
1077  return 4.*(zeta1*dzeta2deta + dzeta1deta*zeta2);
1078 
1079  case 6:
1080  return 4.*(zeta0*dzeta2deta + dzeta0deta*zeta2);
1081 
1082  case 7:
1083  return 4.*(zeta0*dzeta3deta + dzeta0deta*zeta3);
1084 
1085  case 8:
1086  return 4.*(zeta1*dzeta3deta + dzeta1deta*zeta3);
1087 
1088  case 9:
1089  return 4.*(zeta2*dzeta3deta + dzeta2deta*zeta3);
1090 
1091  default:
1092  libmesh_error_msg("Invalid i = " << i);
1093  }
1094  }
1095 
1096  // d()/dzeta
1097  case 2:
1098  {
1099  switch(i)
1100  {
1101  case 0:
1102  return (4.*zeta0 - 1.)*dzeta0dzeta;
1103 
1104  case 1:
1105  return (4.*zeta1 - 1.)*dzeta1dzeta;
1106 
1107  case 2:
1108  return (4.*zeta2 - 1.)*dzeta2dzeta;
1109 
1110  case 3:
1111  return (4.*zeta3 - 1.)*dzeta3dzeta;
1112 
1113  case 4:
1114  return 4.*(zeta0*dzeta1dzeta + dzeta0dzeta*zeta1);
1115 
1116  case 5:
1117  return 4.*(zeta1*dzeta2dzeta + dzeta1dzeta*zeta2);
1118 
1119  case 6:
1120  return 4.*(zeta0*dzeta2dzeta + dzeta0dzeta*zeta2);
1121 
1122  case 7:
1123  return 4.*(zeta0*dzeta3dzeta + dzeta0dzeta*zeta3);
1124 
1125  case 8:
1126  return 4.*(zeta1*dzeta3dzeta + dzeta1dzeta*zeta3);
1127 
1128  case 9:
1129  return 4.*(zeta2*dzeta3dzeta + dzeta2dzeta*zeta3);
1130 
1131  default:
1132  libmesh_error_msg("Invalid i = " << i);
1133  }
1134  }
1135 
1136  default:
1137  libmesh_error_msg("Invalid j = " << j);
1138  }
1139  }
1140 
1141 
1142 
1143  // quadratic prism shape functions
1144  case PRISM18:
1145  {
1146  libmesh_assert_less (i, 18);
1147 
1148  // Compute prism shape functions as a tensor-product
1149  // of a triangle and an edge
1150 
1151  Point p2d(p(0),p(1));
1152  Point p1d(p(2));
1153 
1154  // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
1155  static const unsigned int i0[] = {0, 0, 0, 1, 1, 1, 0, 0, 0, 2, 2, 2, 1, 1, 1, 2, 2, 2};
1156  static const unsigned int i1[] = {0, 1, 2, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 3, 4, 5};
1157 
1158  switch (j)
1159  {
1160  // d()/dxi
1161  case 0:
1162  return (FE<2,L2_LAGRANGE>::shape_deriv(TRI6, SECOND, i1[i], 0, p2d)*
1163  FE<1,L2_LAGRANGE>::shape(EDGE3, SECOND, i0[i], p1d));
1164 
1165  // d()/deta
1166  case 1:
1167  return (FE<2,L2_LAGRANGE>::shape_deriv(TRI6, SECOND, i1[i], 1, p2d)*
1168  FE<1,L2_LAGRANGE>::shape(EDGE3, SECOND, i0[i], p1d));
1169 
1170  // d()/dzeta
1171  case 2:
1172  return (FE<2,L2_LAGRANGE>::shape(TRI6, SECOND, i1[i], p2d)*
1173  FE<1,L2_LAGRANGE>::shape_deriv(EDGE3, SECOND, i0[i], 0, p1d));
1174 
1175  default:
1176  libmesh_error_msg("Invalid shape function derivative j = " << j);
1177  }
1178  }
1179 
1180 
1181  default:
1182  libmesh_error_msg("ERROR: Unsupported 3D element type!: " << type);
1183  }
1184  }
1185 
1186 
1187  // unsupported order
1188  default:
1189  libmesh_error_msg("ERROR: Unsupported 3D FE order!: " << order);
1190  }
1191 
1192 #endif
1193 
1194  libmesh_error_msg("We'll never get here!");
1195  return 0.;
1196 }
1197 
1198 
1199 
1200 template <>
1202  const Order order,
1203  const unsigned int i,
1204  const unsigned int j,
1205  const Point & p)
1206 {
1207  libmesh_assert(elem);
1208 
1209  // call the orientation-independent shape function derivatives
1210  return FE<3,L2_LAGRANGE>::shape_deriv(elem->type(), static_cast<Order>(order + elem->p_level()), i, j, p);
1211 }
1212 
1213 
1214 
1215 template <>
1217  const Order order,
1218  const unsigned int i,
1219  const unsigned int j,
1220  const Point & p)
1221 {
1222 #if LIBMESH_DIM == 3
1223 
1224  libmesh_assert_less (j, 6);
1225 
1226  switch (order)
1227  {
1228  // linear Lagrange shape functions
1229  case FIRST:
1230  {
1231  return 0.;
1232  }
1233 
1234  // quadratic Lagrange shape functions
1235  case SECOND:
1236  {
1237  switch (type)
1238  {
1239 
1240  // serendipity hexahedral quadratic shape functions
1241  case HEX20:
1242  {
1243  static bool warning_given_HEX20 = false;
1244 
1245  if (!warning_given_HEX20)
1246  libMesh::err << "Second derivatives for 3D Lagrangian HEX20"
1247  << " elements are not yet implemented!"
1248  << std::endl;
1249  warning_given_HEX20 = true;
1250  }
1251  libmesh_fallthrough();
1252 
1253  case HEX27:
1254  // triquadratic hexahedral shape functions
1255  {
1256  libmesh_assert_less (i, 27);
1257 
1258  // Compute hex shape functions as a tensor-product
1259  const Real xi = p(0);
1260  const Real eta = p(1);
1261  const Real zeta = p(2);
1262 
1263  // The only way to make any sense of this
1264  // is to look at the mgflo/mg2/mgf documentation
1265  // and make the cut-out cube!
1266  // 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
1267  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};
1268  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};
1269  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};
1270 
1271  switch(j)
1272  {
1273  // d^2()/dxi^2
1274  case 0:
1275  return (FE<1,L2_LAGRANGE>::shape_second_deriv(EDGE3, SECOND, i0[i], 0, xi)*
1276  FE<1,L2_LAGRANGE>::shape (EDGE3, SECOND, i1[i], eta)*
1277  FE<1,L2_LAGRANGE>::shape (EDGE3, SECOND, i2[i], zeta));
1278 
1279  // d^2()/dxideta
1280  case 1:
1281  return (FE<1,L2_LAGRANGE>::shape_deriv(EDGE3, SECOND, i0[i], 0, xi)*
1282  FE<1,L2_LAGRANGE>::shape_deriv(EDGE3, SECOND, i1[i], 0, eta)*
1283  FE<1,L2_LAGRANGE>::shape (EDGE3, SECOND, i2[i], zeta));
1284 
1285  // d^2()/deta^2
1286  case 2:
1287  return (FE<1,L2_LAGRANGE>::shape (EDGE3, SECOND, i0[i], xi)*
1289  FE<1,L2_LAGRANGE>::shape (EDGE3, SECOND, i2[i], zeta));
1290 
1291  // d^2()/dxidzeta
1292  case 3:
1293  return (FE<1,L2_LAGRANGE>::shape_deriv(EDGE3, SECOND, i0[i], 0, xi)*
1294  FE<1,L2_LAGRANGE>::shape (EDGE3, SECOND, i1[i], eta)*
1295  FE<1,L2_LAGRANGE>::shape_deriv(EDGE3, SECOND, i2[i], 0, zeta));
1296 
1297  // d^2()/detadzeta
1298  case 4:
1299  return (FE<1,L2_LAGRANGE>::shape (EDGE3, SECOND, i0[i], xi)*
1300  FE<1,L2_LAGRANGE>::shape_deriv(EDGE3, SECOND, i1[i], 0, eta)*
1301  FE<1,L2_LAGRANGE>::shape_deriv(EDGE3, SECOND, i2[i], 0, zeta));
1302 
1303  // d^2()/dzeta^2
1304  case 5:
1305  return (FE<1,L2_LAGRANGE>::shape (EDGE3, SECOND, i0[i], xi)*
1306  FE<1,L2_LAGRANGE>::shape (EDGE3, SECOND, i1[i], eta)*
1308 
1309  default:
1310  libmesh_error_msg("Invalid j = " << j);
1311  }
1312  }
1313 
1314  // quadratic tetrahedral shape functions
1315  case TET10:
1316  {
1317  // The area coordinates are the same as used for the
1318  // shape() and shape_deriv() functions.
1319  // const Real zeta0 = 1. - zeta1 - zeta2 - zeta3;
1320  // const Real zeta1 = p(0);
1321  // const Real zeta2 = p(1);
1322  // const Real zeta3 = p(2);
1323  static const Real dzetadxi[4][3] =
1324  {
1325  {-1., -1., -1.},
1326  {1., 0., 0.},
1327  {0., 1., 0.},
1328  {0., 0., 1.}
1329  };
1330 
1331  // Convert from j -> (j,k) indices for independent variable
1332  // (0=xi, 1=eta, 2=zeta)
1333  static const unsigned short int independent_var_indices[6][2] =
1334  {
1335  {0, 0}, // d^2 phi / dxi^2
1336  {0, 1}, // d^2 phi / dxi deta
1337  {1, 1}, // d^2 phi / deta^2
1338  {0, 2}, // d^2 phi / dxi dzeta
1339  {1, 2}, // d^2 phi / deta dzeta
1340  {2, 2} // d^2 phi / dzeta^2
1341  };
1342 
1343  // Convert from i -> zeta indices. Each quadratic shape
1344  // function for the Tet10 depends on up to two of the zeta
1345  // area coordinate functions (see the shape() function above).
1346  // This table just tells which two area coords it uses.
1347  static const unsigned short int zeta_indices[10][2] =
1348  {
1349  {0, 0},
1350  {1, 1},
1351  {2, 2},
1352  {3, 3},
1353  {0, 1},
1354  {1, 2},
1355  {2, 0},
1356  {0, 3},
1357  {1, 3},
1358  {2, 3},
1359  };
1360 
1361  // Look up the independent variable indices for this value of j.
1362  const unsigned int my_j = independent_var_indices[j][0];
1363  const unsigned int my_k = independent_var_indices[j][1];
1364 
1365  if (i<4)
1366  {
1367  return 4.*dzetadxi[i][my_j]*dzetadxi[i][my_k];
1368  }
1369 
1370  else if (i<10)
1371  {
1372  const unsigned short int my_m = zeta_indices[i][0];
1373  const unsigned short int my_n = zeta_indices[i][1];
1374 
1375  return 4.*(dzetadxi[my_n][my_j]*dzetadxi[my_m][my_k] +
1376  dzetadxi[my_m][my_j]*dzetadxi[my_n][my_k] );
1377  }
1378  else
1379  libmesh_error_msg("Invalid shape function index " << i);
1380  }
1381 
1382 
1383  // quadratic prism shape functions
1384  case PRISM18:
1385  {
1386  libmesh_assert_less (i, 18);
1387 
1388  // Compute prism shape functions as a tensor-product
1389  // of a triangle and an edge
1390 
1391  Point p2d(p(0),p(1));
1392  Point p1d(p(2));
1393 
1394  // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
1395  static const unsigned int i0[] = {0, 0, 0, 1, 1, 1, 0, 0, 0, 2, 2, 2, 1, 1, 1, 2, 2, 2};
1396  static const unsigned int i1[] = {0, 1, 2, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 3, 4, 5};
1397 
1398  switch (j)
1399  {
1400  // d^2()/dxi^2
1401  case 0:
1402  return (FE<2,L2_LAGRANGE>::shape_second_deriv(TRI6, SECOND, i1[i], 0, p2d)*
1403  FE<1,L2_LAGRANGE>::shape(EDGE3, SECOND, i0[i], p1d));
1404 
1405  // d^2()/dxideta
1406  case 1:
1407  return (FE<2,L2_LAGRANGE>::shape_second_deriv(TRI6, SECOND, i1[i], 1, p2d)*
1408  FE<1,L2_LAGRANGE>::shape(EDGE3, SECOND, i0[i], p1d));
1409 
1410  // d^2()/deta^2
1411  case 2:
1412  return (FE<2,L2_LAGRANGE>::shape_second_deriv(TRI6, SECOND, i1[i], 2, p2d)*
1413  FE<1,L2_LAGRANGE>::shape(EDGE3, SECOND, i0[i], p1d));
1414 
1415  // d^2()/dxidzeta
1416  case 3:
1417  return (FE<2,L2_LAGRANGE>::shape_deriv(TRI6, SECOND, i1[i], 0, p2d)*
1418  FE<1,L2_LAGRANGE>::shape_deriv(EDGE3, SECOND, i0[i], 0, p1d));
1419 
1420  // d^2()/detadzeta
1421  case 4:
1422  return (FE<2,L2_LAGRANGE>::shape_deriv(TRI6, SECOND, i1[i], 1, p2d)*
1423  FE<1,L2_LAGRANGE>::shape_deriv(EDGE3, SECOND, i0[i], 0, p1d));
1424 
1425  // d^2()/dzeta^2
1426  case 5:
1427  return (FE<2,L2_LAGRANGE>::shape(TRI6, SECOND, i1[i], p2d)*
1429 
1430  default:
1431  libmesh_error_msg("Invalid shape function derivative j = " << j);
1432  }
1433  }
1434 
1435 
1436 
1437  default:
1438  libmesh_error_msg("ERROR: Unsupported 3D element type!: " << type);
1439  }
1440  }
1441 
1442 
1443  // unsupported order
1444  default:
1445  libmesh_error_msg("ERROR: Unsupported 3D FE order!: " << order);
1446  }
1447 
1448 #endif
1449 
1450  libmesh_error_msg("We'll never get here!");
1451  return 0.;
1452 }
1453 
1454 
1455 
1456 template <>
1458  const Order order,
1459  const unsigned int i,
1460  const unsigned int j,
1461  const Point & p)
1462 {
1463  libmesh_assert(elem);
1464 
1465  // call the orientation-independent shape function derivatives
1466  return FE<3,L2_LAGRANGE>::shape_second_deriv(elem->type(), static_cast<Order>(order + elem->p_level()), i, j, p);
1467 }
1468 
1469 } // namespace libMesh
OStreamProxy err
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)