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