libMesh
fe_lagrange_shape_3D.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2025 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 // Local includes
20 #include "libmesh/fe.h"
21 #include "libmesh/elem.h"
22 #include "libmesh/fe_lagrange_shape_1D.h"
23 #include "libmesh/enum_to_string.h"
24 #include "libmesh/cell_c0polyhedron.h"
25 #include "libmesh/tensor_value.h"
26 
27 // Anonymous namespace for functions shared by LAGRANGE and
28 // L2_LAGRANGE implementations. Implementations appear at the bottom
29 // of this file.
30 namespace
31 {
32 using namespace libMesh;
33 
34 template <FEFamily T>
35 Real fe_lagrange_3D_shape(const ElemType,
36  const Order order,
37  const Elem * elem,
38  const unsigned int i,
39  const Point & p);
40 
41 template <FEFamily T>
42 Real fe_lagrange_3D_shape_deriv(const ElemType type,
43  const Order order,
44  const Elem * elem,
45  const unsigned int i,
46  const unsigned int j,
47  const Point & p);
48 
49 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
50 
51 template <FEFamily T>
52 Real fe_lagrange_3D_shape_second_deriv(const ElemType type,
53  const Order order,
54  const Elem * elem,
55  const unsigned int i,
56  const unsigned int j,
57  const Point & p);
58 
59 #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
60 
61 } // anonymous namespace
62 
63 namespace libMesh
64 {
65 
66 
67 // TODO: If optimizations for LAGRANGE work well we should do
68 // L2_LAGRANGE too...
69 LIBMESH_DEFAULT_VECTORIZED_FE(3,L2_LAGRANGE)
70 
71 
72 template<>
73 void FE<3,LAGRANGE>::all_shapes
74  (const Elem * elem,
75  const Order o,
76  const std::vector<Point> & p,
77  std::vector<std::vector<OutputShape>> & v,
78  const bool add_p_level)
79 {
80  const ElemType type = elem->type();
81 
82  // Just loop on the harder-to-optimize cases
83  if (type != HEX8 && type != HEX27)
84  {
85  FE<3,LAGRANGE>::default_all_shapes
86  (elem,o,p,v,add_p_level);
87  return;
88  }
89 
90 #if LIBMESH_DIM == 3
91 
92  const unsigned int n_sf = v.size();
93 
94  switch (o)
95  {
96  // linear Lagrange shape functions
97  case FIRST:
98  {
99  switch (type)
100  {
101  // trilinear hexahedral shape functions
102  case HEX8:
103  case HEX20:
104  case HEX27:
105  {
106  libmesh_assert_less_equal (n_sf, 8);
107 
108  // 0 1 2 3 4 5 6 7
109  static const unsigned int i0[] = {0, 1, 1, 0, 0, 1, 1, 0};
110  static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 1, 1};
111  static const unsigned int i2[] = {0, 0, 0, 0, 1, 1, 1, 1};
112 
113  for (auto qp : index_range(p))
114  {
115  const Point & q_point = p[qp];
116  // Compute hex shape functions as a tensor-product
117  const Real xi = q_point(0);
118  const Real eta = q_point(1);
119  const Real zeta = q_point(2);
120 
121  // one_d_shapes[dim][i] = phi_i(p(dim))
122  Real one_d_shapes[3][2] = {
123  {fe_lagrange_1D_linear_shape(0, xi),
124  fe_lagrange_1D_linear_shape(1, xi)},
125  {fe_lagrange_1D_linear_shape(0, eta),
126  fe_lagrange_1D_linear_shape(1, eta)},
127  {fe_lagrange_1D_linear_shape(0, zeta),
128  fe_lagrange_1D_linear_shape(1, zeta)}};
129 
130  for (unsigned int i : make_range(n_sf))
131  v[i][qp] = one_d_shapes[0][i0[i]] *
132  one_d_shapes[1][i1[i]] *
133  one_d_shapes[2][i2[i]];
134  }
135  return;
136  }
137 
138  default:
139  libmesh_error(); // How did we get here?
140  }
141  }
142 
143 
144  // quadratic Lagrange shape functions
145  case SECOND:
146  {
147  switch (type)
148  {
149  // triquadratic hexahedral shape functions
150  case HEX8:
151 // TODO: refactor to optimize this
152 // libmesh_assert_msg(T == L2_LAGRANGE,
153 // "High order on first order elements only supported for L2 families");
154  libmesh_fallthrough();
155  case HEX27:
156  {
157  libmesh_assert_less_equal (n_sf, 27);
158 
159  // The only way to make any sense of this
160  // is to look at the mgflo/mg2/mgf documentation
161  // and make the cut-out cube!
162  // 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
163  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};
164  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};
165  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};
166 
167  for (auto qp : index_range(p))
168  {
169  const Point & q_point = p[qp];
170  // Compute hex shape functions as a tensor-product
171  const Real xi = q_point(0);
172  const Real eta = q_point(1);
173  const Real zeta = q_point(2);
174 
175  // linear_shapes[dim][i] = phi_i(p(dim))
176  Real one_d_shapes[3][3] = {
177  {fe_lagrange_1D_quadratic_shape(0, xi),
178  fe_lagrange_1D_quadratic_shape(1, xi),
179  fe_lagrange_1D_quadratic_shape(2, xi)},
180  {fe_lagrange_1D_quadratic_shape(0, eta),
181  fe_lagrange_1D_quadratic_shape(1, eta),
182  fe_lagrange_1D_quadratic_shape(2, eta)},
183  {fe_lagrange_1D_quadratic_shape(0, zeta),
184  fe_lagrange_1D_quadratic_shape(1, zeta),
185  fe_lagrange_1D_quadratic_shape(2, zeta)}};
186 
187  for (unsigned int i : make_range(n_sf))
188  v[i][qp] = one_d_shapes[0][i0[i]] *
189  one_d_shapes[1][i1[i]] *
190  one_d_shapes[2][i2[i]];
191  }
192  return;
193  }
194 
195  default:
196  libmesh_error(); // How did we get here?
197  }
198  }
199 
200  // unsupported order
201  default:
202  libmesh_error_msg("ERROR: Unsupported 3D FE order on HEX!: " << o);
203  }
204 #else // LIBMESH_DIM != 3
205  libmesh_ignore(elem, o, p, v, add_p_level);
206  libmesh_not_implemented();
207 #endif // LIBMESH_DIM == 3
208 }
209 
210 template<>
211 void FE<3,LAGRANGE>::shapes
212  (const Elem * elem,
213  const Order o,
214  const unsigned int i,
215  const std::vector<Point> & p,
216  std::vector<OutputShape> & v,
217  const bool add_p_level)
218 {
219  FE<3,LAGRANGE>::default_shapes
220  (elem,o,i,p,v,add_p_level);
221 }
222 
223 template<>
224 void FE<3,LAGRANGE>::shape_derivs
225  (const Elem * elem,
226  const Order o,
227  const unsigned int i,
228  const unsigned int j,
229  const std::vector<Point> & p,
230  std::vector<OutputShape> & v,
231  const bool add_p_level)
232 {
233  FE<3,LAGRANGE>::default_shape_derivs
234  (elem,o,i,j,p,v,add_p_level);
235 }
236 
237 template<>
238 void FE<3,LAGRANGE>::all_shape_derivs
239  (const Elem * elem,
240  const Order o,
241  const std::vector<Point> & p,
242  std::vector<std::vector<OutputShape>> * comps[3],
243  const bool add_p_level)
244 {
245  const ElemType type = elem->type();
246 
247  // Just loop on the harder-to-optimize cases
248  if (type != HEX8 && type != HEX27)
249  {
250  FE<3,LAGRANGE>::default_all_shape_derivs
251  (elem,o,p,comps,add_p_level);
252  return;
253  }
254 
255 #if LIBMESH_DIM == 3
256 
257  libmesh_assert(comps[0]);
258  libmesh_assert(comps[1]);
259  libmesh_assert(comps[2]);
260  const unsigned int n_sf = comps[0]->size();
261 
262  switch (o)
263  {
264  // linear Lagrange shape functions
265  case FIRST:
266  {
267  switch (type)
268  {
269  // trilinear hexahedral shape functions
270  case HEX8:
271  case HEX20:
272  case HEX27:
273  {
274  libmesh_assert_equal_to (n_sf, 8);
275 
276  // 0 1 2 3 4 5 6 7
277  static const unsigned int i0[] = {0, 1, 1, 0, 0, 1, 1, 0};
278  static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 1, 1};
279  static const unsigned int i2[] = {0, 0, 0, 0, 1, 1, 1, 1};
280 
281  for (auto qp : index_range(p))
282  {
283  const Point & q_point = p[qp];
284  // Compute hex shape functions as a tensor-product
285  const Real xi = q_point(0);
286  const Real eta = q_point(1);
287  const Real zeta = q_point(2);
288 
289  // one_d_shapes[dim][i] = phi_i(p(dim))
290  Real one_d_shapes[3][2] = {
291  {fe_lagrange_1D_linear_shape(0, xi),
292  fe_lagrange_1D_linear_shape(1, xi)},
293  {fe_lagrange_1D_linear_shape(0, eta),
294  fe_lagrange_1D_linear_shape(1, eta)},
295  {fe_lagrange_1D_linear_shape(0, zeta),
296  fe_lagrange_1D_linear_shape(1, zeta)}};
297 
298  // one_d_derivs[dim][i] = dphi_i/dxi(p(dim))
299  Real one_d_derivs[3][2] = {
300  {fe_lagrange_1D_linear_shape_deriv(0, 0, xi),
301  fe_lagrange_1D_linear_shape_deriv(1, 0, xi)},
302  {fe_lagrange_1D_linear_shape_deriv(0, 0, eta),
303  fe_lagrange_1D_linear_shape_deriv(1, 0, eta)},
304  {fe_lagrange_1D_linear_shape_deriv(0, 0, zeta),
305  fe_lagrange_1D_linear_shape_deriv(1, 0, zeta)}};
306 
307  for (unsigned int i : make_range(n_sf))
308  {
309  (*comps[0])[i][qp] = one_d_derivs[0][i0[i]] *
310  one_d_shapes[1][i1[i]] *
311  one_d_shapes[2][i2[i]];
312  (*comps[1])[i][qp] = one_d_shapes[0][i0[i]] *
313  one_d_derivs[1][i1[i]] *
314  one_d_shapes[2][i2[i]];
315  (*comps[2])[i][qp] = one_d_shapes[0][i0[i]] *
316  one_d_shapes[1][i1[i]] *
317  one_d_derivs[2][i2[i]];
318  }
319  }
320  return;
321  }
322 
323  default:
324  libmesh_error(); // How did we get here?
325  }
326  }
327 
328 
329  // quadratic Lagrange shape functions
330  case SECOND:
331  {
332  switch (type)
333  {
334  // triquadratic hexahedral shape functions
335  case HEX8:
336 // TODO: refactor to optimize this
337 // libmesh_assert_msg(T == L2_LAGRANGE,
338 // "High order on first order elements only supported for L2 families");
339  libmesh_fallthrough();
340  case HEX27:
341  {
342  libmesh_assert_less_equal (n_sf, 27);
343 
344  // The only way to make any sense of this
345  // is to look at the mgflo/mg2/mgf documentation
346  // and make the cut-out cube!
347  // 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
348  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};
349  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};
350  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};
351 
352  for (auto qp : index_range(p))
353  {
354  const Point & q_point = p[qp];
355  // Compute hex shape functions as a tensor-product
356  const Real xi = q_point(0);
357  const Real eta = q_point(1);
358  const Real zeta = q_point(2);
359 
360  // one_d_shapes[dim][i] = phi_i(p(dim))
361  Real one_d_shapes[3][3] = {
362  {fe_lagrange_1D_quadratic_shape(0, xi),
363  fe_lagrange_1D_quadratic_shape(1, xi),
364  fe_lagrange_1D_quadratic_shape(2, xi)},
365  {fe_lagrange_1D_quadratic_shape(0, eta),
366  fe_lagrange_1D_quadratic_shape(1, eta),
367  fe_lagrange_1D_quadratic_shape(2, eta)},
368  {fe_lagrange_1D_quadratic_shape(0, zeta),
369  fe_lagrange_1D_quadratic_shape(1, zeta),
370  fe_lagrange_1D_quadratic_shape(2, zeta)}};
371 
372  // one_d_derivs[dim][i] = dphi_i/dxi(p(dim))
373  Real one_d_derivs[3][3] = {
374  {fe_lagrange_1D_quadratic_shape_deriv(0, 0, xi),
375  fe_lagrange_1D_quadratic_shape_deriv(1, 0, xi),
376  fe_lagrange_1D_quadratic_shape_deriv(2, 0, xi)},
377  {fe_lagrange_1D_quadratic_shape_deriv(0, 0, eta),
378  fe_lagrange_1D_quadratic_shape_deriv(1, 0, eta),
379  fe_lagrange_1D_quadratic_shape_deriv(2, 0, eta)},
380  {fe_lagrange_1D_quadratic_shape_deriv(0, 0, zeta),
381  fe_lagrange_1D_quadratic_shape_deriv(1, 0, zeta),
382  fe_lagrange_1D_quadratic_shape_deriv(2, 0, zeta)}};
383 
384  for (unsigned int i : make_range(n_sf))
385  {
386  (*comps[0])[i][qp] = one_d_derivs[0][i0[i]] *
387  one_d_shapes[1][i1[i]] *
388  one_d_shapes[2][i2[i]];
389  (*comps[1])[i][qp] = one_d_shapes[0][i0[i]] *
390  one_d_derivs[1][i1[i]] *
391  one_d_shapes[2][i2[i]];
392  (*comps[2])[i][qp] = one_d_shapes[0][i0[i]] *
393  one_d_shapes[1][i1[i]] *
394  one_d_derivs[2][i2[i]];
395  }
396  }
397  return;
398  }
399 
400  default:
401  libmesh_error(); // How did we get here?
402  }
403  }
404 
405  // unsupported order
406  default:
407  libmesh_error_msg("ERROR: Unsupported 3D FE order on HEX!: " << o);
408  }
409 #else // LIBMESH_DIM != 3
410  libmesh_ignore(elem, o, p, v, add_p_level);
411  libmesh_not_implemented();
412 #endif // LIBMESH_DIM == 3
413 }
414 
415 
416 
417 
418 template <>
419 Real FE<3,LAGRANGE>::shape(const ElemType type,
420  const Order order,
421  const unsigned int i,
422  const Point & p)
423 {
424  return fe_lagrange_3D_shape<LAGRANGE>(type, order, nullptr, i, p);
425 }
426 
427 
428 
429 template <>
430 Real FE<3,L2_LAGRANGE>::shape(const ElemType type,
431  const Order order,
432  const unsigned int i,
433  const Point & p)
434 {
435  return fe_lagrange_3D_shape<L2_LAGRANGE>(type, order, nullptr, i, p);
436 }
437 
438 
439 
440 template <>
441 Real FE<3,LAGRANGE>::shape(const Elem * elem,
442  const Order order,
443  const unsigned int i,
444  const Point & p,
445  const bool add_p_level)
446 {
447  libmesh_assert(elem);
448 
449  // call the orientation-independent shape functions
450  return fe_lagrange_3D_shape<LAGRANGE>(elem->type(), order + add_p_level*elem->p_level(), elem, i, p);
451 }
452 
453 
454 
455 template <>
456 Real FE<3,L2_LAGRANGE>::shape(const Elem * elem,
457  const Order order,
458  const unsigned int i,
459  const Point & p,
460  const bool add_p_level)
461 {
462  libmesh_assert(elem);
463 
464  // call the orientation-independent shape functions
465  return fe_lagrange_3D_shape<L2_LAGRANGE>(elem->type(), order + add_p_level*elem->p_level(), elem, i, p);
466 }
467 
468 
469 
470 template <>
471 Real FE<3,LAGRANGE>::shape(const FEType fet,
472  const Elem * elem,
473  const unsigned int i,
474  const Point & p,
475  const bool add_p_level)
476 {
477  libmesh_assert(elem);
478  return fe_lagrange_3D_shape<LAGRANGE>(elem->type(), fet.order + add_p_level*elem->p_level(), elem, i, p);
479 }
480 
481 
482 
483 template <>
484 Real FE<3,L2_LAGRANGE>::shape(const FEType fet,
485  const Elem * elem,
486  const unsigned int i,
487  const Point & p,
488  const bool add_p_level)
489 {
490  libmesh_assert(elem);
491  return fe_lagrange_3D_shape<L2_LAGRANGE>(elem->type(), fet.order + add_p_level*elem->p_level(), elem, i, p);
492 }
493 
494 template <>
495 Real FE<3,LAGRANGE>::shape_deriv(const ElemType type,
496  const Order order,
497  const unsigned int i,
498  const unsigned int j,
499  const Point & p)
500 {
501  return fe_lagrange_3D_shape_deriv<LAGRANGE>(type, order, nullptr, i, j, p);
502 }
503 
504 
505 
506 template <>
507 Real FE<3,L2_LAGRANGE>::shape_deriv(const ElemType type,
508  const Order order,
509  const unsigned int i,
510  const unsigned int j,
511  const Point & p)
512 {
513  return fe_lagrange_3D_shape_deriv<L2_LAGRANGE>(type, order, nullptr, i, j, p);
514 }
515 
516 
517 
518 template <>
519 Real FE<3,LAGRANGE>::shape_deriv(const Elem * elem,
520  const Order order,
521  const unsigned int i,
522  const unsigned int j,
523  const Point & p,
524  const bool add_p_level)
525 {
526  libmesh_assert(elem);
527 
528  // call the orientation-independent shape function derivatives
529  return fe_lagrange_3D_shape_deriv<LAGRANGE>(elem->type(), order + add_p_level*elem->p_level(), elem, i, j, p);
530 }
531 
532 
533 template <>
534 Real FE<3,L2_LAGRANGE>::shape_deriv(const Elem * elem,
535  const Order order,
536  const unsigned int i,
537  const unsigned int j,
538  const Point & p,
539  const bool add_p_level)
540 {
541  libmesh_assert(elem);
542 
543  // call the orientation-independent shape function derivatives
544  return fe_lagrange_3D_shape_deriv<L2_LAGRANGE>(elem->type(), order + add_p_level*elem->p_level(), elem, i, j, p);
545 }
546 
547 
548 template <>
549 Real FE<3,LAGRANGE>::shape_deriv(const FEType fet,
550  const Elem * elem,
551  const unsigned int i,
552  const unsigned int j,
553  const Point & p,
554  const bool add_p_level)
555 {
556  libmesh_assert(elem);
557  return fe_lagrange_3D_shape_deriv<LAGRANGE>(elem->type(), fet.order + add_p_level*elem->p_level(), elem, i, j, p);
558 }
559 
560 
561 template <>
562 Real FE<3,L2_LAGRANGE>::shape_deriv(const FEType fet,
563  const Elem * elem,
564  const unsigned int i,
565  const unsigned int j,
566  const Point & p,
567  const bool add_p_level)
568 {
569  libmesh_assert(elem);
570  return fe_lagrange_3D_shape_deriv<L2_LAGRANGE>(elem->type(), fet.order + add_p_level*elem->p_level(), elem, i, j, p);
571 }
572 
573 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
574 
575 template <>
576 Real FE<3,LAGRANGE>::shape_second_deriv(const ElemType type,
577  const Order order,
578  const unsigned int i,
579  const unsigned int j,
580  const Point & p)
581 {
582  return fe_lagrange_3D_shape_second_deriv<LAGRANGE>(type, order, nullptr, i, j, p);
583 }
584 
585 
586 
587 template <>
588 Real FE<3,L2_LAGRANGE>::shape_second_deriv(const ElemType type,
589  const Order order,
590  const unsigned int i,
591  const unsigned int j,
592  const Point & p)
593 {
594  return fe_lagrange_3D_shape_second_deriv<L2_LAGRANGE>(type, order, nullptr, i, j, p);
595 }
596 
597 
598 
599 template <>
600 Real FE<3,LAGRANGE>::shape_second_deriv(const Elem * elem,
601  const Order order,
602  const unsigned int i,
603  const unsigned int j,
604  const Point & p,
605  const bool add_p_level)
606 {
607  libmesh_assert(elem);
608 
609  // call the orientation-independent shape function derivatives
610  return fe_lagrange_3D_shape_second_deriv<LAGRANGE>
611  (elem->type(), order + add_p_level*elem->p_level(), elem, i, j, p);
612 }
613 
614 
615 
616 template <>
617 Real FE<3,L2_LAGRANGE>::shape_second_deriv(const Elem * elem,
618  const Order order,
619  const unsigned int i,
620  const unsigned int j,
621  const Point & p,
622  const bool add_p_level)
623 {
624  libmesh_assert(elem);
625 
626  // call the orientation-independent shape function derivatives
627  return fe_lagrange_3D_shape_second_deriv<L2_LAGRANGE>
628  (elem->type(), order + add_p_level*elem->p_level(), elem, i, j, p);
629 }
630 
631 
632 template <>
633 Real FE<3,LAGRANGE>::shape_second_deriv(const FEType fet,
634  const Elem * elem,
635  const unsigned int i,
636  const unsigned int j,
637  const Point & p,
638  const bool add_p_level)
639 {
640  libmesh_assert(elem);
641  return fe_lagrange_3D_shape_second_deriv<LAGRANGE>
642  (elem->type(), fet.order + add_p_level*elem->p_level(), elem, i, j, p);
643 }
644 
645 
646 
647 template <>
648 Real FE<3,L2_LAGRANGE>::shape_second_deriv(const FEType fet,
649  const Elem * elem,
650  const unsigned int i,
651  const unsigned int j,
652  const Point & p,
653  const bool add_p_level)
654 {
655  libmesh_assert(elem);
656  return fe_lagrange_3D_shape_second_deriv<L2_LAGRANGE>
657  (elem->type(), fet.order + add_p_level*elem->p_level(), elem, i, j, p);
658 }
659 
660 
661 #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
662 
663 } // namespace libMesh
664 
665 
666 
667 namespace
668 {
669 using namespace libMesh;
670 
671 template <FEFamily T>
672 Real fe_lagrange_3D_shape(const ElemType type,
673  const Order order,
674  const Elem * elem,
675  const unsigned int i,
676  const Point & p)
677 {
678 #if LIBMESH_DIM == 3
679 
680  switch (order)
681  {
682  // linear Lagrange shape functions
683  case FIRST:
684  {
685  switch (type)
686  {
687  // trilinear hexahedral shape functions
688  case HEX8:
689  case HEX20:
690  case HEX27:
691  {
692  libmesh_assert_less (i, 8);
693 
694  // Compute hex shape functions as a tensor-product
695  const Real xi = p(0);
696  const Real eta = p(1);
697  const Real zeta = p(2);
698 
699  // 0 1 2 3 4 5 6 7
700  static const unsigned int i0[] = {0, 1, 1, 0, 0, 1, 1, 0};
701  static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 1, 1};
702  static const unsigned int i2[] = {0, 0, 0, 0, 1, 1, 1, 1};
703 
704  return (fe_lagrange_1D_linear_shape(i0[i], xi)*
705  fe_lagrange_1D_linear_shape(i1[i], eta)*
706  fe_lagrange_1D_linear_shape(i2[i], zeta));
707  }
708 
709  // linear tetrahedral shape functions
710  case TET4:
711  case TET10:
712  case TET14:
713  {
714  libmesh_assert_less (i, 4);
715 
716  // Area coordinates, pg. 205, Vol. I, Carey, Oden, Becker FEM
717  const Real zeta1 = p(0);
718  const Real zeta2 = p(1);
719  const Real zeta3 = p(2);
720  const Real zeta0 = 1. - zeta1 - zeta2 - zeta3;
721 
722  switch(i)
723  {
724  case 0:
725  return zeta0;
726 
727  case 1:
728  return zeta1;
729 
730  case 2:
731  return zeta2;
732 
733  case 3:
734  return zeta3;
735 
736  default:
737  libmesh_error_msg("Invalid i = " << i);
738  }
739  }
740 
741  // linear prism shape functions
742  case PRISM6:
743  case PRISM15:
744  case PRISM18:
745  case PRISM20:
746  case PRISM21:
747  {
748  libmesh_assert_less (i, 6);
749 
750  // Compute prism shape functions as a tensor-product
751  // of a triangle and an edge
752 
753  Point p2d(p(0),p(1));
754  Real p1d = p(2);
755 
756  // 0 1 2 3 4 5
757  static const unsigned int i0[] = {0, 0, 0, 1, 1, 1};
758  static const unsigned int i1[] = {0, 1, 2, 0, 1, 2};
759 
760  return (FE<2,LAGRANGE>::shape(TRI3, FIRST, i1[i], p2d)*
761  fe_lagrange_1D_linear_shape(i0[i], p1d));
762  }
763 
764  // linear pyramid shape functions
765  case PYRAMID5:
766  case PYRAMID13:
767  case PYRAMID14:
768  case PYRAMID18:
769  {
770  libmesh_assert_less (i, 5);
771 
772  const Real xi = p(0);
773  const Real eta = p(1);
774  const Real zeta = p(2);
775  const Real eps = 1.e-35;
776 
777  switch(i)
778  {
779  case 0:
780  return .25*(zeta + xi - 1.)*(zeta + eta - 1.)/((1. - zeta) + eps);
781 
782  case 1:
783  return .25*(zeta - xi - 1.)*(zeta + eta - 1.)/((1. - zeta) + eps);
784 
785  case 2:
786  return .25*(zeta - xi - 1.)*(zeta - eta - 1.)/((1. - zeta) + eps);
787 
788  case 3:
789  return .25*(zeta + xi - 1.)*(zeta - eta - 1.)/((1. - zeta) + eps);
790 
791  case 4:
792  return zeta;
793 
794  default:
795  libmesh_error_msg("Invalid i = " << i);
796  }
797  }
798 
799  case C0POLYHEDRON:
800  {
801  // Polyhedra require using newer FE APIs
802  if (!elem)
803  libmesh_error_msg("Code (see stack trace) used an outdated FE function overload.\n"
804  "Shape functions on a polyhedron are not defined by ElemType alone.");
805 
806  libmesh_assert(elem->type() == C0POLYHEDRON);
807 
808  const C0Polyhedron & poly = *cast_ptr<const C0Polyhedron *>(elem);
809 
810  // We can't use a small tolerance here, because in
811  // inverse_map() Newton might hand us intermediate
812  // iterates outside the polyhedron.
813  const auto [s, a, b, c] = poly.subelement_coordinates(p, 100);
814  if (s == invalid_uint)
815  return 0;
816  libmesh_assert_less(s, poly.n_subelements());
817 
818  const auto subtet = poly.subelement(s);
819 
820  // Avoid signed/unsigned comparison warnings
821  const int nodei = i;
822  if (nodei == subtet[0])
823  return 1-a-b-c;
824  if (nodei == subtet[1])
825  return a;
826  if (nodei == subtet[2])
827  return b;
828  if (nodei == subtet[3])
829  return c;
830 
831  // Basis function i is not supported on p's subtet
832  return 0;
833  }
834 
835  default:
836  libmesh_error_msg("ERROR: Unsupported 3D element type!: " << Utility::enum_to_string(type));
837  }
838  }
839 
840 
841  // quadratic Lagrange shape functions
842  case SECOND:
843  {
844  switch (type)
845  {
846 
847  // serendipity hexahedral quadratic shape functions
848  case HEX20:
849  {
850  libmesh_assert_less (i, 20);
851 
852  const Real xi = p(0);
853  const Real eta = p(1);
854  const Real zeta = p(2);
855 
856  // these functions are defined for (x,y,z) in [0,1]^3
857  // so transform the locations
858  const Real x = .5*(xi + 1.);
859  const Real y = .5*(eta + 1.);
860  const Real z = .5*(zeta + 1.);
861 
862  switch (i)
863  {
864  case 0:
865  return (1. - x)*(1. - y)*(1. - z)*(1. - 2.*x - 2.*y - 2.*z);
866 
867  case 1:
868  return x*(1. - y)*(1. - z)*(2.*x - 2.*y - 2.*z - 1.);
869 
870  case 2:
871  return x*y*(1. - z)*(2.*x + 2.*y - 2.*z - 3.);
872 
873  case 3:
874  return (1. - x)*y*(1. - z)*(2.*y - 2.*x - 2.*z - 1.);
875 
876  case 4:
877  return (1. - x)*(1. - y)*z*(2.*z - 2.*x - 2.*y - 1.);
878 
879  case 5:
880  return x*(1. - y)*z*(2.*x - 2.*y + 2.*z - 3.);
881 
882  case 6:
883  return x*y*z*(2.*x + 2.*y + 2.*z - 5.);
884 
885  case 7:
886  return (1. - x)*y*z*(2.*y - 2.*x + 2.*z - 3.);
887 
888  case 8:
889  return 4.*x*(1. - x)*(1. - y)*(1. - z);
890 
891  case 9:
892  return 4.*x*y*(1. - y)*(1. - z);
893 
894  case 10:
895  return 4.*x*(1. - x)*y*(1. - z);
896 
897  case 11:
898  return 4.*(1. - x)*y*(1. - y)*(1. - z);
899 
900  case 12:
901  return 4.*(1. - x)*(1. - y)*z*(1. - z);
902 
903  case 13:
904  return 4.*x*(1. - y)*z*(1. - z);
905 
906  case 14:
907  return 4.*x*y*z*(1. - z);
908 
909  case 15:
910  return 4.*(1. - x)*y*z*(1. - z);
911 
912  case 16:
913  return 4.*x*(1. - x)*(1. - y)*z;
914 
915  case 17:
916  return 4.*x*y*(1. - y)*z;
917 
918  case 18:
919  return 4.*x*(1. - x)*y*z;
920 
921  case 19:
922  return 4.*(1. - x)*y*(1. - y)*z;
923 
924  default:
925  libmesh_error_msg("Invalid i = " << i);
926  }
927  }
928 
929  // triquadratic hexahedral shape functions
930  case HEX8:
931  libmesh_assert_msg(T == L2_LAGRANGE,
932  "High order on first order elements only supported for L2 families");
933  libmesh_fallthrough();
934  case HEX27:
935  {
936  libmesh_assert_less (i, 27);
937 
938  // Compute hex shape functions as a tensor-product
939  const Real xi = p(0);
940  const Real eta = p(1);
941  const Real zeta = p(2);
942 
943  // The only way to make any sense of this
944  // is to look at the mgflo/mg2/mgf documentation
945  // and make the cut-out cube!
946  // 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
947  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};
948  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};
949  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};
950 
951  return (fe_lagrange_1D_quadratic_shape(i0[i], xi)*
952  fe_lagrange_1D_quadratic_shape(i1[i], eta)*
953  fe_lagrange_1D_quadratic_shape(i2[i], zeta));
954  }
955 
956  // quadratic tetrahedral shape functions
957  case TET4:
958  libmesh_assert_msg(T == L2_LAGRANGE,
959  "High order on first order elements only supported for L2 families");
960  libmesh_fallthrough();
961  case TET10:
962  libmesh_assert_less (i, 10);
963  libmesh_fallthrough();
964  case TET14:
965  {
966  libmesh_assert_less (i, 14);
967 
968  // Area coordinates, pg. 205, Vol. I, Carey, Oden, Becker FEM
969  const Real zeta1 = p(0);
970  const Real zeta2 = p(1);
971  const Real zeta3 = p(2);
972  const Real zeta0 = 1. - zeta1 - zeta2 - zeta3;
973 
974  switch(i)
975  {
976  case 0:
977  return zeta0*(2.*zeta0 - 1.);
978 
979  case 1:
980  return zeta1*(2.*zeta1 - 1.);
981 
982  case 2:
983  return zeta2*(2.*zeta2 - 1.);
984 
985  case 3:
986  return zeta3*(2.*zeta3 - 1.);
987 
988  case 4:
989  return 4.*zeta0*zeta1;
990 
991  case 5:
992  return 4.*zeta1*zeta2;
993 
994  case 6:
995  return 4.*zeta2*zeta0;
996 
997  case 7:
998  return 4.*zeta0*zeta3;
999 
1000  case 8:
1001  return 4.*zeta1*zeta3;
1002 
1003  case 9:
1004  return 4.*zeta2*zeta3;
1005 
1006  default:
1007  libmesh_error_msg("Invalid i = " << i);
1008  }
1009  }
1010 
1011  // "serendipity" prism
1012  case PRISM15:
1013  {
1014  libmesh_assert_less (i, 15);
1015 
1016  const Real xi = p(0);
1017  const Real eta = p(1);
1018  const Real zeta = p(2);
1019 
1020  switch(i)
1021  {
1022  case 0:
1023  return (1. - zeta)*(xi + eta - 1.)*(xi + eta + 0.5*zeta);
1024 
1025  case 1:
1026  return (1. - zeta)*xi*(xi - 1. - 0.5*zeta);
1027 
1028  case 2: // phi1 with xi <- eta
1029  return (1. - zeta)*eta*(eta - 1. - 0.5*zeta);
1030 
1031  case 3: // phi0 with zeta <- (-zeta)
1032  return (1. + zeta)*(xi + eta - 1.)*(xi + eta - 0.5*zeta);
1033 
1034  case 4: // phi1 with zeta <- (-zeta)
1035  return (1. + zeta)*xi*(xi - 1. + 0.5*zeta);
1036 
1037  case 5: // phi4 with xi <- eta
1038  return (1. + zeta)*eta*(eta - 1. + 0.5*zeta);
1039 
1040  case 6:
1041  return 2.*(1. - zeta)*xi*(1. - xi - eta);
1042 
1043  case 7:
1044  return 2.*(1. - zeta)*xi*eta;
1045 
1046  case 8:
1047  return 2.*(1. - zeta)*eta*(1. - xi - eta);
1048 
1049  case 9:
1050  return (1. - zeta)*(1. + zeta)*(1. - xi - eta);
1051 
1052  case 10:
1053  return (1. - zeta)*(1. + zeta)*xi;
1054 
1055  case 11: // phi10 with xi <-> eta
1056  return (1. - zeta)*(1. + zeta)*eta;
1057 
1058  case 12: // phi6 with zeta <- (-zeta)
1059  return 2.*(1. + zeta)*xi*(1. - xi - eta);
1060 
1061  case 13: // phi7 with zeta <- (-zeta)
1062  return 2.*(1. + zeta)*xi*eta;
1063 
1064  case 14: // phi8 with zeta <- (-zeta)
1065  return 2.*(1. + zeta)*eta*(1. - xi - eta);
1066 
1067  default:
1068  libmesh_error_msg("Invalid i = " << i);
1069  }
1070  }
1071 
1072  // quadratic prism shape functions
1073  case PRISM6:
1074  libmesh_assert_msg(T == L2_LAGRANGE,
1075  "High order on first order elements only supported for L2 families");
1076  libmesh_fallthrough();
1077  case PRISM18:
1078  case PRISM20:
1079  case PRISM21:
1080  {
1081  libmesh_assert_less (i, 18);
1082 
1083  // Compute prism shape functions as a tensor-product
1084  // of a triangle and an edge
1085 
1086  Point p2d(p(0),p(1));
1087  Real p1d = p(2);
1088 
1089  // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
1090  static const unsigned int i0[] = {0, 0, 0, 1, 1, 1, 0, 0, 0, 2, 2, 2, 1, 1, 1, 2, 2, 2};
1091  static const unsigned int i1[] = {0, 1, 2, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 3, 4, 5};
1092 
1093  return (FE<2,LAGRANGE>::shape(TRI6, SECOND, i1[i], p2d)*
1094  fe_lagrange_1D_quadratic_shape(i0[i], p1d));
1095  }
1096 
1097  // G. Bedrosian, "Shape functions and integration formulas for
1098  // three-dimensional finite element analysis", Int. J. Numerical
1099  // Methods Engineering, vol 35, p. 95-108, 1992.
1100  case PYRAMID13:
1101  {
1102  libmesh_assert_less (i, 13);
1103 
1104  const Real xi = p(0);
1105  const Real eta = p(1);
1106  const Real zeta = p(2);
1107  const Real eps = 1.e-35;
1108 
1109  // Denominators are perturbed by epsilon to avoid
1110  // divide-by-zero issues.
1111  Real den = (1. - zeta + eps);
1112 
1113  switch(i)
1114  {
1115  case 0:
1116  return 0.25*(-xi - eta - 1.)*((1. - xi)*(1. - eta) - zeta + xi*eta*zeta/den);
1117 
1118  case 1:
1119  return 0.25*(-eta + xi - 1.)*((1. + xi)*(1. - eta) - zeta - xi*eta*zeta/den);
1120 
1121  case 2:
1122  return 0.25*(xi + eta - 1.)*((1. + xi)*(1. + eta) - zeta + xi*eta*zeta/den);
1123 
1124  case 3:
1125  return 0.25*(eta - xi - 1.)*((1. - xi)*(1. + eta) - zeta - xi*eta*zeta/den);
1126 
1127  case 4:
1128  return zeta*(2.*zeta - 1.);
1129 
1130  case 5:
1131  return 0.5*(1. + xi - zeta)*(1. - xi - zeta)*(1. - eta - zeta)/den;
1132 
1133  case 6:
1134  return 0.5*(1. + eta - zeta)*(1. - eta - zeta)*(1. + xi - zeta)/den;
1135 
1136  case 7:
1137  return 0.5*(1. + xi - zeta)*(1. - xi - zeta)*(1. + eta - zeta)/den;
1138 
1139  case 8:
1140  return 0.5*(1. + eta - zeta)*(1. - eta - zeta)*(1. - xi - zeta)/den;
1141 
1142  case 9:
1143  return zeta*(1. - xi - zeta)*(1. - eta - zeta)/den;
1144 
1145  case 10:
1146  return zeta*(1. + xi - zeta)*(1. - eta - zeta)/den;
1147 
1148  case 11:
1149  return zeta*(1. + eta - zeta)*(1. + xi - zeta)/den;
1150 
1151  case 12:
1152  return zeta*(1. - xi - zeta)*(1. + eta - zeta)/den;
1153 
1154  default:
1155  libmesh_error_msg("Invalid i = " << i);
1156  }
1157  }
1158 
1159  // Quadratic shape functions, as defined in R. Graglia, "Higher order
1160  // bases on pyramidal elements", IEEE Trans Antennas and Propagation,
1161  // vol 47, no 5, May 1999.
1162  case PYRAMID5:
1163  libmesh_assert_msg(T == L2_LAGRANGE,
1164  "High order on first order elements only supported for L2 families");
1165  libmesh_fallthrough();
1166  case PYRAMID14:
1167  case PYRAMID18:
1168  {
1169  libmesh_assert_less (i, 14);
1170 
1171  const Real xi = p(0);
1172  const Real eta = p(1);
1173  const Real zeta = p(2);
1174  const Real eps = 1.e-35;
1175 
1176  // The "normalized coordinates" defined by Graglia. These are
1177  // the planes which define the faces of the pyramid.
1178  Real
1179  p1 = 0.5*(1. - eta - zeta), // back
1180  p2 = 0.5*(1. + xi - zeta), // left
1181  p3 = 0.5*(1. + eta - zeta), // front
1182  p4 = 0.5*(1. - xi - zeta); // right
1183 
1184  // Denominators are perturbed by epsilon to avoid
1185  // divide-by-zero issues.
1186  Real
1187  den = (-1. + zeta + eps),
1188  den2 = den*den;
1189 
1190  switch(i)
1191  {
1192  case 0:
1193  return p4*p1*(xi*eta - zeta + zeta*zeta)/den2;
1194 
1195  case 1:
1196  return -p1*p2*(xi*eta + zeta - zeta*zeta)/den2;
1197 
1198  case 2:
1199  return p2*p3*(xi*eta - zeta + zeta*zeta)/den2;
1200 
1201  case 3:
1202  return -p3*p4*(xi*eta + zeta - zeta*zeta)/den2;
1203 
1204  case 4:
1205  return zeta*(2.*zeta - 1.);
1206 
1207  case 5:
1208  return -4.*p2*p1*p4*eta/den2;
1209 
1210  case 6:
1211  return 4.*p1*p2*p3*xi/den2;
1212 
1213  case 7:
1214  return 4.*p2*p3*p4*eta/den2;
1215 
1216  case 8:
1217  return -4.*p3*p4*p1*xi/den2;
1218 
1219  case 9:
1220  return -4.*p1*p4*zeta/den;
1221 
1222  case 10:
1223  return -4.*p2*p1*zeta/den;
1224 
1225  case 11:
1226  return -4.*p3*p2*zeta/den;
1227 
1228  case 12:
1229  return -4.*p4*p3*zeta/den;
1230 
1231  case 13:
1232  return 16.*p1*p2*p3*p4/den2;
1233 
1234  default:
1235  libmesh_error_msg("Invalid i = " << i);
1236  }
1237  }
1238 
1239 
1240  default:
1241  libmesh_error_msg("ERROR: Unsupported 3D element type!: " << Utility::enum_to_string(type));
1242  }
1243  }
1244 
1245  case THIRD:
1246  {
1247  switch (type)
1248  {
1249  // quadratic Lagrange shape functions with cubic bubbles
1250  case PRISM20:
1251  {
1252  libmesh_assert_less (i, 20);
1253 
1254  // Compute Prism21 shape functions as a tensor-product
1255  // of a triangle and an edge, then redistribute the
1256  // central bubble function over the other Tri6 nodes
1257  // around it (in a way consistent with the Tri7 shape
1258  // function definitions).
1259  // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
1260  static const unsigned int i0[] = {0, 0, 0, 1, 1, 1, 0, 0, 0, 2, 2, 2, 1, 1, 1, 2, 2, 2, 0, 1};
1261 
1262  Point p2d(p(0),p(1));
1263  Real p1d = p(2);
1264 
1265  const Real mainval = FE<3,LAGRANGE>::shape(PRISM21, THIRD, i, p);
1266 
1267  if (i0[i] != 2)
1268  return mainval;
1269 
1270  const Real bubbleval =
1271  FE<2,LAGRANGE>::shape(TRI7, THIRD, 6, p2d) *
1272  fe_lagrange_1D_quadratic_shape(2, p1d);
1273 
1274  if (i < 12) // vertices
1275  return mainval - bubbleval / 9;
1276 
1277  return mainval + bubbleval * (Real(4) / 9);
1278  }
1279 
1280  // quadratic Lagrange shape functions with cubic bubbles
1281  case PRISM21:
1282  {
1283  libmesh_assert_less (i, 21);
1284 
1285  // Compute prism shape functions as a tensor-product
1286  // of a triangle and an edge
1287 
1288  Point p2d(p(0),p(1));
1289  Real p1d = p(2);
1290 
1291  // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
1292  static const unsigned int i0[] = {0, 0, 0, 1, 1, 1, 0, 0, 0, 2, 2, 2, 1, 1, 1, 2, 2, 2, 0, 1, 2};
1293  static const unsigned int i1[] = {0, 1, 2, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 3, 4, 5, 6, 6, 6};
1294 
1295  return (FE<2,LAGRANGE>::shape(TRI7, THIRD, i1[i], p2d)*
1296  fe_lagrange_1D_quadratic_shape(i0[i], p1d));
1297  }
1298 
1299  // Weird rational shape functions with weirder bubbles...
1300  case PYRAMID18:
1301  {
1302  libmesh_assert_less (i, 18);
1303 
1304  const Real xi = p(0);
1305  const Real eta = p(1);
1306  const Real zeta = p(2);
1307  const Real eps = 1.e-35;
1308 
1309  // The "normalized coordinates" defined by Graglia. These are
1310  // the planes which define the faces of the pyramid.
1311  const Real
1312  p1 = 0.5*(1. - eta - zeta), // back
1313  p2 = 0.5*(1. + xi - zeta), // left
1314  p3 = 0.5*(1. + eta - zeta), // front
1315  p4 = 0.5*(1. - xi - zeta); // right
1316 
1317  // Denominators are perturbed by epsilon to avoid
1318  // divide-by-zero issues.
1319  const Real
1320  den = (-1. + zeta + eps),
1321  den2 = den*den;
1322 
1323  // Bubble functions on triangular sides. We actually
1324  // have a degree of freedom to play with here, and I'm
1325  // not certain how best to use it, so let's leave it
1326  // as a variable in case we figure that out later.
1327  constexpr Real alpha = 0.5;
1328  const Real
1329  bub_f1 = ((1-alpha)*(1-zeta) + alpha*(-eta)),
1330  bub_f2 = ((1-alpha)*(1-zeta) + alpha*(xi)),
1331  bub_f3 = ((1-alpha)*(1-zeta) + alpha*(eta)),
1332  bub_f4 = ((1-alpha)*(1-zeta) + alpha*(-xi));
1333 
1334  const Real
1335  bub1 = bub_f1*p1*p2*p4*zeta/den2,
1336  bub2 = bub_f2*p1*p2*p3*zeta/den2,
1337  bub3 = bub_f3*p2*p3*p4*zeta/den2,
1338  bub4 = bub_f4*p1*p3*p4*zeta/den2;
1339 
1340  switch(i)
1341  {
1342  case 0:
1343  return p4*p1*(xi*eta - zeta + zeta*zeta)/den2 + 3*(bub1+bub4);
1344 
1345  case 1:
1346  return -p1*p2*(xi*eta + zeta - zeta*zeta)/den2 + 3*(bub1+bub2);
1347 
1348  case 2:
1349  return p2*p3*(xi*eta - zeta + zeta*zeta)/den2 + 3*(bub2+bub3);
1350 
1351  case 3:
1352  return -p3*p4*(xi*eta + zeta - zeta*zeta)/den2 + 3*(bub3+bub4);
1353 
1354  case 4:
1355  return zeta*(2.*zeta - 1.) + 3*(bub1+bub2+bub3+bub4);
1356 
1357  case 5:
1358  return -4.*p2*p1*p4*eta/den2 - 12*bub1;
1359 
1360  case 6:
1361  return 4.*p1*p2*p3*xi/den2 - 12*bub2;
1362 
1363  case 7:
1364  return 4.*p2*p3*p4*eta/den2 - 12*bub3;
1365 
1366  case 8:
1367  return -4.*p3*p4*p1*xi/den2 - 12*bub4;
1368 
1369  case 9:
1370  return -4.*p1*p4*zeta/den - 12*(bub1+bub4);
1371 
1372  case 10:
1373  return -4.*p2*p1*zeta/den - 12*(bub1+bub2);
1374 
1375  case 11:
1376  return -4.*p3*p2*zeta/den - 12*(bub2+bub3);
1377 
1378  case 12:
1379  return -4.*p4*p3*zeta/den - 12*(bub3+bub4);
1380 
1381  case 13:
1382  return 16.*p1*p2*p3*p4/den2;
1383 
1384  case 14:
1385  return 27*bub1;
1386 
1387  case 15:
1388  return 27*bub2;
1389 
1390  case 16:
1391  return 27*bub3;
1392 
1393  case 17:
1394  return 27*bub4;
1395 
1396  default:
1397  libmesh_error_msg("Invalid i = " << i);
1398  }
1399  }
1400 
1401  // quadratic Lagrange shape functions with cubic bubbles
1402  case TET14:
1403  {
1404  libmesh_assert_less (i, 14);
1405 
1406  // Area coordinates, pg. 205, Vol. I, Carey, Oden, Becker FEM
1407  const Real zeta1 = p(0);
1408  const Real zeta2 = p(1);
1409  const Real zeta3 = p(2);
1410  const Real zeta0 = 1. - zeta1 - zeta2 - zeta3;
1411 
1412  // Bubble functions (not yet scaled) on side nodes
1413  const Real bubble_012 = zeta0*zeta1*zeta2;
1414  const Real bubble_013 = zeta0*zeta1*zeta3;
1415  const Real bubble_123 = zeta1*zeta2*zeta3;
1416  const Real bubble_023 = zeta0*zeta2*zeta3;
1417 
1418  switch(i)
1419  {
1420  case 0:
1421  return zeta0*(2.*zeta0 - 1.) + 3.*(bubble_012+bubble_013+bubble_023);
1422 
1423  case 1:
1424  return zeta1*(2.*zeta1 - 1.) + 3.*(bubble_012+bubble_013+bubble_123);
1425 
1426  case 2:
1427  return zeta2*(2.*zeta2 - 1.) + 3.*(bubble_012+bubble_023+bubble_123);
1428 
1429  case 3:
1430  return zeta3*(2.*zeta3 - 1.) + 3.*(bubble_013+bubble_023+bubble_123);
1431 
1432  case 4:
1433  return 4.*zeta0*zeta1 - 12.*(bubble_012+bubble_013);
1434 
1435  case 5:
1436  return 4.*zeta1*zeta2 - 12.*(bubble_012+bubble_123);
1437 
1438  case 6:
1439  return 4.*zeta2*zeta0 - 12.*(bubble_012+bubble_023);
1440 
1441  case 7:
1442  return 4.*zeta0*zeta3 - 12.*(bubble_013+bubble_023);
1443 
1444  case 8:
1445  return 4.*zeta1*zeta3 - 12.*(bubble_013+bubble_123);
1446 
1447  case 9:
1448  return 4.*zeta2*zeta3 - 12.*(bubble_023+bubble_123);
1449 
1450  case 10:
1451  return 27.*bubble_012;
1452 
1453  case 11:
1454  return 27.*bubble_013;
1455 
1456  case 12:
1457  return 27.*bubble_123;
1458 
1459  case 13:
1460  return 27.*bubble_023;
1461 
1462  default:
1463  libmesh_error_msg("Invalid i = " << i);
1464  }
1465  }
1466 
1467  default:
1468  libmesh_error_msg("ERROR: Unsupported 3D element type!: " << Utility::enum_to_string(type));
1469  }
1470  }
1471 
1472  // unsupported order
1473  default:
1474  libmesh_error_msg("ERROR: Unsupported 3D FE order!: " << order);
1475  }
1476 
1477 #else // LIBMESH_DIM != 3
1478  libmesh_ignore(type, order, elem, i, p);
1479  libmesh_not_implemented();
1480 #endif
1481 }
1482 
1483 
1484 
1485 template <FEFamily T>
1486 Real fe_lagrange_3D_shape_deriv(const ElemType type,
1487  const Order order,
1488  const Elem * elem,
1489  const unsigned int i,
1490  const unsigned int j,
1491  const Point & p)
1492 {
1493 #if LIBMESH_DIM == 3
1494 
1495  libmesh_assert_less (j, 3);
1496 
1497  switch (order)
1498  {
1499  // linear Lagrange shape functions
1500  case FIRST:
1501  {
1502  switch (type)
1503  {
1504  // trilinear hexahedral shape functions
1505  case HEX8:
1506  case HEX20:
1507  case HEX27:
1508  {
1509  libmesh_assert_less (i, 8);
1510 
1511  // Compute hex shape functions as a tensor-product
1512  const Real xi = p(0);
1513  const Real eta = p(1);
1514  const Real zeta = p(2);
1515 
1516  static const unsigned int i0[] = {0, 1, 1, 0, 0, 1, 1, 0};
1517  static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 1, 1};
1518  static const unsigned int i2[] = {0, 0, 0, 0, 1, 1, 1, 1};
1519 
1520  switch(j)
1521  {
1522  case 0:
1523  return (fe_lagrange_1D_linear_shape_deriv(i0[i], 0, xi)*
1524  fe_lagrange_1D_linear_shape (i1[i], eta)*
1525  fe_lagrange_1D_linear_shape (i2[i], zeta));
1526 
1527  case 1:
1528  return (fe_lagrange_1D_linear_shape (i0[i], xi)*
1529  fe_lagrange_1D_linear_shape_deriv(i1[i], 0, eta)*
1530  fe_lagrange_1D_linear_shape (i2[i], zeta));
1531 
1532  case 2:
1533  return (fe_lagrange_1D_linear_shape (i0[i], xi)*
1534  fe_lagrange_1D_linear_shape (i1[i], eta)*
1535  fe_lagrange_1D_linear_shape_deriv(i2[i], 0, zeta));
1536 
1537  default:
1538  libmesh_error_msg("Invalid j = " << j);
1539  }
1540  }
1541 
1542  // linear tetrahedral shape functions
1543  case TET4:
1544  case TET10:
1545  case TET14:
1546  {
1547  libmesh_assert_less (i, 4);
1548 
1549  // Area coordinates, pg. 205, Vol. I, Carey, Oden, Becker FEM
1550  const Real dzeta0dxi = -1.;
1551  const Real dzeta1dxi = 1.;
1552  const Real dzeta2dxi = 0.;
1553  const Real dzeta3dxi = 0.;
1554 
1555  const Real dzeta0deta = -1.;
1556  const Real dzeta1deta = 0.;
1557  const Real dzeta2deta = 1.;
1558  const Real dzeta3deta = 0.;
1559 
1560  const Real dzeta0dzeta = -1.;
1561  const Real dzeta1dzeta = 0.;
1562  const Real dzeta2dzeta = 0.;
1563  const Real dzeta3dzeta = 1.;
1564 
1565  switch (j)
1566  {
1567  // d()/dxi
1568  case 0:
1569  {
1570  switch(i)
1571  {
1572  case 0:
1573  return dzeta0dxi;
1574 
1575  case 1:
1576  return dzeta1dxi;
1577 
1578  case 2:
1579  return dzeta2dxi;
1580 
1581  case 3:
1582  return dzeta3dxi;
1583 
1584  default:
1585  libmesh_error_msg("Invalid i = " << i);
1586  }
1587  }
1588 
1589  // d()/deta
1590  case 1:
1591  {
1592  switch(i)
1593  {
1594  case 0:
1595  return dzeta0deta;
1596 
1597  case 1:
1598  return dzeta1deta;
1599 
1600  case 2:
1601  return dzeta2deta;
1602 
1603  case 3:
1604  return dzeta3deta;
1605 
1606  default:
1607  libmesh_error_msg("Invalid i = " << i);
1608  }
1609  }
1610 
1611  // d()/dzeta
1612  case 2:
1613  {
1614  switch(i)
1615  {
1616  case 0:
1617  return dzeta0dzeta;
1618 
1619  case 1:
1620  return dzeta1dzeta;
1621 
1622  case 2:
1623  return dzeta2dzeta;
1624 
1625  case 3:
1626  return dzeta3dzeta;
1627 
1628  default:
1629  libmesh_error_msg("Invalid i = " << i);
1630  }
1631  }
1632 
1633  default:
1634  libmesh_error_msg("Invalid shape function derivative j = " << j);
1635  }
1636  }
1637 
1638  // linear prism shape functions
1639  case PRISM6:
1640  case PRISM15:
1641  case PRISM18:
1642  case PRISM20:
1643  case PRISM21:
1644  {
1645  libmesh_assert_less (i, 6);
1646 
1647  // Compute prism shape functions as a tensor-product
1648  // of a triangle and an edge
1649 
1650  Point p2d(p(0),p(1));
1651  Real p1d = p(2);
1652 
1653  // 0 1 2 3 4 5
1654  static const unsigned int i0[] = {0, 0, 0, 1, 1, 1};
1655  static const unsigned int i1[] = {0, 1, 2, 0, 1, 2};
1656 
1657  switch (j)
1658  {
1659  // d()/dxi
1660  case 0:
1661  return (FE<2,LAGRANGE>::shape_deriv(TRI3, FIRST, i1[i], 0, p2d)*
1662  fe_lagrange_1D_linear_shape(i0[i], p1d));
1663 
1664  // d()/deta
1665  case 1:
1666  return (FE<2,LAGRANGE>::shape_deriv(TRI3, FIRST, i1[i], 1, p2d)*
1667  fe_lagrange_1D_linear_shape(i0[i], p1d));
1668 
1669  // d()/dzeta
1670  case 2:
1671  return (FE<2,LAGRANGE>::shape(TRI3, FIRST, i1[i], p2d)*
1672  fe_lagrange_1D_linear_shape_deriv(i0[i], 0, p1d));
1673 
1674  default:
1675  libmesh_error_msg("Invalid shape function derivative j = " << j);
1676  }
1677  }
1678 
1679  // linear pyramid shape functions
1680  case PYRAMID5:
1681  case PYRAMID13:
1682  case PYRAMID14:
1683  case PYRAMID18:
1684  {
1685  libmesh_assert_less (i, 5);
1686 
1687  const Real xi = p(0);
1688  const Real eta = p(1);
1689  const Real zeta = p(2);
1690  const Real eps = 1.e-35;
1691 
1692  switch (j)
1693  {
1694  // d/dxi
1695  case 0:
1696  switch(i)
1697  {
1698  case 0:
1699  return .25*(zeta + eta - 1.)/((1. - zeta) + eps);
1700 
1701  case 1:
1702  return -.25*(zeta + eta - 1.)/((1. - zeta) + eps);
1703 
1704  case 2:
1705  return -.25*(zeta - eta - 1.)/((1. - zeta) + eps);
1706 
1707  case 3:
1708  return .25*(zeta - eta - 1.)/((1. - zeta) + eps);
1709 
1710  case 4:
1711  return 0;
1712 
1713  default:
1714  libmesh_error_msg("Invalid i = " << i);
1715  }
1716 
1717 
1718  // d/deta
1719  case 1:
1720  switch(i)
1721  {
1722  case 0:
1723  return .25*(zeta + xi - 1.)/((1. - zeta) + eps);
1724 
1725  case 1:
1726  return .25*(zeta - xi - 1.)/((1. - zeta) + eps);
1727 
1728  case 2:
1729  return -.25*(zeta - xi - 1.)/((1. - zeta) + eps);
1730 
1731  case 3:
1732  return -.25*(zeta + xi - 1.)/((1. - zeta) + eps);
1733 
1734  case 4:
1735  return 0;
1736 
1737  default:
1738  libmesh_error_msg("Invalid i = " << i);
1739  }
1740 
1741 
1742  // d/dzeta
1743  case 2:
1744  {
1745  // We computed the derivatives with general eps and
1746  // then let eps tend to zero in the numerators...
1747  Real
1748  num = zeta*(2. - zeta) - 1.,
1749  den = (1. - zeta + eps)*(1. - zeta + eps);
1750 
1751  switch(i)
1752  {
1753  case 0:
1754  case 2:
1755  return .25*(num + xi*eta)/den;
1756 
1757  case 1:
1758  case 3:
1759  return .25*(num - xi*eta)/den;
1760 
1761  case 4:
1762  return 1.;
1763 
1764  default:
1765  libmesh_error_msg("Invalid i = " << i);
1766  }
1767  }
1768 
1769  default:
1770  libmesh_error_msg("Invalid j = " << j);
1771  }
1772  }
1773 
1774  case C0POLYHEDRON:
1775  {
1776  // Polyhedra require using newer FE APIs
1777  if (!elem)
1778  libmesh_error_msg("Code (see stack trace) used an outdated FE function overload.\n"
1779  "Shape functions on a polyhedron are not defined by ElemType alone.");
1780 
1781  libmesh_assert(elem->type() == C0POLYHEDRON);
1782 
1783  const C0Polyhedron & poly = *cast_ptr<const C0Polyhedron *>(elem);
1784 
1785  // We can't use a small tolerance here, because in
1786  // inverse_map() Newton might hand us intermediate
1787  // iterates outside the polyhedron.
1788  const auto [s, a, b, c] = poly.subelement_coordinates(p, 100);
1789  if (s == invalid_uint)
1790  return 0;
1791  libmesh_assert_less(s, poly.n_subelements());
1792 
1793  const auto subtet = poly.subelement(s);
1794 
1795  // Find derivatives w.r.t. subtriangle barycentric
1796  // coordinates
1797  Real du_da = 0, du_db = 0, du_dc = 0;
1798 
1799  // Avoid signed/unsigned comparison warnings
1800  const int nodei = i;
1801  if (nodei == subtet[0])
1802  du_da = du_db = du_dc = -1;
1803  else if (nodei == subtet[1])
1804  du_da = 1;
1805  else if (nodei == subtet[2])
1806  du_db = 1;
1807  else if (nodei == subtet[3])
1808  du_dc = 1;
1809  else
1810  // Basis function i is not supported on p's subtet
1811  return 0;
1812 
1813  // We want to return derivatives with respect to
1814  // xi/eta/zeta for the polyhedron, but what we
1815  // calculated above are with respect to xi and eta
1816  // coordinates for a master *triangle*. We need to
1817  // convert from one to the other.
1818 
1819  const auto master_points = poly.master_subelement(s);
1820 
1821  const RealTensor dXi_dA(
1822  master_points[1](0) - master_points[0](0), master_points[2](0) - master_points[0](0), master_points[3](0) - master_points[0](0),
1823  master_points[1](1) - master_points[0](1), master_points[2](1) - master_points[0](1), master_points[3](1) - master_points[0](1),
1824  master_points[1](2) - master_points[0](2), master_points[2](2) - master_points[0](2), master_points[3](2) - master_points[0](2));
1825 
1826  // When we vectorize this we'll want a full inverse, but
1827  // when we're querying one component at a time it's
1828  // cheaper to manually compute a single column.
1829  // const RealTensor dabc_dxietazeta_dabc = dxietazeta_dabc.inverse();
1830  const Real jac = dXi_dA.det();
1831 
1832  switch (j)
1833  {
1834  // d()/dxi
1835  case 0:
1836  {
1837  const Real da_dxi = (dXi_dA(2,2)*dXi_dA(1,1) - dXi_dA(2,1)*dXi_dA(1,2)) / jac;
1838  const Real db_dxi = -(dXi_dA(2,2)*dXi_dA(1,0) - dXi_dA(2,0)*dXi_dA(1,2)) / jac;
1839  const Real dc_dxi = (dXi_dA(2,1)*dXi_dA(1,0) - dXi_dA(2,0)*dXi_dA(1,1)) / jac;
1840  return du_da*da_dxi + du_db*db_dxi + du_dc*dc_dxi;
1841  }
1842  // d()/deta
1843  case 1:
1844  {
1845  const Real da_deta = -(dXi_dA(2,2)*dXi_dA(0,1) - dXi_dA(2,1)*dXi_dA(0,2) ) / jac;
1846  const Real db_deta = (dXi_dA(2,2)*dXi_dA(0,0) - dXi_dA(2,0)*dXi_dA(0,2)) / jac;
1847  const Real dc_deta = -(dXi_dA(2,1)*dXi_dA(0,0) - dXi_dA(2,0)*dXi_dA(0,1)) / jac;
1848  return du_da*da_deta + du_db*db_deta + du_dc*dc_deta;
1849  }
1850  // d()/dzeta
1851  case 2:
1852  {
1853  const Real da_dzeta = (dXi_dA(1,2)*dXi_dA(0,1) - dXi_dA(1,1)*dXi_dA(0,2) ) / jac;
1854  const Real db_dzeta = -(dXi_dA(1,2)*dXi_dA(0,0) - dXi_dA(1,0)*dXi_dA(0,2)) / jac;
1855  const Real dc_dzeta = (dXi_dA(1,1)*dXi_dA(0,0) - dXi_dA(1,0)*dXi_dA(0,1)) / jac;
1856  return du_da*da_dzeta + du_db*db_dzeta + du_dc*dc_dzeta;
1857  }
1858  default:
1859  libmesh_error_msg("ERROR: Invalid derivative index j = " << j);
1860  }
1861  }
1862 
1863  default:
1864  libmesh_error_msg("ERROR: Unsupported 3D element type!: " << Utility::enum_to_string(type));
1865  }
1866  }
1867 
1868 
1869  // quadratic Lagrange shape functions
1870  case SECOND:
1871  {
1872  switch (type)
1873  {
1874 
1875  // serendipity hexahedral quadratic shape functions
1876  case HEX20:
1877  {
1878  libmesh_assert_less (i, 20);
1879 
1880  const Real xi = p(0);
1881  const Real eta = p(1);
1882  const Real zeta = p(2);
1883 
1884  // these functions are defined for (x,y,z) in [0,1]^3
1885  // so transform the locations
1886  const Real x = .5*(xi + 1.);
1887  const Real y = .5*(eta + 1.);
1888  const Real z = .5*(zeta + 1.);
1889 
1890  // and don't forget the chain rule!
1891 
1892  switch (j)
1893  {
1894 
1895  // d/dx*dx/dxi
1896  case 0:
1897  switch (i)
1898  {
1899  case 0:
1900  return .5*(1. - y)*(1. - z)*((1. - x)*(-2.) +
1901  (-1.)*(1. - 2.*x - 2.*y - 2.*z));
1902 
1903  case 1:
1904  return .5*(1. - y)*(1. - z)*(x*(2.) +
1905  (1.)*(2.*x - 2.*y - 2.*z - 1.));
1906 
1907  case 2:
1908  return .5*y*(1. - z)*(x*(2.) +
1909  (1.)*(2.*x + 2.*y - 2.*z - 3.));
1910 
1911  case 3:
1912  return .5*y*(1. - z)*((1. - x)*(-2.) +
1913  (-1.)*(2.*y - 2.*x - 2.*z - 1.));
1914 
1915  case 4:
1916  return .5*(1. - y)*z*((1. - x)*(-2.) +
1917  (-1.)*(2.*z - 2.*x - 2.*y - 1.));
1918 
1919  case 5:
1920  return .5*(1. - y)*z*(x*(2.) +
1921  (1.)*(2.*x - 2.*y + 2.*z - 3.));
1922 
1923  case 6:
1924  return .5*y*z*(x*(2.) +
1925  (1.)*(2.*x + 2.*y + 2.*z - 5.));
1926 
1927  case 7:
1928  return .5*y*z*((1. - x)*(-2.) +
1929  (-1.)*(2.*y - 2.*x + 2.*z - 3.));
1930 
1931  case 8:
1932  return 2.*(1. - y)*(1. - z)*(1. - 2.*x);
1933 
1934  case 9:
1935  return 2.*y*(1. - y)*(1. - z);
1936 
1937  case 10:
1938  return 2.*y*(1. - z)*(1. - 2.*x);
1939 
1940  case 11:
1941  return 2.*y*(1. - y)*(1. - z)*(-1.);
1942 
1943  case 12:
1944  return 2.*(1. - y)*z*(1. - z)*(-1.);
1945 
1946  case 13:
1947  return 2.*(1. - y)*z*(1. - z);
1948 
1949  case 14:
1950  return 2.*y*z*(1. - z);
1951 
1952  case 15:
1953  return 2.*y*z*(1. - z)*(-1.);
1954 
1955  case 16:
1956  return 2.*(1. - y)*z*(1. - 2.*x);
1957 
1958  case 17:
1959  return 2.*y*(1. - y)*z;
1960 
1961  case 18:
1962  return 2.*y*z*(1. - 2.*x);
1963 
1964  case 19:
1965  return 2.*y*(1. - y)*z*(-1.);
1966 
1967  default:
1968  libmesh_error_msg("Invalid i = " << i);
1969  }
1970 
1971 
1972  // d/dy*dy/deta
1973  case 1:
1974  switch (i)
1975  {
1976  case 0:
1977  return .5*(1. - x)*(1. - z)*((1. - y)*(-2.) +
1978  (-1.)*(1. - 2.*x - 2.*y - 2.*z));
1979 
1980  case 1:
1981  return .5*x*(1. - z)*((1. - y)*(-2.) +
1982  (-1.)*(2.*x - 2.*y - 2.*z - 1.));
1983 
1984  case 2:
1985  return .5*x*(1. - z)*(y*(2.) +
1986  (1.)*(2.*x + 2.*y - 2.*z - 3.));
1987 
1988  case 3:
1989  return .5*(1. - x)*(1. - z)*(y*(2.) +
1990  (1.)*(2.*y - 2.*x - 2.*z - 1.));
1991 
1992  case 4:
1993  return .5*(1. - x)*z*((1. - y)*(-2.) +
1994  (-1.)*(2.*z - 2.*x - 2.*y - 1.));
1995 
1996  case 5:
1997  return .5*x*z*((1. - y)*(-2.) +
1998  (-1.)*(2.*x - 2.*y + 2.*z - 3.));
1999 
2000  case 6:
2001  return .5*x*z*(y*(2.) +
2002  (1.)*(2.*x + 2.*y + 2.*z - 5.));
2003 
2004  case 7:
2005  return .5*(1. - x)*z*(y*(2.) +
2006  (1.)*(2.*y - 2.*x + 2.*z - 3.));
2007 
2008  case 8:
2009  return 2.*x*(1. - x)*(1. - z)*(-1.);
2010 
2011  case 9:
2012  return 2.*x*(1. - z)*(1. - 2.*y);
2013 
2014  case 10:
2015  return 2.*x*(1. - x)*(1. - z);
2016 
2017  case 11:
2018  return 2.*(1. - x)*(1. - z)*(1. - 2.*y);
2019 
2020  case 12:
2021  return 2.*(1. - x)*z*(1. - z)*(-1.);
2022 
2023  case 13:
2024  return 2.*x*z*(1. - z)*(-1.);
2025 
2026  case 14:
2027  return 2.*x*z*(1. - z);
2028 
2029  case 15:
2030  return 2.*(1. - x)*z*(1. - z);
2031 
2032  case 16:
2033  return 2.*x*(1. - x)*z*(-1.);
2034 
2035  case 17:
2036  return 2.*x*z*(1. - 2.*y);
2037 
2038  case 18:
2039  return 2.*x*(1. - x)*z;
2040 
2041  case 19:
2042  return 2.*(1. - x)*z*(1. - 2.*y);
2043 
2044  default:
2045  libmesh_error_msg("Invalid i = " << i);
2046  }
2047 
2048 
2049  // d/dz*dz/dzeta
2050  case 2:
2051  switch (i)
2052  {
2053  case 0:
2054  return .5*(1. - x)*(1. - y)*((1. - z)*(-2.) +
2055  (-1.)*(1. - 2.*x - 2.*y - 2.*z));
2056 
2057  case 1:
2058  return .5*x*(1. - y)*((1. - z)*(-2.) +
2059  (-1.)*(2.*x - 2.*y - 2.*z - 1.));
2060 
2061  case 2:
2062  return .5*x*y*((1. - z)*(-2.) +
2063  (-1.)*(2.*x + 2.*y - 2.*z - 3.));
2064 
2065  case 3:
2066  return .5*(1. - x)*y*((1. - z)*(-2.) +
2067  (-1.)*(2.*y - 2.*x - 2.*z - 1.));
2068 
2069  case 4:
2070  return .5*(1. - x)*(1. - y)*(z*(2.) +
2071  (1.)*(2.*z - 2.*x - 2.*y - 1.));
2072 
2073  case 5:
2074  return .5*x*(1. - y)*(z*(2.) +
2075  (1.)*(2.*x - 2.*y + 2.*z - 3.));
2076 
2077  case 6:
2078  return .5*x*y*(z*(2.) +
2079  (1.)*(2.*x + 2.*y + 2.*z - 5.));
2080 
2081  case 7:
2082  return .5*(1. - x)*y*(z*(2.) +
2083  (1.)*(2.*y - 2.*x + 2.*z - 3.));
2084 
2085  case 8:
2086  return 2.*x*(1. - x)*(1. - y)*(-1.);
2087 
2088  case 9:
2089  return 2.*x*y*(1. - y)*(-1.);
2090 
2091  case 10:
2092  return 2.*x*(1. - x)*y*(-1.);
2093 
2094  case 11:
2095  return 2.*(1. - x)*y*(1. - y)*(-1.);
2096 
2097  case 12:
2098  return 2.*(1. - x)*(1. - y)*(1. - 2.*z);
2099 
2100  case 13:
2101  return 2.*x*(1. - y)*(1. - 2.*z);
2102 
2103  case 14:
2104  return 2.*x*y*(1. - 2.*z);
2105 
2106  case 15:
2107  return 2.*(1. - x)*y*(1. - 2.*z);
2108 
2109  case 16:
2110  return 2.*x*(1. - x)*(1. - y);
2111 
2112  case 17:
2113  return 2.*x*y*(1. - y);
2114 
2115  case 18:
2116  return 2.*x*(1. - x)*y;
2117 
2118  case 19:
2119  return 2.*(1. - x)*y*(1. - y);
2120 
2121  default:
2122  libmesh_error_msg("Invalid i = " << i);
2123  }
2124 
2125  default:
2126  libmesh_error_msg("Invalid shape function derivative j = " << j);
2127  }
2128  }
2129 
2130  // triquadratic hexahedral shape functions
2131  case HEX8:
2132  libmesh_assert_msg(T == L2_LAGRANGE,
2133  "High order on first order elements only supported for L2 families");
2134  libmesh_fallthrough();
2135  case HEX27:
2136  {
2137  libmesh_assert_less (i, 27);
2138 
2139  // Compute hex shape functions as a tensor-product
2140  const Real xi = p(0);
2141  const Real eta = p(1);
2142  const Real zeta = p(2);
2143 
2144  // The only way to make any sense of this
2145  // is to look at the mgflo/mg2/mgf documentation
2146  // and make the cut-out cube!
2147  // 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
2148  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};
2149  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};
2150  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};
2151 
2152  switch(j)
2153  {
2154  case 0:
2155  return (fe_lagrange_1D_quadratic_shape_deriv(i0[i], 0, xi)*
2156  fe_lagrange_1D_quadratic_shape (i1[i], eta)*
2157  fe_lagrange_1D_quadratic_shape (i2[i], zeta));
2158 
2159  case 1:
2160  return (fe_lagrange_1D_quadratic_shape (i0[i], xi)*
2161  fe_lagrange_1D_quadratic_shape_deriv(i1[i], 0, eta)*
2162  fe_lagrange_1D_quadratic_shape (i2[i], zeta));
2163 
2164  case 2:
2165  return (fe_lagrange_1D_quadratic_shape (i0[i], xi)*
2166  fe_lagrange_1D_quadratic_shape (i1[i], eta)*
2167  fe_lagrange_1D_quadratic_shape_deriv(i2[i], 0, zeta));
2168 
2169  default:
2170  libmesh_error_msg("Invalid j = " << j);
2171  }
2172  }
2173 
2174  // quadratic tetrahedral shape functions
2175  case TET4:
2176  libmesh_assert_msg(T == L2_LAGRANGE,
2177  "High order on first order elements only supported for L2 families");
2178  libmesh_fallthrough();
2179  case TET10:
2180  case TET14:
2181  {
2182  libmesh_assert_less (i, 10);
2183 
2184  // Area coordinates, pg. 205, Vol. I, Carey, Oden, Becker FEM
2185  const Real zeta1 = p(0);
2186  const Real zeta2 = p(1);
2187  const Real zeta3 = p(2);
2188  const Real zeta0 = 1. - zeta1 - zeta2 - zeta3;
2189 
2190  const Real dzeta0dxi = -1.;
2191  const Real dzeta1dxi = 1.;
2192  const Real dzeta2dxi = 0.;
2193  const Real dzeta3dxi = 0.;
2194 
2195  const Real dzeta0deta = -1.;
2196  const Real dzeta1deta = 0.;
2197  const Real dzeta2deta = 1.;
2198  const Real dzeta3deta = 0.;
2199 
2200  const Real dzeta0dzeta = -1.;
2201  const Real dzeta1dzeta = 0.;
2202  const Real dzeta2dzeta = 0.;
2203  const Real dzeta3dzeta = 1.;
2204 
2205  switch (j)
2206  {
2207  // d()/dxi
2208  case 0:
2209  {
2210  switch(i)
2211  {
2212  case 0:
2213  return (4.*zeta0 - 1.)*dzeta0dxi;
2214 
2215  case 1:
2216  return (4.*zeta1 - 1.)*dzeta1dxi;
2217 
2218  case 2:
2219  return (4.*zeta2 - 1.)*dzeta2dxi;
2220 
2221  case 3:
2222  return (4.*zeta3 - 1.)*dzeta3dxi;
2223 
2224  case 4:
2225  return 4.*(zeta0*dzeta1dxi + dzeta0dxi*zeta1);
2226 
2227  case 5:
2228  return 4.*(zeta1*dzeta2dxi + dzeta1dxi*zeta2);
2229 
2230  case 6:
2231  return 4.*(zeta0*dzeta2dxi + dzeta0dxi*zeta2);
2232 
2233  case 7:
2234  return 4.*(zeta0*dzeta3dxi + dzeta0dxi*zeta3);
2235 
2236  case 8:
2237  return 4.*(zeta1*dzeta3dxi + dzeta1dxi*zeta3);
2238 
2239  case 9:
2240  return 4.*(zeta2*dzeta3dxi + dzeta2dxi*zeta3);
2241 
2242  default:
2243  libmesh_error_msg("Invalid i = " << i);
2244  }
2245  }
2246 
2247  // d()/deta
2248  case 1:
2249  {
2250  switch(i)
2251  {
2252  case 0:
2253  return (4.*zeta0 - 1.)*dzeta0deta;
2254 
2255  case 1:
2256  return (4.*zeta1 - 1.)*dzeta1deta;
2257 
2258  case 2:
2259  return (4.*zeta2 - 1.)*dzeta2deta;
2260 
2261  case 3:
2262  return (4.*zeta3 - 1.)*dzeta3deta;
2263 
2264  case 4:
2265  return 4.*(zeta0*dzeta1deta + dzeta0deta*zeta1);
2266 
2267  case 5:
2268  return 4.*(zeta1*dzeta2deta + dzeta1deta*zeta2);
2269 
2270  case 6:
2271  return 4.*(zeta0*dzeta2deta + dzeta0deta*zeta2);
2272 
2273  case 7:
2274  return 4.*(zeta0*dzeta3deta + dzeta0deta*zeta3);
2275 
2276  case 8:
2277  return 4.*(zeta1*dzeta3deta + dzeta1deta*zeta3);
2278 
2279  case 9:
2280  return 4.*(zeta2*dzeta3deta + dzeta2deta*zeta3);
2281 
2282  default:
2283  libmesh_error_msg("Invalid i = " << i);
2284  }
2285  }
2286 
2287  // d()/dzeta
2288  case 2:
2289  {
2290  switch(i)
2291  {
2292  case 0:
2293  return (4.*zeta0 - 1.)*dzeta0dzeta;
2294 
2295  case 1:
2296  return (4.*zeta1 - 1.)*dzeta1dzeta;
2297 
2298  case 2:
2299  return (4.*zeta2 - 1.)*dzeta2dzeta;
2300 
2301  case 3:
2302  return (4.*zeta3 - 1.)*dzeta3dzeta;
2303 
2304  case 4:
2305  return 4.*(zeta0*dzeta1dzeta + dzeta0dzeta*zeta1);
2306 
2307  case 5:
2308  return 4.*(zeta1*dzeta2dzeta + dzeta1dzeta*zeta2);
2309 
2310  case 6:
2311  return 4.*(zeta0*dzeta2dzeta + dzeta0dzeta*zeta2);
2312 
2313  case 7:
2314  return 4.*(zeta0*dzeta3dzeta + dzeta0dzeta*zeta3);
2315 
2316  case 8:
2317  return 4.*(zeta1*dzeta3dzeta + dzeta1dzeta*zeta3);
2318 
2319  case 9:
2320  return 4.*(zeta2*dzeta3dzeta + dzeta2dzeta*zeta3);
2321 
2322  default:
2323  libmesh_error_msg("Invalid i = " << i);
2324  }
2325  }
2326 
2327  default:
2328  libmesh_error_msg("Invalid j = " << j);
2329  }
2330  }
2331 
2332 
2333  // "serendipity" prism
2334  case PRISM15:
2335  {
2336  libmesh_assert_less (i, 15);
2337 
2338  const Real xi = p(0);
2339  const Real eta = p(1);
2340  const Real zeta = p(2);
2341 
2342  switch (j)
2343  {
2344  // d()/dxi
2345  case 0:
2346  {
2347  switch(i)
2348  {
2349  case 0:
2350  return (2.*xi + 2.*eta + 0.5*zeta - 1.)*(1. - zeta);
2351  case 1:
2352  return (2.*xi - 1. - 0.5*zeta)*(1. - zeta);
2353  case 2:
2354  return 0.;
2355  case 3:
2356  return (2.*xi + 2.*eta - 0.5*zeta - 1.)*(1. + zeta);
2357  case 4:
2358  return (2.*xi - 1. + 0.5*zeta)*(1. + zeta);
2359  case 5:
2360  return 0.;
2361  case 6:
2362  return (4.*xi + 2.*eta - 2.)*(zeta - 1.);
2363  case 7:
2364  return -2.*(zeta - 1.)*eta;
2365  case 8:
2366  return 2.*(zeta - 1.)*eta;
2367  case 9:
2368  return (zeta - 1.)*(1. + zeta);
2369  case 10:
2370  return (1. - zeta)*(1. + zeta);
2371  case 11:
2372  return 0.;
2373  case 12:
2374  return (-4.*xi - 2.*eta + 2.)*(1. + zeta);
2375  case 13:
2376  return 2.*(1. + zeta)*eta;
2377  case 14:
2378  return -2.*(1. + zeta)*eta;
2379  default:
2380  libmesh_error_msg("Invalid i = " << i);
2381  }
2382  }
2383 
2384  // d()/deta
2385  case 1:
2386  {
2387  switch(i)
2388  {
2389  case 0:
2390  return (2.*xi + 2.*eta + 0.5*zeta - 1.)*(1. - zeta);
2391  case 1:
2392  return 0.;
2393  case 2:
2394  return (2.*eta - 1. - 0.5*zeta)*(1. - zeta);
2395  case 3:
2396  return (2.*xi + 2.*eta - 0.5*zeta - 1.)*(1. + zeta);
2397  case 4:
2398  return 0.;
2399  case 5:
2400  return (2.*eta - 1. + 0.5*zeta)*(1. + zeta);
2401  case 6:
2402  return 2.*(zeta - 1.)*xi;
2403  case 7:
2404  return 2.*(1. - zeta)*xi;
2405  case 8:
2406  return (2.*xi + 4.*eta - 2.)*(zeta - 1.);
2407  case 9:
2408  return (zeta - 1.)*(1. + zeta);
2409  case 10:
2410  return 0.;
2411  case 11:
2412  return (1. - zeta)*(1. + zeta);
2413  case 12:
2414  return -2.*(1. + zeta)*xi;
2415  case 13:
2416  return 2.*(1. + zeta)*xi;
2417  case 14:
2418  return (-2.*xi - 4.*eta + 2.)*(1. + zeta);
2419 
2420  default:
2421  libmesh_error_msg("Invalid i = " << i);
2422  }
2423  }
2424 
2425  // d()/dzeta
2426  case 2:
2427  {
2428  switch(i)
2429  {
2430  case 0:
2431  return (-xi - eta - zeta + 0.5)*(xi + eta - 1.);
2432  case 1:
2433  return -0.5*xi*(2.*xi - 1. - 2.*zeta);
2434  case 2:
2435  return -0.5*eta*(2.*eta - 1. - 2.*zeta);
2436  case 3:
2437  return (xi + eta - zeta - 0.5)*(xi + eta - 1.);
2438  case 4:
2439  return 0.5*xi*(2.*xi - 1. + 2.*zeta);
2440  case 5:
2441  return 0.5*eta*(2.*eta - 1. + 2.*zeta);
2442  case 6:
2443  return 2.*xi*(xi + eta - 1.);
2444  case 7:
2445  return -2.*xi*eta;
2446  case 8:
2447  return 2.*eta*(xi + eta - 1.);
2448  case 9:
2449  return 2.*zeta*(xi + eta - 1.);
2450  case 10:
2451  return -2.*xi*zeta;
2452  case 11:
2453  return -2.*eta*zeta;
2454  case 12:
2455  return 2.*xi*(1. - xi - eta);
2456  case 13:
2457  return 2.*xi*eta;
2458  case 14:
2459  return 2.*eta*(1. - xi - eta);
2460 
2461  default:
2462  libmesh_error_msg("Invalid i = " << i);
2463  }
2464  }
2465 
2466  default:
2467  libmesh_error_msg("Invalid j = " << j);
2468  }
2469  }
2470 
2471 
2472 
2473  // quadratic prism shape functions
2474  case PRISM6:
2475  libmesh_assert_msg(T == L2_LAGRANGE,
2476  "High order on first order elements only supported for L2 families");
2477  libmesh_fallthrough();
2478  case PRISM18:
2479  case PRISM20:
2480  case PRISM21:
2481  {
2482  libmesh_assert_less (i, 18);
2483 
2484  // Compute prism shape functions as a tensor-product
2485  // of a triangle and an edge
2486 
2487  Point p2d(p(0),p(1));
2488  Real p1d = p(2);
2489 
2490  // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
2491  static const unsigned int i0[] = {0, 0, 0, 1, 1, 1, 0, 0, 0, 2, 2, 2, 1, 1, 1, 2, 2, 2};
2492  static const unsigned int i1[] = {0, 1, 2, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 3, 4, 5};
2493 
2494  switch (j)
2495  {
2496  // d()/dxi
2497  case 0:
2498  return (FE<2,LAGRANGE>::shape_deriv(TRI6, SECOND, i1[i], 0, p2d)*
2499  fe_lagrange_1D_quadratic_shape(i0[i], p1d));
2500 
2501  // d()/deta
2502  case 1:
2503  return (FE<2,LAGRANGE>::shape_deriv(TRI6, SECOND, i1[i], 1, p2d)*
2504  fe_lagrange_1D_quadratic_shape(i0[i], p1d));
2505 
2506  // d()/dzeta
2507  case 2:
2508  return (FE<2,LAGRANGE>::shape(TRI6, SECOND, i1[i], p2d)*
2509  fe_lagrange_1D_quadratic_shape_deriv(i0[i], 0, p1d));
2510 
2511  default:
2512  libmesh_error_msg("Invalid shape function derivative j = " << j);
2513  }
2514  }
2515 
2516  // G. Bedrosian, "Shape functions and integration formulas for
2517  // three-dimensional finite element analysis", Int. J. Numerical
2518  // Methods Engineering, vol 35, p. 95-108, 1992.
2519  case PYRAMID13:
2520  {
2521  libmesh_assert_less (i, 13);
2522 
2523  const Real xi = p(0);
2524  const Real eta = p(1);
2525  const Real zeta = p(2);
2526  const Real eps = 1.e-35;
2527 
2528  // Denominators are perturbed by epsilon to avoid
2529  // divide-by-zero issues.
2530  Real
2531  den = (-1. + zeta + eps),
2532  den2 = den*den,
2533  xi2 = xi*xi,
2534  eta2 = eta*eta,
2535  zeta2 = zeta*zeta,
2536  zeta3 = zeta2*zeta;
2537 
2538  switch (j)
2539  {
2540  // d/dxi
2541  case 0:
2542  switch(i)
2543  {
2544  case 0:
2545  return 0.25*(-zeta - eta + 2.*eta*zeta - 2.*xi + 2.*zeta*xi + 2.*eta*xi + zeta2 + eta2)/den;
2546 
2547  case 1:
2548  return -0.25*(-zeta - eta + 2.*eta*zeta + 2.*xi - 2.*zeta*xi - 2.*eta*xi + zeta2 + eta2)/den;
2549 
2550  case 2:
2551  return -0.25*(-zeta + eta - 2.*eta*zeta + 2.*xi - 2.*zeta*xi + 2.*eta*xi + zeta2 + eta2)/den;
2552 
2553  case 3:
2554  return 0.25*(-zeta + eta - 2.*eta*zeta - 2.*xi + 2.*zeta*xi - 2.*eta*xi + zeta2 + eta2)/den;
2555 
2556  case 4:
2557  return 0.;
2558 
2559  case 5:
2560  return -(-1. + eta + zeta)*xi/den;
2561 
2562  case 6:
2563  return 0.5*(-1. + eta + zeta)*(1. + eta - zeta)/den;
2564 
2565  case 7:
2566  return (1. + eta - zeta)*xi/den;
2567 
2568  case 8:
2569  return -0.5*(-1. + eta + zeta)*(1. + eta - zeta)/den;
2570 
2571  case 9:
2572  return -(-1. + eta + zeta)*zeta/den;
2573 
2574  case 10:
2575  return (-1. + eta + zeta)*zeta/den;
2576 
2577  case 11:
2578  return -(1. + eta - zeta)*zeta/den;
2579 
2580  case 12:
2581  return (1. + eta - zeta)*zeta/den;
2582 
2583  default:
2584  libmesh_error_msg("Invalid i = " << i);
2585  }
2586 
2587  // d/deta
2588  case 1:
2589  switch(i)
2590  {
2591  case 0:
2592  return 0.25*(-zeta - 2.*eta + 2.*eta*zeta - xi + 2.*zeta*xi + 2.*eta*xi + zeta2 + xi2)/den;
2593 
2594  case 1:
2595  return -0.25*(zeta + 2.*eta - 2.*eta*zeta - xi + 2.*zeta*xi + 2.*eta*xi - zeta2 - xi2)/den;
2596 
2597  case 2:
2598  return -0.25*(-zeta + 2.*eta - 2.*eta*zeta + xi - 2.*zeta*xi + 2.*eta*xi + zeta2 + xi2)/den;
2599 
2600  case 3:
2601  return 0.25*(zeta - 2.*eta + 2.*eta*zeta + xi - 2.*zeta*xi + 2.*eta*xi - zeta2 - xi2)/den;
2602 
2603  case 4:
2604  return 0.;
2605 
2606  case 5:
2607  return -0.5*(-1. + xi + zeta)*(1. + xi - zeta)/den;
2608 
2609  case 6:
2610  return (1. + xi - zeta)*eta/den;
2611 
2612  case 7:
2613  return 0.5*(-1. + xi + zeta)*(1. + xi - zeta)/den;
2614 
2615  case 8:
2616  return -(-1. + xi + zeta)*eta/den;
2617 
2618  case 9:
2619  return -(-1. + xi + zeta)*zeta/den;
2620 
2621  case 10:
2622  return (1. + xi - zeta)*zeta/den;
2623 
2624  case 11:
2625  return -(1. + xi - zeta)*zeta/den;
2626 
2627  case 12:
2628  return (-1. + xi + zeta)*zeta/den;
2629 
2630  default:
2631  libmesh_error_msg("Invalid i = " << i);
2632  }
2633 
2634  // d/dzeta
2635  case 2:
2636  {
2637  switch(i)
2638  {
2639  case 0:
2640  return -0.25*(xi + eta + 1.)*(-1. + 2.*zeta - zeta2 + eta*xi)/den2;
2641 
2642  case 1:
2643  return 0.25*(eta - xi + 1.)*(1. - 2.*zeta + zeta2 + eta*xi)/den2;
2644 
2645  case 2:
2646  return 0.25*(xi + eta - 1.)*(-1. + 2.*zeta - zeta2 + eta*xi)/den2;
2647 
2648  case 3:
2649  return -0.25*(eta - xi - 1.)*(1. - 2.*zeta + zeta2 + eta*xi)/den2;
2650 
2651  case 4:
2652  return 4.*zeta - 1.;
2653 
2654  case 5:
2655  return 0.5*(-2 + eta + 6.*zeta + eta*xi2 + eta*zeta2 - 6.*zeta2 + 2.*zeta3 - 2.*eta*zeta)/den2;
2656 
2657  case 6:
2658  return -0.5*(2 - 6.*zeta + xi + xi*zeta2 + eta2*xi + 6.*zeta2 - 2.*zeta3 - 2.*zeta*xi)/den2;
2659 
2660  case 7:
2661  return -0.5*(2 + eta - 6.*zeta + eta*xi2 + eta*zeta2 + 6.*zeta2 - 2.*zeta3 - 2.*eta*zeta)/den2;
2662 
2663  case 8:
2664  return 0.5*(-2 + 6.*zeta + xi + xi*zeta2 + eta2*xi - 6.*zeta2 + 2.*zeta3 - 2.*zeta*xi)/den2;
2665 
2666  case 9:
2667  return (1. - eta - 4.*zeta - xi - xi*zeta2 - eta*zeta2 + eta*xi + 5.*zeta2 - 2.*zeta3 + 2.*eta*zeta + 2.*zeta*xi)/den2;
2668 
2669  case 10:
2670  return -(-1. + eta + 4.*zeta - xi - xi*zeta2 + eta*zeta2 + eta*xi - 5.*zeta2 + 2.*zeta3 - 2.*eta*zeta + 2.*zeta*xi)/den2;
2671 
2672  case 11:
2673  return (1. + eta - 4.*zeta + xi + xi*zeta2 + eta*zeta2 + eta*xi + 5.*zeta2 - 2.*zeta3 - 2.*eta*zeta - 2.*zeta*xi)/den2;
2674 
2675  case 12:
2676  return -(-1. - eta + 4.*zeta + xi + xi*zeta2 - eta*zeta2 + eta*xi - 5.*zeta2 + 2.*zeta3 + 2.*eta*zeta - 2.*zeta*xi)/den2;
2677 
2678  default:
2679  libmesh_error_msg("Invalid i = " << i);
2680  }
2681  }
2682 
2683  default:
2684  libmesh_error_msg("Invalid j = " << j);
2685  }
2686  }
2687 
2688  // Quadratic shape functions, as defined in R. Graglia, "Higher order
2689  // bases on pyramidal elements", IEEE Trans Antennas and Propagation,
2690  // vol 47, no 5, May 1999.
2691  case PYRAMID5:
2692  libmesh_assert_msg(T == L2_LAGRANGE,
2693  "High order on first order elements only supported for L2 families");
2694  libmesh_fallthrough();
2695  case PYRAMID14:
2696  case PYRAMID18:
2697  {
2698  libmesh_assert_less (i, 14);
2699 
2700  const Real xi = p(0);
2701  const Real eta = p(1);
2702  const Real zeta = p(2);
2703  const Real eps = 1.e-35;
2704 
2705  // The "normalized coordinates" defined by Graglia. These are
2706  // the planes which define the faces of the pyramid.
2707  Real
2708  p1 = 0.5*(1. - eta - zeta), // back
2709  p2 = 0.5*(1. + xi - zeta), // left
2710  p3 = 0.5*(1. + eta - zeta), // front
2711  p4 = 0.5*(1. - xi - zeta); // right
2712 
2713  // Denominators are perturbed by epsilon to avoid
2714  // divide-by-zero issues.
2715  Real
2716  den = (-1. + zeta + eps),
2717  den2 = den*den,
2718  den3 = den2*den;
2719 
2720  switch (j)
2721  {
2722  // d/dxi
2723  case 0:
2724  switch(i)
2725  {
2726  case 0:
2727  return 0.5*p1*(-xi*eta + zeta - zeta*zeta + 2.*p4*eta)/den2;
2728 
2729  case 1:
2730  return -0.5*p1*(xi*eta + zeta - zeta*zeta + 2.*p2*eta)/den2;
2731 
2732  case 2:
2733  return 0.5*p3*(xi*eta - zeta + zeta*zeta + 2.*p2*eta)/den2;
2734 
2735  case 3:
2736  return -0.5*p3*(-xi*eta - zeta + zeta*zeta + 2.*p4*eta)/den2;
2737 
2738  case 4:
2739  return 0.;
2740 
2741  case 5:
2742  return 2.*p1*eta*xi/den2;
2743 
2744  case 6:
2745  return 2.*p1*p3*(xi + 2.*p2)/den2;
2746 
2747  case 7:
2748  return -2.*p3*eta*xi/den2;
2749 
2750  case 8:
2751  return -2.*p1*p3*(-xi + 2.*p4)/den2;
2752 
2753  case 9:
2754  return 2.*p1*zeta/den;
2755 
2756  case 10:
2757  return -2.*p1*zeta/den;
2758 
2759  case 11:
2760  return -2.*p3*zeta/den;
2761 
2762  case 12:
2763  return 2.*p3*zeta/den;
2764 
2765  case 13:
2766  return -8.*p1*p3*xi/den2;
2767 
2768  default:
2769  libmesh_error_msg("Invalid i = " << i);
2770  }
2771 
2772  // d/deta
2773  case 1:
2774  switch(i)
2775  {
2776  case 0:
2777  return -0.5*p4*(xi*eta - zeta + zeta*zeta - 2.*p1*xi)/den2;
2778 
2779  case 1:
2780  return 0.5*p2*(xi*eta + zeta - zeta*zeta - 2.*p1*xi)/den2;
2781 
2782  case 2:
2783  return 0.5*p2*(xi*eta - zeta + zeta*zeta + 2.*p3*xi)/den2;
2784 
2785  case 3:
2786  return -0.5*p4*(xi*eta + zeta - zeta*zeta + 2.*p3*xi)/den2;
2787 
2788  case 4:
2789  return 0.;
2790 
2791  case 5:
2792  return 2.*p2*p4*(eta - 2.*p1)/den2;
2793 
2794  case 6:
2795  return -2.*p2*xi*eta/den2;
2796 
2797  case 7:
2798  return 2.*p2*p4*(eta + 2.*p3)/den2;
2799 
2800  case 8:
2801  return 2.*p4*xi*eta/den2;
2802 
2803  case 9:
2804  return 2.*p4*zeta/den;
2805 
2806  case 10:
2807  return 2.*p2*zeta/den;
2808 
2809  case 11:
2810  return -2.*p2*zeta/den;
2811 
2812  case 12:
2813  return -2.*p4*zeta/den;
2814 
2815  case 13:
2816  return -8.*p2*p4*eta/den2;
2817 
2818  default:
2819  libmesh_error_msg("Invalid i = " << i);
2820  }
2821 
2822 
2823  // d/dzeta
2824  case 2:
2825  {
2826  switch(i)
2827  {
2828  case 0:
2829  return -0.5*p1*(xi*eta - zeta + zeta*zeta)/den2
2830  - 0.5*p4*(xi*eta - zeta + zeta*zeta)/den2
2831  + p4*p1*(2.*zeta - 1)/den2
2832  - 2.*p4*p1*(xi*eta - zeta + zeta*zeta)/den3;
2833 
2834  case 1:
2835  return 0.5*p2*(xi*eta + zeta - zeta*zeta)/den2
2836  + 0.5*p1*(xi*eta + zeta - zeta*zeta)/den2
2837  - p1*p2*(1 - 2.*zeta)/den2
2838  + 2.*p1*p2*(xi*eta + zeta - zeta*zeta)/den3;
2839 
2840  case 2:
2841  return -0.5*p3*(xi*eta - zeta + zeta*zeta)/den2
2842  - 0.5*p2*(xi*eta - zeta + zeta*zeta)/den2
2843  + p2*p3*(2.*zeta - 1)/den2
2844  - 2.*p2*p3*(xi*eta - zeta + zeta*zeta)/den3;
2845 
2846  case 3:
2847  return 0.5*p4*(xi*eta + zeta - zeta*zeta)/den2
2848  + 0.5*p3*(xi*eta + zeta - zeta*zeta)/den2
2849  - p3*p4*(1 - 2.*zeta)/den2
2850  + 2.*p3*p4*(xi*eta + zeta - zeta*zeta)/den3;
2851 
2852  case 4:
2853  return 4.*zeta - 1.;
2854 
2855  case 5:
2856  return 2.*p4*p1*eta/den2
2857  + 2.*p2*p4*eta/den2
2858  + 2.*p1*p2*eta/den2
2859  + 8.*p2*p1*p4*eta/den3;
2860 
2861  case 6:
2862  return -2.*p2*p3*xi/den2
2863  - 2.*p1*p3*xi/den2
2864  - 2.*p1*p2*xi/den2
2865  - 8.*p1*p2*p3*xi/den3;
2866 
2867  case 7:
2868  return -2.*p3*p4*eta/den2
2869  - 2.*p2*p4*eta/den2
2870  - 2.*p2*p3*eta/den2
2871  - 8.*p2*p3*p4*eta/den3;
2872 
2873  case 8:
2874  return 2.*p4*p1*xi/den2
2875  + 2.*p1*p3*xi/den2
2876  + 2.*p3*p4*xi/den2
2877  + 8.*p3*p4*p1*xi/den3;
2878 
2879  case 9:
2880  return 2.*p4*zeta/den
2881  + 2.*p1*zeta/den
2882  - 4.*p1*p4/den
2883  + 4.*p1*p4*zeta/den2;
2884 
2885  case 10:
2886  return 2.*p1*zeta/den
2887  + 2.*p2*zeta/den
2888  - 4.*p2*p1/den
2889  + 4.*p2*p1*zeta/den2;
2890 
2891  case 11:
2892  return 2.*p2*zeta/den
2893  + 2.*p3*zeta/den
2894  - 4.*p3*p2/den
2895  + 4.*p3*p2*zeta/den2;
2896 
2897  case 12:
2898  return 2.*p3*zeta/den
2899  + 2.*p4*zeta/den
2900  - 4.*p4*p3/den
2901  + 4.*p4*p3*zeta/den2;
2902 
2903  case 13:
2904  return -8.*p2*p3*p4/den2
2905  - 8.*p3*p4*p1/den2
2906  - 8.*p2*p1*p4/den2
2907  - 8.*p1*p2*p3/den2
2908  - 32.*p1*p2*p3*p4/den3;
2909 
2910  default:
2911  libmesh_error_msg("Invalid i = " << i);
2912  }
2913  }
2914 
2915  default:
2916  libmesh_error_msg("Invalid j = " << j);
2917  }
2918  }
2919 
2920 
2921  default:
2922  libmesh_error_msg("ERROR: Unsupported 3D element type!: " << Utility::enum_to_string(type));
2923  }
2924  }
2925 
2926  case THIRD:
2927  {
2928  switch (type)
2929  {
2930  // quadratic Lagrange shape functions with a cubic bubble
2931  case TET14:
2932  {
2933  libmesh_assert_less (i, 14);
2934 
2935  // Area coordinates, pg. 205, Vol. I, Carey, Oden, Becker FEM
2936  const Real zeta1 = p(0);
2937  const Real zeta2 = p(1);
2938  const Real zeta3 = p(2);
2939  const Real zeta0 = 1. - zeta1 - zeta2 - zeta3;
2940 
2941  const Real dzeta0dxi = -1.;
2942  const Real dzeta1dxi = 1.;
2943  const Real dzeta2dxi = 0.;
2944  const Real dzeta3dxi = 0.;
2945  const Real dbubble012dxi = (zeta0-zeta1)*zeta2;
2946  const Real dbubble013dxi = (zeta0-zeta1)*zeta3;
2947  const Real dbubble123dxi = zeta2*zeta3;
2948  const Real dbubble023dxi = -zeta2*zeta3;
2949 
2950  const Real dzeta0deta = -1.;
2951  const Real dzeta1deta = 0.;
2952  const Real dzeta2deta = 1.;
2953  const Real dzeta3deta = 0.;
2954  const Real dbubble012deta = (zeta0-zeta2)*zeta1;
2955  const Real dbubble013deta = -zeta1*zeta3;
2956  const Real dbubble123deta = zeta1*zeta3;
2957  const Real dbubble023deta = (zeta0-zeta2)*zeta3;
2958 
2959  const Real dzeta0dzeta = -1.;
2960  const Real dzeta1dzeta = 0.;
2961  const Real dzeta2dzeta = 0.;
2962  const Real dzeta3dzeta = 1.;
2963  const Real dbubble012dzeta = -zeta1*zeta2;
2964  const Real dbubble013dzeta = (zeta0-zeta3)*zeta1;
2965  const Real dbubble123dzeta = zeta1*zeta2;
2966  const Real dbubble023dzeta = (zeta0-zeta3)*zeta2;
2967 
2968  switch (j)
2969  {
2970  // d()/dxi
2971  case 0:
2972  {
2973  switch(i)
2974  {
2975  case 0:
2976  return (4.*zeta0 - 1.)*dzeta0dxi + 3.*(dbubble012dxi+dbubble013dxi+dbubble023dxi);
2977 
2978  case 1:
2979  return (4.*zeta1 - 1.)*dzeta1dxi + 3.*(dbubble012dxi+dbubble013dxi+dbubble123dxi);
2980 
2981  case 2:
2982  return (4.*zeta2 - 1.)*dzeta2dxi + 3.*(dbubble012dxi+dbubble023dxi+dbubble123dxi);
2983 
2984  case 3:
2985  return (4.*zeta3 - 1.)*dzeta3dxi + 3.*(dbubble013dxi+dbubble023dxi+dbubble123dxi);
2986 
2987  case 4:
2988  return 4.*(zeta0*dzeta1dxi + dzeta0dxi*zeta1) - 12.*(dbubble012dxi+dbubble013dxi);
2989 
2990  case 5:
2991  return 4.*(zeta1*dzeta2dxi + dzeta1dxi*zeta2) - 12.*(dbubble012dxi+dbubble123dxi);
2992 
2993  case 6:
2994  return 4.*(zeta0*dzeta2dxi + dzeta0dxi*zeta2) - 12.*(dbubble012dxi+dbubble023dxi);
2995 
2996  case 7:
2997  return 4.*(zeta0*dzeta3dxi + dzeta0dxi*zeta3) - 12.*(dbubble013dxi+dbubble023dxi);
2998 
2999  case 8:
3000  return 4.*(zeta1*dzeta3dxi + dzeta1dxi*zeta3) - 12.*(dbubble013dxi+dbubble123dxi);
3001 
3002  case 9:
3003  return 4.*(zeta2*dzeta3dxi + dzeta2dxi*zeta3) - 12.*(dbubble023dxi+dbubble123dxi);
3004 
3005  case 10:
3006  return 27.*dbubble012dxi;
3007 
3008  case 11:
3009  return 27.*dbubble013dxi;
3010 
3011  case 12:
3012  return 27.*dbubble123dxi;
3013 
3014  case 13:
3015  return 27.*dbubble023dxi;
3016 
3017  default:
3018  libmesh_error_msg("Invalid i = " << i);
3019  }
3020  }
3021 
3022  // d()/deta
3023  case 1:
3024  {
3025  switch(i)
3026  {
3027  case 0:
3028  return (4.*zeta0 - 1.)*dzeta0deta + 3.*(dbubble012deta+dbubble013deta+dbubble023deta);;
3029 
3030  case 1:
3031  return (4.*zeta1 - 1.)*dzeta1deta + 3.*(dbubble012deta+dbubble013deta+dbubble123deta);
3032 
3033  case 2:
3034  return (4.*zeta2 - 1.)*dzeta2deta + 3.*(dbubble012deta+dbubble023deta+dbubble123deta);
3035 
3036  case 3:
3037  return (4.*zeta3 - 1.)*dzeta3deta + 3.*(dbubble013deta+dbubble023deta+dbubble123deta);
3038 
3039  case 4:
3040  return 4.*(zeta0*dzeta1deta + dzeta0deta*zeta1) - 12.*(dbubble012deta+dbubble013deta);
3041 
3042  case 5:
3043  return 4.*(zeta1*dzeta2deta + dzeta1deta*zeta2) - 12.*(dbubble012deta+dbubble123deta);
3044 
3045  case 6:
3046  return 4.*(zeta0*dzeta2deta + dzeta0deta*zeta2) - 12.*(dbubble012deta+dbubble023deta);
3047 
3048  case 7:
3049  return 4.*(zeta0*dzeta3deta + dzeta0deta*zeta3) - 12.*(dbubble013deta+dbubble023deta);
3050 
3051  case 8:
3052  return 4.*(zeta1*dzeta3deta + dzeta1deta*zeta3) - 12.*(dbubble013deta+dbubble123deta);
3053 
3054  case 9:
3055  return 4.*(zeta2*dzeta3deta + dzeta2deta*zeta3) - 12.*(dbubble023deta+dbubble123deta);
3056 
3057  case 10:
3058  return 27.*dbubble012deta;
3059 
3060  case 11:
3061  return 27.*dbubble013deta;
3062 
3063  case 12:
3064  return 27.*dbubble123deta;
3065 
3066  case 13:
3067  return 27.*dbubble023deta;
3068 
3069  default:
3070  libmesh_error_msg("Invalid i = " << i);
3071  }
3072  }
3073 
3074  // d()/dzeta
3075  case 2:
3076  {
3077  switch(i)
3078  {
3079  case 0:
3080  return (4.*zeta0 - 1.)*dzeta0dzeta + 3.*(dbubble012dzeta+dbubble013dzeta+dbubble023dzeta);
3081 
3082  case 1:
3083  return (4.*zeta1 - 1.)*dzeta1dzeta + 3.*(dbubble012dzeta+dbubble013dzeta+dbubble123dzeta);
3084 
3085  case 2:
3086  return (4.*zeta2 - 1.)*dzeta2dzeta + 3.*(dbubble012dzeta+dbubble023dzeta+dbubble123dzeta);
3087 
3088  case 3:
3089  return (4.*zeta3 - 1.)*dzeta3dzeta + 3.*(dbubble013dzeta+dbubble023dzeta+dbubble123dzeta);
3090 
3091  case 4:
3092  return 4.*(zeta0*dzeta1dzeta + dzeta0dzeta*zeta1) - 12.*(dbubble012dzeta+dbubble013dzeta);
3093 
3094  case 5:
3095  return 4.*(zeta1*dzeta2dzeta + dzeta1dzeta*zeta2) - 12.*(dbubble012dzeta+dbubble123dzeta);
3096 
3097  case 6:
3098  return 4.*(zeta0*dzeta2dzeta + dzeta0dzeta*zeta2) - 12.*(dbubble012dzeta+dbubble023dzeta);
3099 
3100  case 7:
3101  return 4.*(zeta0*dzeta3dzeta + dzeta0dzeta*zeta3) - 12.*(dbubble013dzeta+dbubble023dzeta);
3102 
3103  case 8:
3104  return 4.*(zeta1*dzeta3dzeta + dzeta1dzeta*zeta3) - 12.*(dbubble013dzeta+dbubble123dzeta);
3105 
3106  case 9:
3107  return 4.*(zeta2*dzeta3dzeta + dzeta2dzeta*zeta3) - 12.*(dbubble023dzeta+dbubble123dzeta);
3108 
3109  case 10:
3110  return 27.*dbubble012dzeta;
3111 
3112  case 11:
3113  return 27.*dbubble013dzeta;
3114 
3115  case 12:
3116  return 27.*dbubble123dzeta;
3117 
3118  case 13:
3119  return 27.*dbubble023dzeta;
3120 
3121  default:
3122  libmesh_error_msg("Invalid i = " << i);
3123  }
3124  }
3125 
3126  default:
3127  libmesh_error_msg("Invalid j = " << j);
3128  }
3129  }
3130 
3131  case PRISM20:
3132  {
3133  libmesh_assert_less (i, 20);
3134 
3135  // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
3136  static const unsigned int i0[] = {0, 0, 0, 1, 1, 1, 0, 0, 0, 2, 2, 2, 1, 1, 1, 2, 2, 2, 0, 1, 2};
3137 
3138  Point p2d(p(0),p(1));
3139  Real p1d = p(2);
3140 
3141  const Real mainval = FE<3,LAGRANGE>::shape_deriv(PRISM21, THIRD, i, j, p);
3142 
3143  if (i0[i] != 2)
3144  return mainval;
3145 
3146  Real bubbleval = 0;
3147 
3148  switch (j)
3149  {
3150  // d()/dxi
3151  case 0:
3152  bubbleval =
3153  FE<2,LAGRANGE>::shape_deriv(TRI7, THIRD, 6, 0, p2d)*
3154  fe_lagrange_1D_quadratic_shape(2, p1d);
3155  break;
3156 
3157  // d()/deta
3158  case 1:
3159  bubbleval =
3160  FE<2,LAGRANGE>::shape_deriv(TRI7, THIRD, 6, 1, p2d)*
3161  fe_lagrange_1D_quadratic_shape(2, p1d);
3162  break;
3163 
3164  // d()/dzeta
3165  case 2:
3166  bubbleval =
3167  FE<2,LAGRANGE>::shape(TRI7, THIRD, 6, p2d)*
3168  fe_lagrange_1D_quadratic_shape_deriv(2, 0, p1d);
3169  break;
3170 
3171  default:
3172  libmesh_error_msg("Invalid shape function derivative j = " << j);
3173  }
3174 
3175  if (i < 12) // vertices
3176  return mainval - bubbleval / 9;
3177 
3178  return mainval + bubbleval * (Real(4) / 9);
3179  }
3180 
3181  case PRISM21:
3182  {
3183  libmesh_assert_less (i, 21);
3184 
3185  // Compute prism shape functions as a tensor-product
3186  // of a triangle and an edge
3187 
3188  Point p2d(p(0),p(1));
3189  Real p1d = p(2);
3190 
3191  // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
3192  static const unsigned int i0[] = {0, 0, 0, 1, 1, 1, 0, 0, 0, 2, 2, 2, 1, 1, 1, 2, 2, 2, 0, 1, 2};
3193  static const unsigned int i1[] = {0, 1, 2, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 3, 4, 5, 6, 6, 6};
3194 
3195  switch (j)
3196  {
3197  // d()/dxi
3198  case 0:
3199  return (FE<2,LAGRANGE>::shape_deriv(TRI7, THIRD, i1[i], 0, p2d)*
3200  fe_lagrange_1D_quadratic_shape(i0[i], p1d));
3201 
3202  // d()/deta
3203  case 1:
3204  return (FE<2,LAGRANGE>::shape_deriv(TRI7, THIRD, i1[i], 1, p2d)*
3205  fe_lagrange_1D_quadratic_shape(i0[i], p1d));
3206 
3207  // d()/dzeta
3208  case 2:
3209  return (FE<2,LAGRANGE>::shape(TRI7, THIRD, i1[i], p2d)*
3210  fe_lagrange_1D_quadratic_shape_deriv(i0[i], 0, p1d));
3211 
3212  default:
3213  libmesh_error_msg("Invalid shape function derivative j = " << j);
3214  }
3215  }
3216 
3217  case PYRAMID18:
3218  {
3219  return fe_fdm_deriv(type, order, elem, i, j, p, fe_lagrange_3D_shape<T>);
3220  }
3221 
3222  default:
3223  libmesh_error_msg("ERROR: Unsupported 3D element type!: " << Utility::enum_to_string(type));
3224  }
3225  }
3226 
3227  // unsupported order
3228  default:
3229  libmesh_error_msg("ERROR: Unsupported 3D FE order!: " << order);
3230  }
3231 
3232 #else // LIBMESH_DIM != 3
3233  libmesh_ignore(type, order, elem, i, j, p);
3234  libmesh_not_implemented();
3235 #endif
3236 }
3237 
3238 
3239 
3240 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
3241 
3242 template <FEFamily T>
3243 Real fe_lagrange_3D_shape_second_deriv(const ElemType type,
3244  const Order order,
3245  const Elem * elem,
3246  const unsigned int i,
3247  const unsigned int j,
3248  const Point & p)
3249 {
3250 #if LIBMESH_DIM == 3
3251 
3252  libmesh_assert_less (j, 6);
3253 
3254  switch (order)
3255  {
3256  // linear Lagrange shape functions
3257  case FIRST:
3258  {
3259  switch (type)
3260  {
3261  // Linear tets have all second derivatives = 0
3262  case TET4:
3263  case TET10:
3264  case TET14:
3265  {
3266  return 0.;
3267  }
3268 
3269  // The following elements use either tensor product or
3270  // rational basis functions, and therefore probably have
3271  // second derivatives, but we have not implemented them
3272  // yet...
3273  case PRISM6:
3274  case PRISM15:
3275  case PRISM18:
3276  case PRISM20:
3277  case PRISM21:
3278  {
3279  libmesh_assert_less (i, 6);
3280 
3281  // Compute prism shape functions as a tensor-product
3282  // of a triangle and an edge
3283 
3284  Point p2d(p(0),p(1));
3285  Real p1d = p(2);
3286 
3287  // 0 1 2 3 4 5
3288  static const unsigned int i0[] = {0, 0, 0, 1, 1, 1};
3289  static const unsigned int i1[] = {0, 1, 2, 0, 1, 2};
3290 
3291  switch (j)
3292  {
3293  // All repeated second derivatives and the xi-eta derivative are zero on PRISMs
3294  case 0: // d^2()/dxi^2
3295  case 1: // d^2()/dxideta
3296  case 2: // d^2()/deta^2
3297  case 5: // d^2()/dzeta^2
3298  {
3299  return 0.;
3300  }
3301 
3302  case 3: // d^2()/dxidzeta
3303  return (FE<2,LAGRANGE>::shape_deriv(TRI3, FIRST, i1[i], 0, p2d)*
3304  fe_lagrange_1D_linear_shape_deriv(i0[i], 0, p1d));
3305 
3306  case 4: // d^2()/detadzeta
3307  return (FE<2,LAGRANGE>::shape_deriv(TRI3, FIRST, i1[i], 1, p2d)*
3308  fe_lagrange_1D_linear_shape_deriv(i0[i], 0, p1d));
3309 
3310  default:
3311  libmesh_error_msg("Invalid j = " << j);
3312  }
3313  }
3314 
3315  case PYRAMID5:
3316  case PYRAMID13:
3317  case PYRAMID14:
3318  case PYRAMID18:
3319  {
3320  libmesh_assert_less (i, 5);
3321 
3322  const Real xi = p(0);
3323  const Real eta = p(1);
3324  const Real zeta = p(2);
3325  const Real eps = 1.e-35;
3326 
3327  switch (j)
3328  {
3329  // xi-xi and eta-eta derivatives are all zero for PYRAMID5.
3330  case 0: // d^2()/dxi^2
3331  case 2: // d^2()/deta^2
3332  return 0.;
3333 
3334  case 1: // d^2()/dxideta
3335  {
3336  switch (i)
3337  {
3338  case 0:
3339  case 2:
3340  return 0.25/(1. - zeta + eps);
3341  case 1:
3342  case 3:
3343  return -0.25/(1. - zeta + eps);
3344  case 4:
3345  return 0.;
3346  default:
3347  libmesh_error_msg("Invalid i = " << i);
3348  }
3349  }
3350 
3351  case 3: // d^2()/dxidzeta
3352  {
3353  Real den = (1. - zeta + eps)*(1. - zeta + eps);
3354 
3355  switch (i)
3356  {
3357  case 0:
3358  case 2:
3359  return 0.25*eta/den;
3360  case 1:
3361  case 3:
3362  return -0.25*eta/den;
3363  case 4:
3364  return 0.;
3365  default:
3366  libmesh_error_msg("Invalid i = " << i);
3367  }
3368  }
3369 
3370  case 4: // d^2()/detadzeta
3371  {
3372  Real den = (1. - zeta + eps)*(1. - zeta + eps);
3373 
3374  switch (i)
3375  {
3376  case 0:
3377  case 2:
3378  return 0.25*xi/den;
3379  case 1:
3380  case 3:
3381  return -0.25*xi/den;
3382  case 4:
3383  return 0.;
3384  default:
3385  libmesh_error_msg("Invalid i = " << i);
3386  }
3387  }
3388 
3389  case 5: // d^2()/dzeta^2
3390  {
3391  Real den = (1. - zeta + eps)*(1. - zeta + eps)*(1. - zeta + eps);
3392 
3393  switch (i)
3394  {
3395  case 0:
3396  case 2:
3397  return 0.5*xi*eta/den;
3398  case 1:
3399  case 3:
3400  return -0.5*xi*eta/den;
3401  case 4:
3402  return 0.;
3403  default:
3404  libmesh_error_msg("Invalid i = " << i);
3405  }
3406  }
3407 
3408  default:
3409  libmesh_error_msg("Invalid j = " << j);
3410  }
3411  }
3412 
3413  // Trilinear shape functions on HEX8s have nonzero mixed second derivatives
3414  case HEX8:
3415  case HEX20:
3416  case HEX27:
3417  {
3418  libmesh_assert_less (i, 8);
3419 
3420  // Compute hex shape functions as a tensor-product
3421  const Real xi = p(0);
3422  const Real eta = p(1);
3423  const Real zeta = p(2);
3424 
3425  static const unsigned int i0[] = {0, 1, 1, 0, 0, 1, 1, 0};
3426  static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 1, 1};
3427  static const unsigned int i2[] = {0, 0, 0, 0, 1, 1, 1, 1};
3428 
3429  switch (j)
3430  {
3431  // All repeated second derivatives are zero on HEX8
3432  case 0: // d^2()/dxi^2
3433  case 2: // d^2()/deta^2
3434  case 5: // d^2()/dzeta^2
3435  {
3436  return 0.;
3437  }
3438 
3439  case 1: // d^2()/dxideta
3440  return (fe_lagrange_1D_linear_shape_deriv(i0[i], 0, xi)*
3441  fe_lagrange_1D_linear_shape_deriv(i1[i], 0, eta)*
3442  fe_lagrange_1D_linear_shape (i2[i], zeta));
3443 
3444  case 3: // d^2()/dxidzeta
3445  return (fe_lagrange_1D_linear_shape_deriv(i0[i], 0, xi)*
3446  fe_lagrange_1D_linear_shape (i1[i], eta)*
3447  fe_lagrange_1D_linear_shape_deriv(i2[i], 0, zeta));
3448 
3449  case 4: // d^2()/detadzeta
3450  return (fe_lagrange_1D_linear_shape (i0[i], xi)*
3451  fe_lagrange_1D_linear_shape_deriv(i1[i], 0, eta)*
3452  fe_lagrange_1D_linear_shape_deriv(i2[i], 0, zeta));
3453 
3454  default:
3455  libmesh_error_msg("Invalid j = " << j);
3456  }
3457  }
3458 
3459  // All second derivatives for piecewise-linear polyhedra are
3460  // zero or dirac-type distributions, but we can't put the
3461  // latter in a Real, so beware when integrating...
3462  case C0POLYHEDRON:
3463  {
3464  return 0.;
3465  }
3466 
3467  default:
3468  libmesh_error_msg("ERROR: Unsupported 3D element type!: " << Utility::enum_to_string(type));
3469  }
3470 
3471  }
3472 
3473  // quadratic Lagrange shape functions
3474  case SECOND:
3475  {
3476  switch (type)
3477  {
3478 
3479  // serendipity hexahedral quadratic shape functions
3480  case HEX20:
3481  {
3482  libmesh_assert_less (i, 20);
3483 
3484  const Real xi = p(0);
3485  const Real eta = p(1);
3486  const Real zeta = p(2);
3487 
3488  // these functions are defined for (x,y,z) in [0,1]^3
3489  // so transform the locations
3490  const Real x = .5*(xi + 1.);
3491  const Real y = .5*(eta + 1.);
3492  const Real z = .5*(zeta + 1.);
3493 
3494  switch(j)
3495  {
3496  case 0: // d^2()/dxi^2
3497  {
3498  switch(i)
3499  {
3500  case 0:
3501  case 1:
3502  return (1. - y) * (1. - z);
3503  case 2:
3504  case 3:
3505  return y * (1. - z);
3506  case 4:
3507  case 5:
3508  return (1. - y) * z;
3509  case 6:
3510  case 7:
3511  return y * z;
3512  case 8:
3513  return -2. * (1. - y) * (1. - z);
3514  case 10:
3515  return -2. * y * (1. - z);
3516  case 16:
3517  return -2. * (1. - y) * z;
3518  case 18:
3519  return -2. * y * z;
3520  case 9:
3521  case 11:
3522  case 12:
3523  case 13:
3524  case 14:
3525  case 15:
3526  case 17:
3527  case 19:
3528  return 0;
3529  default:
3530  libmesh_error_msg("Invalid i = " << i);
3531  }
3532  }
3533  case 1: // d^2()/dxideta
3534  {
3535  switch(i)
3536  {
3537  case 0:
3538  return (1.25 - x - y - .5*z) * (1. - z);
3539  case 1:
3540  return (-x + y + .5*z - .25) * (1. - z);
3541  case 2:
3542  return (x + y - .5*z - .75) * (1. - z);
3543  case 3:
3544  return (-y + x + .5*z - .25) * (1. - z);
3545  case 4:
3546  return -.25*z * (4.*x + 4.*y - 2.*z - 3);
3547  case 5:
3548  return -.25*z * (-4.*y + 4.*x + 2.*z - 1.);
3549  case 6:
3550  return .25*z * (-5 + 4.*x + 4.*y + 2.*z);
3551  case 7:
3552  return .25*z * (4.*x - 4.*y - 2.*z + 1.);
3553  case 8:
3554  return (-1. + 2.*x) * (1. - z);
3555  case 9:
3556  return (1. - 2.*y) * (1. - z);
3557  case 10:
3558  return (1. - 2.*x) * (1. - z);
3559  case 11:
3560  return (-1. + 2.*y) * (1. - z);
3561  case 12:
3562  return z * (1. - z);
3563  case 13:
3564  return -z * (1. - z);
3565  case 14:
3566  return z * (1. - z);
3567  case 15:
3568  return -z * (1. - z);
3569  case 16:
3570  return (-1. + 2.*x) * z;
3571  case 17:
3572  return (1. - 2.*y) * z;
3573  case 18:
3574  return (1. - 2.*x) * z;
3575  case 19:
3576  return (-1. + 2.*y) * z;
3577  default:
3578  libmesh_error_msg("Invalid i = " << i);
3579  }
3580  }
3581  case 2: // d^2()/deta^2
3582  switch(i)
3583  {
3584  case 0:
3585  case 3:
3586  return (1. - x) * (1. - z);
3587  case 1:
3588  case 2:
3589  return x * (1. - z);
3590  case 4:
3591  case 7:
3592  return (1. - x) * z;
3593  case 5:
3594  case 6:
3595  return x * z;
3596  case 9:
3597  return -2. * x * (1. - z);
3598  case 11:
3599  return -2. * (1. - x) * (1. - z);
3600  case 17:
3601  return -2. * x * z;
3602  case 19:
3603  return -2. * (1. - x) * z;
3604  case 8:
3605  case 10:
3606  case 12:
3607  case 13:
3608  case 14:
3609  case 15:
3610  case 16:
3611  case 18:
3612  return 0.;
3613  default:
3614  libmesh_error_msg("Invalid i = " << i);
3615  }
3616  case 3: // d^2()/dxidzeta
3617  switch(i)
3618  {
3619  case 0:
3620  return (1.25 - x - .5*y - z) * (1. - y);
3621  case 1:
3622  return (-x + .5*y + z - .25) * (1. - y);
3623  case 2:
3624  return -.25*y * (2.*y + 4.*x - 4.*z - 1.);
3625  case 3:
3626  return -.25*y * (-2.*y + 4.*x + 4.*z - 3);
3627  case 4:
3628  return (-z + x + .5*y - .25) * (1. - y);
3629  case 5:
3630  return (x - .5*y + z - .75) * (1. - y);
3631  case 6:
3632  return .25*y * (2.*y + 4.*x + 4.*z - 5);
3633  case 7:
3634  return .25*y * (-2.*y + 4.*x - 4.*z + 1.);
3635  case 8:
3636  return (-1. + 2.*x) * (1. - y);
3637  case 9:
3638  return -y * (1. - y);
3639  case 10:
3640  return (-1. + 2.*x) * y;
3641  case 11:
3642  return y * (1. - y);
3643  case 12:
3644  return (-1. + 2.*z) * (1. - y);
3645  case 13:
3646  return (1. - 2.*z) * (1. - y);
3647  case 14:
3648  return (1. - 2.*z) * y;
3649  case 15:
3650  return (-1. + 2.*z) * y;
3651  case 16:
3652  return (1. - 2.*x) * (1. - y);
3653  case 17:
3654  return y * (1. - y);
3655  case 18:
3656  return (1. - 2.*x) * y;
3657  case 19:
3658  return -y * (1. - y);
3659  default:
3660  libmesh_error_msg("Invalid i = " << i);
3661  }
3662  case 4: // d^2()/detadzeta
3663  switch(i)
3664  {
3665  case 0:
3666  return (1.25 - .5*x - y - z) * (1. - x);
3667  case 1:
3668  return .25*x * (2.*x - 4.*y - 4.*z + 3.);
3669  case 2:
3670  return -.25*x * (2.*x + 4.*y - 4.*z - 1.);
3671  case 3:
3672  return (-y + .5*x + z - .25) * (1. - x);
3673  case 4:
3674  return (-z + .5*x + y - .25) * (1. - x);
3675  case 5:
3676  return -.25*x * (2.*x - 4.*y + 4.*z - 1.);
3677  case 6:
3678  return .25*x * (2.*x + 4.*y + 4.*z - 5);
3679  case 7:
3680  return (y - .5*x + z - .75) * (1. - x);
3681  case 8:
3682  return x * (1. - x);
3683  case 9:
3684  return (-1. + 2.*y) * x;
3685  case 10:
3686  return -x * (1. - x);
3687  case 11:
3688  return (-1. + 2.*y) * (1. - x);
3689  case 12:
3690  return (-1. + 2.*z) * (1. - x);
3691  case 13:
3692  return (-1. + 2.*z) * x;
3693  case 14:
3694  return (1. - 2.*z) * x;
3695  case 15:
3696  return (1. - 2.*z) * (1. - x);
3697  case 16:
3698  return -x * (1. - x);
3699  case 17:
3700  return (1. - 2.*y) * x;
3701  case 18:
3702  return x * (1. - x);
3703  case 19:
3704  return (1. - 2.*y) * (1. - x);
3705  default:
3706  libmesh_error_msg("Invalid i = " << i);
3707  }
3708  case 5: // d^2()/dzeta^2
3709  switch(i)
3710  {
3711  case 0:
3712  case 4:
3713  return (1. - x) * (1. - y);
3714  case 1:
3715  case 5:
3716  return x * (1. - y);
3717  case 2:
3718  case 6:
3719  return x * y;
3720  case 3:
3721  case 7:
3722  return (1. - x) * y;
3723  case 12:
3724  return -2. * (1. - x) * (1. - y);
3725  case 13:
3726  return -2. * x * (1. - y);
3727  case 14:
3728  return -2. * x * y;
3729  case 15:
3730  return -2. * (1. - x) * y;
3731  case 8:
3732  case 9:
3733  case 10:
3734  case 11:
3735  case 16:
3736  case 17:
3737  case 18:
3738  case 19:
3739  return 0.;
3740  default:
3741  libmesh_error_msg("Invalid i = " << i);
3742  }
3743  default:
3744  libmesh_error_msg("Invalid j = " << j);
3745  }
3746  }
3747 
3748  // triquadratic hexahedral shape functions
3749  case HEX8:
3750  libmesh_assert_msg(T == L2_LAGRANGE,
3751  "High order on first order elements only supported for L2 families");
3752  libmesh_fallthrough();
3753  case HEX27:
3754  {
3755  libmesh_assert_less (i, 27);
3756 
3757  // Compute hex shape functions as a tensor-product
3758  const Real xi = p(0);
3759  const Real eta = p(1);
3760  const Real zeta = p(2);
3761 
3762  // The only way to make any sense of this
3763  // is to look at the mgflo/mg2/mgf documentation
3764  // and make the cut-out cube!
3765  // 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
3766  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};
3767  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};
3768  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};
3769 
3770  switch(j)
3771  {
3772  // d^2()/dxi^2
3773  case 0:
3774  return (fe_lagrange_1D_quadratic_shape_second_deriv(i0[i], 0, xi)*
3775  fe_lagrange_1D_quadratic_shape (i1[i], eta)*
3776  fe_lagrange_1D_quadratic_shape (i2[i], zeta));
3777 
3778  // d^2()/dxideta
3779  case 1:
3780  return (fe_lagrange_1D_quadratic_shape_deriv(i0[i], 0, xi)*
3781  fe_lagrange_1D_quadratic_shape_deriv(i1[i], 0, eta)*
3782  fe_lagrange_1D_quadratic_shape (i2[i], zeta));
3783 
3784  // d^2()/deta^2
3785  case 2:
3786  return (fe_lagrange_1D_quadratic_shape (i0[i], xi)*
3787  fe_lagrange_1D_quadratic_shape_second_deriv(i1[i], 0, eta)*
3788  fe_lagrange_1D_quadratic_shape (i2[i], zeta));
3789 
3790  // d^2()/dxidzeta
3791  case 3:
3792  return (fe_lagrange_1D_quadratic_shape_deriv(i0[i], 0, xi)*
3793  fe_lagrange_1D_quadratic_shape (i1[i], eta)*
3794  fe_lagrange_1D_quadratic_shape_deriv(i2[i], 0, zeta));
3795 
3796  // d^2()/detadzeta
3797  case 4:
3798  return (fe_lagrange_1D_quadratic_shape (i0[i], xi)*
3799  fe_lagrange_1D_quadratic_shape_deriv(i1[i], 0, eta)*
3800  fe_lagrange_1D_quadratic_shape_deriv(i2[i], 0, zeta));
3801 
3802  // d^2()/dzeta^2
3803  case 5:
3804  return (fe_lagrange_1D_quadratic_shape (i0[i], xi)*
3805  fe_lagrange_1D_quadratic_shape (i1[i], eta)*
3806  fe_lagrange_1D_quadratic_shape_second_deriv(i2[i], 0, zeta));
3807 
3808  default:
3809  libmesh_error_msg("Invalid j = " << j);
3810  }
3811  }
3812 
3813  // quadratic tetrahedral shape functions
3814  case TET4:
3815  libmesh_assert_msg(T == L2_LAGRANGE,
3816  "High order on first order elements only supported for L2 families");
3817  libmesh_fallthrough();
3818  case TET10:
3819  case TET14:
3820  {
3821  // The area coordinates are the same as used for the
3822  // shape() and shape_deriv() functions.
3823  // const Real zeta0 = 1. - zeta1 - zeta2 - zeta3;
3824  // const Real zeta1 = p(0);
3825  // const Real zeta2 = p(1);
3826  // const Real zeta3 = p(2);
3827  static const Real dzetadxi[4][3] =
3828  {
3829  {-1., -1., -1.},
3830  {1., 0., 0.},
3831  {0., 1., 0.},
3832  {0., 0., 1.}
3833  };
3834 
3835  // Convert from j -> (j,k) indices for independent variable
3836  // (0=xi, 1=eta, 2=zeta)
3837  static const unsigned short int independent_var_indices[6][2] =
3838  {
3839  {0, 0}, // d^2 phi / dxi^2
3840  {0, 1}, // d^2 phi / dxi deta
3841  {1, 1}, // d^2 phi / deta^2
3842  {0, 2}, // d^2 phi / dxi dzeta
3843  {1, 2}, // d^2 phi / deta dzeta
3844  {2, 2} // d^2 phi / dzeta^2
3845  };
3846 
3847  // Convert from i -> zeta indices. Each quadratic shape
3848  // function for the Tet10 depends on up to two of the zeta
3849  // area coordinate functions (see the shape() function above).
3850  // This table just tells which two area coords it uses.
3851  static const unsigned short int zeta_indices[10][2] =
3852  {
3853  {0, 0},
3854  {1, 1},
3855  {2, 2},
3856  {3, 3},
3857  {0, 1},
3858  {1, 2},
3859  {2, 0},
3860  {0, 3},
3861  {1, 3},
3862  {2, 3},
3863  };
3864 
3865  // Look up the independent variable indices for this value of j.
3866  const unsigned int my_j = independent_var_indices[j][0];
3867  const unsigned int my_k = independent_var_indices[j][1];
3868 
3869  if (i<4)
3870  {
3871  return 4.*dzetadxi[i][my_j]*dzetadxi[i][my_k];
3872  }
3873 
3874  else if (i<10)
3875  {
3876  const unsigned short int my_m = zeta_indices[i][0];
3877  const unsigned short int my_n = zeta_indices[i][1];
3878 
3879  return 4.*(dzetadxi[my_n][my_j]*dzetadxi[my_m][my_k] +
3880  dzetadxi[my_m][my_j]*dzetadxi[my_n][my_k] );
3881  }
3882  else
3883  libmesh_error_msg("Invalid shape function index " << i);
3884  }
3885 
3886 
3887 
3888  // "serendipity" prism
3889  case PRISM15:
3890  {
3891  libmesh_assert_less (i, 15);
3892 
3893  const Real xi = p(0);
3894  const Real eta = p(1);
3895  const Real zeta = p(2);
3896 
3897  switch (j)
3898  {
3899  // d^2()/dxi^2
3900  case 0:
3901  {
3902  switch(i)
3903  {
3904  case 0:
3905  case 1:
3906  return 2.*(1. - zeta);
3907  case 2:
3908  case 5:
3909  case 7:
3910  case 8:
3911  case 9:
3912  case 10:
3913  case 11:
3914  case 13:
3915  case 14:
3916  return 0.;
3917  case 3:
3918  case 4:
3919  return 2.*(1. + zeta);
3920  case 6:
3921  return 4.*(zeta - 1);
3922  case 12:
3923  return -4.*(1. + zeta);
3924  default:
3925  libmesh_error_msg("Invalid i = " << i);
3926  }
3927  }
3928 
3929  // d^2()/dxideta
3930  case 1:
3931  {
3932  switch(i)
3933  {
3934  case 0:
3935  case 7:
3936  return 2.*(1. - zeta);
3937  case 1:
3938  case 2:
3939  case 4:
3940  case 5:
3941  case 9:
3942  case 10:
3943  case 11:
3944  return 0.;
3945  case 3:
3946  case 13:
3947  return 2.*(1. + zeta);
3948  case 6:
3949  case 8:
3950  return 2.*(zeta - 1.);
3951  case 12:
3952  case 14:
3953  return -2.*(1. + zeta);
3954  default:
3955  libmesh_error_msg("Invalid i = " << i);
3956  }
3957  }
3958 
3959  // d^2()/deta^2
3960  case 2:
3961  {
3962  switch(i)
3963  {
3964  case 0:
3965  case 2:
3966  return 2.*(1. - zeta);
3967  case 1:
3968  case 4:
3969  case 6:
3970  case 7:
3971  case 9:
3972  case 10:
3973  case 11:
3974  case 12:
3975  case 13:
3976  return 0.;
3977  case 3:
3978  case 5:
3979  return 2.*(1. + zeta);
3980  case 8:
3981  return 4.*(zeta - 1.);
3982  case 14:
3983  return -4.*(1. + zeta);
3984  default:
3985  libmesh_error_msg("Invalid i = " << i);
3986  }
3987  }
3988 
3989  // d^2()/dxidzeta
3990  case 3:
3991  {
3992  switch(i)
3993  {
3994  case 0:
3995  return 1.5 - zeta - 2.*xi - 2.*eta;
3996  case 1:
3997  return 0.5 + zeta - 2.*xi;
3998  case 2:
3999  case 5:
4000  case 11:
4001  return 0.;
4002  case 3:
4003  return -1.5 - zeta + 2.*xi + 2.*eta;
4004  case 4:
4005  return -0.5 + zeta + 2.*xi;
4006  case 6:
4007  return 4.*xi + 2.*eta - 2.;
4008  case 7:
4009  return -2.*eta;
4010  case 8:
4011  return 2.*eta;
4012  case 9:
4013  return 2.*zeta;
4014  case 10:
4015  return -2.*zeta;
4016  case 12:
4017  return -4.*xi - 2.*eta + 2.;
4018  case 13:
4019  return 2.*eta;
4020  case 14:
4021  return -2.*eta;
4022  default:
4023  libmesh_error_msg("Invalid i = " << i);
4024  }
4025  }
4026 
4027  // d^2()/detadzeta
4028  case 4:
4029  {
4030  switch(i)
4031  {
4032  case 0:
4033  return 1.5 - zeta - 2.*xi - 2.*eta;
4034  case 1:
4035  case 4:
4036  case 10:
4037  return 0.;
4038  case 2:
4039  return .5 + zeta - 2.*eta;
4040  case 3:
4041  return -1.5 - zeta + 2.*xi + 2.*eta;
4042  case 5:
4043  return -.5 + zeta + 2.*eta;
4044  case 6:
4045  return 2.*xi;
4046  case 7:
4047  return -2.*xi;
4048  case 8:
4049  return 2.*xi + 4.*eta - 2.;
4050  case 9:
4051  return 2.*zeta;
4052  case 11:
4053  return -2.*zeta;
4054  case 12:
4055  return -2.*xi;
4056  case 13:
4057  return 2.*xi;
4058  case 14:
4059  return -2.*xi - 4.*eta + 2.;
4060  default:
4061  libmesh_error_msg("Invalid i = " << i);
4062  }
4063  }
4064 
4065  // d^2()/dzeta^2
4066  case 5:
4067  {
4068  switch(i)
4069  {
4070  case 0:
4071  case 3:
4072  return 1. - xi - eta;
4073  case 1:
4074  case 4:
4075  return xi;
4076  case 2:
4077  case 5:
4078  return eta;
4079  case 6:
4080  case 7:
4081  case 8:
4082  case 12:
4083  case 13:
4084  case 14:
4085  return 0.;
4086  case 9:
4087  return 2.*xi + 2.*eta - 2.;
4088  case 10:
4089  return -2.*xi;
4090  case 11:
4091  return -2.*eta;
4092  default:
4093  libmesh_error_msg("Invalid i = " << i);
4094  }
4095  }
4096 
4097  default:
4098  libmesh_error_msg("Invalid j = " << j);
4099  }
4100  }
4101 
4102 
4103 
4104  // quadratic prism shape functions
4105  case PRISM6:
4106  libmesh_assert_msg(T == L2_LAGRANGE,
4107  "High order on first order elements only supported for L2 families");
4108  libmesh_fallthrough();
4109  case PRISM18:
4110  case PRISM20:
4111  case PRISM21:
4112  {
4113  libmesh_assert_less (i, 18);
4114 
4115  // Compute prism shape functions as a tensor-product
4116  // of a triangle and an edge
4117 
4118  Point p2d(p(0),p(1));
4119  Real p1d = p(2);
4120 
4121  // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
4122  static const unsigned int i0[] = {0, 0, 0, 1, 1, 1, 0, 0, 0, 2, 2, 2, 1, 1, 1, 2, 2, 2};
4123  static const unsigned int i1[] = {0, 1, 2, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 3, 4, 5};
4124 
4125  switch (j)
4126  {
4127  // d^2()/dxi^2
4128  case 0:
4129  return (FE<2,LAGRANGE>::shape_second_deriv(TRI6, SECOND, i1[i], 0, p2d)*
4130  fe_lagrange_1D_quadratic_shape(i0[i], p1d));
4131 
4132  // d^2()/dxideta
4133  case 1:
4134  return (FE<2,LAGRANGE>::shape_second_deriv(TRI6, SECOND, i1[i], 1, p2d)*
4135  fe_lagrange_1D_quadratic_shape(i0[i], p1d));
4136 
4137  // d^2()/deta^2
4138  case 2:
4139  return (FE<2,LAGRANGE>::shape_second_deriv(TRI6, SECOND, i1[i], 2, p2d)*
4140  fe_lagrange_1D_quadratic_shape(i0[i], p1d));
4141 
4142  // d^2()/dxidzeta
4143  case 3:
4144  return (FE<2,LAGRANGE>::shape_deriv(TRI6, SECOND, i1[i], 0, p2d)*
4145  fe_lagrange_1D_quadratic_shape_deriv(i0[i], 0, p1d));
4146 
4147  // d^2()/detadzeta
4148  case 4:
4149  return (FE<2,LAGRANGE>::shape_deriv(TRI6, SECOND, i1[i], 1, p2d)*
4150  fe_lagrange_1D_quadratic_shape_deriv(i0[i], 0, p1d));
4151 
4152  // d^2()/dzeta^2
4153  case 5:
4154  return (FE<2,LAGRANGE>::shape(TRI6, SECOND, i1[i], p2d)*
4155  fe_lagrange_1D_quadratic_shape_second_deriv(i0[i], 0, p1d));
4156 
4157  default:
4158  libmesh_error_msg("Invalid shape function derivative j = " << j);
4159  }
4160  }
4161 
4162 
4163  // Quadratic shape functions, as defined in R. Graglia, "Higher order
4164  // bases on pyramidal elements", IEEE Trans Antennas and Propagation,
4165  // vol 47, no 5, May 1999.
4166  case PYRAMID5:
4167  libmesh_assert_msg(T == L2_LAGRANGE,
4168  "High order on first order elements only supported for L2 families");
4169  libmesh_fallthrough();
4170  case PYRAMID14:
4171  case PYRAMID18:
4172  {
4173  libmesh_assert_less (i, 14);
4174 
4175  const Real xi = p(0);
4176  const Real eta = p(1);
4177  const Real zeta = p(2);
4178  const Real eps = 1.e-35;
4179 
4180  // The "normalized coordinates" defined by Graglia. These are
4181  // the planes which define the faces of the pyramid.
4182  Real
4183  p1 = 0.5*(1. - eta - zeta), // back
4184  p2 = 0.5*(1. + xi - zeta), // left
4185  p3 = 0.5*(1. + eta - zeta), // front
4186  p4 = 0.5*(1. - xi - zeta); // right
4187 
4188  // Denominators are perturbed by epsilon to avoid
4189  // divide-by-zero issues.
4190  Real
4191  den = (-1. + zeta + eps),
4192  den2 = den*den,
4193  den3 = den2*den,
4194  den4 = den2*den2;
4195 
4196  // These terms are used in several of the derivatives
4197  Real
4198  numer_mp = xi*eta - zeta + zeta*zeta,
4199  numer_pm = xi*eta + zeta - zeta*zeta;
4200 
4201  switch (j)
4202  {
4203  case 0: // d^2()/dxi^2
4204  {
4205  switch(i)
4206  {
4207  case 0:
4208  case 1:
4209  return -p1*eta/den2;
4210  case 2:
4211  case 3:
4212  return p3*eta/den2;
4213  case 4:
4214  case 9:
4215  case 10:
4216  case 11:
4217  case 12:
4218  return 0.;
4219  case 5:
4220  return 2.*p1*eta/den2;
4221  case 6:
4222  case 8:
4223  return 4.*p1*p3/den2;
4224  case 7:
4225  return -2.*p3*eta/den2;
4226  case 13:
4227  return -8.*p1*p3/den2;
4228  default:
4229  libmesh_error_msg("Invalid i = " << i);
4230  }
4231  }
4232 
4233  case 1: // d^2()/dxideta
4234  {
4235  switch(i)
4236  {
4237  case 0:
4238  return 0.25*numer_mp/den2
4239  - 0.5*p1*xi/den2
4240  - 0.5*p4*eta/den2
4241  + p4*p1/den2;
4242 
4243  case 1:
4244  return 0.25*numer_pm/den2
4245  - 0.5*p1*xi/den2
4246  + 0.5*p2*eta/den2
4247  - p1*p2/den2;
4248 
4249  case 2:
4250  return 0.25*numer_mp/den2
4251  + 0.5*p3*xi/den2
4252  + 0.5*p2*eta/den2
4253  + p2*p3/den2;
4254 
4255  case 3:
4256  return 0.25*numer_pm/den2
4257  + 0.5*p3*xi/den2
4258  - 0.5*p4*eta/den2
4259  - p3*p4/den2;
4260 
4261  case 4:
4262  return 0.;
4263 
4264  case 5:
4265  return p4*eta/den2
4266  - 2.*p4*p1/den2
4267  - p2*eta/den2
4268  + 2.*p1*p2/den2;
4269 
4270  case 6:
4271  return -p3*xi/den2
4272  + p1*xi/den2
4273  - 2.*p2*p3/den2
4274  + 2.*p1*p2/den2;
4275 
4276  case 7:
4277  return p4*eta/den2
4278  + 2.*p3*p4/den2
4279  - p2*eta/den2
4280  - 2.*p2*p3/den2;
4281 
4282  case 8:
4283  return -p3*xi/den2
4284  + p1*xi/den2
4285  - 2.*p4*p1/den2
4286  + 2.*p3*p4/den2;
4287 
4288  case 9:
4289  case 11:
4290  return -zeta/den;
4291 
4292  case 10:
4293  case 12:
4294  return zeta/den;
4295 
4296  case 13:
4297  return 4.*p4*p1/den2
4298  - 4.*p3*p4/den2
4299  + 4.*p2*p3/den2
4300  - 4.*p1*p2/den2;
4301 
4302  default:
4303  libmesh_error_msg("Invalid i = " << i);
4304  }
4305  }
4306 
4307 
4308  case 2: // d^2()/deta^2
4309  {
4310  switch(i)
4311  {
4312  case 0:
4313  case 3:
4314  return -p4*xi/den2;
4315  case 1:
4316  case 2:
4317  return p2*xi/den2;
4318  case 4:
4319  case 9:
4320  case 10:
4321  case 11:
4322  case 12:
4323  return 0.;
4324  case 5:
4325  case 7:
4326  return 4.*p2*p4/den2;
4327  case 6:
4328  return -2.*p2*xi/den2;
4329  case 8:
4330  return 2.*p4*xi/den2;
4331  case 13:
4332  return -8.*p2*p4/den2;
4333  default:
4334  libmesh_error_msg("Invalid i = " << i);
4335  }
4336  }
4337 
4338 
4339  case 3: // d^2()/dxidzeta
4340  {
4341  switch(i)
4342  {
4343  case 0:
4344  return 0.25*numer_mp/den2
4345  - 0.5*p1*(2.*zeta - 1.)/den2
4346  + p1*numer_mp/den3
4347  - 0.5*p1*eta/den2
4348  - 0.5*p4*eta/den2
4349  - 2.*p4*p1*eta/den3;
4350 
4351  case 1:
4352  return 0.25*numer_pm/den2
4353  - 0.5*p1*(1 - 2.*zeta)/den2
4354  + p1*numer_pm/den3
4355  + 0.5*p2*eta/den2
4356  + 0.5*p1*eta/den2
4357  + 2.*p1*p2*eta/den3;
4358 
4359  case 2:
4360  return -0.25*numer_mp/den2
4361  + 0.5*p3*(2.*zeta - 1.)/den2
4362  - p3*numer_mp/den3
4363  - 0.5*p3*eta/den2
4364  - 0.5*p2*eta/den2
4365  - 2.*p2*p3*eta/den3;
4366 
4367  case 3:
4368  return -0.25*numer_pm/den2
4369  + 0.5*p3*(1 - 2.*zeta)/den2
4370  - p3*numer_pm/den3
4371  + 0.5*p4*eta/den2
4372  + 0.5*p3*eta/den2
4373  + 2.*p3*p4*eta/den3;
4374 
4375  case 4:
4376  return 0.;
4377 
4378  case 5:
4379  return p4*eta/den2
4380  + 4.*p4*p1*eta/den3
4381  - p2*eta/den2
4382  - 4.*p1*p2*eta/den3;
4383 
4384  case 6:
4385  return -p3*xi/den2
4386  - p1*xi/den2
4387  - 4.*p1*p3*xi/den3
4388  - 2.*p2*p3/den2
4389  - 2.*p1*p3/den2
4390  - 2.*p1*p2/den2
4391  - 8.*p1*p2*p3/den3;
4392 
4393  case 7:
4394  return -p4*eta/den2
4395  - 4.*p3*p4*eta/den3
4396  + p2*eta/den2
4397  + 4.*p2*p3*eta/den3;
4398 
4399  case 8:
4400  return -p3*xi/den2
4401  - p1*xi/den2
4402  - 4.*p1*p3*xi/den3
4403  + 2.*p4*p1/den2
4404  + 2.*p1*p3/den2
4405  + 2.*p3*p4/den2
4406  + 8.*p3*p4*p1/den3;
4407 
4408  case 9:
4409  return -zeta/den
4410  + 2.*p1/den
4411  - 2.*p1*zeta/den2;
4412 
4413  case 10:
4414  return zeta/den
4415  - 2.*p1/den
4416  + 2.*p1*zeta/den2;
4417 
4418  case 11:
4419  return zeta/den
4420  - 2.*p3/den
4421  + 2.*p3*zeta/den2;
4422 
4423  case 12:
4424  return -zeta/den
4425  + 2.*p3/den
4426  - 2.*p3*zeta/den2;
4427 
4428  case 13:
4429  return -4.*p4*p1/den2
4430  - 4.*p3*p4/den2
4431  - 16.*p3*p4*p1/den3
4432  + 4.*p2*p3/den2
4433  + 4.*p1*p2/den2
4434  + 16.*p1*p2*p3/den3;
4435 
4436  default:
4437  libmesh_error_msg("Invalid i = " << i);
4438  }
4439  }
4440 
4441  case 4: // d^2()/detadzeta
4442  {
4443  switch(i)
4444  {
4445  case 0:
4446  return 0.25*numer_mp/den2
4447  - 0.5*p4*(2.*zeta - 1.)/den2
4448  + p4*numer_mp/den3
4449  - 0.5*p1*xi/den2
4450  - 0.5*p4*xi/den2
4451  - 2.*p4*p1*xi/den3;
4452 
4453  case 1:
4454  return -0.25*numer_pm/den2
4455  + 0.5*p2*(1. - 2.*zeta)/den2
4456  - p2*numer_pm/den3
4457  + 0.5*p2*xi/den2
4458  + 0.5*p1*xi/den2
4459  + 2.*p1*p2*xi/den3;
4460 
4461  case 2:
4462  return -0.25*numer_mp/den2
4463  + 0.5*p2*(2.*zeta - 1.)/den2
4464  - p2*numer_mp/den3
4465  - 0.5*p3*xi/den2
4466  - 0.5*p2*xi/den2
4467  - 2.*p2*p3*xi/den3;
4468 
4469  case 3:
4470  return 0.25*numer_pm/den2
4471  - 0.5*p4*(1. - 2.*zeta)/den2
4472  + p4*numer_pm/den3
4473  + 0.5*p4*xi/den2
4474  + 0.5*p3*xi/den2
4475  + 2.*p3*p4*xi/den3;
4476 
4477  case 4:
4478  return 0.;
4479 
4480  case 5:
4481  return -p4*eta/den2
4482  - p2*eta/den2
4483  - 4.*p2*p4*eta/den3
4484  + 2.*p4*p1/den2
4485  + 2.*p2*p4/den2
4486  + 2.*p1*p2/den2
4487  + 8.*p2*p1*p4/den3;
4488 
4489  case 6:
4490  return p3*xi/den2
4491  + 4.*p2*p3*xi/den3
4492  - p1*xi/den2
4493  - 4.*p1*p2*xi/den3;
4494 
4495  case 7:
4496  return -p4*eta/den2
4497  - p2*eta/den2
4498  - 4.*p2*p4*eta/den3
4499  - 2.*p3*p4/den2
4500  - 2.*p2*p4/den2
4501  - 2.*p2*p3/den2
4502  - 8.*p2*p3*p4/den3;
4503 
4504  case 8:
4505  return p1*xi/den2
4506  + 4.*p4*p1*xi/den3
4507  - p3*xi/den2
4508  - 4.*p3*p4*xi/den3;
4509 
4510  case 9:
4511  return -zeta/den
4512  + 2.*p4/den
4513  - 2.*p4*zeta/den2;
4514 
4515  case 10:
4516  return -zeta/den
4517  + 2.*p2/den
4518  - 2.*p2*zeta/den2;
4519 
4520  case 11:
4521  return zeta/den
4522  - 2.*p2/den
4523  + 2.*p2*zeta/den2;
4524 
4525  case 12:
4526  return zeta/den
4527  - 2.*p4/den
4528  + 2.*p4*zeta/den2;
4529 
4530  case 13:
4531  return 4.*p3*p4/den2
4532  + 4.*p2*p3/den2
4533  + 16.*p2*p3*p4/den3
4534  - 4.*p4*p1/den2
4535  - 4.*p1*p2/den2
4536  - 16.*p2*p1*p4/den3;
4537 
4538  default:
4539  libmesh_error_msg("Invalid i = " << i);
4540  }
4541  }
4542 
4543  case 5: // d^2()/dzeta^2
4544  {
4545  switch(i)
4546  {
4547  case 0:
4548  return 0.5*numer_mp/den2
4549  - p1*(2.*zeta - 1.)/den2
4550  + 2.*p1*numer_mp/den3
4551  - p4*(2.*zeta - 1.)/den2
4552  + 2.*p4*numer_mp/den3
4553  + 2.*p4*p1/den2
4554  - 4.*p4*p1*(2.*zeta - 1.)/den3
4555  + 6.*p4*p1*numer_mp/den4;
4556 
4557  case 1:
4558  return -0.5*numer_pm/den2
4559  + p2*(1 - 2.*zeta)/den2
4560  - 2.*p2*numer_pm/den3
4561  + p1*(1 - 2.*zeta)/den2
4562  - 2.*p1*numer_pm/den3
4563  + 2.*p1*p2/den2
4564  + 4.*p1*p2*(1 - 2.*zeta)/den3
4565  - 6.*p1*p2*numer_pm/den4;
4566 
4567  case 2:
4568  return 0.5*numer_mp/den2
4569  - p3*(2.*zeta - 1.)/den2
4570  + 2.*p3*numer_mp/den3
4571  - p2*(2.*zeta - 1.)/den2
4572  + 2.*p2*numer_mp/den3
4573  + 2.*p2*p3/den2
4574  - 4.*p2*p3*(2.*zeta - 1.)/den3
4575  + 6.*p2*p3*numer_mp/den4;
4576 
4577  case 3:
4578  return -0.5*numer_pm/den2
4579  + p4*(1 - 2.*zeta)/den2
4580  - 2.*p4*numer_pm/den3
4581  + p3*(1 - 2.*zeta)/den2
4582  - 2.*p3*numer_pm/den3
4583  + 2.*p3*p4/den2
4584  + 4.*p3*p4*(1 - 2.*zeta)/den3
4585  - 6.*p3*p4*numer_pm/den4;
4586 
4587  case 4:
4588  return 4.;
4589 
4590  case 5:
4591  return -2.*p1*eta/den2
4592  - 2.*p4*eta/den2
4593  - 8.*p4*p1*eta/den3
4594  - 2.*p2*eta/den2
4595  - 8.*p2*p4*eta/den3
4596  - 8.*p1*p2*eta/den3
4597  - 24.*p2*p1*p4*eta/den4;
4598 
4599  case 6:
4600  return 2.*p3*xi/den2
4601  + 2.*p2*xi/den2
4602  + 8.*p2*p3*xi/den3
4603  + 2.*p1*xi/den2
4604  + 8.*p1*p3*xi/den3
4605  + 8.*p1*p2*xi/den3
4606  + 24.*p1*p2*p3*xi/den4;
4607 
4608  case 7:
4609  return 2.*p4*eta/den2
4610  + 2.*p3*eta/den2
4611  + 8.*p3*p4*eta/den3
4612  + 2.*p2*eta/den2
4613  + 8.*p2*p4*eta/den3
4614  + 8.*p2*p3*eta/den3
4615  + 24.*p2*p3*p4*eta/den4;
4616 
4617  case 8:
4618  return -2.*p1*xi/den2
4619  - 2.*p4*xi/den2
4620  - 8.*p4*p1*xi/den3
4621  - 2.*p3*xi/den2
4622  - 8.*p1*p3*xi/den3
4623  - 8.*p3*p4*xi/den3
4624  - 24.*p3*p4*p1*xi/den4;
4625 
4626  case 9:
4627  return -2.*zeta/den
4628  + 4.*p4/den
4629  - 4.*p4*zeta/den2
4630  + 4.*p1/den
4631  - 4.*p1*zeta/den2
4632  + 8.*p4*p1/den2
4633  - 8.*p1*p4*zeta/den3;
4634 
4635  case 10:
4636  return -2.*zeta/den
4637  + 4.*p1/den
4638  - 4.*p1*zeta/den2
4639  + 4.*p2/den
4640  - 4.*p2*zeta/den2
4641  + 8.*p1*p2/den2
4642  - 8.*p2*p1*zeta/den3;
4643 
4644  case 11:
4645  return -2.*zeta/den
4646  + 4.*p2/den
4647  - 4.*p2*zeta/den2
4648  + 4.*p3/den
4649  - 4.*p3*zeta/den2
4650  + 8.*p2*p3/den2
4651  - 8.*p3*p2*zeta/den3;
4652 
4653  case 12:
4654  return -2.*zeta/den
4655  + 4.*p3/den
4656  - 4.*p3*zeta/den2
4657  + 4.*p4/den
4658  - 4.*p4*zeta/den2
4659  + 8.*p3*p4/den2
4660  - 8.*p4*p3*zeta/den3;
4661 
4662  case 13:
4663  return 8.*p3*p4/den2
4664  + 8.*p2*p4/den2
4665  + 8.*p2*p3/den2
4666  + 32.*p2*p3*p4/den3
4667  + 8.*p4*p1/den2
4668  + 8.*p1*p3/den2
4669  + 32.*p3*p4*p1/den3
4670  + 8.*p1*p2/den2
4671  + 32.*p2*p1*p4/den3
4672  + 32.*p1*p2*p3/den3
4673  + 96.*p1*p2*p3*p4/den4;
4674 
4675  default:
4676  libmesh_error_msg("Invalid i = " << i);
4677  }
4678  }
4679 
4680  default:
4681  libmesh_error_msg("Invalid j = " << j);
4682  }
4683  }
4684 
4685  // G. Bedrosian, "Shape functions and integration formulas for
4686  // three-dimensional finite element analysis", Int. J. Numerical
4687  // Methods Engineering, vol 35, p. 95-108, 1992.
4688  case PYRAMID13:
4689  {
4690  libmesh_assert_less (i, 13);
4691 
4692  const Real xi = p(0);
4693  const Real eta = p(1);
4694  const Real zeta = p(2);
4695  const Real eps = 1.e-35;
4696 
4697  // Denominators are perturbed by epsilon to avoid
4698  // divide-by-zero issues.
4699  Real
4700  den = (-1. + zeta + eps),
4701  den2 = den*den,
4702  den3 = den2*den,
4703  xi2 = xi*xi,
4704  eta2 = eta*eta,
4705  zeta2 = zeta*zeta,
4706  zeta3 = zeta2*zeta;
4707 
4708  switch (j)
4709  {
4710  case 0: // d^2()/dxi^2
4711  {
4712  switch(i)
4713  {
4714  case 0:
4715  case 1:
4716  return 0.5*(-1. + zeta + eta)/den;
4717 
4718  case 2:
4719  case 3:
4720  return 0.5*(-1. + zeta - eta)/den;
4721 
4722  case 4:
4723  case 6:
4724  case 8:
4725  case 9:
4726  case 10:
4727  case 11:
4728  case 12:
4729  return 0.;
4730 
4731  case 5:
4732  return (1. - eta - zeta)/den;
4733 
4734  case 7:
4735  return (1. + eta - zeta)/den;
4736 
4737  default:
4738  libmesh_error_msg("Invalid i = " << i);
4739  }
4740  }
4741 
4742  case 1: // d^2()/dxideta
4743  {
4744  switch(i)
4745  {
4746  case 0:
4747  return 0.25*(-1. + 2.*zeta + 2.*xi + 2.*eta)/den;
4748 
4749  case 1:
4750  return -0.25*(-1. + 2.*zeta - 2.*xi + 2.*eta)/den;
4751 
4752  case 2:
4753  return -0.25*(1. - 2.*zeta + 2.*xi + 2.*eta)/den;
4754 
4755  case 3:
4756  return 0.25*(1. - 2.*zeta - 2.*xi + 2.*eta)/den;
4757 
4758  case 4:
4759  return 0.;
4760 
4761  case 5:
4762  return -xi/den;
4763 
4764  case 6:
4765  return eta/den;
4766 
4767  case 7:
4768  return xi/den;
4769 
4770  case 8:
4771  return -eta/den;
4772 
4773  case 9:
4774  return -zeta/den;
4775 
4776  case 10:
4777  return zeta/den;
4778 
4779  case 11:
4780  return -zeta/den;
4781 
4782  case 12:
4783  return zeta/den;
4784 
4785  default:
4786  libmesh_error_msg("Invalid i = " << i);
4787  }
4788  }
4789 
4790 
4791  case 2: // d^2()/deta^2
4792  {
4793  switch(i)
4794  {
4795  case 0:
4796  case 3:
4797  return 0.5*(-1. + zeta + xi)/den;
4798 
4799  case 1:
4800  case 2:
4801  return 0.5*(-1. + zeta - xi)/den;
4802 
4803  case 4:
4804  case 5:
4805  case 7:
4806  case 9:
4807  case 10:
4808  case 11:
4809  case 12:
4810  return 0.;
4811 
4812  case 6:
4813  return (1. + xi - zeta)/den;
4814 
4815  case 8:
4816  return (1. - xi - zeta)/den;
4817 
4818  default:
4819  libmesh_error_msg("Invalid i = " << i);
4820  }
4821  }
4822 
4823 
4824  case 3: // d^2()/dxidzeta
4825  {
4826  switch(i)
4827  {
4828  case 0:
4829  return -0.25*(-1. + 2.*zeta - zeta2 + eta + 2.*eta*xi + eta2)/den2;
4830 
4831  case 1:
4832  return 0.25*(-1. + 2.*zeta - zeta2 + eta - 2.*eta*xi + eta2)/den2;
4833 
4834  case 2:
4835  return 0.25*(-1. + 2.*zeta - zeta2 - eta + 2.*eta*xi + eta2)/den2;
4836 
4837  case 3:
4838  return -0.25*(-1. + 2.*zeta - zeta2 - eta - 2.*eta*xi + eta2)/den2;
4839 
4840  case 4:
4841  return 0.;
4842 
4843  case 5:
4844  return eta*xi/den2;
4845 
4846  case 6:
4847  return -0.5*(1. + zeta2 + eta2 - 2.*zeta)/den2;
4848 
4849  case 7:
4850  return -eta*xi/den2;
4851 
4852  case 8:
4853  return 0.5*(1. + zeta2 + eta2 - 2.*zeta)/den2;
4854 
4855  case 9:
4856  return (-1. - zeta2 + eta + 2.*zeta)/den2;
4857 
4858  case 10:
4859  return -(-1. - zeta2 + eta + 2.*zeta)/den2;
4860 
4861  case 11:
4862  return (1. + zeta2 + eta - 2.*zeta)/den2;
4863 
4864  case 12:
4865  return -(1. + zeta2 + eta - 2.*zeta)/den2;
4866 
4867  default:
4868  libmesh_error_msg("Invalid i = " << i);
4869  }
4870  }
4871 
4872  case 4: // d^2()/detadzeta
4873  {
4874  switch(i)
4875  {
4876  case 0:
4877  return -0.25*(-1. + 2.*zeta - zeta2 + xi + 2.*eta*xi + xi2)/den2;
4878 
4879  case 1:
4880  return 0.25*(1. - 2.*zeta + zeta2 + xi + 2.*eta*xi - xi2)/den2;
4881 
4882  case 2:
4883  return 0.25*(-1. + 2.*zeta - zeta2 - xi + 2.*eta*xi + xi2)/den2;
4884 
4885  case 3:
4886  return -0.25*(1. - 2.*zeta + zeta2 - xi + 2.*eta*xi - xi2)/den2;
4887 
4888  case 4:
4889  return 0.;
4890 
4891  case 5:
4892  return 0.5*(1. + xi2 + zeta2 - 2.*zeta)/den2;
4893 
4894  case 6:
4895  return -eta*xi/den2;
4896 
4897  case 7:
4898  return -0.5*(1. + xi2 + zeta2 - 2.*zeta)/den2;
4899 
4900  case 8:
4901  return eta*xi/den2;
4902 
4903  case 9:
4904  return (-1. - zeta2 + xi + 2.*zeta)/den2;
4905 
4906  case 10:
4907  return -(1. + zeta2 + xi - 2.*zeta)/den2;
4908 
4909  case 11:
4910  return (1. + zeta2 + xi - 2.*zeta)/den2;
4911 
4912  case 12:
4913  return -(-1. - zeta2 + xi + 2.*zeta)/den2;
4914 
4915  default:
4916  libmesh_error_msg("Invalid i = " << i);
4917  }
4918  }
4919 
4920  case 5: // d^2()/dzeta^2
4921  {
4922  switch(i)
4923  {
4924  case 0:
4925  return 0.5*(xi + eta + 1.)*eta*xi/den3;
4926 
4927  case 1:
4928  return -0.5*(eta - xi + 1.)*eta*xi/den3;
4929 
4930  case 2:
4931  return -0.5*(xi + eta - 1.)*eta*xi/den3;
4932 
4933  case 3:
4934  return 0.5*(eta - xi - 1.)*eta*xi/den3;
4935 
4936  case 4:
4937  return 4.;
4938 
4939  case 5:
4940  return -(1. - 3.*zeta + 3.*zeta2 - zeta3 + eta*xi2)/den3;
4941 
4942  case 6:
4943  return (-1. + 3.*zeta - 3.*zeta2 + zeta3 + eta2*xi)/den3;
4944 
4945  case 7:
4946  return (-1. + 3.*zeta - 3.*zeta2 + zeta3 + eta*xi2)/den3;
4947 
4948  case 8:
4949  return -(1. - 3.*zeta + 3.*zeta2 - zeta3 + eta2*xi)/den3;
4950 
4951  case 9:
4952  return -2.*(-1. + 3.*zeta - 3.*zeta2 + zeta3 + eta*xi)/den3;
4953 
4954  case 10:
4955  return 2.*(1. - 3.*zeta + 3.*zeta2 - zeta3 + eta*xi)/den3;
4956 
4957  case 11:
4958  return -2.*(-1. + 3.*zeta - 3.*zeta2 + zeta3 + eta*xi)/den3;
4959 
4960  case 12:
4961  return 2.*(1. - 3.*zeta + 3.*zeta2 - zeta3 + eta*xi)/den3;
4962 
4963  default:
4964  libmesh_error_msg("Invalid i = " << i);
4965  }
4966  }
4967 
4968  default:
4969  libmesh_error_msg("Invalid j = " << j);
4970  }
4971  }
4972 
4973  default:
4974  libmesh_error_msg("ERROR: Unsupported 3D element type!: " << Utility::enum_to_string(type));
4975  }
4976  }
4977 
4978  case THIRD:
4979  {
4980  switch (type)
4981  {
4982  // quadratic Lagrange shape functions with a cubic bubble
4983  case TET14:
4984  {
4985  libmesh_assert_less (i, 14);
4986 
4987  // The area coordinates are the same as used for the
4988  // shape() and shape_deriv() functions.
4989  // const Real zeta0 = 1. - zeta1 - zeta2 - zeta3;
4990  // const Real zeta1 = p(0);
4991  // const Real zeta2 = p(1);
4992  // const Real zeta3 = p(2);
4993  static const Real dzetadxi[4][3] =
4994  {
4995  {-1., -1., -1.},
4996  {1., 0., 0.},
4997  {0., 1., 0.},
4998  {0., 0., 1.}
4999  };
5000 
5001  // Convert from j -> (j,k) indices for independent variable
5002  // (0=xi, 1=eta, 2=zeta)
5003  static const unsigned short int independent_var_indices[6][2] =
5004  {
5005  {0, 0}, // d^2 phi / dxi^2
5006  {0, 1}, // d^2 phi / dxi deta
5007  {1, 1}, // d^2 phi / deta^2
5008  {0, 2}, // d^2 phi / dxi dzeta
5009  {1, 2}, // d^2 phi / deta dzeta
5010  {2, 2} // d^2 phi / dzeta^2
5011  };
5012 
5013  // Convert from i -> zeta indices. Each quadratic shape
5014  // function for the Tet10 depends on up to two of the zeta
5015  // area coordinate functions (see the shape() function above).
5016  // This table just tells which two area coords it uses.
5017  static const unsigned short int zeta_indices[10][2] =
5018  {
5019  {0, 0},
5020  {1, 1},
5021  {2, 2},
5022  {3, 3},
5023  {0, 1},
5024  {1, 2},
5025  {2, 0},
5026  {0, 3},
5027  {1, 3},
5028  {2, 3},
5029  };
5030 
5031  // Look up the independent variable indices for this value of j.
5032  const unsigned int my_j = independent_var_indices[j][0];
5033  const unsigned int my_k = independent_var_indices[j][1];
5034 
5035  Real returnval = 0;
5036  if (i<4)
5037  returnval = 4.*dzetadxi[i][my_j]*dzetadxi[i][my_k];
5038 
5039  else if (i<10)
5040  {
5041  const unsigned short int my_m = zeta_indices[i][0];
5042  const unsigned short int my_n = zeta_indices[i][1];
5043 
5044  returnval =
5045  4.*(dzetadxi[my_n][my_j]*dzetadxi[my_m][my_k] +
5046  dzetadxi[my_m][my_j]*dzetadxi[my_n][my_k] );
5047  }
5048 
5049  const Real zeta1 = p(0);
5050  const Real zeta2 = p(1);
5051  const Real zeta3 = p(2);
5052  const Real zeta0 = 1. - zeta1 - zeta2 - zeta3;
5053 
5054  // Fill these with whichever derivative we're concerned
5055  // with
5056  Real d2bubble012, d2bubble013, d2bubble023, d2bubble123;
5057  switch (j)
5058  {
5059  // d^2()/dxi^2
5060  case 0:
5061  {
5062  d2bubble012 = -2.*zeta2;
5063  d2bubble013 = -2.*zeta3;
5064  d2bubble023 = 0.;
5065  d2bubble123 = 0.;
5066  break;
5067  }
5068 
5069  // d^2()/dxideta
5070  case 1:
5071  {
5072  d2bubble012 = (zeta0-zeta1)-zeta2;
5073  d2bubble013 = -zeta3;
5074  d2bubble123 = zeta3;
5075  d2bubble023 = -zeta3;
5076  break;
5077  }
5078 
5079  // d^2()/deta^2
5080  case 2:
5081  {
5082  d2bubble012 = -2.*zeta1;
5083  d2bubble013 = 0.;
5084  d2bubble123 = 0.;
5085  d2bubble023 = -2.*zeta3;
5086  break;
5087  }
5088 
5089  // d^2()/dxi dzeta
5090  case 3:
5091  {
5092  d2bubble012 = -zeta2;
5093  d2bubble013 = (zeta0-zeta3)-zeta1;
5094  d2bubble123 = zeta2;
5095  d2bubble023 = -zeta2;
5096  break;
5097  }
5098 
5099  // d^2()/deta dzeta
5100  case 4:
5101  {
5102  d2bubble012 = -zeta1;
5103  d2bubble013 = -zeta1;
5104  d2bubble123 = zeta1;
5105  d2bubble023 = (zeta0-zeta3)-zeta2;
5106  break;
5107  }
5108 
5109  // d^2()/dzeta^2
5110  case 5:
5111  {
5112  d2bubble012 = 0.;
5113  d2bubble013 = -2.*zeta1;
5114  d2bubble123 = 0.;
5115  d2bubble023 = -2.*zeta2;
5116  break;
5117  }
5118 
5119  default:
5120  libmesh_error_msg("Invalid j = " << j);
5121  }
5122 
5123  switch (i)
5124  {
5125  case 0:
5126  return returnval + 3.*(d2bubble012+d2bubble013+d2bubble023);
5127 
5128  case 1:
5129  return returnval + 3.*(d2bubble012+d2bubble013+d2bubble123);
5130 
5131  case 2:
5132  return returnval + 3.*(d2bubble012+d2bubble023+d2bubble123);
5133 
5134  case 3:
5135  return returnval + 3.*(d2bubble013+d2bubble023+d2bubble123);
5136 
5137  case 4:
5138  return returnval - 12.*(d2bubble012+d2bubble013);
5139 
5140  case 5:
5141  return returnval - 12.*(d2bubble012+d2bubble123);
5142 
5143  case 6:
5144  return returnval - 12.*(d2bubble012+d2bubble023);
5145 
5146  case 7:
5147  return returnval - 12.*(d2bubble013+d2bubble023);
5148 
5149  case 8:
5150  return returnval - 12.*(d2bubble013+d2bubble123);
5151 
5152  case 9:
5153  return returnval - 12.*(d2bubble023+d2bubble123);
5154 
5155  case 10:
5156  return 27.*d2bubble012;
5157 
5158  case 11:
5159  return 27.*d2bubble013;
5160 
5161  case 12:
5162  return 27.*d2bubble123;
5163 
5164  case 13:
5165  return 27.*d2bubble023;
5166 
5167  default:
5168  libmesh_error_msg("Invalid i = " << i);
5169  }
5170  }
5171 
5172  case PRISM20:
5173  {
5174  libmesh_assert_less (i, 20);
5175 
5176  // Compute prism shape functions as a tensor-product
5177  // of a triangle and an edge
5178 
5179  Point p2d(p(0),p(1));
5180  Real p1d = p(2);
5181 
5182  // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
5183  static const unsigned int i0[] = {0, 0, 0, 1, 1, 1, 0, 0, 0, 2, 2, 2, 1, 1, 1, 2, 2, 2, 0, 1};
5184 
5185  const Real mainval = FE<3,LAGRANGE>::shape_second_deriv(PRISM21, THIRD, i, j, p);
5186 
5187  if (i0[i] != 2)
5188  return mainval;
5189 
5190  Real bubbleval = 0;
5191 
5192  switch (j)
5193  {
5194  // d^2()/dxi^2
5195  case 0:
5196  bubbleval =
5197  FE<2,LAGRANGE>::shape_second_deriv(TRI7, THIRD, 6, 0, p2d)*
5198  fe_lagrange_1D_quadratic_shape(2, p1d);
5199  break;
5200 
5201  // d^2()/dxideta
5202  case 1:
5203  bubbleval =
5204  FE<2,LAGRANGE>::shape_second_deriv(TRI7, THIRD, 6, 1, p2d)*
5205  fe_lagrange_1D_quadratic_shape(2, p1d);
5206  break;
5207 
5208  // d^2()/deta^2
5209  case 2:
5210  bubbleval =
5211  FE<2,LAGRANGE>::shape_second_deriv(TRI7, THIRD, 6, 2, p2d)*
5212  fe_lagrange_1D_quadratic_shape(2, p1d);
5213  break;
5214 
5215  // d^2()/dxidzeta
5216  case 3:
5217  bubbleval =
5218  FE<2,LAGRANGE>::shape_deriv(TRI7, THIRD, 6, 0, p2d)*
5219  fe_lagrange_1D_quadratic_shape_deriv(2, 0, p1d);
5220  break;
5221 
5222  // d^2()/detadzeta
5223  case 4:
5224  bubbleval =
5225  FE<2,LAGRANGE>::shape_deriv(TRI7, THIRD, 6, 1, p2d)*
5226  fe_lagrange_1D_quadratic_shape_deriv(2, 0, p1d);
5227  break;
5228 
5229  // d^2()/dzeta^2
5230  case 5:
5231  bubbleval =
5232  FE<2,LAGRANGE>::shape(TRI7, THIRD, 6, p2d)*
5233  fe_lagrange_1D_quadratic_shape_second_deriv(2, 0, p1d);
5234  break;
5235 
5236  default:
5237  libmesh_error_msg("Invalid shape function derivative j = " << j);
5238  }
5239 
5240  if (i < 12) // vertices
5241  return mainval - bubbleval / 9;
5242 
5243  return mainval + bubbleval * (Real(4) / 9);
5244  }
5245 
5246  case PRISM21:
5247  {
5248  libmesh_assert_less (i, 21);
5249 
5250  // Compute prism shape functions as a tensor-product
5251  // of a triangle and an edge
5252 
5253  Point p2d(p(0),p(1));
5254  Real p1d = p(2);
5255 
5256  // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
5257  static const unsigned int i0[] = {0, 0, 0, 1, 1, 1, 0, 0, 0, 2, 2, 2, 1, 1, 1, 2, 2, 2, 0, 1, 2};
5258  static const unsigned int i1[] = {0, 1, 2, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 3, 4, 5, 6, 6, 6};
5259 
5260  switch (j)
5261  {
5262  // d^2()/dxi^2
5263  case 0:
5264  return (FE<2,LAGRANGE>::shape_second_deriv(TRI7, THIRD, i1[i], 0, p2d)*
5265  fe_lagrange_1D_quadratic_shape(i0[i], p1d));
5266 
5267  // d^2()/dxideta
5268  case 1:
5269  return (FE<2,LAGRANGE>::shape_second_deriv(TRI7, THIRD, i1[i], 1, p2d)*
5270  fe_lagrange_1D_quadratic_shape(i0[i], p1d));
5271 
5272  // d^2()/deta^2
5273  case 2:
5274  return (FE<2,LAGRANGE>::shape_second_deriv(TRI7, THIRD, i1[i], 2, p2d)*
5275  fe_lagrange_1D_quadratic_shape(i0[i], p1d));
5276 
5277  // d^2()/dxidzeta
5278  case 3:
5279  return (FE<2,LAGRANGE>::shape_deriv(TRI7, THIRD, i1[i], 0, p2d)*
5280  fe_lagrange_1D_quadratic_shape_deriv(i0[i], 0, p1d));
5281 
5282  // d^2()/detadzeta
5283  case 4:
5284  return (FE<2,LAGRANGE>::shape_deriv(TRI7, THIRD, i1[i], 1, p2d)*
5285  fe_lagrange_1D_quadratic_shape_deriv(i0[i], 0, p1d));
5286 
5287  // d^2()/dzeta^2
5288  case 5:
5289  return (FE<2,LAGRANGE>::shape(TRI7, THIRD, i1[i], p2d)*
5290  fe_lagrange_1D_quadratic_shape_second_deriv(i0[i], 0, p1d));
5291 
5292  default:
5293  libmesh_error_msg("Invalid shape function derivative j = " << j);
5294  }
5295  }
5296 
5297  case PYRAMID18:
5298  {
5299  return fe_fdm_second_deriv(type, order, elem, i, j, p,
5300  fe_lagrange_3D_shape_deriv<T>);
5301  }
5302 
5303  default:
5304  libmesh_error_msg("ERROR: Unsupported 3D element type!: " << Utility::enum_to_string(type));
5305  }
5306  }
5307 
5308  // unsupported order
5309  default:
5310  libmesh_error_msg("ERROR: Unsupported 3D FE order!: " << order);
5311  }
5312 
5313 #else // LIBMESH_DIM != 3
5314  libmesh_ignore(type, order, elem, i, j, p);
5315  libmesh_not_implemented();
5316 #endif
5317 }
5318 
5319 #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
5320 
5321 
5322 } // anonymous namespace
libMesh::fe_lagrange_1D_quadratic_shape_second_deriv
Real fe_lagrange_1D_quadratic_shape_second_deriv(const unsigned int i, const unsigned int libmesh_dbg_var(j), const Real)
Definition: fe_lagrange_shape_1D.h:228
libMesh::FEType
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:196
libMesh::fe_lagrange_1D_linear_shape_deriv
Real fe_lagrange_1D_linear_shape_deriv(const unsigned int i, const unsigned int libmesh_dbg_var(j), const Real)
Definition: fe_lagrange_shape_1D.h:124
libMesh::ElemType
ElemType
Defines an enum for geometric element types.
Definition: enum_elem_type.h:33
libMesh::Order
Order
defines an enum for polynomial orders.
Definition: enum_order.h:40
libMesh::fe_lagrange_1D_quadratic_shape_deriv
Real fe_lagrange_1D_quadratic_shape_deriv(const unsigned int i, const unsigned int libmesh_dbg_var(j), const Real xi)
Definition: fe_lagrange_shape_1D.h:146
libMesh::invalid_uint
const unsigned int invalid_uint
A number which is used quite often to represent an invalid or uninitialized value for an unsigned int...
Definition: libmesh.h:310
libMesh::FE::shape
static OutputShape shape(const ElemType t, const Order o, const unsigned int i, const Point &p)
libMesh::fe_lagrange_1D_quadratic_shape
Real fe_lagrange_1D_quadratic_shape(const unsigned int i, const Real xi)
Definition: fe_lagrange_shape_1D.h:51
libMesh::FE::all_shape_derivs
static void all_shape_derivs(const Elem *elem, const Order o, const std::vector< Point > &p, std::vector< std::vector< OutputShape >> *comps[3], const bool add_p_level=true)
Fills comps with dphidxi (and in higher dimensions, eta/zeta) derivative component values for all sha...
libMesh::FE::default_all_shape_derivs
static void default_all_shape_derivs(const Elem *elem, const Order o, const std::vector< Point > &p, std::vector< std::vector< OutputShape >> *comps[3], const bool add_p_level=true)
A default implementation for all_shape_derivs.
Definition: fe.C:736
libMesh::Elem
This is the base class from which all geometric element types are derived.
Definition: elem.h:94
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:3109
libMesh::FEType::order
OrderWrapper order
The approximation order of the element.
Definition: fe_type.h:215
libMesh
The libMesh namespace provides an interface to certain functionality in the library.
libMesh::FE::default_all_shapes
static void default_all_shapes(const Elem *elem, const Order o, const std::vector< Point > &p, std::vector< std::vector< OutputShape >> &v, const bool add_p_level=true)
A default implementation for all_shapes.
Definition: fe.h:753
libMesh::fe_fdm_deriv
OutputShape fe_fdm_deriv(const Elem *elem, const Order order, const unsigned int i, const unsigned int j, const Point &p, const bool add_p_level, OutputShape(*shape_func)(const Elem *, const Order, const unsigned int, const Point &, const bool))
Helper functions for finite differenced derivatives in cases where analytical calculations haven&#39;t be...
Definition: fe.C:911
libMesh::LIBMESH_DEFAULT_VECTORIZED_FE
LIBMESH_DEFAULT_VECTORIZED_FE(template<>Real FE< 0, BERNSTEIN)
Definition: fe_bernstein_shape_0D.C:30
libMesh::FE
A specific instantiation of the FEBase class.
Definition: fe.h:127
libMesh::libmesh_ignore
void libmesh_ignore(const Args &...)
Definition: libmesh_common.h:515
libMesh::Polyhedron::subelement_coordinates
std::tuple< unsigned int, Real, Real, Real > subelement_coordinates(const Point &p, Real tol=TOLERANCE *TOLERANCE) const
Definition: cell_polyhedron.C:688
libMesh::libmesh_assert
libmesh_assert(ctx)
libMesh::fe_fdm_second_deriv
OutputShape fe_fdm_second_deriv(const Elem *elem, const Order order, const unsigned int i, const unsigned int j, const Point &p, const bool add_p_level, OutputShape(*deriv_func)(const Elem *, const Order, const unsigned int, const unsigned int, const Point &, const bool))
Definition: fe.C:969
libMesh::FE::default_shape_derivs
static void default_shape_derivs(const Elem *elem, const Order o, const unsigned int i, const unsigned int j, const std::vector< Point > &p, std::vector< OutputShape > &v, const bool add_p_level=true)
A default implementation for shape_derivs.
Definition: fe.h:769
libMesh::Utility::enum_to_string
std::string enum_to_string(const T e)
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:144
libMesh::FE::shape_derivs
static void shape_derivs(const Elem *elem, const Order o, const unsigned int i, const unsigned int j, const std::vector< Point > &p, std::vector< OutputShape > &v, const bool add_p_level=true)
Fills v with the derivative of the shape function, evaluated at all points p.
libMesh::make_range
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
Definition: int_range.h:140
libMesh::FE::shapes
static void shapes(const Elem *elem, const Order o, const unsigned int i, const std::vector< Point > &p, std::vector< OutputShape > &v, const bool add_p_level=true)
Fills v with the values of the shape function, evaluated at all points p.
libMesh::C0Polyhedron
The C0Polyhedron is an element in 3D with an arbitrary (but fixed) number of polygonal first-order (C...
Definition: cell_c0polyhedron.h:45
libMesh::Elem::type
virtual ElemType type() const =0
libMesh::Point
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
libMesh::index_range
auto index_range(const T &sizable)
Helper function that returns an IntRange<std::size_t> representing all the indices of the passed-in v...
Definition: int_range.h:117
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)
libMesh::TensorValue
This class defines a tensor in LIBMESH_DIM dimensional Real or Complex space.
Definition: exact_solution.h:46
libMesh::FE::default_shapes
static void default_shapes(const Elem *elem, const Order o, const unsigned int i, const std::vector< Point > &p, std::vector< OutputShape > &v, const bool add_p_level=true)
A default implementation for shapes.
Definition: fe.h:738
libMesh::fe_lagrange_1D_linear_shape
Real fe_lagrange_1D_linear_shape(const unsigned int i, const Real xi)
Definition: fe_lagrange_shape_1D.h:32
libMesh::FE::all_shapes
static void all_shapes(const Elem *elem, const Order o, const std::vector< Point > &p, std::vector< std::vector< OutputShape > > &v, const bool add_p_level=true)
Fills v[i][qp] with the values of the shape functions, evaluated at all points in p...
Generated on Thu Oct 30 2025 06:31:48 for libMesh by   doxygen 1.8.14