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 15334774 : LIBMESH_DEFAULT_VECTORIZED_FE(3,L2_LAGRANGE)
70 :
71 :
72 : template<>
73 6120913 : 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 6120913 : const ElemType type = elem->type();
81 :
82 : // Just loop on the harder-to-optimize cases
83 6120913 : if (type != HEX8 && type != HEX27)
84 : {
85 : FE<3,LAGRANGE>::default_all_shapes
86 4948764 : (elem,o,p,v,add_p_level);
87 4948764 : return;
88 : }
89 :
90 : #if LIBMESH_DIM == 3
91 :
92 1172149 : const unsigned int n_sf = v.size();
93 :
94 1172149 : switch (o)
95 : {
96 : // linear Lagrange shape functions
97 936714 : case FIRST:
98 : {
99 492464 : switch (type)
100 : {
101 : // trilinear hexahedral shape functions
102 246232 : case HEX8:
103 : case HEX20:
104 : case HEX27:
105 : {
106 246232 : 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 7110833 : for (auto qp : index_range(p))
114 : {
115 1630453 : const Point & q_point = p[qp];
116 : // Compute hex shape functions as a tensor-product
117 6174119 : const Real xi = q_point(0);
118 6174119 : const Real eta = q_point(1);
119 6174119 : 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 1630453 : {fe_lagrange_1D_linear_shape(0, xi),
124 3260906 : fe_lagrange_1D_linear_shape(1, xi)},
125 3260906 : {fe_lagrange_1D_linear_shape(0, eta),
126 3260906 : fe_lagrange_1D_linear_shape(1, eta)},
127 3260906 : {fe_lagrange_1D_linear_shape(0, zeta),
128 14326384 : fe_lagrange_1D_linear_shape(1, zeta)}};
129 :
130 55567071 : for (unsigned int i : make_range(n_sf))
131 62436576 : v[i][qp] = one_d_shapes[0][i0[i]] *
132 75480200 : one_d_shapes[1][i1[i]] *
133 49392952 : one_d_shapes[2][i2[i]];
134 : }
135 246232 : 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 235435 : case SECOND:
146 : {
147 235435 : switch (type)
148 : {
149 : // triquadratic hexahedral shape functions
150 57351 : 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 57351 : 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 2351002 : for (auto qp : index_range(p))
168 : {
169 503098 : const Point & q_point = p[qp];
170 : // Compute hex shape functions as a tensor-product
171 2115567 : const Real xi = q_point(0);
172 2115567 : const Real eta = q_point(1);
173 2115567 : const Real zeta = q_point(2);
174 :
175 : // linear_shapes[dim][i] = phi_i(p(dim))
176 : Real one_d_shapes[3][3] = {
177 503098 : {fe_lagrange_1D_quadratic_shape(0, xi),
178 1006196 : fe_lagrange_1D_quadratic_shape(1, xi),
179 1006196 : fe_lagrange_1D_quadratic_shape(2, xi)},
180 1006196 : {fe_lagrange_1D_quadratic_shape(0, eta),
181 1006196 : fe_lagrange_1D_quadratic_shape(1, eta),
182 1006196 : fe_lagrange_1D_quadratic_shape(2, eta)},
183 1006196 : {fe_lagrange_1D_quadratic_shape(0, zeta),
184 1006196 : fe_lagrange_1D_quadratic_shape(1, zeta),
185 6140351 : fe_lagrange_1D_quadratic_shape(2, zeta)}};
186 :
187 59235876 : for (unsigned int i : make_range(n_sf))
188 70703955 : v[i][qp] = one_d_shapes[0][i0[i]] *
189 84287601 : one_d_shapes[1][i1[i]] *
190 57120309 : one_d_shapes[2][i2[i]];
191 : }
192 57351 : 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 69348122 : 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 69348122 : (elem,o,i,p,v,add_p_level);
221 69348122 : }
222 :
223 : template<>
224 260918799 : 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 260918799 : (elem,o,i,j,p,v,add_p_level);
235 260918799 : }
236 :
237 : template<>
238 8801512 : 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 8801512 : const ElemType type = elem->type();
246 :
247 : // Just loop on the harder-to-optimize cases
248 8801512 : if (type != HEX8 && type != HEX27)
249 : {
250 : FE<3,LAGRANGE>::default_all_shape_derivs
251 6039042 : (elem,o,p,comps,add_p_level);
252 6039042 : return;
253 : }
254 :
255 : #if LIBMESH_DIM == 3
256 :
257 691867 : libmesh_assert(comps[0]);
258 691867 : libmesh_assert(comps[1]);
259 691867 : libmesh_assert(comps[2]);
260 2762470 : const unsigned int n_sf = comps[0]->size();
261 :
262 2762470 : switch (o)
263 : {
264 : // linear Lagrange shape functions
265 2160273 : case FIRST:
266 : {
267 1087328 : switch (type)
268 : {
269 : // trilinear hexahedral shape functions
270 543664 : case HEX8:
271 : case HEX20:
272 : case HEX27:
273 : {
274 543664 : 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 5800200 : for (auto qp : index_range(p))
282 : {
283 928077 : const Point & q_point = p[qp];
284 : // Compute hex shape functions as a tensor-product
285 3639927 : const Real xi = q_point(0);
286 3639927 : const Real eta = q_point(1);
287 3639927 : 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 928077 : {fe_lagrange_1D_linear_shape(0, xi),
292 1856154 : fe_lagrange_1D_linear_shape(1, xi)},
293 1856154 : {fe_lagrange_1D_linear_shape(0, eta),
294 1856154 : fe_lagrange_1D_linear_shape(1, eta)},
295 1856154 : {fe_lagrange_1D_linear_shape(0, zeta),
296 7352235 : 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 928077 : {fe_lagrange_1D_linear_shape_deriv(0, 0, xi),
301 1856154 : fe_lagrange_1D_linear_shape_deriv(1, 0, xi)},
302 1856154 : {fe_lagrange_1D_linear_shape_deriv(0, 0, eta),
303 1856154 : fe_lagrange_1D_linear_shape_deriv(1, 0, eta)},
304 1856154 : {fe_lagrange_1D_linear_shape_deriv(0, 0, zeta),
305 3639927 : fe_lagrange_1D_linear_shape_deriv(1, 0, zeta)}};
306 :
307 32759343 : for (unsigned int i : make_range(n_sf))
308 : {
309 29119416 : (*comps[0])[i][qp] = one_d_derivs[0][i0[i]] *
310 43968648 : one_d_shapes[1][i1[i]] *
311 29119416 : one_d_shapes[2][i2[i]];
312 29119416 : (*comps[1])[i][qp] = one_d_shapes[0][i0[i]] *
313 36544032 : one_d_derivs[1][i1[i]] *
314 7424616 : one_d_shapes[2][i2[i]];
315 36544032 : (*comps[2])[i][qp] = one_d_shapes[0][i0[i]] *
316 43968648 : one_d_shapes[1][i1[i]] *
317 29119416 : one_d_derivs[2][i2[i]];
318 : }
319 : }
320 543664 : 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 602197 : case SECOND:
331 : {
332 602197 : switch (type)
333 : {
334 : // triquadratic hexahedral shape functions
335 148203 : 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 148203 : 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 3847197 : for (auto qp : index_range(p))
353 : {
354 763214 : const Point & q_point = p[qp];
355 : // Compute hex shape functions as a tensor-product
356 3245000 : const Real xi = q_point(0);
357 3245000 : const Real eta = q_point(1);
358 3245000 : 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 763214 : {fe_lagrange_1D_quadratic_shape(0, xi),
363 1526428 : fe_lagrange_1D_quadratic_shape(1, xi),
364 1526428 : fe_lagrange_1D_quadratic_shape(2, xi)},
365 1526428 : {fe_lagrange_1D_quadratic_shape(0, eta),
366 1526428 : fe_lagrange_1D_quadratic_shape(1, eta),
367 1526428 : fe_lagrange_1D_quadratic_shape(2, eta)},
368 1526428 : {fe_lagrange_1D_quadratic_shape(0, zeta),
369 1526428 : fe_lagrange_1D_quadratic_shape(1, zeta),
370 9350712 : 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 763214 : {fe_lagrange_1D_quadratic_shape_deriv(0, 0, xi),
375 1526428 : fe_lagrange_1D_quadratic_shape_deriv(1, 0, xi),
376 1526428 : fe_lagrange_1D_quadratic_shape_deriv(2, 0, xi)},
377 1526428 : {fe_lagrange_1D_quadratic_shape_deriv(0, 0, eta),
378 1526428 : fe_lagrange_1D_quadratic_shape_deriv(1, 0, eta),
379 1526428 : fe_lagrange_1D_quadratic_shape_deriv(2, 0, eta)},
380 1526428 : {fe_lagrange_1D_quadratic_shape_deriv(0, 0, zeta),
381 1526428 : fe_lagrange_1D_quadratic_shape_deriv(1, 0, zeta),
382 9350712 : fe_lagrange_1D_quadratic_shape_deriv(2, 0, zeta)}};
383 :
384 90860000 : for (unsigned int i : make_range(n_sf))
385 : {
386 87615000 : (*comps[0])[i][qp] = one_d_derivs[0][i0[i]] *
387 128828556 : one_d_shapes[1][i1[i]] *
388 87615000 : one_d_shapes[2][i2[i]];
389 87615000 : (*comps[1])[i][qp] = one_d_shapes[0][i0[i]] *
390 108221778 : one_d_derivs[1][i1[i]] *
391 20606778 : one_d_shapes[2][i2[i]];
392 108221778 : (*comps[2])[i][qp] = one_d_shapes[0][i0[i]] *
393 128828556 : one_d_shapes[1][i1[i]] *
394 87615000 : one_d_derivs[2][i2[i]];
395 : }
396 : }
397 148203 : 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 77920816 : Real FE<3,LAGRANGE>::shape(const ElemType type,
420 : const Order order,
421 : const unsigned int i,
422 : const Point & p)
423 : {
424 77920816 : 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 350250285 : 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 88345395 : libmesh_assert(elem);
448 :
449 : // call the orientation-independent shape functions
450 526941075 : 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 22269341 : 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 5618364 : libmesh_assert(elem);
463 :
464 : // call the orientation-independent shape functions
465 33506069 : 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 444785385 : 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 130786109 : libmesh_assert(elem);
478 674213145 : 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 83943852 : 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 83943852 : 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 1163344362 : 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 299933946 : libmesh_assert(elem);
527 :
528 : // call the orientation-independent shape function derivatives
529 1742812494 : 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 13769610 : 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 3593712 : libmesh_assert(elem);
542 :
543 : // call the orientation-independent shape function derivatives
544 20957034 : 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 3593311602 : 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 921489978 : libmesh_assert(elem);
557 5434515138 : 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 30320760 : 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 30320760 : 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 36188688 : 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 9392010 : libmesh_assert(elem);
608 :
609 : // call the orientation-independent shape function derivatives
610 : return fe_lagrange_3D_shape_second_deriv<LAGRANGE>
611 54972708 : (elem->type(), order + add_p_level*elem->p_level(), elem, i, j, p);
612 : }
613 :
614 :
615 :
616 : template <>
617 16229172 : 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 4359504 : libmesh_assert(elem);
625 :
626 : // call the orientation-independent shape function derivatives
627 : return fe_lagrange_3D_shape_second_deriv<L2_LAGRANGE>
628 24948180 : (elem->type(), order + add_p_level*elem->p_level(), elem, i, j, p);
629 : }
630 :
631 :
632 : template <>
633 629042586 : 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 152996850 : libmesh_assert(elem);
641 : return fe_lagrange_3D_shape_second_deriv<LAGRANGE>
642 935036286 : (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 974100495 : 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 974100495 : switch (order)
681 : {
682 : // linear Lagrange shape functions
683 137700205 : case FIRST:
684 : {
685 137700205 : switch (type)
686 : {
687 : // trilinear hexahedral shape functions
688 65024160 : case HEX8:
689 : case HEX20:
690 : case HEX27:
691 : {
692 19146576 : libmesh_assert_less (i, 8);
693 :
694 : // Compute hex shape functions as a tensor-product
695 65024160 : const Real xi = p(0);
696 65024160 : const Real eta = p(1);
697 65024160 : 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 65024160 : return (fe_lagrange_1D_linear_shape(i0[i], xi)*
705 65024160 : fe_lagrange_1D_linear_shape(i1[i], eta)*
706 110901744 : fe_lagrange_1D_linear_shape(i2[i], zeta));
707 : }
708 :
709 : // linear tetrahedral shape functions
710 62040564 : case TET4:
711 : case TET10:
712 : case TET14:
713 : {
714 15725072 : libmesh_assert_less (i, 4);
715 :
716 : // Area coordinates, pg. 205, Vol. I, Carey, Oden, Becker FEM
717 62040564 : const Real zeta1 = p(0);
718 62040564 : const Real zeta2 = p(1);
719 62040564 : const Real zeta3 = p(2);
720 62040564 : const Real zeta0 = 1. - zeta1 - zeta2 - zeta3;
721 :
722 62040564 : switch(i)
723 : {
724 3931268 : case 0:
725 3931268 : return zeta0;
726 :
727 15510141 : case 1:
728 15510141 : return zeta1;
729 :
730 15510141 : case 2:
731 15510141 : return zeta2;
732 :
733 15510141 : case 3:
734 15510141 : return zeta3;
735 :
736 0 : default:
737 0 : libmesh_error_msg("Invalid i = " << i);
738 : }
739 : }
740 :
741 : // linear prism shape functions
742 5438940 : case PRISM6:
743 : case PRISM15:
744 : case PRISM18:
745 : case PRISM20:
746 : case PRISM21:
747 : {
748 1642458 : libmesh_assert_less (i, 6);
749 :
750 : // Compute prism shape functions as a tensor-product
751 : // of a triangle and an edge
752 :
753 5438940 : Point p2d(p(0),p(1));
754 5438940 : 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 5438940 : return (FE<2,LAGRANGE>::shape(TRI3, FIRST, i1[i], p2d)*
761 9235422 : fe_lagrange_1D_linear_shape(i0[i], p1d));
762 : }
763 :
764 : // linear pyramid shape functions
765 330335 : case PYRAMID5:
766 : case PYRAMID13:
767 : case PYRAMID14:
768 : case PYRAMID18:
769 : {
770 93140 : libmesh_assert_less (i, 5);
771 :
772 330335 : const Real xi = p(0);
773 330335 : const Real eta = p(1);
774 330335 : const Real zeta = p(2);
775 93140 : const Real eps = 1.e-35;
776 :
777 330335 : switch(i)
778 : {
779 66067 : case 0:
780 66067 : return .25*(zeta + xi - 1.)*(zeta + eta - 1.)/((1. - zeta) + eps);
781 :
782 66067 : case 1:
783 66067 : return .25*(zeta - xi - 1.)*(zeta + eta - 1.)/((1. - zeta) + eps);
784 :
785 66067 : case 2:
786 66067 : return .25*(zeta - xi - 1.)*(zeta - eta - 1.)/((1. - zeta) + eps);
787 :
788 66067 : case 3:
789 66067 : return .25*(zeta + xi - 1.)*(zeta - eta - 1.)/((1. - zeta) + eps);
790 :
791 18628 : case 4:
792 18628 : return zeta;
793 :
794 0 : default:
795 0 : libmesh_error_msg("Invalid i = " << i);
796 : }
797 : }
798 :
799 4866206 : case C0POLYHEDRON:
800 : {
801 : // Polyhedra require using newer FE APIs
802 4866206 : 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 1361249 : libmesh_assert(elem->type() == C0POLYHEDRON);
807 :
808 1361249 : 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 4866206 : const auto [s, a, b, c] = poly.subelement_coordinates(p, 100);
814 4866206 : if (s == invalid_uint)
815 0 : return 0;
816 1361249 : libmesh_assert_less(s, poly.n_subelements());
817 :
818 4866206 : const auto subtet = poly.subelement(s);
819 :
820 : // Avoid signed/unsigned comparison warnings
821 4866206 : const int nodei = i;
822 4866206 : if (nodei == subtet[0])
823 576766 : return 1-a-b-c;
824 4289440 : if (nodei == subtet[1])
825 576766 : return a;
826 3712674 : if (nodei == subtet[2])
827 576766 : return b;
828 3135908 : if (nodei == subtet[3])
829 576766 : return c;
830 :
831 : // Basis function i is not supported on p's subtet
832 715813 : 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 276831562 : case SECOND:
843 : {
844 276831562 : switch (type)
845 : {
846 :
847 : // serendipity hexahedral quadratic shape functions
848 22160180 : case HEX20:
849 : {
850 6182820 : libmesh_assert_less (i, 20);
851 :
852 22160180 : const Real xi = p(0);
853 22160180 : const Real eta = p(1);
854 22160180 : 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 22160180 : const Real x = .5*(xi + 1.);
859 22160180 : const Real y = .5*(eta + 1.);
860 22160180 : const Real z = .5*(zeta + 1.);
861 :
862 22160180 : switch (i)
863 : {
864 1108009 : case 0:
865 1108009 : return (1. - x)*(1. - y)*(1. - z)*(1. - 2.*x - 2.*y - 2.*z);
866 :
867 1108009 : case 1:
868 1108009 : return x*(1. - y)*(1. - z)*(2.*x - 2.*y - 2.*z - 1.);
869 :
870 1108009 : case 2:
871 1108009 : return x*y*(1. - z)*(2.*x + 2.*y - 2.*z - 3.);
872 :
873 1108009 : case 3:
874 1108009 : return (1. - x)*y*(1. - z)*(2.*y - 2.*x - 2.*z - 1.);
875 :
876 1108009 : case 4:
877 1108009 : return (1. - x)*(1. - y)*z*(2.*z - 2.*x - 2.*y - 1.);
878 :
879 1108009 : case 5:
880 1108009 : return x*(1. - y)*z*(2.*x - 2.*y + 2.*z - 3.);
881 :
882 1108009 : case 6:
883 1108009 : return x*y*z*(2.*x + 2.*y + 2.*z - 5.);
884 :
885 1108009 : case 7:
886 1108009 : return (1. - x)*y*z*(2.*y - 2.*x + 2.*z - 3.);
887 :
888 1108009 : case 8:
889 1108009 : return 4.*x*(1. - x)*(1. - y)*(1. - z);
890 :
891 1108009 : case 9:
892 1108009 : return 4.*x*y*(1. - y)*(1. - z);
893 :
894 1108009 : case 10:
895 1108009 : return 4.*x*(1. - x)*y*(1. - z);
896 :
897 1108009 : case 11:
898 1108009 : return 4.*(1. - x)*y*(1. - y)*(1. - z);
899 :
900 1108009 : case 12:
901 1108009 : return 4.*(1. - x)*(1. - y)*z*(1. - z);
902 :
903 1108009 : case 13:
904 1108009 : return 4.*x*(1. - y)*z*(1. - z);
905 :
906 1108009 : case 14:
907 1108009 : return 4.*x*y*z*(1. - z);
908 :
909 1108009 : case 15:
910 1108009 : return 4.*(1. - x)*y*z*(1. - z);
911 :
912 1108009 : case 16:
913 1108009 : return 4.*x*(1. - x)*(1. - y)*z;
914 :
915 1108009 : case 17:
916 1108009 : return 4.*x*y*(1. - y)*z;
917 :
918 1108009 : case 18:
919 1108009 : return 4.*x*(1. - x)*y*z;
920 :
921 1108009 : case 19:
922 1108009 : 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 102248514 : case HEX8:
931 859032 : libmesh_assert_msg(T == L2_LAGRANGE,
932 : "High order on first order elements only supported for L2 families");
933 : libmesh_fallthrough();
934 38835099 : case HEX27:
935 : {
936 39805911 : libmesh_assert_less (i, 27);
937 :
938 : // Compute hex shape functions as a tensor-product
939 141195393 : const Real xi = p(0);
940 141195393 : const Real eta = p(1);
941 141195393 : 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 141195393 : return (fe_lagrange_1D_quadratic_shape(i0[i], xi)*
952 141195393 : fe_lagrange_1D_quadratic_shape(i1[i], eta)*
953 175844925 : fe_lagrange_1D_quadratic_shape(i2[i], zeta));
954 : }
955 :
956 : // quadratic tetrahedral shape functions
957 34943340 : case TET4:
958 51000 : libmesh_assert_msg(T == L2_LAGRANGE,
959 : "High order on first order elements only supported for L2 families");
960 : libmesh_fallthrough();
961 12001770 : case TET10:
962 12103740 : libmesh_assert_less (i, 10);
963 : libmesh_fallthrough();
964 : case TET14:
965 : {
966 13493510 : libmesh_assert_less (i, 14);
967 :
968 : // Area coordinates, pg. 205, Vol. I, Carey, Oden, Becker FEM
969 48385850 : const Real zeta1 = p(0);
970 48385850 : const Real zeta2 = p(1);
971 48385850 : const Real zeta3 = p(2);
972 48385850 : const Real zeta0 = 1. - zeta1 - zeta2 - zeta3;
973 :
974 48385850 : switch(i)
975 : {
976 4838585 : case 0:
977 4838585 : return zeta0*(2.*zeta0 - 1.);
978 :
979 4838585 : case 1:
980 4838585 : return zeta1*(2.*zeta1 - 1.);
981 :
982 4838585 : case 2:
983 4838585 : return zeta2*(2.*zeta2 - 1.);
984 :
985 4838585 : case 3:
986 4838585 : return zeta3*(2.*zeta3 - 1.);
987 :
988 4838585 : case 4:
989 4838585 : return 4.*zeta0*zeta1;
990 :
991 4838585 : case 5:
992 4838585 : return 4.*zeta1*zeta2;
993 :
994 4838585 : case 6:
995 4838585 : return 4.*zeta2*zeta0;
996 :
997 4838585 : case 7:
998 4838585 : return 4.*zeta0*zeta3;
999 :
1000 4838585 : case 8:
1001 4838585 : return 4.*zeta1*zeta3;
1002 :
1003 4838585 : case 9:
1004 4838585 : return 4.*zeta2*zeta3;
1005 :
1006 0 : default:
1007 0 : libmesh_error_msg("Invalid i = " << i);
1008 : }
1009 : }
1010 :
1011 : // "serendipity" prism
1012 17041830 : case PRISM15:
1013 : {
1014 5164680 : libmesh_assert_less (i, 15);
1015 :
1016 17041830 : const Real xi = p(0);
1017 17041830 : const Real eta = p(1);
1018 17041830 : const Real zeta = p(2);
1019 :
1020 17041830 : switch(i)
1021 : {
1022 1136122 : case 0:
1023 1136122 : return (1. - zeta)*(xi + eta - 1.)*(xi + eta + 0.5*zeta);
1024 :
1025 1136122 : case 1:
1026 1136122 : return (1. - zeta)*xi*(xi - 1. - 0.5*zeta);
1027 :
1028 1136122 : case 2: // phi1 with xi <- eta
1029 1136122 : return (1. - zeta)*eta*(eta - 1. - 0.5*zeta);
1030 :
1031 1136122 : case 3: // phi0 with zeta <- (-zeta)
1032 1136122 : return (1. + zeta)*(xi + eta - 1.)*(xi + eta - 0.5*zeta);
1033 :
1034 1136122 : case 4: // phi1 with zeta <- (-zeta)
1035 1136122 : return (1. + zeta)*xi*(xi - 1. + 0.5*zeta);
1036 :
1037 1136122 : case 5: // phi4 with xi <- eta
1038 1136122 : return (1. + zeta)*eta*(eta - 1. + 0.5*zeta);
1039 :
1040 1136122 : case 6:
1041 1136122 : return 2.*(1. - zeta)*xi*(1. - xi - eta);
1042 :
1043 1136122 : case 7:
1044 1136122 : return 2.*(1. - zeta)*xi*eta;
1045 :
1046 1136122 : case 8:
1047 1136122 : return 2.*(1. - zeta)*eta*(1. - xi - eta);
1048 :
1049 1136122 : case 9:
1050 1136122 : return (1. - zeta)*(1. + zeta)*(1. - xi - eta);
1051 :
1052 1136122 : case 10:
1053 1136122 : return (1. - zeta)*(1. + zeta)*xi;
1054 :
1055 1136122 : case 11: // phi10 with xi <-> eta
1056 1136122 : return (1. - zeta)*(1. + zeta)*eta;
1057 :
1058 1136122 : case 12: // phi6 with zeta <- (-zeta)
1059 1136122 : return 2.*(1. + zeta)*xi*(1. - xi - eta);
1060 :
1061 1136122 : case 13: // phi7 with zeta <- (-zeta)
1062 1136122 : return 2.*(1. + zeta)*xi*eta;
1063 :
1064 1136122 : case 14: // phi8 with zeta <- (-zeta)
1065 1136122 : 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 27220698 : case PRISM6:
1074 79704 : libmesh_assert_msg(T == L2_LAGRANGE,
1075 : "High order on first order elements only supported for L2 families");
1076 : libmesh_fallthrough();
1077 11970720 : case PRISM18:
1078 : case PRISM20:
1079 : case PRISM21:
1080 : {
1081 12291012 : libmesh_assert_less (i, 18);
1082 :
1083 : // Compute prism shape functions as a tensor-product
1084 : // of a triangle and an edge
1085 :
1086 39432006 : Point p2d(p(0),p(1));
1087 39432006 : 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 39432006 : return (FE<2,LAGRANGE>::shape(TRI6, SECOND, i1[i], p2d)*
1094 50707818 : 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 4044911 : case PYRAMID13:
1101 : {
1102 1146548 : libmesh_assert_less (i, 13);
1103 :
1104 4044911 : const Real xi = p(0);
1105 4044911 : const Real eta = p(1);
1106 4044911 : const Real zeta = p(2);
1107 1146548 : const Real eps = 1.e-35;
1108 :
1109 : // Denominators are perturbed by epsilon to avoid
1110 : // divide-by-zero issues.
1111 4044911 : Real den = (1. - zeta + eps);
1112 :
1113 4044911 : switch(i)
1114 : {
1115 311147 : case 0:
1116 311147 : return 0.25*(-xi - eta - 1.)*((1. - xi)*(1. - eta) - zeta + xi*eta*zeta/den);
1117 :
1118 311147 : case 1:
1119 311147 : return 0.25*(-eta + xi - 1.)*((1. + xi)*(1. - eta) - zeta - xi*eta*zeta/den);
1120 :
1121 311147 : case 2:
1122 311147 : return 0.25*(xi + eta - 1.)*((1. + xi)*(1. + eta) - zeta + xi*eta*zeta/den);
1123 :
1124 311147 : case 3:
1125 311147 : return 0.25*(eta - xi - 1.)*((1. - xi)*(1. + eta) - zeta - xi*eta*zeta/den);
1126 :
1127 311147 : case 4:
1128 311147 : return zeta*(2.*zeta - 1.);
1129 :
1130 311147 : case 5:
1131 311147 : return 0.5*(1. + xi - zeta)*(1. - xi - zeta)*(1. - eta - zeta)/den;
1132 :
1133 311147 : case 6:
1134 311147 : return 0.5*(1. + eta - zeta)*(1. - eta - zeta)*(1. + xi - zeta)/den;
1135 :
1136 311147 : case 7:
1137 311147 : return 0.5*(1. + xi - zeta)*(1. - xi - zeta)*(1. + eta - zeta)/den;
1138 :
1139 311147 : case 8:
1140 311147 : return 0.5*(1. + eta - zeta)*(1. - eta - zeta)*(1. - xi - zeta)/den;
1141 :
1142 311147 : case 9:
1143 311147 : return zeta*(1. - xi - zeta)*(1. - eta - zeta)/den;
1144 :
1145 311147 : case 10:
1146 311147 : return zeta*(1. + xi - zeta)*(1. - eta - zeta)/den;
1147 :
1148 311147 : case 11:
1149 311147 : return zeta*(1. + eta - zeta)*(1. + xi - zeta)/den;
1150 :
1151 311147 : case 12:
1152 311147 : 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 3276336 : 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 1295056 : case PYRAMID14:
1167 : case PYRAMID18:
1168 : {
1169 1295056 : libmesh_assert_less (i, 14);
1170 :
1171 4571392 : const Real xi = p(0);
1172 4571392 : const Real eta = p(1);
1173 4571392 : const Real zeta = p(2);
1174 1295056 : 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 4571392 : p1 = 0.5*(1. - eta - zeta), // back
1180 4571392 : p2 = 0.5*(1. + xi - zeta), // left
1181 4571392 : p3 = 0.5*(1. + eta - zeta), // front
1182 4571392 : p4 = 0.5*(1. - xi - zeta); // right
1183 :
1184 : // Denominators are perturbed by epsilon to avoid
1185 : // divide-by-zero issues.
1186 : Real
1187 4571392 : den = (-1. + zeta + eps),
1188 4571392 : den2 = den*den;
1189 :
1190 4571392 : switch(i)
1191 : {
1192 326528 : case 0:
1193 326528 : return p4*p1*(xi*eta - zeta + zeta*zeta)/den2;
1194 :
1195 326528 : case 1:
1196 326528 : return -p1*p2*(xi*eta + zeta - zeta*zeta)/den2;
1197 :
1198 326528 : case 2:
1199 326528 : return p2*p3*(xi*eta - zeta + zeta*zeta)/den2;
1200 :
1201 326528 : case 3:
1202 326528 : return -p3*p4*(xi*eta + zeta - zeta*zeta)/den2;
1203 :
1204 326528 : case 4:
1205 326528 : return zeta*(2.*zeta - 1.);
1206 :
1207 326528 : case 5:
1208 326528 : return -4.*p2*p1*p4*eta/den2;
1209 :
1210 326528 : case 6:
1211 326528 : return 4.*p1*p2*p3*xi/den2;
1212 :
1213 326528 : case 7:
1214 326528 : return 4.*p2*p3*p4*eta/den2;
1215 :
1216 326528 : case 8:
1217 326528 : return -4.*p3*p4*p1*xi/den2;
1218 :
1219 326528 : case 9:
1220 326528 : return -4.*p1*p4*zeta/den;
1221 :
1222 326528 : case 10:
1223 326528 : return -4.*p2*p1*zeta/den;
1224 :
1225 326528 : case 11:
1226 326528 : return -4.*p3*p2*zeta/den;
1227 :
1228 326528 : case 12:
1229 326528 : return -4.*p4*p3*zeta/den;
1230 :
1231 326528 : case 13:
1232 326528 : 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 559568728 : case THIRD:
1246 : {
1247 559568728 : switch (type)
1248 : {
1249 : // quadratic Lagrange shape functions with cubic bubbles
1250 25870300 : case PRISM20:
1251 : {
1252 7657960 : 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 25870300 : Point p2d(p(0),p(1));
1263 25870300 : Real p1d = p(2);
1264 :
1265 25870300 : const Real mainval = FE<3,LAGRANGE>::shape(PRISM21, THIRD, i, p);
1266 :
1267 25870300 : if (i0[i] != 2)
1268 5360572 : return mainval;
1269 :
1270 7761090 : const Real bubbleval =
1271 7761090 : FE<2,LAGRANGE>::shape(TRI7, THIRD, 6, p2d) *
1272 2297388 : fe_lagrange_1D_quadratic_shape(2, p1d);
1273 :
1274 7761090 : if (i < 12) // vertices
1275 3880545 : return mainval - bubbleval / 9;
1276 :
1277 3880545 : return mainval + bubbleval * (Real(4) / 9);
1278 : }
1279 :
1280 : // quadratic Lagrange shape functions with cubic bubbles
1281 67635352 : case PRISM21:
1282 : {
1283 20003566 : libmesh_assert_less (i, 21);
1284 :
1285 : // Compute prism shape functions as a tensor-product
1286 : // of a triangle and an edge
1287 :
1288 67635352 : Point p2d(p(0),p(1));
1289 67635352 : 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 67635352 : return (FE<2,LAGRANGE>::shape(TRI7, THIRD, i1[i], p2d)*
1296 85745078 : fe_lagrange_1D_quadratic_shape(i0[i], p1d));
1297 : }
1298 :
1299 : // Weird rational shape functions with weirder bubbles...
1300 90004050 : case PYRAMID18:
1301 : {
1302 25210692 : libmesh_assert_less (i, 18);
1303 :
1304 90004050 : const Real xi = p(0);
1305 90004050 : const Real eta = p(1);
1306 90004050 : const Real zeta = p(2);
1307 25210692 : 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 90004050 : p1 = 0.5*(1. - eta - zeta), // back
1313 90004050 : p2 = 0.5*(1. + xi - zeta), // left
1314 90004050 : p3 = 0.5*(1. + eta - zeta), // front
1315 90004050 : 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 90004050 : den = (-1. + zeta + eps),
1321 90004050 : 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 25210692 : constexpr Real alpha = 0.5;
1328 : const Real
1329 90004050 : bub_f1 = ((1-alpha)*(1-zeta) + alpha*(-eta)),
1330 90004050 : bub_f2 = ((1-alpha)*(1-zeta) + alpha*(xi)),
1331 90004050 : bub_f3 = ((1-alpha)*(1-zeta) + alpha*(eta)),
1332 90004050 : bub_f4 = ((1-alpha)*(1-zeta) + alpha*(-xi));
1333 :
1334 : const Real
1335 90004050 : bub1 = bub_f1*p1*p2*p4*zeta/den2,
1336 90004050 : bub2 = bub_f2*p1*p2*p3*zeta/den2,
1337 90004050 : bub3 = bub_f3*p2*p3*p4*zeta/den2,
1338 90004050 : bub4 = bub_f4*p1*p3*p4*zeta/den2;
1339 :
1340 90004050 : switch(i)
1341 : {
1342 5000225 : case 0:
1343 5000225 : return p4*p1*(xi*eta - zeta + zeta*zeta)/den2 + 3*(bub1+bub4);
1344 :
1345 5000225 : case 1:
1346 5000225 : return -p1*p2*(xi*eta + zeta - zeta*zeta)/den2 + 3*(bub1+bub2);
1347 :
1348 5000225 : case 2:
1349 5000225 : return p2*p3*(xi*eta - zeta + zeta*zeta)/den2 + 3*(bub2+bub3);
1350 :
1351 5000225 : case 3:
1352 5000225 : return -p3*p4*(xi*eta + zeta - zeta*zeta)/den2 + 3*(bub3+bub4);
1353 :
1354 5000225 : case 4:
1355 5000225 : return zeta*(2.*zeta - 1.) + 3*(bub1+bub2+bub3+bub4);
1356 :
1357 5000225 : case 5:
1358 5000225 : return -4.*p2*p1*p4*eta/den2 - 12*bub1;
1359 :
1360 5000225 : case 6:
1361 5000225 : return 4.*p1*p2*p3*xi/den2 - 12*bub2;
1362 :
1363 5000225 : case 7:
1364 5000225 : return 4.*p2*p3*p4*eta/den2 - 12*bub3;
1365 :
1366 5000225 : case 8:
1367 5000225 : return -4.*p3*p4*p1*xi/den2 - 12*bub4;
1368 :
1369 5000225 : case 9:
1370 5000225 : return -4.*p1*p4*zeta/den - 12*(bub1+bub4);
1371 :
1372 5000225 : case 10:
1373 5000225 : return -4.*p2*p1*zeta/den - 12*(bub1+bub2);
1374 :
1375 5000225 : case 11:
1376 5000225 : return -4.*p3*p2*zeta/den - 12*(bub2+bub3);
1377 :
1378 5000225 : case 12:
1379 5000225 : return -4.*p4*p3*zeta/den - 12*(bub3+bub4);
1380 :
1381 5000225 : case 13:
1382 5000225 : return 16.*p1*p2*p3*p4/den2;
1383 :
1384 5000225 : case 14:
1385 5000225 : return 27*bub1;
1386 :
1387 5000225 : case 15:
1388 5000225 : return 27*bub2;
1389 :
1390 5000225 : case 16:
1391 5000225 : return 27*bub3;
1392 :
1393 5000225 : case 17:
1394 5000225 : 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 376059026 : case TET14:
1403 : {
1404 97263824 : libmesh_assert_less (i, 14);
1405 :
1406 : // Area coordinates, pg. 205, Vol. I, Carey, Oden, Becker FEM
1407 376059026 : const Real zeta1 = p(0);
1408 376059026 : const Real zeta2 = p(1);
1409 376059026 : const Real zeta3 = p(2);
1410 376059026 : const Real zeta0 = 1. - zeta1 - zeta2 - zeta3;
1411 :
1412 : // Bubble functions (not yet scaled) on side nodes
1413 376059026 : const Real bubble_012 = zeta0*zeta1*zeta2;
1414 376059026 : const Real bubble_013 = zeta0*zeta1*zeta3;
1415 376059026 : const Real bubble_123 = zeta1*zeta2*zeta3;
1416 376059026 : const Real bubble_023 = zeta0*zeta2*zeta3;
1417 :
1418 376059026 : switch(i)
1419 : {
1420 26861359 : case 0:
1421 26861359 : return zeta0*(2.*zeta0 - 1.) + 3.*(bubble_012+bubble_013+bubble_023);
1422 :
1423 26861359 : case 1:
1424 26861359 : return zeta1*(2.*zeta1 - 1.) + 3.*(bubble_012+bubble_013+bubble_123);
1425 :
1426 26861359 : case 2:
1427 26861359 : return zeta2*(2.*zeta2 - 1.) + 3.*(bubble_012+bubble_023+bubble_123);
1428 :
1429 26861359 : case 3:
1430 26861359 : return zeta3*(2.*zeta3 - 1.) + 3.*(bubble_013+bubble_023+bubble_123);
1431 :
1432 26861359 : case 4:
1433 26861359 : return 4.*zeta0*zeta1 - 12.*(bubble_012+bubble_013);
1434 :
1435 26861359 : case 5:
1436 26861359 : return 4.*zeta1*zeta2 - 12.*(bubble_012+bubble_123);
1437 :
1438 26861359 : case 6:
1439 26861359 : return 4.*zeta2*zeta0 - 12.*(bubble_012+bubble_023);
1440 :
1441 26861359 : case 7:
1442 26861359 : return 4.*zeta0*zeta3 - 12.*(bubble_013+bubble_023);
1443 :
1444 26861359 : case 8:
1445 26861359 : return 4.*zeta1*zeta3 - 12.*(bubble_013+bubble_123);
1446 :
1447 26861359 : case 9:
1448 26861359 : return 4.*zeta2*zeta3 - 12.*(bubble_023+bubble_123);
1449 :
1450 26861359 : case 10:
1451 26861359 : return 27.*bubble_012;
1452 :
1453 26861359 : case 11:
1454 26861359 : return 27.*bubble_013;
1455 :
1456 26861359 : case 12:
1457 26861359 : return 27.*bubble_123;
1458 :
1459 26861359 : case 13:
1460 26861359 : 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 4857872514 : 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 1249283916 : libmesh_assert_less (j, 3);
1496 :
1497 4857872514 : switch (order)
1498 : {
1499 : // linear Lagrange shape functions
1500 202301184 : case FIRST:
1501 : {
1502 202301184 : switch (type)
1503 : {
1504 : // trilinear hexahedral shape functions
1505 168184800 : case HEX8:
1506 : case HEX20:
1507 : case HEX27:
1508 : {
1509 42325584 : libmesh_assert_less (i, 8);
1510 :
1511 : // Compute hex shape functions as a tensor-product
1512 168184800 : const Real xi = p(0);
1513 168184800 : const Real eta = p(1);
1514 168184800 : 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 168184800 : switch(j)
1521 : {
1522 56061600 : case 0:
1523 56061600 : return (fe_lagrange_1D_linear_shape_deriv(i0[i], 0, xi)*
1524 56061600 : fe_lagrange_1D_linear_shape (i1[i], eta)*
1525 98014672 : fe_lagrange_1D_linear_shape (i2[i], zeta));
1526 :
1527 56061600 : case 1:
1528 56061600 : return (fe_lagrange_1D_linear_shape (i0[i], xi)*
1529 56061600 : fe_lagrange_1D_linear_shape_deriv(i1[i], 0, eta)*
1530 98014672 : fe_lagrange_1D_linear_shape (i2[i], zeta));
1531 :
1532 56061600 : case 2:
1533 56061600 : return (fe_lagrange_1D_linear_shape (i0[i], xi)*
1534 56061600 : fe_lagrange_1D_linear_shape (i1[i], eta)*
1535 77038136 : 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 1052616 : case TET4:
1544 : case TET10:
1545 : case TET14:
1546 : {
1547 340884 : libmesh_assert_less (i, 4);
1548 :
1549 : // Area coordinates, pg. 205, Vol. I, Carey, Oden, Becker FEM
1550 340884 : const Real dzeta0dxi = -1.;
1551 340884 : const Real dzeta1dxi = 1.;
1552 340884 : const Real dzeta2dxi = 0.;
1553 340884 : const Real dzeta3dxi = 0.;
1554 :
1555 340884 : const Real dzeta0deta = -1.;
1556 340884 : const Real dzeta1deta = 0.;
1557 340884 : const Real dzeta2deta = 1.;
1558 340884 : const Real dzeta3deta = 0.;
1559 :
1560 340884 : const Real dzeta0dzeta = -1.;
1561 340884 : const Real dzeta1dzeta = 0.;
1562 340884 : const Real dzeta2dzeta = 0.;
1563 340884 : const Real dzeta3dzeta = 1.;
1564 :
1565 1052616 : switch (j)
1566 : {
1567 : // d()/dxi
1568 350872 : case 0:
1569 : {
1570 237204 : switch(i)
1571 : {
1572 28407 : case 0:
1573 28407 : return dzeta0dxi;
1574 :
1575 28407 : case 1:
1576 28407 : return dzeta1dxi;
1577 :
1578 28407 : case 2:
1579 28407 : return dzeta2dxi;
1580 :
1581 28407 : case 3:
1582 28407 : return dzeta3dxi;
1583 :
1584 0 : default:
1585 0 : libmesh_error_msg("Invalid i = " << i);
1586 : }
1587 : }
1588 :
1589 : // d()/deta
1590 350872 : case 1:
1591 : {
1592 237204 : switch(i)
1593 : {
1594 28407 : case 0:
1595 28407 : return dzeta0deta;
1596 :
1597 28407 : case 1:
1598 28407 : return dzeta1deta;
1599 :
1600 28407 : case 2:
1601 28407 : return dzeta2deta;
1602 :
1603 28407 : case 3:
1604 28407 : return dzeta3deta;
1605 :
1606 0 : default:
1607 0 : libmesh_error_msg("Invalid i = " << i);
1608 : }
1609 : }
1610 :
1611 : // d()/dzeta
1612 350872 : case 2:
1613 : {
1614 237204 : switch(i)
1615 : {
1616 28407 : case 0:
1617 28407 : return dzeta0dzeta;
1618 :
1619 28407 : case 1:
1620 28407 : return dzeta1dzeta;
1621 :
1622 28407 : case 2:
1623 28407 : return dzeta2dzeta;
1624 :
1625 28407 : case 3:
1626 28407 : 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 16866756 : case PRISM6:
1640 : case PRISM15:
1641 : case PRISM18:
1642 : case PRISM20:
1643 : case PRISM21:
1644 : {
1645 4703472 : libmesh_assert_less (i, 6);
1646 :
1647 : // Compute prism shape functions as a tensor-product
1648 : // of a triangle and an edge
1649 :
1650 16866756 : Point p2d(p(0),p(1));
1651 16866756 : 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 16866756 : switch (j)
1658 : {
1659 : // d()/dxi
1660 5622252 : case 0:
1661 5622252 : return (FE<2,LAGRANGE>::shape_deriv(TRI3, FIRST, i1[i], 0, p2d)*
1662 9676680 : fe_lagrange_1D_linear_shape(i0[i], p1d));
1663 :
1664 : // d()/deta
1665 5622252 : case 1:
1666 5622252 : return (FE<2,LAGRANGE>::shape_deriv(TRI3, FIRST, i1[i], 1, p2d)*
1667 9676680 : fe_lagrange_1D_linear_shape(i0[i], p1d));
1668 :
1669 : // d()/dzeta
1670 5622252 : case 2:
1671 5622252 : return (FE<2,LAGRANGE>::shape(TRI3, FIRST, i1[i], p2d)*
1672 7649466 : 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 1134600 : case PYRAMID5:
1681 : case PYRAMID13:
1682 : case PYRAMID14:
1683 : case PYRAMID18:
1684 : {
1685 308040 : libmesh_assert_less (i, 5);
1686 :
1687 1134600 : const Real xi = p(0);
1688 1134600 : const Real eta = p(1);
1689 1134600 : const Real zeta = p(2);
1690 308040 : const Real eps = 1.e-35;
1691 :
1692 1134600 : switch (j)
1693 : {
1694 : // d/dxi
1695 378200 : case 0:
1696 378200 : switch(i)
1697 : {
1698 75640 : case 0:
1699 75640 : return .25*(zeta + eta - 1.)/((1. - zeta) + eps);
1700 :
1701 75640 : case 1:
1702 75640 : return -.25*(zeta + eta - 1.)/((1. - zeta) + eps);
1703 :
1704 75640 : case 2:
1705 75640 : return -.25*(zeta - eta - 1.)/((1. - zeta) + eps);
1706 :
1707 75640 : case 3:
1708 75640 : return .25*(zeta - eta - 1.)/((1. - zeta) + eps);
1709 :
1710 20536 : case 4:
1711 20536 : return 0;
1712 :
1713 0 : default:
1714 0 : libmesh_error_msg("Invalid i = " << i);
1715 : }
1716 :
1717 :
1718 : // d/deta
1719 378200 : case 1:
1720 378200 : switch(i)
1721 : {
1722 75640 : case 0:
1723 75640 : return .25*(zeta + xi - 1.)/((1. - zeta) + eps);
1724 :
1725 75640 : case 1:
1726 75640 : return .25*(zeta - xi - 1.)/((1. - zeta) + eps);
1727 :
1728 75640 : case 2:
1729 75640 : return -.25*(zeta - xi - 1.)/((1. - zeta) + eps);
1730 :
1731 75640 : case 3:
1732 75640 : return -.25*(zeta + xi - 1.)/((1. - zeta) + eps);
1733 :
1734 20536 : case 4:
1735 20536 : return 0;
1736 :
1737 0 : default:
1738 0 : libmesh_error_msg("Invalid i = " << i);
1739 : }
1740 :
1741 :
1742 : // d/dzeta
1743 378200 : 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 378200 : num = zeta*(2. - zeta) - 1.,
1749 378200 : den = (1. - zeta + eps)*(1. - zeta + eps);
1750 :
1751 378200 : switch(i)
1752 : {
1753 151280 : case 0:
1754 : case 2:
1755 151280 : return .25*(num + xi*eta)/den;
1756 :
1757 151280 : case 1:
1758 : case 3:
1759 151280 : return .25*(num - xi*eta)/den;
1760 :
1761 20536 : case 4:
1762 20536 : 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 15062412 : case C0POLYHEDRON:
1775 : {
1776 : // Polyhedra require using newer FE APIs
1777 15062412 : 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 3766728 : libmesh_assert(elem->type() == C0POLYHEDRON);
1782 :
1783 3766728 : 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 15062412 : const auto [s, a, b, c] = poly.subelement_coordinates(p, 100);
1789 15062412 : if (s == invalid_uint)
1790 0 : return 0;
1791 3766728 : libmesh_assert_less(s, poly.n_subelements());
1792 :
1793 15062412 : const auto subtet = poly.subelement(s);
1794 :
1795 : // Find derivatives w.r.t. subtriangle barycentric
1796 : // coordinates
1797 3766728 : Real du_da = 0, du_db = 0, du_dc = 0;
1798 :
1799 : // Avoid signed/unsigned comparison warnings
1800 15062412 : const int nodei = i;
1801 15062412 : if (nodei == subtet[0])
1802 446463 : du_da = du_db = du_dc = -1;
1803 13276935 : else if (nodei == subtet[1])
1804 446463 : du_da = 1;
1805 11491458 : else if (nodei == subtet[2])
1806 446463 : du_db = 1;
1807 9705981 : else if (nodei == subtet[3])
1808 446463 : du_dc = 1;
1809 : else
1810 : // Basis function i is not supported on p's subtet
1811 1980876 : 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 7141908 : const auto master_points = poly.master_subelement(s);
1820 :
1821 1785852 : const RealTensor dXi_dA(
1822 7141908 : 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 7141908 : 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 7141908 : 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 7141908 : const Real jac = dXi_dA.det();
1831 :
1832 7141908 : switch (j)
1833 : {
1834 : // d()/dxi
1835 595284 : case 0:
1836 : {
1837 2380636 : const Real da_dxi = (dXi_dA(2,2)*dXi_dA(1,1) - dXi_dA(2,1)*dXi_dA(1,2)) / jac;
1838 2380636 : const Real db_dxi = -(dXi_dA(2,2)*dXi_dA(1,0) - dXi_dA(2,0)*dXi_dA(1,2)) / jac;
1839 2380636 : const Real dc_dxi = (dXi_dA(2,1)*dXi_dA(1,0) - dXi_dA(2,0)*dXi_dA(1,1)) / jac;
1840 2380636 : return du_da*da_dxi + du_db*db_dxi + du_dc*dc_dxi;
1841 : }
1842 : // d()/deta
1843 595284 : case 1:
1844 : {
1845 2380636 : const Real da_deta = -(dXi_dA(2,2)*dXi_dA(0,1) - dXi_dA(2,1)*dXi_dA(0,2) ) / jac;
1846 2380636 : const Real db_deta = (dXi_dA(2,2)*dXi_dA(0,0) - dXi_dA(2,0)*dXi_dA(0,2)) / jac;
1847 2380636 : const Real dc_deta = -(dXi_dA(2,1)*dXi_dA(0,0) - dXi_dA(2,0)*dXi_dA(0,1)) / jac;
1848 2380636 : return du_da*da_deta + du_db*db_deta + du_dc*dc_deta;
1849 : }
1850 : // d()/dzeta
1851 595284 : case 2:
1852 : {
1853 2380636 : const Real da_dzeta = (dXi_dA(1,2)*dXi_dA(0,1) - dXi_dA(1,1)*dXi_dA(0,2) ) / jac;
1854 2380636 : const Real db_dzeta = -(dXi_dA(1,2)*dXi_dA(0,0) - dXi_dA(1,0)*dXi_dA(0,2)) / jac;
1855 2380636 : const Real dc_dzeta = (dXi_dA(1,1)*dXi_dA(0,0) - dXi_dA(1,0)*dXi_dA(0,1)) / jac;
1856 2380636 : 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 3108636069 : case SECOND:
1871 : {
1872 3108636069 : switch (type)
1873 : {
1874 :
1875 : // serendipity hexahedral quadratic shape functions
1876 72459120 : case HEX20:
1877 : {
1878 18239580 : libmesh_assert_less (i, 20);
1879 :
1880 72459120 : const Real xi = p(0);
1881 72459120 : const Real eta = p(1);
1882 72459120 : 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 72459120 : const Real x = .5*(xi + 1.);
1887 72459120 : const Real y = .5*(eta + 1.);
1888 72459120 : const Real z = .5*(zeta + 1.);
1889 :
1890 : // and don't forget the chain rule!
1891 :
1892 72459120 : switch (j)
1893 : {
1894 :
1895 : // d/dx*dx/dxi
1896 24153040 : case 0:
1897 24153040 : switch (i)
1898 : {
1899 1207652 : case 0:
1900 1815582 : return .5*(1. - y)*(1. - z)*((1. - x)*(-2.) +
1901 1207652 : (-1.)*(1. - 2.*x - 2.*y - 2.*z));
1902 :
1903 1207652 : case 1:
1904 1815582 : return .5*(1. - y)*(1. - z)*(x*(2.) +
1905 1207652 : (1.)*(2.*x - 2.*y - 2.*z - 1.));
1906 :
1907 1207652 : case 2:
1908 1815582 : return .5*y*(1. - z)*(x*(2.) +
1909 1207652 : (1.)*(2.*x + 2.*y - 2.*z - 3.));
1910 :
1911 1207652 : case 3:
1912 1815582 : return .5*y*(1. - z)*((1. - x)*(-2.) +
1913 1207652 : (-1.)*(2.*y - 2.*x - 2.*z - 1.));
1914 :
1915 1207652 : case 4:
1916 1815582 : return .5*(1. - y)*z*((1. - x)*(-2.) +
1917 1207652 : (-1.)*(2.*z - 2.*x - 2.*y - 1.));
1918 :
1919 1207652 : case 5:
1920 1815582 : return .5*(1. - y)*z*(x*(2.) +
1921 1207652 : (1.)*(2.*x - 2.*y + 2.*z - 3.));
1922 :
1923 1207652 : case 6:
1924 1815582 : return .5*y*z*(x*(2.) +
1925 1207652 : (1.)*(2.*x + 2.*y + 2.*z - 5.));
1926 :
1927 1207652 : case 7:
1928 1815582 : return .5*y*z*((1. - x)*(-2.) +
1929 1207652 : (-1.)*(2.*y - 2.*x + 2.*z - 3.));
1930 :
1931 1207652 : case 8:
1932 1207652 : return 2.*(1. - y)*(1. - z)*(1. - 2.*x);
1933 :
1934 1207652 : case 9:
1935 1207652 : return 2.*y*(1. - y)*(1. - z);
1936 :
1937 1207652 : case 10:
1938 1207652 : return 2.*y*(1. - z)*(1. - 2.*x);
1939 :
1940 1207652 : case 11:
1941 1207652 : return 2.*y*(1. - y)*(1. - z)*(-1.);
1942 :
1943 1207652 : case 12:
1944 1207652 : return 2.*(1. - y)*z*(1. - z)*(-1.);
1945 :
1946 1207652 : case 13:
1947 1207652 : return 2.*(1. - y)*z*(1. - z);
1948 :
1949 1207652 : case 14:
1950 1207652 : return 2.*y*z*(1. - z);
1951 :
1952 1207652 : case 15:
1953 1207652 : return 2.*y*z*(1. - z)*(-1.);
1954 :
1955 1207652 : case 16:
1956 1207652 : return 2.*(1. - y)*z*(1. - 2.*x);
1957 :
1958 1207652 : case 17:
1959 1207652 : return 2.*y*(1. - y)*z;
1960 :
1961 1207652 : case 18:
1962 1207652 : return 2.*y*z*(1. - 2.*x);
1963 :
1964 1207652 : case 19:
1965 1207652 : 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 24153040 : case 1:
1974 24153040 : switch (i)
1975 : {
1976 1207652 : case 0:
1977 1815582 : return .5*(1. - x)*(1. - z)*((1. - y)*(-2.) +
1978 1207652 : (-1.)*(1. - 2.*x - 2.*y - 2.*z));
1979 :
1980 1207652 : case 1:
1981 1815582 : return .5*x*(1. - z)*((1. - y)*(-2.) +
1982 1207652 : (-1.)*(2.*x - 2.*y - 2.*z - 1.));
1983 :
1984 1207652 : case 2:
1985 1815582 : return .5*x*(1. - z)*(y*(2.) +
1986 1207652 : (1.)*(2.*x + 2.*y - 2.*z - 3.));
1987 :
1988 1207652 : case 3:
1989 1815582 : return .5*(1. - x)*(1. - z)*(y*(2.) +
1990 1207652 : (1.)*(2.*y - 2.*x - 2.*z - 1.));
1991 :
1992 1207652 : case 4:
1993 1815582 : return .5*(1. - x)*z*((1. - y)*(-2.) +
1994 1207652 : (-1.)*(2.*z - 2.*x - 2.*y - 1.));
1995 :
1996 1207652 : case 5:
1997 1815582 : return .5*x*z*((1. - y)*(-2.) +
1998 1207652 : (-1.)*(2.*x - 2.*y + 2.*z - 3.));
1999 :
2000 1207652 : case 6:
2001 1815582 : return .5*x*z*(y*(2.) +
2002 1207652 : (1.)*(2.*x + 2.*y + 2.*z - 5.));
2003 :
2004 1207652 : case 7:
2005 1815582 : return .5*(1. - x)*z*(y*(2.) +
2006 1207652 : (1.)*(2.*y - 2.*x + 2.*z - 3.));
2007 :
2008 1207652 : case 8:
2009 1207652 : return 2.*x*(1. - x)*(1. - z)*(-1.);
2010 :
2011 1207652 : case 9:
2012 1207652 : return 2.*x*(1. - z)*(1. - 2.*y);
2013 :
2014 1207652 : case 10:
2015 1207652 : return 2.*x*(1. - x)*(1. - z);
2016 :
2017 1207652 : case 11:
2018 1207652 : return 2.*(1. - x)*(1. - z)*(1. - 2.*y);
2019 :
2020 1207652 : case 12:
2021 1207652 : return 2.*(1. - x)*z*(1. - z)*(-1.);
2022 :
2023 1207652 : case 13:
2024 1207652 : return 2.*x*z*(1. - z)*(-1.);
2025 :
2026 1207652 : case 14:
2027 1207652 : return 2.*x*z*(1. - z);
2028 :
2029 1207652 : case 15:
2030 1207652 : return 2.*(1. - x)*z*(1. - z);
2031 :
2032 1207652 : case 16:
2033 1207652 : return 2.*x*(1. - x)*z*(-1.);
2034 :
2035 1207652 : case 17:
2036 1207652 : return 2.*x*z*(1. - 2.*y);
2037 :
2038 1207652 : case 18:
2039 1207652 : return 2.*x*(1. - x)*z;
2040 :
2041 1207652 : case 19:
2042 1207652 : 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 24153040 : case 2:
2051 24153040 : switch (i)
2052 : {
2053 1207652 : case 0:
2054 1815582 : return .5*(1. - x)*(1. - y)*((1. - z)*(-2.) +
2055 1207652 : (-1.)*(1. - 2.*x - 2.*y - 2.*z));
2056 :
2057 1207652 : case 1:
2058 1815582 : return .5*x*(1. - y)*((1. - z)*(-2.) +
2059 1207652 : (-1.)*(2.*x - 2.*y - 2.*z - 1.));
2060 :
2061 1207652 : case 2:
2062 1815582 : return .5*x*y*((1. - z)*(-2.) +
2063 1207652 : (-1.)*(2.*x + 2.*y - 2.*z - 3.));
2064 :
2065 1207652 : case 3:
2066 1815582 : return .5*(1. - x)*y*((1. - z)*(-2.) +
2067 1207652 : (-1.)*(2.*y - 2.*x - 2.*z - 1.));
2068 :
2069 1207652 : case 4:
2070 1815582 : return .5*(1. - x)*(1. - y)*(z*(2.) +
2071 1207652 : (1.)*(2.*z - 2.*x - 2.*y - 1.));
2072 :
2073 1207652 : case 5:
2074 1815582 : return .5*x*(1. - y)*(z*(2.) +
2075 1207652 : (1.)*(2.*x - 2.*y + 2.*z - 3.));
2076 :
2077 1207652 : case 6:
2078 1815582 : return .5*x*y*(z*(2.) +
2079 1207652 : (1.)*(2.*x + 2.*y + 2.*z - 5.));
2080 :
2081 1207652 : case 7:
2082 1815582 : return .5*(1. - x)*y*(z*(2.) +
2083 1207652 : (1.)*(2.*y - 2.*x + 2.*z - 3.));
2084 :
2085 1207652 : case 8:
2086 1207652 : return 2.*x*(1. - x)*(1. - y)*(-1.);
2087 :
2088 1207652 : case 9:
2089 1207652 : return 2.*x*y*(1. - y)*(-1.);
2090 :
2091 1207652 : case 10:
2092 1207652 : return 2.*x*(1. - x)*y*(-1.);
2093 :
2094 1207652 : case 11:
2095 1207652 : return 2.*(1. - x)*y*(1. - y)*(-1.);
2096 :
2097 1207652 : case 12:
2098 1207652 : return 2.*(1. - x)*(1. - y)*(1. - 2.*z);
2099 :
2100 1207652 : case 13:
2101 1207652 : return 2.*x*(1. - y)*(1. - 2.*z);
2102 :
2103 1207652 : case 14:
2104 1207652 : return 2.*x*y*(1. - 2.*z);
2105 :
2106 1207652 : case 15:
2107 1207652 : return 2.*(1. - x)*y*(1. - 2.*z);
2108 :
2109 1207652 : case 16:
2110 1207652 : return 2.*x*(1. - x)*(1. - y);
2111 :
2112 1207652 : case 17:
2113 1207652 : return 2.*x*y*(1. - y);
2114 :
2115 1207652 : case 18:
2116 1207652 : return 2.*x*(1. - x)*y;
2117 :
2118 1207652 : case 19:
2119 1207652 : 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 2034727209 : case HEX8:
2132 1152711 : libmesh_assert_msg(T == L2_LAGRANGE,
2133 : "High order on first order elements only supported for L2 families");
2134 : libmesh_fallthrough();
2135 679535406 : case HEX27:
2136 : {
2137 680990652 : libmesh_assert_less (i, 27);
2138 :
2139 : // Compute hex shape functions as a tensor-product
2140 2714565150 : const Real xi = p(0);
2141 2714565150 : const Real eta = p(1);
2142 2714565150 : 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 2714565150 : switch(j)
2153 : {
2154 904855050 : case 0:
2155 904855050 : return (fe_lagrange_1D_quadratic_shape_deriv(i0[i], 0, xi)*
2156 904855050 : fe_lagrange_1D_quadratic_shape (i1[i], eta)*
2157 1131660072 : fe_lagrange_1D_quadratic_shape (i2[i], zeta));
2158 :
2159 904855050 : case 1:
2160 1355908194 : return (fe_lagrange_1D_quadratic_shape (i0[i], xi)*
2161 904855050 : fe_lagrange_1D_quadratic_shape_deriv(i1[i], 0, eta)*
2162 1131660072 : fe_lagrange_1D_quadratic_shape (i2[i], zeta));
2163 :
2164 904855050 : case 2:
2165 904855050 : return (fe_lagrange_1D_quadratic_shape (i0[i], xi)*
2166 904855050 : fe_lagrange_1D_quadratic_shape (i1[i], eta)*
2167 1582713216 : 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 96321900 : case TET4:
2176 20880 : libmesh_assert_msg(T == L2_LAGRANGE,
2177 : "High order on first order elements only supported for L2 families");
2178 : libmesh_fallthrough();
2179 37050960 : case TET10:
2180 : case TET14:
2181 : {
2182 37555950 : libmesh_assert_less (i, 10);
2183 :
2184 : // Area coordinates, pg. 205, Vol. I, Carey, Oden, Becker FEM
2185 133856970 : const Real zeta1 = p(0);
2186 133856970 : const Real zeta2 = p(1);
2187 133856970 : const Real zeta3 = p(2);
2188 133856970 : const Real zeta0 = 1. - zeta1 - zeta2 - zeta3;
2189 :
2190 37555950 : const Real dzeta0dxi = -1.;
2191 37555950 : const Real dzeta1dxi = 1.;
2192 37555950 : const Real dzeta2dxi = 0.;
2193 37555950 : const Real dzeta3dxi = 0.;
2194 :
2195 37555950 : const Real dzeta0deta = -1.;
2196 37555950 : const Real dzeta1deta = 0.;
2197 37555950 : const Real dzeta2deta = 1.;
2198 37555950 : const Real dzeta3deta = 0.;
2199 :
2200 37555950 : const Real dzeta0dzeta = -1.;
2201 37555950 : const Real dzeta1dzeta = 0.;
2202 37555950 : const Real dzeta2dzeta = 0.;
2203 37555950 : const Real dzeta3dzeta = 1.;
2204 :
2205 133856970 : switch (j)
2206 : {
2207 : // d()/dxi
2208 44618990 : case 0:
2209 : {
2210 44618990 : switch(i)
2211 : {
2212 4461899 : case 0:
2213 4461899 : return (4.*zeta0 - 1.)*dzeta0dxi;
2214 :
2215 4461899 : case 1:
2216 4461899 : return (4.*zeta1 - 1.)*dzeta1dxi;
2217 :
2218 4461899 : case 2:
2219 4461899 : return (4.*zeta2 - 1.)*dzeta2dxi;
2220 :
2221 4461899 : case 3:
2222 4461899 : return (4.*zeta3 - 1.)*dzeta3dxi;
2223 :
2224 4461899 : case 4:
2225 4461899 : return 4.*(zeta0*dzeta1dxi + dzeta0dxi*zeta1);
2226 :
2227 4461899 : case 5:
2228 4461899 : return 4.*(zeta1*dzeta2dxi + dzeta1dxi*zeta2);
2229 :
2230 4461899 : case 6:
2231 4461899 : return 4.*(zeta0*dzeta2dxi + dzeta0dxi*zeta2);
2232 :
2233 4461899 : case 7:
2234 4461899 : return 4.*(zeta0*dzeta3dxi + dzeta0dxi*zeta3);
2235 :
2236 4461899 : case 8:
2237 4461899 : return 4.*(zeta1*dzeta3dxi + dzeta1dxi*zeta3);
2238 :
2239 4461899 : case 9:
2240 4461899 : 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 44618990 : case 1:
2249 : {
2250 44618990 : switch(i)
2251 : {
2252 4461899 : case 0:
2253 4461899 : return (4.*zeta0 - 1.)*dzeta0deta;
2254 :
2255 4461899 : case 1:
2256 4461899 : return (4.*zeta1 - 1.)*dzeta1deta;
2257 :
2258 4461899 : case 2:
2259 4461899 : return (4.*zeta2 - 1.)*dzeta2deta;
2260 :
2261 4461899 : case 3:
2262 4461899 : return (4.*zeta3 - 1.)*dzeta3deta;
2263 :
2264 4461899 : case 4:
2265 4461899 : return 4.*(zeta0*dzeta1deta + dzeta0deta*zeta1);
2266 :
2267 4461899 : case 5:
2268 4461899 : return 4.*(zeta1*dzeta2deta + dzeta1deta*zeta2);
2269 :
2270 4461899 : case 6:
2271 4461899 : return 4.*(zeta0*dzeta2deta + dzeta0deta*zeta2);
2272 :
2273 4461899 : case 7:
2274 4461899 : return 4.*(zeta0*dzeta3deta + dzeta0deta*zeta3);
2275 :
2276 4461899 : case 8:
2277 4461899 : return 4.*(zeta1*dzeta3deta + dzeta1deta*zeta3);
2278 :
2279 4461899 : case 9:
2280 4461899 : 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 44618990 : case 2:
2289 : {
2290 44618990 : switch(i)
2291 : {
2292 4461899 : case 0:
2293 4461899 : return (4.*zeta0 - 1.)*dzeta0dzeta;
2294 :
2295 4461899 : case 1:
2296 4461899 : return (4.*zeta1 - 1.)*dzeta1dzeta;
2297 :
2298 4461899 : case 2:
2299 4461899 : return (4.*zeta2 - 1.)*dzeta2dzeta;
2300 :
2301 4461899 : case 3:
2302 4461899 : return (4.*zeta3 - 1.)*dzeta3dzeta;
2303 :
2304 4461899 : case 4:
2305 4461899 : return 4.*(zeta0*dzeta1dzeta + dzeta0dzeta*zeta1);
2306 :
2307 4461899 : case 5:
2308 4461899 : return 4.*(zeta1*dzeta2dzeta + dzeta1dzeta*zeta2);
2309 :
2310 4461899 : case 6:
2311 4461899 : return 4.*(zeta0*dzeta2dzeta + dzeta0dzeta*zeta2);
2312 :
2313 4461899 : case 7:
2314 4461899 : return 4.*(zeta0*dzeta3dzeta + dzeta0dzeta*zeta3);
2315 :
2316 4461899 : case 8:
2317 4461899 : return 4.*(zeta1*dzeta3dzeta + dzeta1dzeta*zeta3);
2318 :
2319 4461899 : case 9:
2320 4461899 : 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 53932455 : case PRISM15:
2335 : {
2336 15202530 : libmesh_assert_less (i, 15);
2337 :
2338 53932455 : const Real xi = p(0);
2339 53932455 : const Real eta = p(1);
2340 53932455 : const Real zeta = p(2);
2341 :
2342 53932455 : switch (j)
2343 : {
2344 : // d()/dxi
2345 17977485 : case 0:
2346 : {
2347 17977485 : switch(i)
2348 : {
2349 1198499 : case 0:
2350 1198499 : return (2.*xi + 2.*eta + 0.5*zeta - 1.)*(1. - zeta);
2351 1198499 : case 1:
2352 1198499 : return (2.*xi - 1. - 0.5*zeta)*(1. - zeta);
2353 337834 : case 2:
2354 337834 : return 0.;
2355 1198499 : case 3:
2356 1198499 : return (2.*xi + 2.*eta - 0.5*zeta - 1.)*(1. + zeta);
2357 1198499 : case 4:
2358 1198499 : return (2.*xi - 1. + 0.5*zeta)*(1. + zeta);
2359 337834 : case 5:
2360 337834 : return 0.;
2361 1198499 : case 6:
2362 1198499 : return (4.*xi + 2.*eta - 2.)*(zeta - 1.);
2363 1198499 : case 7:
2364 1198499 : return -2.*(zeta - 1.)*eta;
2365 1198499 : case 8:
2366 1198499 : return 2.*(zeta - 1.)*eta;
2367 1198499 : case 9:
2368 1198499 : return (zeta - 1.)*(1. + zeta);
2369 1198499 : case 10:
2370 1198499 : return (1. - zeta)*(1. + zeta);
2371 337834 : case 11:
2372 337834 : return 0.;
2373 1198499 : case 12:
2374 1198499 : return (-4.*xi - 2.*eta + 2.)*(1. + zeta);
2375 1198499 : case 13:
2376 1198499 : return 2.*(1. + zeta)*eta;
2377 1198499 : case 14:
2378 1198499 : return -2.*(1. + zeta)*eta;
2379 0 : default:
2380 0 : libmesh_error_msg("Invalid i = " << i);
2381 : }
2382 : }
2383 :
2384 : // d()/deta
2385 17977485 : case 1:
2386 : {
2387 17977485 : switch(i)
2388 : {
2389 1198499 : case 0:
2390 1198499 : return (2.*xi + 2.*eta + 0.5*zeta - 1.)*(1. - zeta);
2391 337834 : case 1:
2392 337834 : return 0.;
2393 1198499 : case 2:
2394 1198499 : return (2.*eta - 1. - 0.5*zeta)*(1. - zeta);
2395 1198499 : case 3:
2396 1198499 : return (2.*xi + 2.*eta - 0.5*zeta - 1.)*(1. + zeta);
2397 337834 : case 4:
2398 337834 : return 0.;
2399 1198499 : case 5:
2400 1198499 : return (2.*eta - 1. + 0.5*zeta)*(1. + zeta);
2401 1198499 : case 6:
2402 1198499 : return 2.*(zeta - 1.)*xi;
2403 1198499 : case 7:
2404 1198499 : return 2.*(1. - zeta)*xi;
2405 1198499 : case 8:
2406 1198499 : return (2.*xi + 4.*eta - 2.)*(zeta - 1.);
2407 1198499 : case 9:
2408 1198499 : return (zeta - 1.)*(1. + zeta);
2409 337834 : case 10:
2410 337834 : return 0.;
2411 1198499 : case 11:
2412 1198499 : return (1. - zeta)*(1. + zeta);
2413 1198499 : case 12:
2414 1198499 : return -2.*(1. + zeta)*xi;
2415 1198499 : case 13:
2416 1198499 : return 2.*(1. + zeta)*xi;
2417 1198499 : case 14:
2418 1198499 : 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 17977485 : case 2:
2427 : {
2428 17977485 : switch(i)
2429 : {
2430 1198499 : case 0:
2431 1198499 : return (-xi - eta - zeta + 0.5)*(xi + eta - 1.);
2432 1198499 : case 1:
2433 1198499 : return -0.5*xi*(2.*xi - 1. - 2.*zeta);
2434 1198499 : case 2:
2435 1198499 : return -0.5*eta*(2.*eta - 1. - 2.*zeta);
2436 1198499 : case 3:
2437 1198499 : return (xi + eta - zeta - 0.5)*(xi + eta - 1.);
2438 1198499 : case 4:
2439 1198499 : return 0.5*xi*(2.*xi - 1. + 2.*zeta);
2440 1198499 : case 5:
2441 1198499 : return 0.5*eta*(2.*eta - 1. + 2.*zeta);
2442 1198499 : case 6:
2443 1198499 : return 2.*xi*(xi + eta - 1.);
2444 1198499 : case 7:
2445 1198499 : return -2.*xi*eta;
2446 1198499 : case 8:
2447 1198499 : return 2.*eta*(xi + eta - 1.);
2448 1198499 : case 9:
2449 1198499 : return 2.*zeta*(xi + eta - 1.);
2450 1198499 : case 10:
2451 1198499 : return -2.*xi*zeta;
2452 1198499 : case 11:
2453 1198499 : return -2.*eta*zeta;
2454 1198499 : case 12:
2455 1198499 : return 2.*xi*(1. - xi - eta);
2456 1198499 : case 13:
2457 1198499 : return 2.*xi*eta;
2458 1198499 : case 14:
2459 1198499 : 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 75477744 : case PRISM6:
2475 206712 : libmesh_assert_msg(T == L2_LAGRANGE,
2476 : "High order on first order elements only supported for L2 families");
2477 : libmesh_fallthrough();
2478 30511512 : case PRISM18:
2479 : case PRISM20:
2480 : case PRISM21:
2481 : {
2482 31338360 : libmesh_assert_less (i, 18);
2483 :
2484 : // Compute prism shape functions as a tensor-product
2485 : // of a triangle and an edge
2486 :
2487 106609392 : Point p2d(p(0),p(1));
2488 106609392 : 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 106609392 : switch (j)
2495 : {
2496 : // d()/dxi
2497 35536464 : case 0:
2498 35536464 : return (FE<2,LAGRANGE>::shape_deriv(TRI6, SECOND, i1[i], 0, p2d)*
2499 45984420 : fe_lagrange_1D_quadratic_shape(i0[i], p1d));
2500 :
2501 : // d()/deta
2502 35536464 : case 1:
2503 35536464 : return (FE<2,LAGRANGE>::shape_deriv(TRI6, SECOND, i1[i], 1, p2d)*
2504 45984420 : fe_lagrange_1D_quadratic_shape(i0[i], p1d));
2505 :
2506 : // d()/dzeta
2507 35536464 : case 2:
2508 35536464 : return (FE<2,LAGRANGE>::shape(TRI6, SECOND, i1[i], p2d)*
2509 60626808 : 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 12798162 : case PYRAMID13:
2520 : {
2521 3588468 : libmesh_assert_less (i, 13);
2522 :
2523 12798162 : const Real xi = p(0);
2524 12798162 : const Real eta = p(1);
2525 12798162 : const Real zeta = p(2);
2526 3588468 : const Real eps = 1.e-35;
2527 :
2528 : // Denominators are perturbed by epsilon to avoid
2529 : // divide-by-zero issues.
2530 : Real
2531 12798162 : den = (-1. + zeta + eps),
2532 12798162 : den2 = den*den,
2533 12798162 : xi2 = xi*xi,
2534 12798162 : eta2 = eta*eta,
2535 12798162 : zeta2 = zeta*zeta,
2536 12798162 : zeta3 = zeta2*zeta;
2537 :
2538 12798162 : switch (j)
2539 : {
2540 : // d/dxi
2541 4266054 : case 0:
2542 4266054 : switch(i)
2543 : {
2544 328158 : case 0:
2545 328158 : return 0.25*(-zeta - eta + 2.*eta*zeta - 2.*xi + 2.*zeta*xi + 2.*eta*xi + zeta2 + eta2)/den;
2546 :
2547 328158 : case 1:
2548 328158 : return -0.25*(-zeta - eta + 2.*eta*zeta + 2.*xi - 2.*zeta*xi - 2.*eta*xi + zeta2 + eta2)/den;
2549 :
2550 328158 : case 2:
2551 328158 : return -0.25*(-zeta + eta - 2.*eta*zeta + 2.*xi - 2.*zeta*xi + 2.*eta*xi + zeta2 + eta2)/den;
2552 :
2553 328158 : case 3:
2554 328158 : return 0.25*(-zeta + eta - 2.*eta*zeta - 2.*xi + 2.*zeta*xi - 2.*eta*xi + zeta2 + eta2)/den;
2555 :
2556 92012 : case 4:
2557 92012 : return 0.;
2558 :
2559 328158 : case 5:
2560 328158 : return -(-1. + eta + zeta)*xi/den;
2561 :
2562 328158 : case 6:
2563 328158 : return 0.5*(-1. + eta + zeta)*(1. + eta - zeta)/den;
2564 :
2565 328158 : case 7:
2566 328158 : return (1. + eta - zeta)*xi/den;
2567 :
2568 328158 : case 8:
2569 328158 : return -0.5*(-1. + eta + zeta)*(1. + eta - zeta)/den;
2570 :
2571 328158 : case 9:
2572 328158 : return -(-1. + eta + zeta)*zeta/den;
2573 :
2574 328158 : case 10:
2575 328158 : return (-1. + eta + zeta)*zeta/den;
2576 :
2577 328158 : case 11:
2578 328158 : return -(1. + eta - zeta)*zeta/den;
2579 :
2580 328158 : case 12:
2581 328158 : return (1. + eta - zeta)*zeta/den;
2582 :
2583 0 : default:
2584 0 : libmesh_error_msg("Invalid i = " << i);
2585 : }
2586 :
2587 : // d/deta
2588 4266054 : case 1:
2589 4266054 : switch(i)
2590 : {
2591 328158 : case 0:
2592 328158 : return 0.25*(-zeta - 2.*eta + 2.*eta*zeta - xi + 2.*zeta*xi + 2.*eta*xi + zeta2 + xi2)/den;
2593 :
2594 328158 : case 1:
2595 328158 : return -0.25*(zeta + 2.*eta - 2.*eta*zeta - xi + 2.*zeta*xi + 2.*eta*xi - zeta2 - xi2)/den;
2596 :
2597 328158 : case 2:
2598 328158 : return -0.25*(-zeta + 2.*eta - 2.*eta*zeta + xi - 2.*zeta*xi + 2.*eta*xi + zeta2 + xi2)/den;
2599 :
2600 328158 : case 3:
2601 328158 : return 0.25*(zeta - 2.*eta + 2.*eta*zeta + xi - 2.*zeta*xi + 2.*eta*xi - zeta2 - xi2)/den;
2602 :
2603 92012 : case 4:
2604 92012 : return 0.;
2605 :
2606 328158 : case 5:
2607 328158 : return -0.5*(-1. + xi + zeta)*(1. + xi - zeta)/den;
2608 :
2609 328158 : case 6:
2610 328158 : return (1. + xi - zeta)*eta/den;
2611 :
2612 328158 : case 7:
2613 328158 : return 0.5*(-1. + xi + zeta)*(1. + xi - zeta)/den;
2614 :
2615 328158 : case 8:
2616 328158 : return -(-1. + xi + zeta)*eta/den;
2617 :
2618 328158 : case 9:
2619 328158 : return -(-1. + xi + zeta)*zeta/den;
2620 :
2621 328158 : case 10:
2622 328158 : return (1. + xi - zeta)*zeta/den;
2623 :
2624 328158 : case 11:
2625 328158 : return -(1. + xi - zeta)*zeta/den;
2626 :
2627 328158 : case 12:
2628 328158 : return (-1. + xi + zeta)*zeta/den;
2629 :
2630 0 : default:
2631 0 : libmesh_error_msg("Invalid i = " << i);
2632 : }
2633 :
2634 : // d/dzeta
2635 4266054 : case 2:
2636 : {
2637 4266054 : switch(i)
2638 : {
2639 328158 : case 0:
2640 328158 : return -0.25*(xi + eta + 1.)*(-1. + 2.*zeta - zeta2 + eta*xi)/den2;
2641 :
2642 328158 : case 1:
2643 328158 : return 0.25*(eta - xi + 1.)*(1. - 2.*zeta + zeta2 + eta*xi)/den2;
2644 :
2645 328158 : case 2:
2646 328158 : return 0.25*(xi + eta - 1.)*(-1. + 2.*zeta - zeta2 + eta*xi)/den2;
2647 :
2648 328158 : case 3:
2649 328158 : return -0.25*(eta - xi - 1.)*(1. - 2.*zeta + zeta2 + eta*xi)/den2;
2650 :
2651 328158 : case 4:
2652 328158 : return 4.*zeta - 1.;
2653 :
2654 328158 : case 5:
2655 328158 : return 0.5*(-2 + eta + 6.*zeta + eta*xi2 + eta*zeta2 - 6.*zeta2 + 2.*zeta3 - 2.*eta*zeta)/den2;
2656 :
2657 328158 : case 6:
2658 328158 : return -0.5*(2 - 6.*zeta + xi + xi*zeta2 + eta2*xi + 6.*zeta2 - 2.*zeta3 - 2.*zeta*xi)/den2;
2659 :
2660 328158 : case 7:
2661 328158 : return -0.5*(2 + eta - 6.*zeta + eta*xi2 + eta*zeta2 + 6.*zeta2 - 2.*zeta3 - 2.*eta*zeta)/den2;
2662 :
2663 328158 : case 8:
2664 328158 : return 0.5*(-2 + 6.*zeta + xi + xi*zeta2 + eta2*xi - 6.*zeta2 + 2.*zeta3 - 2.*zeta*xi)/den2;
2665 :
2666 328158 : case 9:
2667 328158 : 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 328158 : case 10:
2670 328158 : 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 328158 : case 11:
2673 328158 : 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 328158 : case 12:
2676 328158 : 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 10374252 : 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 4040568 : case PYRAMID14:
2696 : case PYRAMID18:
2697 : {
2698 4040568 : libmesh_assert_less (i, 14);
2699 :
2700 14414820 : const Real xi = p(0);
2701 14414820 : const Real eta = p(1);
2702 14414820 : const Real zeta = p(2);
2703 4040568 : 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 14414820 : p1 = 0.5*(1. - eta - zeta), // back
2709 14414820 : p2 = 0.5*(1. + xi - zeta), // left
2710 14414820 : p3 = 0.5*(1. + eta - zeta), // front
2711 14414820 : p4 = 0.5*(1. - xi - zeta); // right
2712 :
2713 : // Denominators are perturbed by epsilon to avoid
2714 : // divide-by-zero issues.
2715 : Real
2716 14414820 : den = (-1. + zeta + eps),
2717 14414820 : den2 = den*den,
2718 14414820 : den3 = den2*den;
2719 :
2720 14414820 : switch (j)
2721 : {
2722 : // d/dxi
2723 4804940 : case 0:
2724 4804940 : switch(i)
2725 : {
2726 343210 : case 0:
2727 343210 : return 0.5*p1*(-xi*eta + zeta - zeta*zeta + 2.*p4*eta)/den2;
2728 :
2729 343210 : case 1:
2730 343210 : return -0.5*p1*(xi*eta + zeta - zeta*zeta + 2.*p2*eta)/den2;
2731 :
2732 343210 : case 2:
2733 343210 : return 0.5*p3*(xi*eta - zeta + zeta*zeta + 2.*p2*eta)/den2;
2734 :
2735 343210 : case 3:
2736 343210 : return -0.5*p3*(-xi*eta - zeta + zeta*zeta + 2.*p4*eta)/den2;
2737 :
2738 96204 : case 4:
2739 96204 : return 0.;
2740 :
2741 343210 : case 5:
2742 343210 : return 2.*p1*eta*xi/den2;
2743 :
2744 343210 : case 6:
2745 343210 : return 2.*p1*p3*(xi + 2.*p2)/den2;
2746 :
2747 343210 : case 7:
2748 343210 : return -2.*p3*eta*xi/den2;
2749 :
2750 343210 : case 8:
2751 343210 : return -2.*p1*p3*(-xi + 2.*p4)/den2;
2752 :
2753 343210 : case 9:
2754 343210 : return 2.*p1*zeta/den;
2755 :
2756 343210 : case 10:
2757 343210 : return -2.*p1*zeta/den;
2758 :
2759 343210 : case 11:
2760 343210 : return -2.*p3*zeta/den;
2761 :
2762 343210 : case 12:
2763 343210 : return 2.*p3*zeta/den;
2764 :
2765 343210 : case 13:
2766 343210 : return -8.*p1*p3*xi/den2;
2767 :
2768 0 : default:
2769 0 : libmesh_error_msg("Invalid i = " << i);
2770 : }
2771 :
2772 : // d/deta
2773 4804940 : case 1:
2774 4804940 : switch(i)
2775 : {
2776 343210 : case 0:
2777 343210 : return -0.5*p4*(xi*eta - zeta + zeta*zeta - 2.*p1*xi)/den2;
2778 :
2779 343210 : case 1:
2780 343210 : return 0.5*p2*(xi*eta + zeta - zeta*zeta - 2.*p1*xi)/den2;
2781 :
2782 343210 : case 2:
2783 343210 : return 0.5*p2*(xi*eta - zeta + zeta*zeta + 2.*p3*xi)/den2;
2784 :
2785 343210 : case 3:
2786 343210 : return -0.5*p4*(xi*eta + zeta - zeta*zeta + 2.*p3*xi)/den2;
2787 :
2788 96204 : case 4:
2789 96204 : return 0.;
2790 :
2791 343210 : case 5:
2792 343210 : return 2.*p2*p4*(eta - 2.*p1)/den2;
2793 :
2794 343210 : case 6:
2795 343210 : return -2.*p2*xi*eta/den2;
2796 :
2797 343210 : case 7:
2798 343210 : return 2.*p2*p4*(eta + 2.*p3)/den2;
2799 :
2800 343210 : case 8:
2801 343210 : return 2.*p4*xi*eta/den2;
2802 :
2803 343210 : case 9:
2804 343210 : return 2.*p4*zeta/den;
2805 :
2806 343210 : case 10:
2807 343210 : return 2.*p2*zeta/den;
2808 :
2809 343210 : case 11:
2810 343210 : return -2.*p2*zeta/den;
2811 :
2812 343210 : case 12:
2813 343210 : return -2.*p4*zeta/den;
2814 :
2815 343210 : case 13:
2816 343210 : 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 4804940 : case 2:
2825 : {
2826 4804940 : switch(i)
2827 : {
2828 343210 : case 0:
2829 343210 : return -0.5*p1*(xi*eta - zeta + zeta*zeta)/den2
2830 343210 : - 0.5*p4*(xi*eta - zeta + zeta*zeta)/den2
2831 343210 : + p4*p1*(2.*zeta - 1)/den2
2832 343210 : - 2.*p4*p1*(xi*eta - zeta + zeta*zeta)/den3;
2833 :
2834 343210 : case 1:
2835 343210 : return 0.5*p2*(xi*eta + zeta - zeta*zeta)/den2
2836 343210 : + 0.5*p1*(xi*eta + zeta - zeta*zeta)/den2
2837 343210 : - p1*p2*(1 - 2.*zeta)/den2
2838 343210 : + 2.*p1*p2*(xi*eta + zeta - zeta*zeta)/den3;
2839 :
2840 343210 : case 2:
2841 343210 : return -0.5*p3*(xi*eta - zeta + zeta*zeta)/den2
2842 343210 : - 0.5*p2*(xi*eta - zeta + zeta*zeta)/den2
2843 343210 : + p2*p3*(2.*zeta - 1)/den2
2844 343210 : - 2.*p2*p3*(xi*eta - zeta + zeta*zeta)/den3;
2845 :
2846 343210 : case 3:
2847 343210 : return 0.5*p4*(xi*eta + zeta - zeta*zeta)/den2
2848 343210 : + 0.5*p3*(xi*eta + zeta - zeta*zeta)/den2
2849 343210 : - p3*p4*(1 - 2.*zeta)/den2
2850 343210 : + 2.*p3*p4*(xi*eta + zeta - zeta*zeta)/den3;
2851 :
2852 343210 : case 4:
2853 343210 : return 4.*zeta - 1.;
2854 :
2855 343210 : case 5:
2856 343210 : return 2.*p4*p1*eta/den2
2857 343210 : + 2.*p2*p4*eta/den2
2858 343210 : + 2.*p1*p2*eta/den2
2859 343210 : + 8.*p2*p1*p4*eta/den3;
2860 :
2861 343210 : case 6:
2862 343210 : return -2.*p2*p3*xi/den2
2863 343210 : - 2.*p1*p3*xi/den2
2864 343210 : - 2.*p1*p2*xi/den2
2865 343210 : - 8.*p1*p2*p3*xi/den3;
2866 :
2867 343210 : case 7:
2868 343210 : return -2.*p3*p4*eta/den2
2869 343210 : - 2.*p2*p4*eta/den2
2870 343210 : - 2.*p2*p3*eta/den2
2871 343210 : - 8.*p2*p3*p4*eta/den3;
2872 :
2873 343210 : case 8:
2874 343210 : return 2.*p4*p1*xi/den2
2875 343210 : + 2.*p1*p3*xi/den2
2876 343210 : + 2.*p3*p4*xi/den2
2877 343210 : + 8.*p3*p4*p1*xi/den3;
2878 :
2879 343210 : case 9:
2880 343210 : return 2.*p4*zeta/den
2881 343210 : + 2.*p1*zeta/den
2882 343210 : - 4.*p1*p4/den
2883 343210 : + 4.*p1*p4*zeta/den2;
2884 :
2885 343210 : case 10:
2886 343210 : return 2.*p1*zeta/den
2887 343210 : + 2.*p2*zeta/den
2888 343210 : - 4.*p2*p1/den
2889 343210 : + 4.*p2*p1*zeta/den2;
2890 :
2891 343210 : case 11:
2892 343210 : return 2.*p2*zeta/den
2893 343210 : + 2.*p3*zeta/den
2894 343210 : - 4.*p3*p2/den
2895 343210 : + 4.*p3*p2*zeta/den2;
2896 :
2897 343210 : case 12:
2898 343210 : return 2.*p3*zeta/den
2899 343210 : + 2.*p4*zeta/den
2900 343210 : - 4.*p4*p3/den
2901 343210 : + 4.*p4*p3*zeta/den2;
2902 :
2903 343210 : case 13:
2904 343210 : return -8.*p2*p3*p4/den2
2905 343210 : - 8.*p3*p4*p1/den2
2906 343210 : - 8.*p2*p1*p4/den2
2907 343210 : - 8.*p1*p2*p3/den2
2908 343210 : - 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 1546935261 : case THIRD:
2927 : {
2928 1546935261 : switch (type)
2929 : {
2930 : // quadratic Lagrange shape functions with a cubic bubble
2931 1215501294 : case TET14:
2932 : {
2933 314983284 : libmesh_assert_less (i, 14);
2934 :
2935 : // Area coordinates, pg. 205, Vol. I, Carey, Oden, Becker FEM
2936 1215501294 : const Real zeta1 = p(0);
2937 1215501294 : const Real zeta2 = p(1);
2938 1215501294 : const Real zeta3 = p(2);
2939 1215501294 : const Real zeta0 = 1. - zeta1 - zeta2 - zeta3;
2940 :
2941 314983284 : const Real dzeta0dxi = -1.;
2942 314983284 : const Real dzeta1dxi = 1.;
2943 314983284 : const Real dzeta2dxi = 0.;
2944 314983284 : const Real dzeta3dxi = 0.;
2945 1215501294 : const Real dbubble012dxi = (zeta0-zeta1)*zeta2;
2946 1215501294 : const Real dbubble013dxi = (zeta0-zeta1)*zeta3;
2947 1215501294 : const Real dbubble123dxi = zeta2*zeta3;
2948 619799628 : const Real dbubble023dxi = -zeta2*zeta3;
2949 :
2950 314983284 : const Real dzeta0deta = -1.;
2951 314983284 : const Real dzeta1deta = 0.;
2952 314983284 : const Real dzeta2deta = 1.;
2953 314983284 : const Real dzeta3deta = 0.;
2954 1215501294 : const Real dbubble012deta = (zeta0-zeta2)*zeta1;
2955 1215501294 : const Real dbubble013deta = -zeta1*zeta3;
2956 619799628 : const Real dbubble123deta = zeta1*zeta3;
2957 1215501294 : const Real dbubble023deta = (zeta0-zeta2)*zeta3;
2958 :
2959 314983284 : const Real dzeta0dzeta = -1.;
2960 314983284 : const Real dzeta1dzeta = 0.;
2961 314983284 : const Real dzeta2dzeta = 0.;
2962 314983284 : const Real dzeta3dzeta = 1.;
2963 1215501294 : const Real dbubble012dzeta = -zeta1*zeta2;
2964 1215501294 : const Real dbubble013dzeta = (zeta0-zeta3)*zeta1;
2965 619799628 : const Real dbubble123dzeta = zeta1*zeta2;
2966 1215501294 : const Real dbubble023dzeta = (zeta0-zeta3)*zeta2;
2967 :
2968 1215501294 : switch (j)
2969 : {
2970 : // d()/dxi
2971 405167098 : case 0:
2972 : {
2973 405167098 : switch(i)
2974 : {
2975 28940507 : case 0:
2976 28940507 : return (4.*zeta0 - 1.)*dzeta0dxi + 3.*(dbubble012dxi+dbubble013dxi+dbubble023dxi);
2977 :
2978 28940507 : case 1:
2979 28940507 : return (4.*zeta1 - 1.)*dzeta1dxi + 3.*(dbubble012dxi+dbubble013dxi+dbubble123dxi);
2980 :
2981 28940507 : case 2:
2982 28940507 : return (4.*zeta2 - 1.)*dzeta2dxi + 3.*(dbubble012dxi+dbubble023dxi+dbubble123dxi);
2983 :
2984 28940507 : case 3:
2985 28940507 : return (4.*zeta3 - 1.)*dzeta3dxi + 3.*(dbubble013dxi+dbubble023dxi+dbubble123dxi);
2986 :
2987 28940507 : case 4:
2988 28940507 : return 4.*(zeta0*dzeta1dxi + dzeta0dxi*zeta1) - 12.*(dbubble012dxi+dbubble013dxi);
2989 :
2990 28940507 : case 5:
2991 28940507 : return 4.*(zeta1*dzeta2dxi + dzeta1dxi*zeta2) - 12.*(dbubble012dxi+dbubble123dxi);
2992 :
2993 28940507 : case 6:
2994 28940507 : return 4.*(zeta0*dzeta2dxi + dzeta0dxi*zeta2) - 12.*(dbubble012dxi+dbubble023dxi);
2995 :
2996 28940507 : case 7:
2997 28940507 : return 4.*(zeta0*dzeta3dxi + dzeta0dxi*zeta3) - 12.*(dbubble013dxi+dbubble023dxi);
2998 :
2999 28940507 : case 8:
3000 28940507 : return 4.*(zeta1*dzeta3dxi + dzeta1dxi*zeta3) - 12.*(dbubble013dxi+dbubble123dxi);
3001 :
3002 28940507 : case 9:
3003 28940507 : return 4.*(zeta2*dzeta3dxi + dzeta2dxi*zeta3) - 12.*(dbubble023dxi+dbubble123dxi);
3004 :
3005 28940507 : case 10:
3006 28940507 : return 27.*dbubble012dxi;
3007 :
3008 28940507 : case 11:
3009 28940507 : return 27.*dbubble013dxi;
3010 :
3011 28940507 : case 12:
3012 28940507 : return 27.*dbubble123dxi;
3013 :
3014 28940507 : case 13:
3015 28940507 : return 27.*dbubble023dxi;
3016 :
3017 0 : default:
3018 0 : libmesh_error_msg("Invalid i = " << i);
3019 : }
3020 : }
3021 :
3022 : // d()/deta
3023 405167098 : case 1:
3024 : {
3025 405167098 : switch(i)
3026 : {
3027 28940507 : case 0:
3028 28940507 : return (4.*zeta0 - 1.)*dzeta0deta + 3.*(dbubble012deta+dbubble013deta+dbubble023deta);;
3029 :
3030 28940507 : case 1:
3031 28940507 : return (4.*zeta1 - 1.)*dzeta1deta + 3.*(dbubble012deta+dbubble013deta+dbubble123deta);
3032 :
3033 28940507 : case 2:
3034 28940507 : return (4.*zeta2 - 1.)*dzeta2deta + 3.*(dbubble012deta+dbubble023deta+dbubble123deta);
3035 :
3036 28940507 : case 3:
3037 28940507 : return (4.*zeta3 - 1.)*dzeta3deta + 3.*(dbubble013deta+dbubble023deta+dbubble123deta);
3038 :
3039 28940507 : case 4:
3040 28940507 : return 4.*(zeta0*dzeta1deta + dzeta0deta*zeta1) - 12.*(dbubble012deta+dbubble013deta);
3041 :
3042 28940507 : case 5:
3043 28940507 : return 4.*(zeta1*dzeta2deta + dzeta1deta*zeta2) - 12.*(dbubble012deta+dbubble123deta);
3044 :
3045 28940507 : case 6:
3046 28940507 : return 4.*(zeta0*dzeta2deta + dzeta0deta*zeta2) - 12.*(dbubble012deta+dbubble023deta);
3047 :
3048 28940507 : case 7:
3049 28940507 : return 4.*(zeta0*dzeta3deta + dzeta0deta*zeta3) - 12.*(dbubble013deta+dbubble023deta);
3050 :
3051 28940507 : case 8:
3052 28940507 : return 4.*(zeta1*dzeta3deta + dzeta1deta*zeta3) - 12.*(dbubble013deta+dbubble123deta);
3053 :
3054 28940507 : case 9:
3055 28940507 : return 4.*(zeta2*dzeta3deta + dzeta2deta*zeta3) - 12.*(dbubble023deta+dbubble123deta);
3056 :
3057 28940507 : case 10:
3058 28940507 : return 27.*dbubble012deta;
3059 :
3060 28940507 : case 11:
3061 28940507 : return 27.*dbubble013deta;
3062 :
3063 28940507 : case 12:
3064 28940507 : return 27.*dbubble123deta;
3065 :
3066 28940507 : case 13:
3067 28940507 : return 27.*dbubble023deta;
3068 :
3069 0 : default:
3070 0 : libmesh_error_msg("Invalid i = " << i);
3071 : }
3072 : }
3073 :
3074 : // d()/dzeta
3075 405167098 : case 2:
3076 : {
3077 405167098 : switch(i)
3078 : {
3079 28940507 : case 0:
3080 28940507 : return (4.*zeta0 - 1.)*dzeta0dzeta + 3.*(dbubble012dzeta+dbubble013dzeta+dbubble023dzeta);
3081 :
3082 28940507 : case 1:
3083 28940507 : return (4.*zeta1 - 1.)*dzeta1dzeta + 3.*(dbubble012dzeta+dbubble013dzeta+dbubble123dzeta);
3084 :
3085 28940507 : case 2:
3086 28940507 : return (4.*zeta2 - 1.)*dzeta2dzeta + 3.*(dbubble012dzeta+dbubble023dzeta+dbubble123dzeta);
3087 :
3088 28940507 : case 3:
3089 28940507 : return (4.*zeta3 - 1.)*dzeta3dzeta + 3.*(dbubble013dzeta+dbubble023dzeta+dbubble123dzeta);
3090 :
3091 28940507 : case 4:
3092 28940507 : return 4.*(zeta0*dzeta1dzeta + dzeta0dzeta*zeta1) - 12.*(dbubble012dzeta+dbubble013dzeta);
3093 :
3094 28940507 : case 5:
3095 28940507 : return 4.*(zeta1*dzeta2dzeta + dzeta1dzeta*zeta2) - 12.*(dbubble012dzeta+dbubble123dzeta);
3096 :
3097 28940507 : case 6:
3098 28940507 : return 4.*(zeta0*dzeta2dzeta + dzeta0dzeta*zeta2) - 12.*(dbubble012dzeta+dbubble023dzeta);
3099 :
3100 28940507 : case 7:
3101 28940507 : return 4.*(zeta0*dzeta3dzeta + dzeta0dzeta*zeta3) - 12.*(dbubble013dzeta+dbubble023dzeta);
3102 :
3103 28940507 : case 8:
3104 28940507 : return 4.*(zeta1*dzeta3dzeta + dzeta1dzeta*zeta3) - 12.*(dbubble013dzeta+dbubble123dzeta);
3105 :
3106 28940507 : case 9:
3107 28940507 : return 4.*(zeta2*dzeta3dzeta + dzeta2dzeta*zeta3) - 12.*(dbubble023dzeta+dbubble123dzeta);
3108 :
3109 28940507 : case 10:
3110 28940507 : return 27.*dbubble012dzeta;
3111 :
3112 28940507 : case 11:
3113 28940507 : return 27.*dbubble013dzeta;
3114 :
3115 28940507 : case 12:
3116 28940507 : return 27.*dbubble123dzeta;
3117 :
3118 28940507 : case 13:
3119 28940507 : 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 83879340 : case PRISM20:
3132 : {
3133 23297880 : 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 83879340 : Point p2d(p(0),p(1));
3139 83879340 : Real p1d = p(2);
3140 :
3141 83879340 : const Real mainval = FE<3,LAGRANGE>::shape_deriv(PRISM21, THIRD, i, j, p);
3142 :
3143 83879340 : if (i0[i] != 2)
3144 16308516 : return mainval;
3145 :
3146 6989364 : Real bubbleval = 0;
3147 :
3148 6989364 : switch (j)
3149 : {
3150 : // d()/dxi
3151 8387934 : case 0:
3152 8387934 : bubbleval =
3153 8387934 : FE<2,LAGRANGE>::shape_deriv(TRI7, THIRD, 6, 0, p2d)*
3154 2329788 : fe_lagrange_1D_quadratic_shape(2, p1d);
3155 8387934 : break;
3156 :
3157 : // d()/deta
3158 8387934 : case 1:
3159 8387934 : bubbleval =
3160 8387934 : FE<2,LAGRANGE>::shape_deriv(TRI7, THIRD, 6, 1, p2d)*
3161 2329788 : fe_lagrange_1D_quadratic_shape(2, p1d);
3162 8387934 : break;
3163 :
3164 : // d()/dzeta
3165 8387934 : case 2:
3166 8387934 : bubbleval =
3167 8387934 : FE<2,LAGRANGE>::shape(TRI7, THIRD, 6, p2d)*
3168 2329788 : fe_lagrange_1D_quadratic_shape_deriv(2, 0, p1d);
3169 8387934 : break;
3170 :
3171 0 : default:
3172 0 : libmesh_error_msg("Invalid shape function derivative j = " << j);
3173 : }
3174 :
3175 25163802 : if (i < 12) // vertices
3176 12581901 : return mainval - bubbleval / 9;
3177 :
3178 12581901 : return mainval + bubbleval * (Real(4) / 9);
3179 : }
3180 :
3181 208117293 : case PRISM21:
3182 : {
3183 57572652 : libmesh_assert_less (i, 21);
3184 :
3185 : // Compute prism shape functions as a tensor-product
3186 : // of a triangle and an edge
3187 :
3188 208117293 : Point p2d(p(0),p(1));
3189 208117293 : 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 208117293 : switch (j)
3196 : {
3197 : // d()/dxi
3198 69372431 : case 0:
3199 69372431 : return (FE<2,LAGRANGE>::shape_deriv(TRI7, THIRD, i1[i], 0, p2d)*
3200 88556061 : fe_lagrange_1D_quadratic_shape(i0[i], p1d));
3201 :
3202 : // d()/deta
3203 69372431 : case 1:
3204 69372431 : return (FE<2,LAGRANGE>::shape_deriv(TRI7, THIRD, i1[i], 1, p2d)*
3205 88556061 : fe_lagrange_1D_quadratic_shape(i0[i], p1d));
3206 :
3207 : // d()/dzeta
3208 69372431 : case 2:
3209 69372431 : return (FE<2,LAGRANGE>::shape(TRI7, THIRD, i1[i], p2d)*
3210 119553978 : 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 39437334 : case PYRAMID18:
3218 : {
3219 39437334 : 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 711781206 : 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 174912924 : libmesh_assert_less (j, 6);
3253 :
3254 711781206 : switch (order)
3255 : {
3256 : // linear Lagrange shape functions
3257 59996400 : case FIRST:
3258 : {
3259 59996400 : switch (type)
3260 : {
3261 : // Linear tets have all second derivatives = 0
3262 191280 : case TET4:
3263 : case TET10:
3264 : case TET14:
3265 : {
3266 191280 : 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 5447556 : case PRISM6:
3274 : case PRISM15:
3275 : case PRISM18:
3276 : case PRISM20:
3277 : case PRISM21:
3278 : {
3279 1514016 : libmesh_assert_less (i, 6);
3280 :
3281 : // Compute prism shape functions as a tensor-product
3282 : // of a triangle and an edge
3283 :
3284 5447556 : Point p2d(p(0),p(1));
3285 1514016 : 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 5447556 : switch (j)
3292 : {
3293 : // All repeated second derivatives and the xi-eta derivative are zero on PRISMs
3294 1009344 : 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 1009344 : return 0.;
3300 : }
3301 :
3302 907926 : case 3: // d^2()/dxidzeta
3303 907926 : return (FE<2,LAGRANGE>::shape_deriv(TRI3, FIRST, i1[i], 0, p2d)*
3304 1235721 : fe_lagrange_1D_linear_shape_deriv(i0[i], 0, p1d));
3305 :
3306 907926 : case 4: // d^2()/detadzeta
3307 907926 : return (FE<2,LAGRANGE>::shape_deriv(TRI3, FIRST, i1[i], 1, p2d)*
3308 1235721 : 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 269340 : case PYRAMID5:
3316 : case PYRAMID13:
3317 : case PYRAMID14:
3318 : case PYRAMID18:
3319 : {
3320 72240 : libmesh_assert_less (i, 5);
3321 :
3322 269340 : const Real xi = p(0);
3323 269340 : const Real eta = p(1);
3324 269340 : const Real zeta = p(2);
3325 72240 : const Real eps = 1.e-35;
3326 :
3327 269340 : switch (j)
3328 : {
3329 : // xi-xi and eta-eta derivatives are all zero for PYRAMID5.
3330 24080 : case 0: // d^2()/dxi^2
3331 : case 2: // d^2()/deta^2
3332 24080 : return 0.;
3333 :
3334 44890 : case 1: // d^2()/dxideta
3335 : {
3336 44890 : switch (i)
3337 : {
3338 17956 : case 0:
3339 : case 2:
3340 17956 : return 0.25/(1. - zeta + eps);
3341 17956 : case 1:
3342 : case 3:
3343 17956 : return -0.25/(1. - zeta + eps);
3344 2408 : case 4:
3345 2408 : return 0.;
3346 0 : default:
3347 0 : libmesh_error_msg("Invalid i = " << i);
3348 : }
3349 : }
3350 :
3351 44890 : case 3: // d^2()/dxidzeta
3352 : {
3353 44890 : Real den = (1. - zeta + eps)*(1. - zeta + eps);
3354 :
3355 44890 : switch (i)
3356 : {
3357 17956 : case 0:
3358 : case 2:
3359 17956 : return 0.25*eta/den;
3360 17956 : case 1:
3361 : case 3:
3362 17956 : return -0.25*eta/den;
3363 2408 : case 4:
3364 2408 : return 0.;
3365 0 : default:
3366 0 : libmesh_error_msg("Invalid i = " << i);
3367 : }
3368 : }
3369 :
3370 44890 : case 4: // d^2()/detadzeta
3371 : {
3372 44890 : Real den = (1. - zeta + eps)*(1. - zeta + eps);
3373 :
3374 44890 : switch (i)
3375 : {
3376 17956 : case 0:
3377 : case 2:
3378 17956 : return 0.25*xi/den;
3379 17956 : case 1:
3380 : case 3:
3381 17956 : return -0.25*xi/den;
3382 2408 : case 4:
3383 2408 : return 0.;
3384 0 : default:
3385 0 : libmesh_error_msg("Invalid i = " << i);
3386 : }
3387 : }
3388 :
3389 44890 : case 5: // d^2()/dzeta^2
3390 : {
3391 44890 : Real den = (1. - zeta + eps)*(1. - zeta + eps)*(1. - zeta + eps);
3392 :
3393 44890 : switch (i)
3394 : {
3395 17956 : case 0:
3396 : case 2:
3397 17956 : return 0.5*xi*eta/den;
3398 17956 : case 1:
3399 : case 3:
3400 17956 : return -0.5*xi*eta/den;
3401 2408 : case 4:
3402 2408 : 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 48069600 : case HEX8:
3415 : case HEX20:
3416 : case HEX27:
3417 : {
3418 12020592 : libmesh_assert_less (i, 8);
3419 :
3420 : // Compute hex shape functions as a tensor-product
3421 48069600 : const Real xi = p(0);
3422 48069600 : const Real eta = p(1);
3423 48069600 : 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 48069600 : switch (j)
3430 : {
3431 : // All repeated second derivatives are zero on HEX8
3432 6010296 : case 0: // d^2()/dxi^2
3433 : case 2: // d^2()/deta^2
3434 : case 5: // d^2()/dzeta^2
3435 : {
3436 6010296 : return 0.;
3437 : }
3438 :
3439 8011600 : case 1: // d^2()/dxideta
3440 8011600 : return (fe_lagrange_1D_linear_shape_deriv(i0[i], 0, xi)*
3441 8011600 : fe_lagrange_1D_linear_shape_deriv(i1[i], 0, eta)*
3442 14019768 : fe_lagrange_1D_linear_shape (i2[i], zeta));
3443 :
3444 8011600 : case 3: // d^2()/dxidzeta
3445 8011600 : return (fe_lagrange_1D_linear_shape_deriv(i0[i], 0, xi)*
3446 8011600 : fe_lagrange_1D_linear_shape (i1[i], eta)*
3447 11015684 : fe_lagrange_1D_linear_shape_deriv(i2[i], 0, zeta));
3448 :
3449 8011600 : case 4: // d^2()/detadzeta
3450 8011600 : return (fe_lagrange_1D_linear_shape (i0[i], xi)*
3451 8011600 : fe_lagrange_1D_linear_shape_deriv(i1[i], 0, eta)*
3452 11015684 : 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 317275128 : case SECOND:
3475 : {
3476 317275128 : switch (type)
3477 : {
3478 :
3479 : // serendipity hexahedral quadratic shape functions
3480 24134400 : case HEX20:
3481 : {
3482 6033600 : libmesh_assert_less (i, 20);
3483 :
3484 24134400 : const Real xi = p(0);
3485 24134400 : const Real eta = p(1);
3486 24134400 : 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 24134400 : const Real x = .5*(xi + 1.);
3491 24134400 : const Real y = .5*(eta + 1.);
3492 24134400 : const Real z = .5*(zeta + 1.);
3493 :
3494 24134400 : switch(j)
3495 : {
3496 4022400 : case 0: // d^2()/dxi^2
3497 : {
3498 4022400 : switch(i)
3499 : {
3500 402240 : case 0:
3501 : case 1:
3502 402240 : return (1. - y) * (1. - z);
3503 402240 : case 2:
3504 : case 3:
3505 402240 : return y * (1. - z);
3506 402240 : case 4:
3507 : case 5:
3508 402240 : return (1. - y) * z;
3509 402240 : case 6:
3510 : case 7:
3511 402240 : return y * z;
3512 201120 : case 8:
3513 201120 : return -2. * (1. - y) * (1. - z);
3514 201120 : case 10:
3515 201120 : return -2. * y * (1. - z);
3516 201120 : case 16:
3517 201120 : return -2. * (1. - y) * z;
3518 201120 : case 18:
3519 201120 : return -2. * y * z;
3520 402240 : case 9:
3521 : case 11:
3522 : case 12:
3523 : case 13:
3524 : case 14:
3525 : case 15:
3526 : case 17:
3527 : case 19:
3528 402240 : return 0;
3529 0 : default:
3530 0 : libmesh_error_msg("Invalid i = " << i);
3531 : }
3532 : }
3533 4022400 : case 1: // d^2()/dxideta
3534 : {
3535 4022400 : switch(i)
3536 : {
3537 201120 : case 0:
3538 201120 : return (1.25 - x - y - .5*z) * (1. - z);
3539 201120 : case 1:
3540 201120 : return (-x + y + .5*z - .25) * (1. - z);
3541 201120 : case 2:
3542 201120 : return (x + y - .5*z - .75) * (1. - z);
3543 201120 : case 3:
3544 201120 : return (-y + x + .5*z - .25) * (1. - z);
3545 201120 : case 4:
3546 201120 : return -.25*z * (4.*x + 4.*y - 2.*z - 3);
3547 201120 : case 5:
3548 201120 : return -.25*z * (-4.*y + 4.*x + 2.*z - 1.);
3549 201120 : case 6:
3550 201120 : return .25*z * (-5 + 4.*x + 4.*y + 2.*z);
3551 201120 : case 7:
3552 201120 : return .25*z * (4.*x - 4.*y - 2.*z + 1.);
3553 201120 : case 8:
3554 201120 : return (-1. + 2.*x) * (1. - z);
3555 201120 : case 9:
3556 201120 : return (1. - 2.*y) * (1. - z);
3557 201120 : case 10:
3558 201120 : return (1. - 2.*x) * (1. - z);
3559 201120 : case 11:
3560 201120 : return (-1. + 2.*y) * (1. - z);
3561 201120 : case 12:
3562 201120 : return z * (1. - z);
3563 201120 : case 13:
3564 201120 : return -z * (1. - z);
3565 201120 : case 14:
3566 201120 : return z * (1. - z);
3567 201120 : case 15:
3568 201120 : return -z * (1. - z);
3569 201120 : case 16:
3570 201120 : return (-1. + 2.*x) * z;
3571 201120 : case 17:
3572 201120 : return (1. - 2.*y) * z;
3573 201120 : case 18:
3574 201120 : return (1. - 2.*x) * z;
3575 201120 : case 19:
3576 201120 : return (-1. + 2.*y) * z;
3577 0 : default:
3578 0 : libmesh_error_msg("Invalid i = " << i);
3579 : }
3580 : }
3581 4022400 : case 2: // d^2()/deta^2
3582 4022400 : switch(i)
3583 : {
3584 402240 : case 0:
3585 : case 3:
3586 402240 : return (1. - x) * (1. - z);
3587 402240 : case 1:
3588 : case 2:
3589 402240 : return x * (1. - z);
3590 402240 : case 4:
3591 : case 7:
3592 402240 : return (1. - x) * z;
3593 402240 : case 5:
3594 : case 6:
3595 402240 : return x * z;
3596 201120 : case 9:
3597 201120 : return -2. * x * (1. - z);
3598 201120 : case 11:
3599 201120 : return -2. * (1. - x) * (1. - z);
3600 201120 : case 17:
3601 201120 : return -2. * x * z;
3602 201120 : case 19:
3603 201120 : return -2. * (1. - x) * z;
3604 402240 : case 8:
3605 : case 10:
3606 : case 12:
3607 : case 13:
3608 : case 14:
3609 : case 15:
3610 : case 16:
3611 : case 18:
3612 402240 : return 0.;
3613 0 : default:
3614 0 : libmesh_error_msg("Invalid i = " << i);
3615 : }
3616 4022400 : case 3: // d^2()/dxidzeta
3617 4022400 : switch(i)
3618 : {
3619 201120 : case 0:
3620 201120 : return (1.25 - x - .5*y - z) * (1. - y);
3621 201120 : case 1:
3622 201120 : return (-x + .5*y + z - .25) * (1. - y);
3623 201120 : case 2:
3624 201120 : return -.25*y * (2.*y + 4.*x - 4.*z - 1.);
3625 201120 : case 3:
3626 201120 : return -.25*y * (-2.*y + 4.*x + 4.*z - 3);
3627 201120 : case 4:
3628 201120 : return (-z + x + .5*y - .25) * (1. - y);
3629 201120 : case 5:
3630 201120 : return (x - .5*y + z - .75) * (1. - y);
3631 201120 : case 6:
3632 201120 : return .25*y * (2.*y + 4.*x + 4.*z - 5);
3633 201120 : case 7:
3634 201120 : return .25*y * (-2.*y + 4.*x - 4.*z + 1.);
3635 201120 : case 8:
3636 201120 : return (-1. + 2.*x) * (1. - y);
3637 201120 : case 9:
3638 201120 : return -y * (1. - y);
3639 201120 : case 10:
3640 201120 : return (-1. + 2.*x) * y;
3641 201120 : case 11:
3642 201120 : return y * (1. - y);
3643 201120 : case 12:
3644 201120 : return (-1. + 2.*z) * (1. - y);
3645 201120 : case 13:
3646 201120 : return (1. - 2.*z) * (1. - y);
3647 201120 : case 14:
3648 201120 : return (1. - 2.*z) * y;
3649 201120 : case 15:
3650 201120 : return (-1. + 2.*z) * y;
3651 201120 : case 16:
3652 201120 : return (1. - 2.*x) * (1. - y);
3653 201120 : case 17:
3654 201120 : return y * (1. - y);
3655 201120 : case 18:
3656 201120 : return (1. - 2.*x) * y;
3657 201120 : case 19:
3658 201120 : return -y * (1. - y);
3659 0 : default:
3660 0 : libmesh_error_msg("Invalid i = " << i);
3661 : }
3662 4022400 : case 4: // d^2()/detadzeta
3663 4022400 : switch(i)
3664 : {
3665 201120 : case 0:
3666 201120 : return (1.25 - .5*x - y - z) * (1. - x);
3667 201120 : case 1:
3668 201120 : return .25*x * (2.*x - 4.*y - 4.*z + 3.);
3669 201120 : case 2:
3670 201120 : return -.25*x * (2.*x + 4.*y - 4.*z - 1.);
3671 201120 : case 3:
3672 201120 : return (-y + .5*x + z - .25) * (1. - x);
3673 201120 : case 4:
3674 201120 : return (-z + .5*x + y - .25) * (1. - x);
3675 201120 : case 5:
3676 201120 : return -.25*x * (2.*x - 4.*y + 4.*z - 1.);
3677 201120 : case 6:
3678 201120 : return .25*x * (2.*x + 4.*y + 4.*z - 5);
3679 201120 : case 7:
3680 201120 : return (y - .5*x + z - .75) * (1. - x);
3681 201120 : case 8:
3682 201120 : return x * (1. - x);
3683 201120 : case 9:
3684 201120 : return (-1. + 2.*y) * x;
3685 201120 : case 10:
3686 201120 : return -x * (1. - x);
3687 201120 : case 11:
3688 201120 : return (-1. + 2.*y) * (1. - x);
3689 201120 : case 12:
3690 201120 : return (-1. + 2.*z) * (1. - x);
3691 201120 : case 13:
3692 201120 : return (-1. + 2.*z) * x;
3693 201120 : case 14:
3694 201120 : return (1. - 2.*z) * x;
3695 201120 : case 15:
3696 201120 : return (1. - 2.*z) * (1. - x);
3697 201120 : case 16:
3698 201120 : return -x * (1. - x);
3699 201120 : case 17:
3700 201120 : return (1. - 2.*y) * x;
3701 201120 : case 18:
3702 201120 : return x * (1. - x);
3703 201120 : case 19:
3704 201120 : return (1. - 2.*y) * (1. - x);
3705 0 : default:
3706 0 : libmesh_error_msg("Invalid i = " << i);
3707 : }
3708 4022400 : case 5: // d^2()/dzeta^2
3709 4022400 : switch(i)
3710 : {
3711 402240 : case 0:
3712 : case 4:
3713 402240 : return (1. - x) * (1. - y);
3714 402240 : case 1:
3715 : case 5:
3716 402240 : return x * (1. - y);
3717 402240 : case 2:
3718 : case 6:
3719 402240 : return x * y;
3720 402240 : case 3:
3721 : case 7:
3722 402240 : return (1. - x) * y;
3723 201120 : case 12:
3724 201120 : return -2. * (1. - x) * (1. - y);
3725 201120 : case 13:
3726 201120 : return -2. * x * (1. - y);
3727 201120 : case 14:
3728 201120 : return -2. * x * y;
3729 201120 : case 15:
3730 201120 : return -2. * (1. - x) * y;
3731 402240 : case 8:
3732 : case 9:
3733 : case 10:
3734 : case 11:
3735 : case 16:
3736 : case 17:
3737 : case 18:
3738 : case 19:
3739 402240 : 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 190427436 : case HEX8:
3750 625806 : libmesh_assert_msg(T == L2_LAGRANGE,
3751 : "High order on first order elements only supported for L2 families");
3752 : libmesh_fallthrough();
3753 55906200 : case HEX27:
3754 : {
3755 57137076 : libmesh_assert_less (i, 27);
3756 :
3757 : // Compute hex shape functions as a tensor-product
3758 246938706 : const Real xi = p(0);
3759 246938706 : const Real eta = p(1);
3760 246938706 : 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 246938706 : switch(j)
3771 : {
3772 : // d^2()/dxi^2
3773 41156451 : case 0:
3774 41156451 : return (fe_lagrange_1D_quadratic_shape_second_deriv(i0[i], 0, xi)*
3775 41156451 : fe_lagrange_1D_quadratic_shape (i1[i], eta)*
3776 50679297 : fe_lagrange_1D_quadratic_shape (i2[i], zeta));
3777 :
3778 : // d^2()/dxideta
3779 41156451 : case 1:
3780 41156451 : return (fe_lagrange_1D_quadratic_shape_deriv(i0[i], 0, xi)*
3781 41156451 : fe_lagrange_1D_quadratic_shape_deriv(i1[i], 0, eta)*
3782 50679297 : fe_lagrange_1D_quadratic_shape (i2[i], zeta));
3783 :
3784 : // d^2()/deta^2
3785 41156451 : case 2:
3786 63267210 : return (fe_lagrange_1D_quadratic_shape (i0[i], xi)*
3787 41156451 : fe_lagrange_1D_quadratic_shape_second_deriv(i1[i], 0, eta)*
3788 50679297 : fe_lagrange_1D_quadratic_shape (i2[i], zeta));
3789 :
3790 : // d^2()/dxidzeta
3791 41156451 : case 3:
3792 41156451 : return (fe_lagrange_1D_quadratic_shape_deriv(i0[i], 0, xi)*
3793 41156451 : fe_lagrange_1D_quadratic_shape (i1[i], eta)*
3794 72790056 : fe_lagrange_1D_quadratic_shape_deriv(i2[i], 0, zeta));
3795 :
3796 : // d^2()/detadzeta
3797 41156451 : case 4:
3798 63267210 : return (fe_lagrange_1D_quadratic_shape (i0[i], xi)*
3799 41156451 : fe_lagrange_1D_quadratic_shape_deriv(i1[i], 0, eta)*
3800 72790056 : fe_lagrange_1D_quadratic_shape_deriv(i2[i], 0, zeta));
3801 :
3802 : // d^2()/dzeta^2
3803 41156451 : case 5:
3804 41156451 : return (fe_lagrange_1D_quadratic_shape (i0[i], xi)*
3805 41156451 : fe_lagrange_1D_quadratic_shape (i1[i], eta)*
3806 51700986 : 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 1584060 : case TET4:
3815 41760 : libmesh_assert_msg(T == L2_LAGRANGE,
3816 : "High order on first order elements only supported for L2 families");
3817 : libmesh_fallthrough();
3818 529560 : 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 2160240 : const unsigned int my_j = independent_var_indices[j][0];
3867 2160240 : const unsigned int my_k = independent_var_indices[j][1];
3868 :
3869 2160240 : if (i<4)
3870 : {
3871 864096 : return 4.*dzetadxi[i][my_j]*dzetadxi[i][my_k];
3872 : }
3873 :
3874 1296144 : else if (i<10)
3875 : {
3876 1296144 : const unsigned short int my_m = zeta_indices[i][0];
3877 1296144 : const unsigned short int my_n = zeta_indices[i][1];
3878 :
3879 2037672 : return 4.*(dzetadxi[my_n][my_j]*dzetadxi[my_m][my_k] +
3880 1296144 : 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 14292990 : case PRISM15:
3890 : {
3891 3977640 : libmesh_assert_less (i, 15);
3892 :
3893 14292990 : const Real xi = p(0);
3894 14292990 : const Real eta = p(1);
3895 14292990 : const Real zeta = p(2);
3896 :
3897 14292990 : switch (j)
3898 : {
3899 : // d^2()/dxi^2
3900 2382165 : case 0:
3901 : {
3902 2382165 : switch(i)
3903 : {
3904 317622 : case 0:
3905 : case 1:
3906 317622 : return 2.*(1. - zeta);
3907 397764 : 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 397764 : return 0.;
3917 317622 : case 3:
3918 : case 4:
3919 317622 : return 2.*(1. + zeta);
3920 158811 : case 6:
3921 158811 : return 4.*(zeta - 1);
3922 158811 : case 12:
3923 158811 : return -4.*(1. + zeta);
3924 0 : default:
3925 0 : libmesh_error_msg("Invalid i = " << i);
3926 : }
3927 : }
3928 :
3929 : // d^2()/dxideta
3930 2382165 : case 1:
3931 : {
3932 2382165 : switch(i)
3933 : {
3934 317622 : case 0:
3935 : case 7:
3936 317622 : return 2.*(1. - zeta);
3937 309372 : case 1:
3938 : case 2:
3939 : case 4:
3940 : case 5:
3941 : case 9:
3942 : case 10:
3943 : case 11:
3944 309372 : return 0.;
3945 317622 : case 3:
3946 : case 13:
3947 317622 : return 2.*(1. + zeta);
3948 317622 : case 6:
3949 : case 8:
3950 317622 : return 2.*(zeta - 1.);
3951 317622 : case 12:
3952 : case 14:
3953 317622 : return -2.*(1. + zeta);
3954 0 : default:
3955 0 : libmesh_error_msg("Invalid i = " << i);
3956 : }
3957 : }
3958 :
3959 : // d^2()/deta^2
3960 2382165 : case 2:
3961 : {
3962 2382165 : switch(i)
3963 : {
3964 317622 : case 0:
3965 : case 2:
3966 317622 : return 2.*(1. - zeta);
3967 397764 : 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 397764 : return 0.;
3977 317622 : case 3:
3978 : case 5:
3979 317622 : return 2.*(1. + zeta);
3980 158811 : case 8:
3981 158811 : return 4.*(zeta - 1.);
3982 158811 : case 14:
3983 158811 : return -4.*(1. + zeta);
3984 0 : default:
3985 0 : libmesh_error_msg("Invalid i = " << i);
3986 : }
3987 : }
3988 :
3989 : // d^2()/dxidzeta
3990 2382165 : case 3:
3991 : {
3992 2382165 : switch(i)
3993 : {
3994 158811 : case 0:
3995 158811 : return 1.5 - zeta - 2.*xi - 2.*eta;
3996 158811 : case 1:
3997 158811 : return 0.5 + zeta - 2.*xi;
3998 132588 : case 2:
3999 : case 5:
4000 : case 11:
4001 132588 : return 0.;
4002 158811 : case 3:
4003 158811 : return -1.5 - zeta + 2.*xi + 2.*eta;
4004 158811 : case 4:
4005 158811 : return -0.5 + zeta + 2.*xi;
4006 158811 : case 6:
4007 158811 : return 4.*xi + 2.*eta - 2.;
4008 158811 : case 7:
4009 158811 : return -2.*eta;
4010 158811 : case 8:
4011 158811 : return 2.*eta;
4012 158811 : case 9:
4013 158811 : return 2.*zeta;
4014 158811 : case 10:
4015 158811 : return -2.*zeta;
4016 158811 : case 12:
4017 158811 : return -4.*xi - 2.*eta + 2.;
4018 158811 : case 13:
4019 158811 : return 2.*eta;
4020 158811 : case 14:
4021 158811 : return -2.*eta;
4022 0 : default:
4023 0 : libmesh_error_msg("Invalid i = " << i);
4024 : }
4025 : }
4026 :
4027 : // d^2()/detadzeta
4028 2382165 : case 4:
4029 : {
4030 2382165 : switch(i)
4031 : {
4032 158811 : case 0:
4033 158811 : return 1.5 - zeta - 2.*xi - 2.*eta;
4034 132588 : case 1:
4035 : case 4:
4036 : case 10:
4037 132588 : return 0.;
4038 158811 : case 2:
4039 158811 : return .5 + zeta - 2.*eta;
4040 158811 : case 3:
4041 158811 : return -1.5 - zeta + 2.*xi + 2.*eta;
4042 158811 : case 5:
4043 158811 : return -.5 + zeta + 2.*eta;
4044 158811 : case 6:
4045 158811 : return 2.*xi;
4046 158811 : case 7:
4047 158811 : return -2.*xi;
4048 158811 : case 8:
4049 158811 : return 2.*xi + 4.*eta - 2.;
4050 158811 : case 9:
4051 158811 : return 2.*zeta;
4052 158811 : case 11:
4053 158811 : return -2.*zeta;
4054 158811 : case 12:
4055 158811 : return -2.*xi;
4056 158811 : case 13:
4057 158811 : return 2.*xi;
4058 158811 : case 14:
4059 158811 : 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 2382165 : case 5:
4067 : {
4068 2382165 : switch(i)
4069 : {
4070 317622 : case 0:
4071 : case 3:
4072 317622 : return 1. - xi - eta;
4073 88392 : case 1:
4074 : case 4:
4075 88392 : return xi;
4076 317622 : case 2:
4077 : case 5:
4078 317622 : return eta;
4079 530352 : case 6:
4080 : case 7:
4081 : case 8:
4082 : case 12:
4083 : case 13:
4084 : case 14:
4085 530352 : return 0.;
4086 158811 : case 9:
4087 158811 : return 2.*xi + 2.*eta - 2.;
4088 158811 : case 10:
4089 158811 : return -2.*xi;
4090 158811 : case 11:
4091 158811 : 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 20446344 : case PRISM6:
4106 413424 : libmesh_assert_msg(T == L2_LAGRANGE,
4107 : "High order on first order elements only supported for L2 families");
4108 : libmesh_fallthrough();
4109 6067872 : case PRISM18:
4110 : case PRISM20:
4111 : case PRISM21:
4112 : {
4113 7721568 : libmesh_assert_less (i, 18);
4114 :
4115 : // Compute prism shape functions as a tensor-product
4116 : // of a triangle and an edge
4117 :
4118 27754488 : Point p2d(p(0),p(1));
4119 27754488 : 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 27754488 : switch (j)
4126 : {
4127 : // d^2()/dxi^2
4128 4625748 : case 0:
4129 4625748 : return (FE<2,LAGRANGE>::shape_second_deriv(TRI6, SECOND, i1[i], 0, p2d)*
4130 5912676 : fe_lagrange_1D_quadratic_shape(i0[i], p1d));
4131 :
4132 : // d^2()/dxideta
4133 4625748 : case 1:
4134 4625748 : return (FE<2,LAGRANGE>::shape_second_deriv(TRI6, SECOND, i1[i], 1, p2d)*
4135 5912676 : fe_lagrange_1D_quadratic_shape(i0[i], p1d));
4136 :
4137 : // d^2()/deta^2
4138 4625748 : case 2:
4139 4625748 : return (FE<2,LAGRANGE>::shape_second_deriv(TRI6, SECOND, i1[i], 2, p2d)*
4140 5912676 : fe_lagrange_1D_quadratic_shape(i0[i], p1d));
4141 :
4142 : // d^2()/dxidzeta
4143 4625748 : case 3:
4144 4625748 : return (FE<2,LAGRANGE>::shape_deriv(TRI6, SECOND, i1[i], 0, p2d)*
4145 7964568 : fe_lagrange_1D_quadratic_shape_deriv(i0[i], 0, p1d));
4146 :
4147 : // d^2()/detadzeta
4148 4625748 : case 4:
4149 4625748 : return (FE<2,LAGRANGE>::shape_deriv(TRI6, SECOND, i1[i], 1, p2d)*
4150 7964568 : fe_lagrange_1D_quadratic_shape_deriv(i0[i], 0, p1d));
4151 :
4152 : // d^2()/dzeta^2
4153 4625748 : case 5:
4154 4625748 : return (FE<2,LAGRANGE>::shape(TRI6, SECOND, i1[i], p2d)*
4155 5738688 : 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 899640 : 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 332136 : case PYRAMID14:
4171 : case PYRAMID18:
4172 : {
4173 332136 : libmesh_assert_less (i, 14);
4174 :
4175 1231776 : const Real xi = p(0);
4176 1231776 : const Real eta = p(1);
4177 1231776 : const Real zeta = p(2);
4178 332136 : 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 1231776 : p1 = 0.5*(1. - eta - zeta), // back
4184 1231776 : p2 = 0.5*(1. + xi - zeta), // left
4185 1231776 : p3 = 0.5*(1. + eta - zeta), // front
4186 1231776 : p4 = 0.5*(1. - xi - zeta); // right
4187 :
4188 : // Denominators are perturbed by epsilon to avoid
4189 : // divide-by-zero issues.
4190 : Real
4191 1231776 : den = (-1. + zeta + eps),
4192 1231776 : den2 = den*den,
4193 1231776 : den3 = den2*den,
4194 1231776 : den4 = den2*den2;
4195 :
4196 : // These terms are used in several of the derivatives
4197 : Real
4198 1231776 : numer_mp = xi*eta - zeta + zeta*zeta,
4199 1231776 : numer_pm = xi*eta + zeta - zeta*zeta;
4200 :
4201 1231776 : switch (j)
4202 : {
4203 205296 : case 0: // d^2()/dxi^2
4204 : {
4205 205296 : switch(i)
4206 : {
4207 29328 : case 0:
4208 : case 1:
4209 29328 : return -p1*eta/den2;
4210 29328 : case 2:
4211 : case 3:
4212 29328 : return p3*eta/den2;
4213 19770 : case 4:
4214 : case 9:
4215 : case 10:
4216 : case 11:
4217 : case 12:
4218 19770 : return 0.;
4219 14664 : case 5:
4220 14664 : return 2.*p1*eta/den2;
4221 29328 : case 6:
4222 : case 8:
4223 29328 : return 4.*p1*p3/den2;
4224 14664 : case 7:
4225 14664 : return -2.*p3*eta/den2;
4226 14664 : case 13:
4227 14664 : return -8.*p1*p3/den2;
4228 0 : default:
4229 0 : libmesh_error_msg("Invalid i = " << i);
4230 : }
4231 : }
4232 :
4233 205296 : case 1: // d^2()/dxideta
4234 : {
4235 205296 : switch(i)
4236 : {
4237 14664 : case 0:
4238 14664 : return 0.25*numer_mp/den2
4239 14664 : - 0.5*p1*xi/den2
4240 14664 : - 0.5*p4*eta/den2
4241 14664 : + p4*p1/den2;
4242 :
4243 14664 : case 1:
4244 14664 : return 0.25*numer_pm/den2
4245 14664 : - 0.5*p1*xi/den2
4246 14664 : + 0.5*p2*eta/den2
4247 14664 : - p1*p2/den2;
4248 :
4249 14664 : case 2:
4250 14664 : return 0.25*numer_mp/den2
4251 14664 : + 0.5*p3*xi/den2
4252 14664 : + 0.5*p2*eta/den2
4253 14664 : + p2*p3/den2;
4254 :
4255 14664 : case 3:
4256 14664 : return 0.25*numer_pm/den2
4257 14664 : + 0.5*p3*xi/den2
4258 14664 : - 0.5*p4*eta/den2
4259 14664 : - p3*p4/den2;
4260 :
4261 3954 : case 4:
4262 3954 : return 0.;
4263 :
4264 14664 : case 5:
4265 14664 : return p4*eta/den2
4266 14664 : - 2.*p4*p1/den2
4267 14664 : - p2*eta/den2
4268 14664 : + 2.*p1*p2/den2;
4269 :
4270 14664 : case 6:
4271 14664 : return -p3*xi/den2
4272 14664 : + p1*xi/den2
4273 14664 : - 2.*p2*p3/den2
4274 14664 : + 2.*p1*p2/den2;
4275 :
4276 14664 : case 7:
4277 14664 : return p4*eta/den2
4278 14664 : + 2.*p3*p4/den2
4279 14664 : - p2*eta/den2
4280 14664 : - 2.*p2*p3/den2;
4281 :
4282 14664 : case 8:
4283 14664 : return -p3*xi/den2
4284 14664 : + p1*xi/den2
4285 14664 : - 2.*p4*p1/den2
4286 14664 : + 2.*p3*p4/den2;
4287 :
4288 29328 : case 9:
4289 : case 11:
4290 29328 : return -zeta/den;
4291 :
4292 29328 : case 10:
4293 : case 12:
4294 29328 : return zeta/den;
4295 :
4296 14664 : case 13:
4297 14664 : return 4.*p4*p1/den2
4298 14664 : - 4.*p3*p4/den2
4299 14664 : + 4.*p2*p3/den2
4300 14664 : - 4.*p1*p2/den2;
4301 :
4302 0 : default:
4303 0 : libmesh_error_msg("Invalid i = " << i);
4304 : }
4305 : }
4306 :
4307 :
4308 205296 : case 2: // d^2()/deta^2
4309 : {
4310 205296 : switch(i)
4311 : {
4312 29328 : case 0:
4313 : case 3:
4314 29328 : return -p4*xi/den2;
4315 29328 : case 1:
4316 : case 2:
4317 29328 : return p2*xi/den2;
4318 19770 : case 4:
4319 : case 9:
4320 : case 10:
4321 : case 11:
4322 : case 12:
4323 19770 : return 0.;
4324 29328 : case 5:
4325 : case 7:
4326 29328 : return 4.*p2*p4/den2;
4327 14664 : case 6:
4328 14664 : return -2.*p2*xi/den2;
4329 14664 : case 8:
4330 14664 : return 2.*p4*xi/den2;
4331 14664 : case 13:
4332 14664 : return -8.*p2*p4/den2;
4333 0 : default:
4334 0 : libmesh_error_msg("Invalid i = " << i);
4335 : }
4336 : }
4337 :
4338 :
4339 205296 : case 3: // d^2()/dxidzeta
4340 : {
4341 205296 : switch(i)
4342 : {
4343 14664 : case 0:
4344 14664 : return 0.25*numer_mp/den2
4345 14664 : - 0.5*p1*(2.*zeta - 1.)/den2
4346 14664 : + p1*numer_mp/den3
4347 14664 : - 0.5*p1*eta/den2
4348 14664 : - 0.5*p4*eta/den2
4349 14664 : - 2.*p4*p1*eta/den3;
4350 :
4351 14664 : case 1:
4352 14664 : return 0.25*numer_pm/den2
4353 14664 : - 0.5*p1*(1 - 2.*zeta)/den2
4354 14664 : + p1*numer_pm/den3
4355 14664 : + 0.5*p2*eta/den2
4356 14664 : + 0.5*p1*eta/den2
4357 14664 : + 2.*p1*p2*eta/den3;
4358 :
4359 14664 : case 2:
4360 14664 : return -0.25*numer_mp/den2
4361 14664 : + 0.5*p3*(2.*zeta - 1.)/den2
4362 14664 : - p3*numer_mp/den3
4363 14664 : - 0.5*p3*eta/den2
4364 14664 : - 0.5*p2*eta/den2
4365 14664 : - 2.*p2*p3*eta/den3;
4366 :
4367 14664 : case 3:
4368 14664 : return -0.25*numer_pm/den2
4369 14664 : + 0.5*p3*(1 - 2.*zeta)/den2
4370 14664 : - p3*numer_pm/den3
4371 14664 : + 0.5*p4*eta/den2
4372 14664 : + 0.5*p3*eta/den2
4373 14664 : + 2.*p3*p4*eta/den3;
4374 :
4375 3954 : case 4:
4376 3954 : return 0.;
4377 :
4378 14664 : case 5:
4379 14664 : return p4*eta/den2
4380 14664 : + 4.*p4*p1*eta/den3
4381 14664 : - p2*eta/den2
4382 14664 : - 4.*p1*p2*eta/den3;
4383 :
4384 14664 : case 6:
4385 14664 : return -p3*xi/den2
4386 14664 : - p1*xi/den2
4387 14664 : - 4.*p1*p3*xi/den3
4388 14664 : - 2.*p2*p3/den2
4389 14664 : - 2.*p1*p3/den2
4390 14664 : - 2.*p1*p2/den2
4391 14664 : - 8.*p1*p2*p3/den3;
4392 :
4393 14664 : case 7:
4394 14664 : return -p4*eta/den2
4395 14664 : - 4.*p3*p4*eta/den3
4396 14664 : + p2*eta/den2
4397 14664 : + 4.*p2*p3*eta/den3;
4398 :
4399 14664 : case 8:
4400 14664 : return -p3*xi/den2
4401 14664 : - p1*xi/den2
4402 14664 : - 4.*p1*p3*xi/den3
4403 14664 : + 2.*p4*p1/den2
4404 14664 : + 2.*p1*p3/den2
4405 14664 : + 2.*p3*p4/den2
4406 14664 : + 8.*p3*p4*p1/den3;
4407 :
4408 14664 : case 9:
4409 14664 : return -zeta/den
4410 14664 : + 2.*p1/den
4411 14664 : - 2.*p1*zeta/den2;
4412 :
4413 14664 : case 10:
4414 14664 : return zeta/den
4415 14664 : - 2.*p1/den
4416 14664 : + 2.*p1*zeta/den2;
4417 :
4418 14664 : case 11:
4419 14664 : return zeta/den
4420 14664 : - 2.*p3/den
4421 14664 : + 2.*p3*zeta/den2;
4422 :
4423 14664 : case 12:
4424 14664 : return -zeta/den
4425 14664 : + 2.*p3/den
4426 14664 : - 2.*p3*zeta/den2;
4427 :
4428 14664 : case 13:
4429 14664 : return -4.*p4*p1/den2
4430 14664 : - 4.*p3*p4/den2
4431 14664 : - 16.*p3*p4*p1/den3
4432 14664 : + 4.*p2*p3/den2
4433 14664 : + 4.*p1*p2/den2
4434 14664 : + 16.*p1*p2*p3/den3;
4435 :
4436 0 : default:
4437 0 : libmesh_error_msg("Invalid i = " << i);
4438 : }
4439 : }
4440 :
4441 205296 : case 4: // d^2()/detadzeta
4442 : {
4443 205296 : switch(i)
4444 : {
4445 14664 : case 0:
4446 14664 : return 0.25*numer_mp/den2
4447 14664 : - 0.5*p4*(2.*zeta - 1.)/den2
4448 14664 : + p4*numer_mp/den3
4449 14664 : - 0.5*p1*xi/den2
4450 14664 : - 0.5*p4*xi/den2
4451 14664 : - 2.*p4*p1*xi/den3;
4452 :
4453 14664 : case 1:
4454 14664 : return -0.25*numer_pm/den2
4455 14664 : + 0.5*p2*(1. - 2.*zeta)/den2
4456 14664 : - p2*numer_pm/den3
4457 14664 : + 0.5*p2*xi/den2
4458 14664 : + 0.5*p1*xi/den2
4459 14664 : + 2.*p1*p2*xi/den3;
4460 :
4461 14664 : case 2:
4462 14664 : return -0.25*numer_mp/den2
4463 14664 : + 0.5*p2*(2.*zeta - 1.)/den2
4464 14664 : - p2*numer_mp/den3
4465 14664 : - 0.5*p3*xi/den2
4466 14664 : - 0.5*p2*xi/den2
4467 14664 : - 2.*p2*p3*xi/den3;
4468 :
4469 14664 : case 3:
4470 14664 : return 0.25*numer_pm/den2
4471 14664 : - 0.5*p4*(1. - 2.*zeta)/den2
4472 14664 : + p4*numer_pm/den3
4473 14664 : + 0.5*p4*xi/den2
4474 14664 : + 0.5*p3*xi/den2
4475 14664 : + 2.*p3*p4*xi/den3;
4476 :
4477 3954 : case 4:
4478 3954 : return 0.;
4479 :
4480 14664 : case 5:
4481 14664 : return -p4*eta/den2
4482 14664 : - p2*eta/den2
4483 14664 : - 4.*p2*p4*eta/den3
4484 14664 : + 2.*p4*p1/den2
4485 14664 : + 2.*p2*p4/den2
4486 14664 : + 2.*p1*p2/den2
4487 14664 : + 8.*p2*p1*p4/den3;
4488 :
4489 14664 : case 6:
4490 14664 : return p3*xi/den2
4491 14664 : + 4.*p2*p3*xi/den3
4492 14664 : - p1*xi/den2
4493 14664 : - 4.*p1*p2*xi/den3;
4494 :
4495 14664 : case 7:
4496 14664 : return -p4*eta/den2
4497 14664 : - p2*eta/den2
4498 14664 : - 4.*p2*p4*eta/den3
4499 14664 : - 2.*p3*p4/den2
4500 14664 : - 2.*p2*p4/den2
4501 14664 : - 2.*p2*p3/den2
4502 14664 : - 8.*p2*p3*p4/den3;
4503 :
4504 14664 : case 8:
4505 14664 : return p1*xi/den2
4506 14664 : + 4.*p4*p1*xi/den3
4507 14664 : - p3*xi/den2
4508 14664 : - 4.*p3*p4*xi/den3;
4509 :
4510 14664 : case 9:
4511 14664 : return -zeta/den
4512 14664 : + 2.*p4/den
4513 14664 : - 2.*p4*zeta/den2;
4514 :
4515 14664 : case 10:
4516 14664 : return -zeta/den
4517 14664 : + 2.*p2/den
4518 14664 : - 2.*p2*zeta/den2;
4519 :
4520 14664 : case 11:
4521 14664 : return zeta/den
4522 14664 : - 2.*p2/den
4523 14664 : + 2.*p2*zeta/den2;
4524 :
4525 14664 : case 12:
4526 14664 : return zeta/den
4527 14664 : - 2.*p4/den
4528 14664 : + 2.*p4*zeta/den2;
4529 :
4530 14664 : case 13:
4531 14664 : return 4.*p3*p4/den2
4532 14664 : + 4.*p2*p3/den2
4533 14664 : + 16.*p2*p3*p4/den3
4534 14664 : - 4.*p4*p1/den2
4535 14664 : - 4.*p1*p2/den2
4536 14664 : - 16.*p2*p1*p4/den3;
4537 :
4538 0 : default:
4539 0 : libmesh_error_msg("Invalid i = " << i);
4540 : }
4541 : }
4542 :
4543 205296 : case 5: // d^2()/dzeta^2
4544 : {
4545 205296 : switch(i)
4546 : {
4547 14664 : case 0:
4548 14664 : return 0.5*numer_mp/den2
4549 14664 : - p1*(2.*zeta - 1.)/den2
4550 14664 : + 2.*p1*numer_mp/den3
4551 14664 : - p4*(2.*zeta - 1.)/den2
4552 14664 : + 2.*p4*numer_mp/den3
4553 14664 : + 2.*p4*p1/den2
4554 14664 : - 4.*p4*p1*(2.*zeta - 1.)/den3
4555 14664 : + 6.*p4*p1*numer_mp/den4;
4556 :
4557 14664 : case 1:
4558 14664 : return -0.5*numer_pm/den2
4559 14664 : + p2*(1 - 2.*zeta)/den2
4560 14664 : - 2.*p2*numer_pm/den3
4561 14664 : + p1*(1 - 2.*zeta)/den2
4562 14664 : - 2.*p1*numer_pm/den3
4563 14664 : + 2.*p1*p2/den2
4564 14664 : + 4.*p1*p2*(1 - 2.*zeta)/den3
4565 14664 : - 6.*p1*p2*numer_pm/den4;
4566 :
4567 14664 : case 2:
4568 14664 : return 0.5*numer_mp/den2
4569 14664 : - p3*(2.*zeta - 1.)/den2
4570 14664 : + 2.*p3*numer_mp/den3
4571 14664 : - p2*(2.*zeta - 1.)/den2
4572 14664 : + 2.*p2*numer_mp/den3
4573 14664 : + 2.*p2*p3/den2
4574 14664 : - 4.*p2*p3*(2.*zeta - 1.)/den3
4575 14664 : + 6.*p2*p3*numer_mp/den4;
4576 :
4577 14664 : case 3:
4578 14664 : return -0.5*numer_pm/den2
4579 14664 : + p4*(1 - 2.*zeta)/den2
4580 14664 : - 2.*p4*numer_pm/den3
4581 14664 : + p3*(1 - 2.*zeta)/den2
4582 14664 : - 2.*p3*numer_pm/den3
4583 14664 : + 2.*p3*p4/den2
4584 14664 : + 4.*p3*p4*(1 - 2.*zeta)/den3
4585 14664 : - 6.*p3*p4*numer_pm/den4;
4586 :
4587 3954 : case 4:
4588 3954 : return 4.;
4589 :
4590 14664 : case 5:
4591 14664 : return -2.*p1*eta/den2
4592 14664 : - 2.*p4*eta/den2
4593 14664 : - 8.*p4*p1*eta/den3
4594 14664 : - 2.*p2*eta/den2
4595 14664 : - 8.*p2*p4*eta/den3
4596 14664 : - 8.*p1*p2*eta/den3
4597 14664 : - 24.*p2*p1*p4*eta/den4;
4598 :
4599 14664 : case 6:
4600 14664 : return 2.*p3*xi/den2
4601 14664 : + 2.*p2*xi/den2
4602 14664 : + 8.*p2*p3*xi/den3
4603 14664 : + 2.*p1*xi/den2
4604 14664 : + 8.*p1*p3*xi/den3
4605 14664 : + 8.*p1*p2*xi/den3
4606 14664 : + 24.*p1*p2*p3*xi/den4;
4607 :
4608 14664 : case 7:
4609 14664 : return 2.*p4*eta/den2
4610 14664 : + 2.*p3*eta/den2
4611 14664 : + 8.*p3*p4*eta/den3
4612 14664 : + 2.*p2*eta/den2
4613 14664 : + 8.*p2*p4*eta/den3
4614 14664 : + 8.*p2*p3*eta/den3
4615 14664 : + 24.*p2*p3*p4*eta/den4;
4616 :
4617 14664 : case 8:
4618 14664 : return -2.*p1*xi/den2
4619 14664 : - 2.*p4*xi/den2
4620 14664 : - 8.*p4*p1*xi/den3
4621 14664 : - 2.*p3*xi/den2
4622 14664 : - 8.*p1*p3*xi/den3
4623 14664 : - 8.*p3*p4*xi/den3
4624 14664 : - 24.*p3*p4*p1*xi/den4;
4625 :
4626 14664 : case 9:
4627 14664 : return -2.*zeta/den
4628 14664 : + 4.*p4/den
4629 14664 : - 4.*p4*zeta/den2
4630 14664 : + 4.*p1/den
4631 14664 : - 4.*p1*zeta/den2
4632 14664 : + 8.*p4*p1/den2
4633 14664 : - 8.*p1*p4*zeta/den3;
4634 :
4635 14664 : case 10:
4636 14664 : return -2.*zeta/den
4637 14664 : + 4.*p1/den
4638 14664 : - 4.*p1*zeta/den2
4639 14664 : + 4.*p2/den
4640 14664 : - 4.*p2*zeta/den2
4641 14664 : + 8.*p1*p2/den2
4642 14664 : - 8.*p2*p1*zeta/den3;
4643 :
4644 14664 : case 11:
4645 14664 : return -2.*zeta/den
4646 14664 : + 4.*p2/den
4647 14664 : - 4.*p2*zeta/den2
4648 14664 : + 4.*p3/den
4649 14664 : - 4.*p3*zeta/den2
4650 14664 : + 8.*p2*p3/den2
4651 14664 : - 8.*p3*p2*zeta/den3;
4652 :
4653 14664 : case 12:
4654 14664 : return -2.*zeta/den
4655 14664 : + 4.*p3/den
4656 14664 : - 4.*p3*zeta/den2
4657 14664 : + 4.*p4/den
4658 14664 : - 4.*p4*zeta/den2
4659 14664 : + 8.*p3*p4/den2
4660 14664 : - 8.*p4*p3*zeta/den3;
4661 :
4662 14664 : case 13:
4663 14664 : return 8.*p3*p4/den2
4664 14664 : + 8.*p2*p4/den2
4665 14664 : + 8.*p2*p3/den2
4666 14664 : + 32.*p2*p3*p4/den3
4667 14664 : + 8.*p4*p1/den2
4668 14664 : + 8.*p1*p3/den2
4669 14664 : + 32.*p3*p4*p1/den3
4670 14664 : + 8.*p1*p2/den2
4671 14664 : + 32.*p2*p1*p4/den3
4672 14664 : + 32.*p1*p2*p3/den3
4673 14664 : + 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 762528 : case PYRAMID13:
4689 : {
4690 205608 : libmesh_assert_less (i, 13);
4691 :
4692 762528 : const Real xi = p(0);
4693 762528 : const Real eta = p(1);
4694 762528 : const Real zeta = p(2);
4695 205608 : const Real eps = 1.e-35;
4696 :
4697 : // Denominators are perturbed by epsilon to avoid
4698 : // divide-by-zero issues.
4699 : Real
4700 762528 : den = (-1. + zeta + eps),
4701 762528 : den2 = den*den,
4702 762528 : den3 = den2*den,
4703 762528 : xi2 = xi*xi,
4704 762528 : eta2 = eta*eta,
4705 762528 : zeta2 = zeta*zeta,
4706 762528 : zeta3 = zeta2*zeta;
4707 :
4708 762528 : switch (j)
4709 : {
4710 127088 : case 0: // d^2()/dxi^2
4711 : {
4712 127088 : switch(i)
4713 : {
4714 19552 : case 0:
4715 : case 1:
4716 19552 : return 0.5*(-1. + zeta + eta)/den;
4717 :
4718 19552 : case 2:
4719 : case 3:
4720 19552 : return 0.5*(-1. + zeta - eta)/den;
4721 :
4722 18452 : case 4:
4723 : case 6:
4724 : case 8:
4725 : case 9:
4726 : case 10:
4727 : case 11:
4728 : case 12:
4729 18452 : return 0.;
4730 :
4731 9776 : case 5:
4732 9776 : return (1. - eta - zeta)/den;
4733 :
4734 9776 : case 7:
4735 9776 : return (1. + eta - zeta)/den;
4736 :
4737 0 : default:
4738 0 : libmesh_error_msg("Invalid i = " << i);
4739 : }
4740 : }
4741 :
4742 127088 : case 1: // d^2()/dxideta
4743 : {
4744 127088 : switch(i)
4745 : {
4746 9776 : case 0:
4747 9776 : return 0.25*(-1. + 2.*zeta + 2.*xi + 2.*eta)/den;
4748 :
4749 9776 : case 1:
4750 9776 : return -0.25*(-1. + 2.*zeta - 2.*xi + 2.*eta)/den;
4751 :
4752 9776 : case 2:
4753 9776 : return -0.25*(1. - 2.*zeta + 2.*xi + 2.*eta)/den;
4754 :
4755 9776 : case 3:
4756 9776 : return 0.25*(1. - 2.*zeta - 2.*xi + 2.*eta)/den;
4757 :
4758 2636 : case 4:
4759 2636 : return 0.;
4760 :
4761 9776 : case 5:
4762 9776 : return -xi/den;
4763 :
4764 9776 : case 6:
4765 9776 : return eta/den;
4766 :
4767 9776 : case 7:
4768 9776 : return xi/den;
4769 :
4770 9776 : case 8:
4771 9776 : return -eta/den;
4772 :
4773 9776 : case 9:
4774 9776 : return -zeta/den;
4775 :
4776 9776 : case 10:
4777 9776 : return zeta/den;
4778 :
4779 9776 : case 11:
4780 9776 : return -zeta/den;
4781 :
4782 9776 : case 12:
4783 9776 : return zeta/den;
4784 :
4785 0 : default:
4786 0 : libmesh_error_msg("Invalid i = " << i);
4787 : }
4788 : }
4789 :
4790 :
4791 127088 : case 2: // d^2()/deta^2
4792 : {
4793 127088 : switch(i)
4794 : {
4795 19552 : case 0:
4796 : case 3:
4797 19552 : return 0.5*(-1. + zeta + xi)/den;
4798 :
4799 19552 : case 1:
4800 : case 2:
4801 19552 : return 0.5*(-1. + zeta - xi)/den;
4802 :
4803 18452 : case 4:
4804 : case 5:
4805 : case 7:
4806 : case 9:
4807 : case 10:
4808 : case 11:
4809 : case 12:
4810 18452 : return 0.;
4811 :
4812 9776 : case 6:
4813 9776 : return (1. + xi - zeta)/den;
4814 :
4815 9776 : case 8:
4816 9776 : return (1. - xi - zeta)/den;
4817 :
4818 0 : default:
4819 0 : libmesh_error_msg("Invalid i = " << i);
4820 : }
4821 : }
4822 :
4823 :
4824 127088 : case 3: // d^2()/dxidzeta
4825 : {
4826 127088 : switch(i)
4827 : {
4828 9776 : case 0:
4829 9776 : return -0.25*(-1. + 2.*zeta - zeta2 + eta + 2.*eta*xi + eta2)/den2;
4830 :
4831 9776 : case 1:
4832 9776 : return 0.25*(-1. + 2.*zeta - zeta2 + eta - 2.*eta*xi + eta2)/den2;
4833 :
4834 9776 : case 2:
4835 9776 : return 0.25*(-1. + 2.*zeta - zeta2 - eta + 2.*eta*xi + eta2)/den2;
4836 :
4837 9776 : case 3:
4838 9776 : return -0.25*(-1. + 2.*zeta - zeta2 - eta - 2.*eta*xi + eta2)/den2;
4839 :
4840 2636 : case 4:
4841 2636 : return 0.;
4842 :
4843 9776 : case 5:
4844 9776 : return eta*xi/den2;
4845 :
4846 9776 : case 6:
4847 9776 : return -0.5*(1. + zeta2 + eta2 - 2.*zeta)/den2;
4848 :
4849 9776 : case 7:
4850 9776 : return -eta*xi/den2;
4851 :
4852 9776 : case 8:
4853 9776 : return 0.5*(1. + zeta2 + eta2 - 2.*zeta)/den2;
4854 :
4855 9776 : case 9:
4856 9776 : return (-1. - zeta2 + eta + 2.*zeta)/den2;
4857 :
4858 9776 : case 10:
4859 9776 : return -(-1. - zeta2 + eta + 2.*zeta)/den2;
4860 :
4861 9776 : case 11:
4862 9776 : return (1. + zeta2 + eta - 2.*zeta)/den2;
4863 :
4864 9776 : case 12:
4865 9776 : return -(1. + zeta2 + eta - 2.*zeta)/den2;
4866 :
4867 0 : default:
4868 0 : libmesh_error_msg("Invalid i = " << i);
4869 : }
4870 : }
4871 :
4872 127088 : case 4: // d^2()/detadzeta
4873 : {
4874 127088 : switch(i)
4875 : {
4876 9776 : case 0:
4877 9776 : return -0.25*(-1. + 2.*zeta - zeta2 + xi + 2.*eta*xi + xi2)/den2;
4878 :
4879 9776 : case 1:
4880 9776 : return 0.25*(1. - 2.*zeta + zeta2 + xi + 2.*eta*xi - xi2)/den2;
4881 :
4882 9776 : case 2:
4883 9776 : return 0.25*(-1. + 2.*zeta - zeta2 - xi + 2.*eta*xi + xi2)/den2;
4884 :
4885 9776 : case 3:
4886 9776 : return -0.25*(1. - 2.*zeta + zeta2 - xi + 2.*eta*xi - xi2)/den2;
4887 :
4888 2636 : case 4:
4889 2636 : return 0.;
4890 :
4891 9776 : case 5:
4892 9776 : return 0.5*(1. + xi2 + zeta2 - 2.*zeta)/den2;
4893 :
4894 9776 : case 6:
4895 9776 : return -eta*xi/den2;
4896 :
4897 9776 : case 7:
4898 9776 : return -0.5*(1. + xi2 + zeta2 - 2.*zeta)/den2;
4899 :
4900 9776 : case 8:
4901 9776 : return eta*xi/den2;
4902 :
4903 9776 : case 9:
4904 9776 : return (-1. - zeta2 + xi + 2.*zeta)/den2;
4905 :
4906 9776 : case 10:
4907 9776 : return -(1. + zeta2 + xi - 2.*zeta)/den2;
4908 :
4909 9776 : case 11:
4910 9776 : return (1. + zeta2 + xi - 2.*zeta)/den2;
4911 :
4912 9776 : case 12:
4913 9776 : return -(-1. - zeta2 + xi + 2.*zeta)/den2;
4914 :
4915 0 : default:
4916 0 : libmesh_error_msg("Invalid i = " << i);
4917 : }
4918 : }
4919 :
4920 127088 : case 5: // d^2()/dzeta^2
4921 : {
4922 127088 : switch(i)
4923 : {
4924 9776 : case 0:
4925 9776 : return 0.5*(xi + eta + 1.)*eta*xi/den3;
4926 :
4927 9776 : case 1:
4928 9776 : return -0.5*(eta - xi + 1.)*eta*xi/den3;
4929 :
4930 9776 : case 2:
4931 9776 : return -0.5*(xi + eta - 1.)*eta*xi/den3;
4932 :
4933 9776 : case 3:
4934 9776 : return 0.5*(eta - xi - 1.)*eta*xi/den3;
4935 :
4936 2636 : case 4:
4937 2636 : return 4.;
4938 :
4939 9776 : case 5:
4940 9776 : return -(1. - 3.*zeta + 3.*zeta2 - zeta3 + eta*xi2)/den3;
4941 :
4942 9776 : case 6:
4943 9776 : return (-1. + 3.*zeta - 3.*zeta2 + zeta3 + eta2*xi)/den3;
4944 :
4945 9776 : case 7:
4946 9776 : return (-1. + 3.*zeta - 3.*zeta2 + zeta3 + eta*xi2)/den3;
4947 :
4948 9776 : case 8:
4949 9776 : return -(1. - 3.*zeta + 3.*zeta2 - zeta3 + eta2*xi)/den3;
4950 :
4951 9776 : case 9:
4952 9776 : return -2.*(-1. + 3.*zeta - 3.*zeta2 + zeta3 + eta*xi)/den3;
4953 :
4954 9776 : case 10:
4955 9776 : return 2.*(1. - 3.*zeta + 3.*zeta2 - zeta3 + eta*xi)/den3;
4956 :
4957 9776 : case 11:
4958 9776 : return -2.*(-1. + 3.*zeta - 3.*zeta2 + zeta3 + eta*xi)/den3;
4959 :
4960 9776 : case 12:
4961 9776 : 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 334509678 : case THIRD:
4979 : {
4980 334509678 : switch (type)
4981 : {
4982 : // quadratic Lagrange shape functions with a cubic bubble
4983 233848524 : case TET14:
4984 : {
4985 56529648 : 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 233848524 : const unsigned int my_j = independent_var_indices[j][0];
5033 233848524 : const unsigned int my_k = independent_var_indices[j][1];
5034 :
5035 56529648 : Real returnval = 0;
5036 233848524 : if (i<4)
5037 66813864 : returnval = 4.*dzetadxi[i][my_j]*dzetadxi[i][my_k];
5038 :
5039 167034660 : else if (i<10)
5040 : {
5041 100220796 : const unsigned short int my_m = zeta_indices[i][0];
5042 100220796 : const unsigned short int my_n = zeta_indices[i][1];
5043 :
5044 100220796 : returnval =
5045 148674780 : 4.*(dzetadxi[my_n][my_j]*dzetadxi[my_m][my_k] +
5046 100220796 : dzetadxi[my_m][my_j]*dzetadxi[my_n][my_k] );
5047 : }
5048 :
5049 233848524 : const Real zeta1 = p(0);
5050 233848524 : const Real zeta2 = p(1);
5051 233848524 : const Real zeta3 = p(2);
5052 233848524 : 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 56529648 : switch (j)
5058 : {
5059 : // d^2()/dxi^2
5060 38974754 : case 0:
5061 : {
5062 38974754 : d2bubble012 = -2.*zeta2;
5063 38974754 : d2bubble013 = -2.*zeta3;
5064 9421608 : d2bubble023 = 0.;
5065 9421608 : d2bubble123 = 0.;
5066 38974754 : break;
5067 : }
5068 :
5069 : // d^2()/dxideta
5070 38974754 : case 1:
5071 : {
5072 38974754 : d2bubble012 = (zeta0-zeta1)-zeta2;
5073 38974754 : d2bubble013 = -zeta3;
5074 9421608 : d2bubble123 = zeta3;
5075 9421608 : d2bubble023 = -zeta3;
5076 38974754 : break;
5077 : }
5078 :
5079 : // d^2()/deta^2
5080 38974754 : case 2:
5081 : {
5082 38974754 : d2bubble012 = -2.*zeta1;
5083 9421608 : d2bubble013 = 0.;
5084 9421608 : d2bubble123 = 0.;
5085 38974754 : d2bubble023 = -2.*zeta3;
5086 38974754 : break;
5087 : }
5088 :
5089 : // d^2()/dxi dzeta
5090 38974754 : case 3:
5091 : {
5092 38974754 : d2bubble012 = -zeta2;
5093 38974754 : d2bubble013 = (zeta0-zeta3)-zeta1;
5094 9421608 : d2bubble123 = zeta2;
5095 9421608 : d2bubble023 = -zeta2;
5096 38974754 : break;
5097 : }
5098 :
5099 : // d^2()/deta dzeta
5100 38974754 : case 4:
5101 : {
5102 38974754 : d2bubble012 = -zeta1;
5103 9421608 : d2bubble013 = -zeta1;
5104 9421608 : d2bubble123 = zeta1;
5105 38974754 : d2bubble023 = (zeta0-zeta3)-zeta2;
5106 38974754 : break;
5107 : }
5108 :
5109 : // d^2()/dzeta^2
5110 38974754 : case 5:
5111 : {
5112 9421608 : d2bubble012 = 0.;
5113 38974754 : d2bubble013 = -2.*zeta1;
5114 9421608 : d2bubble123 = 0.;
5115 38974754 : d2bubble023 = -2.*zeta2;
5116 38974754 : break;
5117 : }
5118 :
5119 0 : default:
5120 0 : libmesh_error_msg("Invalid j = " << j);
5121 : }
5122 :
5123 233848524 : switch (i)
5124 : {
5125 16703466 : case 0:
5126 16703466 : return returnval + 3.*(d2bubble012+d2bubble013+d2bubble023);
5127 :
5128 16703466 : case 1:
5129 16703466 : return returnval + 3.*(d2bubble012+d2bubble013+d2bubble123);
5130 :
5131 16703466 : case 2:
5132 16703466 : return returnval + 3.*(d2bubble012+d2bubble023+d2bubble123);
5133 :
5134 16703466 : case 3:
5135 16703466 : return returnval + 3.*(d2bubble013+d2bubble023+d2bubble123);
5136 :
5137 16703466 : case 4:
5138 16703466 : return returnval - 12.*(d2bubble012+d2bubble013);
5139 :
5140 16703466 : case 5:
5141 16703466 : return returnval - 12.*(d2bubble012+d2bubble123);
5142 :
5143 16703466 : case 6:
5144 16703466 : return returnval - 12.*(d2bubble012+d2bubble023);
5145 :
5146 16703466 : case 7:
5147 16703466 : return returnval - 12.*(d2bubble013+d2bubble023);
5148 :
5149 16703466 : case 8:
5150 16703466 : return returnval - 12.*(d2bubble013+d2bubble123);
5151 :
5152 16703466 : case 9:
5153 16703466 : return returnval - 12.*(d2bubble023+d2bubble123);
5154 :
5155 16703466 : case 10:
5156 16703466 : return 27.*d2bubble012;
5157 :
5158 16703466 : case 11:
5159 16703466 : return 27.*d2bubble013;
5160 :
5161 16703466 : case 12:
5162 16703466 : return 27.*d2bubble123;
5163 :
5164 16703466 : case 13:
5165 16703466 : return 27.*d2bubble023;
5166 :
5167 0 : default:
5168 0 : libmesh_error_msg("Invalid i = " << i);
5169 : }
5170 : }
5171 :
5172 30320760 : case PRISM20:
5173 : {
5174 8164560 : libmesh_assert_less (i, 20);
5175 :
5176 : // Compute prism shape functions as a tensor-product
5177 : // of a triangle and an edge
5178 :
5179 30320760 : Point p2d(p(0),p(1));
5180 30320760 : 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 30320760 : const Real mainval = FE<3,LAGRANGE>::shape_second_deriv(PRISM21, THIRD, i, j, p);
5186 :
5187 30320760 : if (i0[i] != 2)
5188 5715192 : return mainval;
5189 :
5190 2449368 : Real bubbleval = 0;
5191 :
5192 2449368 : switch (j)
5193 : {
5194 : // d^2()/dxi^2
5195 1516038 : case 0:
5196 1516038 : bubbleval =
5197 1516038 : FE<2,LAGRANGE>::shape_second_deriv(TRI7, THIRD, 6, 0, p2d)*
5198 408228 : fe_lagrange_1D_quadratic_shape(2, p1d);
5199 1516038 : break;
5200 :
5201 : // d^2()/dxideta
5202 1516038 : case 1:
5203 1516038 : bubbleval =
5204 1516038 : FE<2,LAGRANGE>::shape_second_deriv(TRI7, THIRD, 6, 1, p2d)*
5205 408228 : fe_lagrange_1D_quadratic_shape(2, p1d);
5206 1516038 : break;
5207 :
5208 : // d^2()/deta^2
5209 1516038 : case 2:
5210 1516038 : bubbleval =
5211 1516038 : FE<2,LAGRANGE>::shape_second_deriv(TRI7, THIRD, 6, 2, p2d)*
5212 408228 : fe_lagrange_1D_quadratic_shape(2, p1d);
5213 1516038 : break;
5214 :
5215 : // d^2()/dxidzeta
5216 1516038 : case 3:
5217 1516038 : bubbleval =
5218 1516038 : FE<2,LAGRANGE>::shape_deriv(TRI7, THIRD, 6, 0, p2d)*
5219 408228 : fe_lagrange_1D_quadratic_shape_deriv(2, 0, p1d);
5220 1516038 : break;
5221 :
5222 : // d^2()/detadzeta
5223 1516038 : case 4:
5224 1516038 : bubbleval =
5225 1516038 : FE<2,LAGRANGE>::shape_deriv(TRI7, THIRD, 6, 1, p2d)*
5226 408228 : fe_lagrange_1D_quadratic_shape_deriv(2, 0, p1d);
5227 1516038 : break;
5228 :
5229 : // d^2()/dzeta^2
5230 1516038 : case 5:
5231 1516038 : bubbleval =
5232 1516038 : FE<2,LAGRANGE>::shape(TRI7, THIRD, 6, p2d)*
5233 408228 : fe_lagrange_1D_quadratic_shape_second_deriv(2, 0, p1d);
5234 1516038 : break;
5235 :
5236 0 : default:
5237 0 : libmesh_error_msg("Invalid shape function derivative j = " << j);
5238 : }
5239 :
5240 9096228 : if (i < 12) // vertices
5241 4548114 : return mainval - bubbleval / 9;
5242 :
5243 4548114 : return mainval + bubbleval * (Real(4) / 9);
5244 : }
5245 :
5246 68588850 : case PRISM21:
5247 : {
5248 18534360 : libmesh_assert_less (i, 21);
5249 :
5250 : // Compute prism shape functions as a tensor-product
5251 : // of a triangle and an edge
5252 :
5253 68588850 : Point p2d(p(0),p(1));
5254 68588850 : 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 68588850 : switch (j)
5261 : {
5262 : // d^2()/dxi^2
5263 11431475 : case 0:
5264 11431475 : return (FE<2,LAGRANGE>::shape_second_deriv(TRI7, THIRD, i1[i], 0, p2d)*
5265 14520535 : fe_lagrange_1D_quadratic_shape(i0[i], p1d));
5266 :
5267 : // d^2()/dxideta
5268 11431475 : case 1:
5269 11431475 : return (FE<2,LAGRANGE>::shape_second_deriv(TRI7, THIRD, i1[i], 1, p2d)*
5270 14520535 : fe_lagrange_1D_quadratic_shape(i0[i], p1d));
5271 :
5272 : // d^2()/deta^2
5273 11431475 : case 2:
5274 11431475 : return (FE<2,LAGRANGE>::shape_second_deriv(TRI7, THIRD, i1[i], 2, p2d)*
5275 14520535 : fe_lagrange_1D_quadratic_shape(i0[i], p1d));
5276 :
5277 : // d^2()/dxidzeta
5278 11431475 : case 3:
5279 11431475 : return (FE<2,LAGRANGE>::shape_deriv(TRI7, THIRD, i1[i], 0, p2d)*
5280 19773890 : fe_lagrange_1D_quadratic_shape_deriv(i0[i], 0, p1d));
5281 :
5282 : // d^2()/detadzeta
5283 11431475 : case 4:
5284 11431475 : return (FE<2,LAGRANGE>::shape_deriv(TRI7, THIRD, i1[i], 1, p2d)*
5285 19773890 : fe_lagrange_1D_quadratic_shape_deriv(i0[i], 0, p1d));
5286 :
5287 : // d^2()/dzeta^2
5288 11431475 : case 5:
5289 11431475 : return (FE<2,LAGRANGE>::shape(TRI7, THIRD, i1[i], p2d)*
5290 14089190 : 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 1751544 : case PYRAMID18:
5298 : {
5299 1751544 : return fe_fdm_second_deriv(type, order, elem, i, j, p,
5300 1751544 : 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
|
