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

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

Hosted By:
SourceForge.net Logo