fe_l2_lagrange_shape_3D.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2018 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 // C++ includes
20 
21 // Local includes
22 #include "libmesh/fe.h"
23 #include "libmesh/elem.h"
24 
25 namespace libMesh
26 {
27 
28 
29 
30 
31 template <>
32 Real FE<3,L2_LAGRANGE>::shape(const ElemType type,
33  const Order order,
34  const unsigned int i,
35  const Point & p)
36 {
37 #if LIBMESH_DIM == 3
38 
39 
40  switch (order)
41  {
42  // linear Lagrange shape functions
43  case FIRST:
44  {
45  switch (type)
46  {
47  // trilinear hexahedral shape functions
48  case HEX8:
49  case HEX20:
50  case HEX27:
51  {
52  libmesh_assert_less (i, 8);
53 
54  // Compute hex shape functions as a tensor-product
55  const Real xi = p(0);
56  const Real eta = p(1);
57  const Real zeta = p(2);
58 
59  // 0 1 2 3 4 5 6 7
60  static const unsigned int i0[] = {0, 1, 1, 0, 0, 1, 1, 0};
61  static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 1, 1};
62  static const unsigned int i2[] = {0, 0, 0, 0, 1, 1, 1, 1};
63 
64  return (FE<1,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 
353 
354 
355 template <>
356 Real FE<3,L2_LAGRANGE>::shape(const Elem * elem,
357  const Order order,
358  const unsigned int i,
359  const Point & p)
360 {
361  libmesh_assert(elem);
362 
363  // call the orientation-independent shape functions
364  return FE<3,L2_LAGRANGE>::shape(elem->type(), static_cast<Order>(order + elem->p_level()), i, p);
365 }
366 
367 
368 
369 
370 template <>
371 Real FE<3,L2_LAGRANGE>::shape_deriv(const ElemType type,
372  const Order order,
373  const unsigned int i,
374  const unsigned int j,
375  const Point & p)
376 {
377 #if LIBMESH_DIM == 3
378 
379  libmesh_assert_less (j, 3);
380 
381  switch (order)
382  {
383  // linear Lagrange shape functions
384  case FIRST:
385  {
386  switch (type)
387  {
388  // trilinear hexahedral shape functions
389  case HEX8:
390  case HEX20:
391  case HEX27:
392  {
393  libmesh_assert_less (i, 8);
394 
395  // Compute hex shape functions as a tensor-product
396  const Real xi = p(0);
397  const Real eta = p(1);
398  const Real zeta = p(2);
399 
400  static const unsigned int i0[] = {0, 1, 1, 0, 0, 1, 1, 0};
401  static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 1, 1};
402  static const unsigned int i2[] = {0, 0, 0, 0, 1, 1, 1, 1};
403 
404  switch(j)
405  {
406  case 0:
407  return (FE<1,L2_LAGRANGE>::shape_deriv(EDGE2, FIRST, i0[i], 0, xi)*
408  FE<1,L2_LAGRANGE>::shape (EDGE2, FIRST, i1[i], eta)*
409  FE<1,L2_LAGRANGE>::shape (EDGE2, FIRST, i2[i], zeta));
410 
411  case 1:
412  return (FE<1,L2_LAGRANGE>::shape (EDGE2, FIRST, i0[i], xi)*
413  FE<1,L2_LAGRANGE>::shape_deriv(EDGE2, FIRST, i1[i], 0, eta)*
414  FE<1,L2_LAGRANGE>::shape (EDGE2, FIRST, i2[i], zeta));
415 
416  case 2:
417  return (FE<1,L2_LAGRANGE>::shape (EDGE2, FIRST, i0[i], xi)*
418  FE<1,L2_LAGRANGE>::shape (EDGE2, FIRST, i1[i], eta)*
419  FE<1,L2_LAGRANGE>::shape_deriv(EDGE2, FIRST, i2[i], 0, zeta));
420 
421  default:
422  libmesh_error_msg("Invalid j = " << j);
423  }
424  }
425 
426  // linear tetrahedral shape functions
427  case TET4:
428  case TET10:
429  {
430  libmesh_assert_less (i, 4);
431 
432  // Area coordinates, pg. 205, Vol. I, Carey, Oden, Becker FEM
433  const Real dzeta0dxi = -1.;
434  const Real dzeta1dxi = 1.;
435  const Real dzeta2dxi = 0.;
436  const Real dzeta3dxi = 0.;
437 
438  const Real dzeta0deta = -1.;
439  const Real dzeta1deta = 0.;
440  const Real dzeta2deta = 1.;
441  const Real dzeta3deta = 0.;
442 
443  const Real dzeta0dzeta = -1.;
444  const Real dzeta1dzeta = 0.;
445  const Real dzeta2dzeta = 0.;
446  const Real dzeta3dzeta = 1.;
447 
448  switch (j)
449  {
450  // d()/dxi
451  case 0:
452  {
453  switch(i)
454  {
455  case 0:
456  return dzeta0dxi;
457 
458  case 1:
459  return dzeta1dxi;
460 
461  case 2:
462  return dzeta2dxi;
463 
464  case 3:
465  return dzeta3dxi;
466 
467  default:
468  libmesh_error_msg("Invalid i = " << i);
469  }
470  }
471 
472  // d()/deta
473  case 1:
474  {
475  switch(i)
476  {
477  case 0:
478  return dzeta0deta;
479 
480  case 1:
481  return dzeta1deta;
482 
483  case 2:
484  return dzeta2deta;
485 
486  case 3:
487  return dzeta3deta;
488 
489  default:
490  libmesh_error_msg("Invalid i = " << i);
491  }
492  }
493 
494  // d()/dzeta
495  case 2:
496  {
497  switch(i)
498  {
499  case 0:
500  return dzeta0dzeta;
501 
502  case 1:
503  return dzeta1dzeta;
504 
505  case 2:
506  return dzeta2dzeta;
507 
508  case 3:
509  return dzeta3dzeta;
510 
511  default:
512  libmesh_error_msg("Invalid i = " << i);
513  }
514  }
515 
516  default:
517  libmesh_error_msg("Invalid shape function derivative j = " << j);
518  }
519  }
520 
521  // linear prism shape functions
522  case PRISM6:
523  case PRISM15:
524  case PRISM18:
525  {
526  libmesh_assert_less (i, 6);
527 
528  // Compute prism shape functions as a tensor-product
529  // of a triangle and an edge
530 
531  Point p2d(p(0),p(1));
532  Point p1d(p(2));
533 
534  // 0 1 2 3 4 5
535  static const unsigned int i0[] = {0, 0, 0, 1, 1, 1};
536  static const unsigned int i1[] = {0, 1, 2, 0, 1, 2};
537 
538  switch (j)
539  {
540  // d()/dxi
541  case 0:
542  return (FE<2,L2_LAGRANGE>::shape_deriv(TRI3, FIRST, i1[i], 0, p2d)*
543  FE<1,L2_LAGRANGE>::shape(EDGE2, FIRST, i0[i], p1d));
544 
545  // d()/deta
546  case 1:
547  return (FE<2,L2_LAGRANGE>::shape_deriv(TRI3, FIRST, i1[i], 1, p2d)*
548  FE<1,L2_LAGRANGE>::shape(EDGE2, FIRST, i0[i], p1d));
549 
550  // d()/dzeta
551  case 2:
552  return (FE<2,L2_LAGRANGE>::shape(TRI3, FIRST, i1[i], p2d)*
553  FE<1,L2_LAGRANGE>::shape_deriv(EDGE2, FIRST, i0[i], 0, p1d));
554 
555  default:
556  libmesh_error_msg("Invalid shape function derivative j = " << j);
557  }
558  }
559 
560  // linear pyramid shape functions
561  case PYRAMID5:
562  {
563  libmesh_assert_less (i, 5);
564 
565  const Real xi = p(0);
566  const Real eta = p(1);
567  const Real zeta = p(2);
568  const Real eps = 1.e-35;
569 
570  switch (j)
571  {
572  // d/dxi
573  case 0:
574  switch(i)
575  {
576  case 0:
577  return .25*(zeta + eta - 1.)/((1. - zeta) + eps);
578 
579  case 1:
580  return -.25*(zeta + eta - 1.)/((1. - zeta) + eps);
581 
582  case 2:
583  return -.25*(zeta - eta - 1.)/((1. - zeta) + eps);
584 
585  case 3:
586  return .25*(zeta - eta - 1.)/((1. - zeta) + eps);
587 
588  case 4:
589  return 0;
590 
591  default:
592  libmesh_error_msg("Invalid i = " << i);
593  }
594 
595 
596  // d/deta
597  case 1:
598  switch(i)
599  {
600  case 0:
601  return .25*(zeta + xi - 1.)/((1. - zeta) + eps);
602 
603  case 1:
604  return .25*(zeta - xi - 1.)/((1. - zeta) + eps);
605 
606  case 2:
607  return -.25*(zeta - xi - 1.)/((1. - zeta) + eps);
608 
609  case 3:
610  return -.25*(zeta + xi - 1.)/((1. - zeta) + eps);
611 
612  case 4:
613  return 0;
614 
615  default:
616  libmesh_error_msg("Invalid i = " << i);
617  }
618 
619 
620  // d/dzeta
621  case 2:
622  switch(i)
623  {
624  case 0:
625  {
626  const Real a=1.;
627  const Real b=1.;
628 
629  return .25*(((zeta + a*xi - 1.)*(zeta + b*eta - 1.) +
630  (1. - zeta)*((zeta + a*xi -1.) + (zeta + b*eta - 1.)))/
631  ((1. - zeta)*(1. - zeta) + eps));
632  }
633 
634  case 1:
635  {
636  const Real a=-1.;
637  const Real b=1.;
638 
639  return .25*(((zeta + a*xi - 1.)*(zeta + b*eta - 1.) +
640  (1. - zeta)*((zeta + a*xi -1.) + (zeta + b*eta - 1.)))/
641  ((1. - zeta)*(1. - zeta) + eps));
642  }
643 
644  case 2:
645  {
646  const Real a=-1.;
647  const Real b=-1.;
648 
649  return .25*(((zeta + a*xi - 1.)*(zeta + b*eta - 1.) +
650  (1. - zeta)*((zeta + a*xi -1.) + (zeta + b*eta - 1.)))/
651  ((1. - zeta)*(1. - zeta) + eps));
652  }
653 
654  case 3:
655  {
656  const Real a=1.;
657  const Real b=-1.;
658 
659  return .25*(((zeta + a*xi - 1.)*(zeta + b*eta - 1.) +
660  (1. - zeta)*((zeta + a*xi -1.) + (zeta + b*eta - 1.)))/
661  ((1. - zeta)*(1. - zeta) + eps));
662  }
663 
664  case 4:
665  return 1.;
666 
667  default:
668  libmesh_error_msg("Invalid i = " << i);
669  }
670 
671 
672  default:
673  libmesh_error_msg("Invalid j = " << j);
674  }
675  }
676 
677 
678  default:
679  libmesh_error_msg("ERROR: Unsupported 3D element type!: " << type);
680  }
681  }
682 
683 
684  // quadratic Lagrange shape functions
685  case SECOND:
686  {
687  switch (type)
688  {
689 
690  // serendipity hexahedral quadratic shape functions
691  case HEX20:
692  {
693  libmesh_assert_less (i, 20);
694 
695  const Real xi = p(0);
696  const Real eta = p(1);
697  const Real zeta = p(2);
698 
699  // these functions are defined for (x,y,z) in [0,1]^3
700  // so transform the locations
701  const Real x = .5*(xi + 1.);
702  const Real y = .5*(eta + 1.);
703  const Real z = .5*(zeta + 1.);
704 
705  // and don't forget the chain rule!
706 
707  switch (j)
708  {
709 
710  // d/dx*dx/dxi
711  case 0:
712  switch (i)
713  {
714  case 0:
715  return .5*(1. - y)*(1. - z)*((1. - x)*(-2.) +
716  (-1.)*(1. - 2.*x - 2.*y - 2.*z));
717 
718  case 1:
719  return .5*(1. - y)*(1. - z)*(x*(2.) +
720  (1.)*(2.*x - 2.*y - 2.*z - 1.));
721 
722  case 2:
723  return .5*y*(1. - z)*(x*(2.) +
724  (1.)*(2.*x + 2.*y - 2.*z - 3.));
725 
726  case 3:
727  return .5*y*(1. - z)*((1. - x)*(-2.) +
728  (-1.)*(2.*y - 2.*x - 2.*z - 1.));
729 
730  case 4:
731  return .5*(1. - y)*z*((1. - x)*(-2.) +
732  (-1.)*(2.*z - 2.*x - 2.*y - 1.));
733 
734  case 5:
735  return .5*(1. - y)*z*(x*(2.) +
736  (1.)*(2.*x - 2.*y + 2.*z - 3.));
737 
738  case 6:
739  return .5*y*z*(x*(2.) +
740  (1.)*(2.*x + 2.*y + 2.*z - 5.));
741 
742  case 7:
743  return .5*y*z*((1. - x)*(-2.) +
744  (-1.)*(2.*y - 2.*x + 2.*z - 3.));
745 
746  case 8:
747  return 2.*(1. - y)*(1. - z)*(1. - 2.*x);
748 
749  case 9:
750  return 2.*y*(1. - y)*(1. - z);
751 
752  case 10:
753  return 2.*y*(1. - z)*(1. - 2.*x);
754 
755  case 11:
756  return 2.*y*(1. - y)*(1. - z)*(-1.);
757 
758  case 12:
759  return 2.*(1. - y)*z*(1. - z)*(-1.);
760 
761  case 13:
762  return 2.*(1. - y)*z*(1. - z);
763 
764  case 14:
765  return 2.*y*z*(1. - z);
766 
767  case 15:
768  return 2.*y*z*(1. - z)*(-1.);
769 
770  case 16:
771  return 2.*(1. - y)*z*(1. - 2.*x);
772 
773  case 17:
774  return 2.*y*(1. - y)*z;
775 
776  case 18:
777  return 2.*y*z*(1. - 2.*x);
778 
779  case 19:
780  return 2.*y*(1. - y)*z*(-1.);
781 
782  default:
783  libmesh_error_msg("Invalid i = " << i);
784  }
785 
786 
787  // d/dy*dy/deta
788  case 1:
789  switch (i)
790  {
791  case 0:
792  return .5*(1. - x)*(1. - z)*((1. - y)*(-2.) +
793  (-1.)*(1. - 2.*x - 2.*y - 2.*z));
794 
795  case 1:
796  return .5*x*(1. - z)*((1. - y)*(-2.) +
797  (-1.)*(2.*x - 2.*y - 2.*z - 1.));
798 
799  case 2:
800  return .5*x*(1. - z)*(y*(2.) +
801  (1.)*(2.*x + 2.*y - 2.*z - 3.));
802 
803  case 3:
804  return .5*(1. - x)*(1. - z)*(y*(2.) +
805  (1.)*(2.*y - 2.*x - 2.*z - 1.));
806 
807  case 4:
808  return .5*(1. - x)*z*((1. - y)*(-2.) +
809  (-1.)*(2.*z - 2.*x - 2.*y - 1.));
810 
811  case 5:
812  return .5*x*z*((1. - y)*(-2.) +
813  (-1.)*(2.*x - 2.*y + 2.*z - 3.));
814 
815  case 6:
816  return .5*x*z*(y*(2.) +
817  (1.)*(2.*x + 2.*y + 2.*z - 5.));
818 
819  case 7:
820  return .5*(1. - x)*z*(y*(2.) +
821  (1.)*(2.*y - 2.*x + 2.*z - 3.));
822 
823  case 8:
824  return 2.*x*(1. - x)*(1. - z)*(-1.);
825 
826  case 9:
827  return 2.*x*(1. - z)*(1. - 2.*y);
828 
829  case 10:
830  return 2.*x*(1. - x)*(1. - z);
831 
832  case 11:
833  return 2.*(1. - x)*(1. - z)*(1. - 2.*y);
834 
835  case 12:
836  return 2.*(1. - x)*z*(1. - z)*(-1.);
837 
838  case 13:
839  return 2.*x*z*(1. - z)*(-1.);
840 
841  case 14:
842  return 2.*x*z*(1. - z);
843 
844  case 15:
845  return 2.*(1. - x)*z*(1. - z);
846 
847  case 16:
848  return 2.*x*(1. - x)*z*(-1.);
849 
850  case 17:
851  return 2.*x*z*(1. - 2.*y);
852 
853  case 18:
854  return 2.*x*(1. - x)*z;
855 
856  case 19:
857  return 2.*(1. - x)*z*(1. - 2.*y);
858 
859  default:
860  libmesh_error_msg("Invalid i = " << i);
861  }
862 
863 
864  // d/dz*dz/dzeta
865  case 2:
866  switch (i)
867  {
868  case 0:
869  return .5*(1. - x)*(1. - y)*((1. - z)*(-2.) +
870  (-1.)*(1. - 2.*x - 2.*y - 2.*z));
871 
872  case 1:
873  return .5*x*(1. - y)*((1. - z)*(-2.) +
874  (-1.)*(2.*x - 2.*y - 2.*z - 1.));
875 
876  case 2:
877  return .5*x*y*((1. - z)*(-2.) +
878  (-1.)*(2.*x + 2.*y - 2.*z - 3.));
879 
880  case 3:
881  return .5*(1. - x)*y*((1. - z)*(-2.) +
882  (-1.)*(2.*y - 2.*x - 2.*z - 1.));
883 
884  case 4:
885  return .5*(1. - x)*(1. - y)*(z*(2.) +
886  (1.)*(2.*z - 2.*x - 2.*y - 1.));
887 
888  case 5:
889  return .5*x*(1. - y)*(z*(2.) +
890  (1.)*(2.*x - 2.*y + 2.*z - 3.));
891 
892  case 6:
893  return .5*x*y*(z*(2.) +
894  (1.)*(2.*x + 2.*y + 2.*z - 5.));
895 
896  case 7:
897  return .5*(1. - x)*y*(z*(2.) +
898  (1.)*(2.*y - 2.*x + 2.*z - 3.));
899 
900  case 8:
901  return 2.*x*(1. - x)*(1. - y)*(-1.);
902 
903  case 9:
904  return 2.*x*y*(1. - y)*(-1.);
905 
906  case 10:
907  return 2.*x*(1. - x)*y*(-1.);
908 
909  case 11:
910  return 2.*(1. - x)*y*(1. - y)*(-1.);
911 
912  case 12:
913  return 2.*(1. - x)*(1. - y)*(1. - 2.*z);
914 
915  case 13:
916  return 2.*x*(1. - y)*(1. - 2.*z);
917 
918  case 14:
919  return 2.*x*y*(1. - 2.*z);
920 
921  case 15:
922  return 2.*(1. - x)*y*(1. - 2.*z);
923 
924  case 16:
925  return 2.*x*(1. - x)*(1. - y);
926 
927  case 17:
928  return 2.*x*y*(1. - y);
929 
930  case 18:
931  return 2.*x*(1. - x)*y;
932 
933  case 19:
934  return 2.*(1. - x)*y*(1. - y);
935 
936  default:
937  libmesh_error_msg("Invalid i = " << i);
938  }
939 
940  default:
941  libmesh_error_msg("Invalid shape function derivative j = " << j);
942  }
943  }
944 
945  // triquadratic hexahedral shape functions
946  case HEX27:
947  {
948  libmesh_assert_less (i, 27);
949 
950  // Compute hex shape functions as a tensor-product
951  const Real xi = p(0);
952  const Real eta = p(1);
953  const Real zeta = p(2);
954 
955  // The only way to make any sense of this
956  // is to look at the mgflo/mg2/mgf documentation
957  // and make the cut-out cube!
958  // 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
959  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};
960  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};
961  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};
962 
963  switch(j)
964  {
965  case 0:
966  return (FE<1,L2_LAGRANGE>::shape_deriv(EDGE3, SECOND, i0[i], 0, xi)*
967  FE<1,L2_LAGRANGE>::shape (EDGE3, SECOND, i1[i], eta)*
968  FE<1,L2_LAGRANGE>::shape (EDGE3, SECOND, i2[i], zeta));
969 
970  case 1:
971  return (FE<1,L2_LAGRANGE>::shape (EDGE3, SECOND, i0[i], xi)*
972  FE<1,L2_LAGRANGE>::shape_deriv(EDGE3, SECOND, i1[i], 0, eta)*
973  FE<1,L2_LAGRANGE>::shape (EDGE3, SECOND, i2[i], zeta));
974 
975  case 2:
976  return (FE<1,L2_LAGRANGE>::shape (EDGE3, SECOND, i0[i], xi)*
977  FE<1,L2_LAGRANGE>::shape (EDGE3, SECOND, i1[i], eta)*
978  FE<1,L2_LAGRANGE>::shape_deriv(EDGE3, SECOND, i2[i], 0, zeta));
979 
980  default:
981  libmesh_error_msg("Invalid j = " << j);
982  }
983  }
984 
985  // quadratic tetrahedral shape functions
986  case TET10:
987  {
988  libmesh_assert_less (i, 10);
989 
990  // Area coordinates, pg. 205, Vol. I, Carey, Oden, Becker FEM
991  const Real zeta1 = p(0);
992  const Real zeta2 = p(1);
993  const Real zeta3 = p(2);
994  const Real zeta0 = 1. - zeta1 - zeta2 - zeta3;
995 
996  const Real dzeta0dxi = -1.;
997  const Real dzeta1dxi = 1.;
998  const Real dzeta2dxi = 0.;
999  const Real dzeta3dxi = 0.;
1000 
1001  const Real dzeta0deta = -1.;
1002  const Real dzeta1deta = 0.;
1003  const Real dzeta2deta = 1.;
1004  const Real dzeta3deta = 0.;
1005 
1006  const Real dzeta0dzeta = -1.;
1007  const Real dzeta1dzeta = 0.;
1008  const Real dzeta2dzeta = 0.;
1009  const Real dzeta3dzeta = 1.;
1010 
1011  switch (j)
1012  {
1013  // d()/dxi
1014  case 0:
1015  {
1016  switch(i)
1017  {
1018  case 0:
1019  return (4.*zeta0 - 1.)*dzeta0dxi;
1020 
1021  case 1:
1022  return (4.*zeta1 - 1.)*dzeta1dxi;
1023 
1024  case 2:
1025  return (4.*zeta2 - 1.)*dzeta2dxi;
1026 
1027  case 3:
1028  return (4.*zeta3 - 1.)*dzeta3dxi;
1029 
1030  case 4:
1031  return 4.*(zeta0*dzeta1dxi + dzeta0dxi*zeta1);
1032 
1033  case 5:
1034  return 4.*(zeta1*dzeta2dxi + dzeta1dxi*zeta2);
1035 
1036  case 6:
1037  return 4.*(zeta0*dzeta2dxi + dzeta0dxi*zeta2);
1038 
1039  case 7:
1040  return 4.*(zeta0*dzeta3dxi + dzeta0dxi*zeta3);
1041 
1042  case 8:
1043  return 4.*(zeta1*dzeta3dxi + dzeta1dxi*zeta3);
1044 
1045  case 9:
1046  return 4.*(zeta2*dzeta3dxi + dzeta2dxi*zeta3);
1047 
1048  default:
1049  libmesh_error_msg("Invalid i = " << i);
1050  }
1051  }
1052 
1053  // d()/deta
1054  case 1:
1055  {
1056  switch(i)
1057  {
1058  case 0:
1059  return (4.*zeta0 - 1.)*dzeta0deta;
1060 
1061  case 1:
1062  return (4.*zeta1 - 1.)*dzeta1deta;
1063 
1064  case 2:
1065  return (4.*zeta2 - 1.)*dzeta2deta;
1066 
1067  case 3:
1068  return (4.*zeta3 - 1.)*dzeta3deta;
1069 
1070  case 4:
1071  return 4.*(zeta0*dzeta1deta + dzeta0deta*zeta1);
1072 
1073  case 5:
1074  return 4.*(zeta1*dzeta2deta + dzeta1deta*zeta2);
1075 
1076  case 6:
1077  return 4.*(zeta0*dzeta2deta + dzeta0deta*zeta2);
1078 
1079  case 7:
1080  return 4.*(zeta0*dzeta3deta + dzeta0deta*zeta3);
1081 
1082  case 8:
1083  return 4.*(zeta1*dzeta3deta + dzeta1deta*zeta3);
1084 
1085  case 9:
1086  return 4.*(zeta2*dzeta3deta + dzeta2deta*zeta3);
1087 
1088  default:
1089  libmesh_error_msg("Invalid i = " << i);
1090  }
1091  }
1092 
1093  // d()/dzeta
1094  case 2:
1095  {
1096  switch(i)
1097  {
1098  case 0:
1099  return (4.*zeta0 - 1.)*dzeta0dzeta;
1100 
1101  case 1:
1102  return (4.*zeta1 - 1.)*dzeta1dzeta;
1103 
1104  case 2:
1105  return (4.*zeta2 - 1.)*dzeta2dzeta;
1106 
1107  case 3:
1108  return (4.*zeta3 - 1.)*dzeta3dzeta;
1109 
1110  case 4:
1111  return 4.*(zeta0*dzeta1dzeta + dzeta0dzeta*zeta1);
1112 
1113  case 5:
1114  return 4.*(zeta1*dzeta2dzeta + dzeta1dzeta*zeta2);
1115 
1116  case 6:
1117  return 4.*(zeta0*dzeta2dzeta + dzeta0dzeta*zeta2);
1118 
1119  case 7:
1120  return 4.*(zeta0*dzeta3dzeta + dzeta0dzeta*zeta3);
1121 
1122  case 8:
1123  return 4.*(zeta1*dzeta3dzeta + dzeta1dzeta*zeta3);
1124 
1125  case 9:
1126  return 4.*(zeta2*dzeta3dzeta + dzeta2dzeta*zeta3);
1127 
1128  default:
1129  libmesh_error_msg("Invalid i = " << i);
1130  }
1131  }
1132 
1133  default:
1134  libmesh_error_msg("Invalid j = " << j);
1135  }
1136  }
1137 
1138 
1139 
1140  // quadratic prism shape functions
1141  case PRISM18:
1142  {
1143  libmesh_assert_less (i, 18);
1144 
1145  // Compute prism shape functions as a tensor-product
1146  // of a triangle and an edge
1147 
1148  Point p2d(p(0),p(1));
1149  Point p1d(p(2));
1150 
1151  // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
1152  static const unsigned int i0[] = {0, 0, 0, 1, 1, 1, 0, 0, 0, 2, 2, 2, 1, 1, 1, 2, 2, 2};
1153  static const unsigned int i1[] = {0, 1, 2, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 3, 4, 5};
1154 
1155  switch (j)
1156  {
1157  // d()/dxi
1158  case 0:
1159  return (FE<2,L2_LAGRANGE>::shape_deriv(TRI6, SECOND, i1[i], 0, p2d)*
1160  FE<1,L2_LAGRANGE>::shape(EDGE3, SECOND, i0[i], p1d));
1161 
1162  // d()/deta
1163  case 1:
1164  return (FE<2,L2_LAGRANGE>::shape_deriv(TRI6, SECOND, i1[i], 1, p2d)*
1165  FE<1,L2_LAGRANGE>::shape(EDGE3, SECOND, i0[i], p1d));
1166 
1167  // d()/dzeta
1168  case 2:
1169  return (FE<2,L2_LAGRANGE>::shape(TRI6, SECOND, i1[i], p2d)*
1170  FE<1,L2_LAGRANGE>::shape_deriv(EDGE3, SECOND, i0[i], 0, p1d));
1171 
1172  default:
1173  libmesh_error_msg("Invalid shape function derivative j = " << j);
1174  }
1175  }
1176 
1177 
1178  default:
1179  libmesh_error_msg("ERROR: Unsupported 3D element type!: " << type);
1180  }
1181  }
1182 
1183 
1184  // unsupported order
1185  default:
1186  libmesh_error_msg("ERROR: Unsupported 3D FE order!: " << order);
1187  }
1188 
1189 #endif
1190 }
1191 
1192 
1193 
1194 template <>
1195 Real FE<3,L2_LAGRANGE>::shape_deriv(const Elem * elem,
1196  const Order order,
1197  const unsigned int i,
1198  const unsigned int j,
1199  const Point & p)
1200 {
1201  libmesh_assert(elem);
1202 
1203  // call the orientation-independent shape function derivatives
1204  return FE<3,L2_LAGRANGE>::shape_deriv(elem->type(), static_cast<Order>(order + elem->p_level()), i, j, p);
1205 }
1206 
1207 
1208 
1209 template <>
1210 Real FE<3,L2_LAGRANGE>::shape_second_deriv(const ElemType type,
1211  const Order order,
1212  const unsigned int i,
1213  const unsigned int j,
1214  const Point & p)
1215 {
1216 #if LIBMESH_DIM == 3
1217 
1218  libmesh_assert_less (j, 6);
1219 
1220  switch (order)
1221  {
1222  // linear Lagrange shape functions
1223  case FIRST:
1224  {
1225  return 0.;
1226  }
1227 
1228  // quadratic Lagrange shape functions
1229  case SECOND:
1230  {
1231  switch (type)
1232  {
1233 
1234  // serendipity hexahedral quadratic shape functions
1235  case HEX20:
1236  {
1237  static bool warning_given_HEX20 = false;
1238 
1239  if (!warning_given_HEX20)
1240  libMesh::err << "Second derivatives for 3D Lagrangian HEX20"
1241  << " elements are not yet implemented!"
1242  << std::endl;
1243  warning_given_HEX20 = true;
1244  }
1245  libmesh_fallthrough();
1246 
1247  case HEX27:
1248  // triquadratic hexahedral shape functions
1249  {
1250  libmesh_assert_less (i, 27);
1251 
1252  // Compute hex shape functions as a tensor-product
1253  const Real xi = p(0);
1254  const Real eta = p(1);
1255  const Real zeta = p(2);
1256 
1257  // The only way to make any sense of this
1258  // is to look at the mgflo/mg2/mgf documentation
1259  // and make the cut-out cube!
1260  // 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
1261  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};
1262  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};
1263  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};
1264 
1265  switch(j)
1266  {
1267  // d^2()/dxi^2
1268  case 0:
1269  return (FE<1,L2_LAGRANGE>::shape_second_deriv(EDGE3, SECOND, i0[i], 0, xi)*
1270  FE<1,L2_LAGRANGE>::shape (EDGE3, SECOND, i1[i], eta)*
1271  FE<1,L2_LAGRANGE>::shape (EDGE3, SECOND, i2[i], zeta));
1272 
1273  // d^2()/dxideta
1274  case 1:
1275  return (FE<1,L2_LAGRANGE>::shape_deriv(EDGE3, SECOND, i0[i], 0, xi)*
1276  FE<1,L2_LAGRANGE>::shape_deriv(EDGE3, SECOND, i1[i], 0, eta)*
1277  FE<1,L2_LAGRANGE>::shape (EDGE3, SECOND, i2[i], zeta));
1278 
1279  // d^2()/deta^2
1280  case 2:
1281  return (FE<1,L2_LAGRANGE>::shape (EDGE3, SECOND, i0[i], xi)*
1282  FE<1,L2_LAGRANGE>::shape_second_deriv(EDGE3, SECOND, i1[i], 0, eta)*
1283  FE<1,L2_LAGRANGE>::shape (EDGE3, SECOND, i2[i], zeta));
1284 
1285  // d^2()/dxidzeta
1286  case 3:
1287  return (FE<1,L2_LAGRANGE>::shape_deriv(EDGE3, SECOND, i0[i], 0, xi)*
1288  FE<1,L2_LAGRANGE>::shape (EDGE3, SECOND, i1[i], eta)*
1289  FE<1,L2_LAGRANGE>::shape_deriv(EDGE3, SECOND, i2[i], 0, zeta));
1290 
1291  // d^2()/detadzeta
1292  case 4:
1293  return (FE<1,L2_LAGRANGE>::shape (EDGE3, SECOND, i0[i], xi)*
1294  FE<1,L2_LAGRANGE>::shape_deriv(EDGE3, SECOND, i1[i], 0, eta)*
1295  FE<1,L2_LAGRANGE>::shape_deriv(EDGE3, SECOND, i2[i], 0, zeta));
1296 
1297  // d^2()/dzeta^2
1298  case 5:
1299  return (FE<1,L2_LAGRANGE>::shape (EDGE3, SECOND, i0[i], xi)*
1300  FE<1,L2_LAGRANGE>::shape (EDGE3, SECOND, i1[i], eta)*
1301  FE<1,L2_LAGRANGE>::shape_second_deriv(EDGE3, SECOND, i2[i], 0, zeta));
1302 
1303  default:
1304  libmesh_error_msg("Invalid j = " << j);
1305  }
1306  }
1307 
1308  // quadratic tetrahedral shape functions
1309  case TET10:
1310  {
1311  // The area coordinates are the same as used for the
1312  // shape() and shape_deriv() functions.
1313  // const Real zeta0 = 1. - zeta1 - zeta2 - zeta3;
1314  // const Real zeta1 = p(0);
1315  // const Real zeta2 = p(1);
1316  // const Real zeta3 = p(2);
1317  static const Real dzetadxi[4][3] =
1318  {
1319  {-1., -1., -1.},
1320  {1., 0., 0.},
1321  {0., 1., 0.},
1322  {0., 0., 1.}
1323  };
1324 
1325  // Convert from j -> (j,k) indices for independent variable
1326  // (0=xi, 1=eta, 2=zeta)
1327  static const unsigned short int independent_var_indices[6][2] =
1328  {
1329  {0, 0}, // d^2 phi / dxi^2
1330  {0, 1}, // d^2 phi / dxi deta
1331  {1, 1}, // d^2 phi / deta^2
1332  {0, 2}, // d^2 phi / dxi dzeta
1333  {1, 2}, // d^2 phi / deta dzeta
1334  {2, 2} // d^2 phi / dzeta^2
1335  };
1336 
1337  // Convert from i -> zeta indices. Each quadratic shape
1338  // function for the Tet10 depends on up to two of the zeta
1339  // area coordinate functions (see the shape() function above).
1340  // This table just tells which two area coords it uses.
1341  static const unsigned short int zeta_indices[10][2] =
1342  {
1343  {0, 0},
1344  {1, 1},
1345  {2, 2},
1346  {3, 3},
1347  {0, 1},
1348  {1, 2},
1349  {2, 0},
1350  {0, 3},
1351  {1, 3},
1352  {2, 3},
1353  };
1354 
1355  // Look up the independent variable indices for this value of j.
1356  const unsigned int my_j = independent_var_indices[j][0];
1357  const unsigned int my_k = independent_var_indices[j][1];
1358 
1359  if (i<4)
1360  {
1361  return 4.*dzetadxi[i][my_j]*dzetadxi[i][my_k];
1362  }
1363 
1364  else if (i<10)
1365  {
1366  const unsigned short int my_m = zeta_indices[i][0];
1367  const unsigned short int my_n = zeta_indices[i][1];
1368 
1369  return 4.*(dzetadxi[my_n][my_j]*dzetadxi[my_m][my_k] +
1370  dzetadxi[my_m][my_j]*dzetadxi[my_n][my_k] );
1371  }
1372  else
1373  libmesh_error_msg("Invalid shape function index " << i);
1374  }
1375 
1376 
1377  // quadratic prism shape functions
1378  case PRISM18:
1379  {
1380  libmesh_assert_less (i, 18);
1381 
1382  // Compute prism shape functions as a tensor-product
1383  // of a triangle and an edge
1384 
1385  Point p2d(p(0),p(1));
1386  Point p1d(p(2));
1387 
1388  // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
1389  static const unsigned int i0[] = {0, 0, 0, 1, 1, 1, 0, 0, 0, 2, 2, 2, 1, 1, 1, 2, 2, 2};
1390  static const unsigned int i1[] = {0, 1, 2, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 3, 4, 5};
1391 
1392  switch (j)
1393  {
1394  // d^2()/dxi^2
1395  case 0:
1396  return (FE<2,L2_LAGRANGE>::shape_second_deriv(TRI6, SECOND, i1[i], 0, p2d)*
1397  FE<1,L2_LAGRANGE>::shape(EDGE3, SECOND, i0[i], p1d));
1398 
1399  // d^2()/dxideta
1400  case 1:
1401  return (FE<2,L2_LAGRANGE>::shape_second_deriv(TRI6, SECOND, i1[i], 1, p2d)*
1402  FE<1,L2_LAGRANGE>::shape(EDGE3, SECOND, i0[i], p1d));
1403 
1404  // d^2()/deta^2
1405  case 2:
1406  return (FE<2,L2_LAGRANGE>::shape_second_deriv(TRI6, SECOND, i1[i], 2, p2d)*
1407  FE<1,L2_LAGRANGE>::shape(EDGE3, SECOND, i0[i], p1d));
1408 
1409  // d^2()/dxidzeta
1410  case 3:
1411  return (FE<2,L2_LAGRANGE>::shape_deriv(TRI6, SECOND, i1[i], 0, p2d)*
1412  FE<1,L2_LAGRANGE>::shape_deriv(EDGE3, SECOND, i0[i], 0, p1d));
1413 
1414  // d^2()/detadzeta
1415  case 4:
1416  return (FE<2,L2_LAGRANGE>::shape_deriv(TRI6, SECOND, i1[i], 1, p2d)*
1417  FE<1,L2_LAGRANGE>::shape_deriv(EDGE3, SECOND, i0[i], 0, p1d));
1418 
1419  // d^2()/dzeta^2
1420  case 5:
1421  return (FE<2,L2_LAGRANGE>::shape(TRI6, SECOND, i1[i], p2d)*
1422  FE<1,L2_LAGRANGE>::shape_second_deriv(EDGE3, SECOND, i0[i], 0, p1d));
1423 
1424  default:
1425  libmesh_error_msg("Invalid shape function derivative j = " << j);
1426  }
1427  }
1428 
1429 
1430 
1431  default:
1432  libmesh_error_msg("ERROR: Unsupported 3D element type!: " << type);
1433  }
1434  }
1435 
1436 
1437  // unsupported order
1438  default:
1439  libmesh_error_msg("ERROR: Unsupported 3D FE order!: " << order);
1440  }
1441 
1442 #endif
1443 }
1444 
1445 
1446 
1447 template <>
1448 Real FE<3,L2_LAGRANGE>::shape_second_deriv(const Elem * elem,
1449  const Order order,
1450  const unsigned int i,
1451  const unsigned int j,
1452  const Point & p)
1453 {
1454  libmesh_assert(elem);
1455 
1456  // call the orientation-independent shape function derivatives
1457  return FE<3,L2_LAGRANGE>::shape_second_deriv(elem->type(), static_cast<Order>(order + elem->p_level()), i, j, p);
1458 }
1459 
1460 } // namespace libMesh
libMesh::FE::shape
static OutputShape shape(const ElemType t, const Order o, const unsigned int i, const Point &p)
libMesh::Elem
The base class for all geometric element types.
Definition: elem.h:100
libMesh::FE::shape_deriv
static OutputShape shape_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)
libMesh::Elem::p_level
unsigned int p_level() const
Definition: elem.h:2555
libMesh::FE
Template class which generates the different FE families and orders.
Definition: fe.h:89
libMesh::err
OStreamProxy err(std::cerr)
Definition: libmesh_common.h:229
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:131
libMesh::Elem::type
virtual ElemType type() const =0
libMesh::Point
A geometric point in (x,y,z) space.
Definition: point.h:38
libMesh::FE::shape_second_deriv
static OutputShape shape_second_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)
generated by doxygen