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