libMesh
Loading...
Searching...
No Matches
fe_lagrange_shape_3D.C
Go to the documentation of this file.
1// The libMesh Finite Element Library.
2// Copyright (C) 2002-2026 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.
30namespace
31{
32using namespace libMesh;
33
34template <FEFamily T>
35Real fe_lagrange_3D_shape(const ElemType,
36 const Order order,
37 const Elem * elem,
38 const unsigned int i,
39 const Point & p);
40
41template <FEFamily T>
42Real 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
51template <FEFamily T>
52Real 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
63namespace libMesh
64{
65
66
67// TODO: If optimizations for LAGRANGE work well we should do
68// L2_LAGRANGE too...
69LIBMESH_DEFAULT_VECTORIZED_FE(3,L2_LAGRANGE)
70
71
72template<>
73void 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
210template<>
211void 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
223template<>
224void 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
237template<>
238void 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
418template <>
419Real 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
429template <>
430Real 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
440template <>
441Real 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
455template <>
456Real 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
470template <>
471Real 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
483template <>
484Real 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
494template <>
495Real 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
506template <>
507Real 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
518template <>
519Real 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
533template <>
534Real 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
548template <>
549Real 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
561template <>
562Real 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
575template <>
576Real 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
587template <>
588Real 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
599template <>
600Real 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
616template <>
617Real 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
632template <>
633Real 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
647template <>
648Real 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
667namespace
668{
669using namespace libMesh;
670
671template <FEFamily T>
672Real 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
1485template <FEFamily T>
1486Real 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
3242template <FEFamily T>
3243Real 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::C0Polyhedron
The C0Polyhedron is an element in 3D with an arbitrary (but fixed) number of polygonal first-order (C...
Definition cell_c0polyhedron.h:48
libMesh::Elem
This is the base class from which all geometric element types are derived.
Definition elem.h:96
libMesh::Elem::type
virtual ElemType type() const =0
libMesh::Elem::p_level
unsigned int p_level() const
Definition elem.h:3125
libMesh::FEType
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition fe_type.h:197
libMesh::FEType::order
OrderWrapper order
The approximation order of the element (at 0 p-refinement level).
Definition fe_type.h:203
libMesh::FE
A specific instantiation of the FEBase class.
Definition fe.h:128
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::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.
libMesh::FE::shape
static OutputShape shape(const ElemType t, const Order o, const unsigned int i, const Point &p)
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::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::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::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::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::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::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::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::Point
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition point.h:40
libMesh::Polyhedron::master_subelement
virtual std::array< Point, 4 > master_subelement(unsigned int i) const
Definition cell_polyhedron.C:655
libMesh::Polyhedron::subelement
virtual std::array< int, 4 > subelement(unsigned int i) const
Definition cell_polyhedron.h:289
libMesh::Polyhedron::n_subelements
unsigned int n_subelements() const
Definition cell_polyhedron.h:281
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:670
libMesh::TensorValue
This class defines a tensor in LIBMESH_DIM dimensional Real or Complex space.
Definition tensor_value.h:45
b
static const Real b
Definition disjoint_neighbor_test.C:24
libMesh::Utility::enum_to_string
std::string enum_to_string(const T e)
libMesh
The libMesh namespace provides an interface to certain functionality in the library.
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::LIBMESH_DEFAULT_VECTORIZED_FE
LIBMESH_DEFAULT_VECTORIZED_FE(template<>Real FE< 0, BERNSTEIN)
Definition fe_bernstein_shape_0D.C:30
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::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:153
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::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't be...
Definition fe.C:911
libMesh::ElemType
ElemType
Defines an enum for geometric element types.
Definition enum_elem_type.h:33
libMesh::libmesh_ignore
void libmesh_ignore(const Args &...)
Definition libmesh_common.h:569
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::libmesh_assert
libmesh_assert(ctx)
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::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:303
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::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition libmesh_common.h:144
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:176
Generated on Fri Aug 21 2026 13:32:28 for libMesh by doxygen 1.9.8