Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2026 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 :
19 : // Local includes
20 : #include "libmesh/fe.h"
21 : #include "libmesh/elem.h"
22 : #include "libmesh/fe_lagrange_shape_1D.h"
23 : #include "libmesh/enum_to_string.h"
24 : #include "libmesh/face_c0polygon.h"
25 :
26 : // Anonymous namespace for functions shared by LAGRANGE and
27 : // L2_LAGRANGE implementations. Implementations appear at the bottom
28 : // of this file.
29 : namespace
30 : {
31 : using namespace libMesh;
32 :
33 : template <FEFamily T>
34 : Real fe_lagrange_2D_shape(const ElemType,
35 : const Elem * elem,
36 : const Order order,
37 : const unsigned int i,
38 : const Point & p);
39 :
40 : template <FEFamily T>
41 : Real fe_lagrange_2D_shape_deriv(const ElemType type,
42 : const Elem * elem,
43 : const Order order,
44 : const unsigned int i,
45 : const unsigned int j,
46 : const Point & p);
47 :
48 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
49 :
50 : template <FEFamily T>
51 : Real fe_lagrange_2D_shape_second_deriv(const ElemType type,
52 : const Elem * elem,
53 : const Order order,
54 : const unsigned int i,
55 : const unsigned int j,
56 : const Point & p);
57 :
58 : #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
59 :
60 : } // anonymous namespace
61 :
62 :
63 :
64 : namespace libMesh
65 : {
66 :
67 :
68 242076497 : LIBMESH_DEFAULT_VECTORIZED_FE(2,LAGRANGE)
69 1333199 : LIBMESH_DEFAULT_VECTORIZED_FE(2,L2_LAGRANGE)
70 :
71 :
72 : template <>
73 249476816 : Real FE<2,LAGRANGE>::shape(const ElemType type,
74 : const Order order,
75 : const unsigned int i,
76 : const Point & p)
77 : {
78 249476816 : return fe_lagrange_2D_shape<LAGRANGE>(type, nullptr, order, i, p);
79 : }
80 :
81 :
82 :
83 : template <>
84 0 : Real FE<2,L2_LAGRANGE>::shape(const ElemType type,
85 : const Order order,
86 : const unsigned int i,
87 : const Point & p)
88 : {
89 0 : return fe_lagrange_2D_shape<L2_LAGRANGE>(type, nullptr, order, i, p);
90 : }
91 :
92 :
93 : template <>
94 540199804 : Real FE<2,LAGRANGE>::shape(const Elem * elem,
95 : const Order order,
96 : const unsigned int i,
97 : const Point & p,
98 : const bool add_p_level)
99 : {
100 167105469 : libmesh_assert(elem);
101 :
102 : // call the orientation-independent shape functions
103 853104758 : return fe_lagrange_2D_shape<LAGRANGE>(elem->type(), elem, order + add_p_level*elem->p_level(), i, p);
104 : }
105 :
106 :
107 :
108 : template <>
109 2363318 : Real FE<2,L2_LAGRANGE>::shape(const Elem * elem,
110 : const Order order,
111 : const unsigned int i,
112 : const Point & p,
113 : const bool add_p_level)
114 : {
115 606904 : libmesh_assert(elem);
116 :
117 : // call the orientation-independent shape functions
118 3577126 : return fe_lagrange_2D_shape<L2_LAGRANGE>(elem->type(), elem, order + add_p_level*elem->p_level(), i, p);
119 : }
120 :
121 :
122 : template <>
123 441910137 : Real FE<2,LAGRANGE>::shape(const FEType fet,
124 : const Elem * elem,
125 : const unsigned int i,
126 : const Point & p,
127 : const bool add_p_level)
128 : {
129 138754920 : libmesh_assert(elem);
130 671558515 : return fe_lagrange_2D_shape<LAGRANGE>(elem->type(), elem, fet.order + add_p_level*elem->p_level(), i, p);
131 : }
132 :
133 :
134 :
135 : template <>
136 0 : Real FE<2,L2_LAGRANGE>::shape(const FEType fet,
137 : const Elem * elem,
138 : const unsigned int i,
139 : const Point & p,
140 : const bool add_p_level)
141 : {
142 0 : libmesh_assert(elem);
143 0 : return fe_lagrange_2D_shape<L2_LAGRANGE>(elem->type(), elem, fet.order + add_p_level*elem->p_level(), i, p);
144 : }
145 :
146 :
147 : template <>
148 262325710 : Real FE<2,LAGRANGE>::shape_deriv(const ElemType type,
149 : const Order order,
150 : const unsigned int i,
151 : const unsigned int j,
152 : const Point & p)
153 : {
154 262325710 : return fe_lagrange_2D_shape_deriv<LAGRANGE>(type, nullptr, order, i, j, p);
155 : }
156 :
157 :
158 :
159 : template <>
160 0 : Real FE<2,L2_LAGRANGE>::shape_deriv(const ElemType type,
161 : const Order order,
162 : const unsigned int i,
163 : const unsigned int j,
164 : const Point & p)
165 : {
166 0 : return fe_lagrange_2D_shape_deriv<L2_LAGRANGE>(type, nullptr, order, i, j, p);
167 : }
168 :
169 :
170 :
171 : template <>
172 228495712 : Real FE<2,LAGRANGE>::shape_deriv(const Elem * elem,
173 : const Order order,
174 : const unsigned int i,
175 : const unsigned int j,
176 : const Point & p,
177 : const bool add_p_level)
178 : {
179 64125526 : libmesh_assert(elem);
180 :
181 : // call the orientation-independent shape functions
182 356722700 : return fe_lagrange_2D_shape_deriv<LAGRANGE>(elem->type(), elem, order + add_p_level*elem->p_level(), i, j, p);
183 : }
184 :
185 :
186 :
187 : template <>
188 543656 : Real FE<2,L2_LAGRANGE>::shape_deriv(const Elem * elem,
189 : const Order order,
190 : const unsigned int i,
191 : const unsigned int j,
192 : const Point & p,
193 : const bool add_p_level)
194 : {
195 140040 : libmesh_assert(elem);
196 :
197 : // call the orientation-independent shape functions
198 823736 : return fe_lagrange_2D_shape_deriv<L2_LAGRANGE>(elem->type(), elem, order + add_p_level*elem->p_level(), i, j, p);
199 : }
200 :
201 :
202 : template <>
203 4578152226 : Real FE<2,LAGRANGE>::shape_deriv(const FEType fet,
204 : const Elem * elem,
205 : const unsigned int i,
206 : const unsigned int j,
207 : const Point & p,
208 : const bool add_p_level)
209 : {
210 1170050916 : libmesh_assert(elem);
211 6913111582 : return fe_lagrange_2D_shape_deriv<LAGRANGE>(elem->type(), elem, fet.order + add_p_level*elem->p_level(), i, j, p);
212 : }
213 :
214 :
215 :
216 : template <>
217 0 : Real FE<2,L2_LAGRANGE>::shape_deriv(const FEType fet,
218 : const Elem * elem,
219 : const unsigned int i,
220 : const unsigned int j,
221 : const Point & p,
222 : const bool add_p_level)
223 : {
224 0 : libmesh_assert(elem);
225 0 : return fe_lagrange_2D_shape_deriv<L2_LAGRANGE>(elem->type(), elem, fet.order + add_p_level*elem->p_level(), i, j, p);
226 : }
227 :
228 :
229 :
230 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
231 :
232 : template <>
233 48406047 : Real FE<2,LAGRANGE>::shape_second_deriv(const ElemType type,
234 : const Order order,
235 : const unsigned int i,
236 : const unsigned int j,
237 : const Point & p)
238 : {
239 48406047 : return fe_lagrange_2D_shape_second_deriv<LAGRANGE>(type, nullptr, order, i, j, p);
240 : }
241 :
242 :
243 :
244 : template <>
245 0 : Real FE<2,L2_LAGRANGE>::shape_second_deriv(const ElemType type,
246 : const Order order,
247 : const unsigned int i,
248 : const unsigned int j,
249 : const Point & p)
250 : {
251 0 : return fe_lagrange_2D_shape_second_deriv<L2_LAGRANGE>(type, nullptr, order, i, j, p);
252 : }
253 :
254 :
255 :
256 : template <>
257 1066743 : Real FE<2,LAGRANGE>::shape_second_deriv(const Elem * elem,
258 : const Order order,
259 : const unsigned int i,
260 : const unsigned int j,
261 : const Point & p,
262 : const bool add_p_level)
263 : {
264 336408 : libmesh_assert(elem);
265 :
266 : // call the orientation-independent shape functions
267 1739559 : return fe_lagrange_2D_shape_second_deriv<LAGRANGE>(elem->type(), elem, order + add_p_level*elem->p_level(), i, j, p);
268 : }
269 :
270 :
271 :
272 : template <>
273 292044 : Real FE<2,L2_LAGRANGE>::shape_second_deriv(const Elem * elem,
274 : const Order order,
275 : const unsigned int i,
276 : const unsigned int j,
277 : const Point & p,
278 : const bool add_p_level)
279 : {
280 78876 : libmesh_assert(elem);
281 :
282 : // call the orientation-independent shape functions
283 449796 : return fe_lagrange_2D_shape_second_deriv<L2_LAGRANGE>(elem->type(), elem, order + add_p_level*elem->p_level(), i, j, p);
284 : }
285 :
286 :
287 : template <>
288 47606409 : Real FE<2,LAGRANGE>::shape_second_deriv(const FEType fet,
289 : const Elem * elem,
290 : const unsigned int i,
291 : const unsigned int j,
292 : const Point & p,
293 : const bool add_p_level)
294 : {
295 12044136 : libmesh_assert(elem);
296 71694681 : return fe_lagrange_2D_shape_second_deriv<LAGRANGE>(elem->type(), elem, fet.order + add_p_level*elem->p_level(), i, j, p);
297 : }
298 :
299 :
300 :
301 : template <>
302 0 : Real FE<2,L2_LAGRANGE>::shape_second_deriv(const FEType fet,
303 : const Elem * elem,
304 : const unsigned int i,
305 : const unsigned int j,
306 : const Point & p,
307 : const bool add_p_level)
308 : {
309 0 : libmesh_assert(elem);
310 0 : return fe_lagrange_2D_shape_second_deriv<L2_LAGRANGE>(elem->type(), elem, fet.order + add_p_level*elem->p_level(), i, j, p);
311 : }
312 :
313 : #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
314 :
315 : } // namespace libMesh
316 :
317 :
318 :
319 :
320 : // Anonymous namespace function definitions
321 : namespace
322 : {
323 : using namespace libMesh;
324 :
325 : template <FEFamily T>
326 1233950075 : Real fe_lagrange_2D_shape(const ElemType type,
327 : const Elem * elem,
328 : const Order order,
329 : const unsigned int i,
330 : const Point & p)
331 : {
332 : #if LIBMESH_DIM > 1
333 :
334 1233950075 : switch (order)
335 : {
336 : // linear Lagrange shape functions
337 371756947 : case FIRST:
338 : {
339 258489506 : switch (type)
340 : {
341 316277660 : case QUAD4:
342 : case QUADSHELL4:
343 : case QUAD8:
344 : case QUADSHELL8:
345 : case QUAD9:
346 : case QUADSHELL9:
347 : {
348 : // Compute quad shape functions as a tensor-product
349 316277660 : const Real xi = p(0);
350 316277660 : const Real eta = p(1);
351 :
352 102979408 : libmesh_assert_less (i, 4);
353 :
354 : // 0 1 2 3
355 : static const unsigned int i0[] = {0, 1, 1, 0};
356 : static const unsigned int i1[] = {0, 0, 1, 1};
357 :
358 316277660 : return (fe_lagrange_1D_linear_shape(i0[i], xi)*
359 529575912 : fe_lagrange_1D_linear_shape(i1[i], eta));
360 : }
361 :
362 55238487 : case TRI3:
363 : case TRISHELL3:
364 : case TRI6:
365 : case TRI7:
366 : {
367 55238487 : const Real zeta1 = p(0);
368 55238487 : const Real zeta2 = p(1);
369 55238487 : const Real zeta0 = 1. - zeta1 - zeta2;
370 :
371 14917638 : libmesh_assert_less (i, 3);
372 :
373 55238487 : switch(i)
374 : {
375 4972546 : case 0:
376 4972546 : return zeta0;
377 :
378 18412829 : case 1:
379 18412829 : return zeta1;
380 :
381 18412829 : case 2:
382 18412829 : return zeta2;
383 :
384 0 : default:
385 0 : libmesh_error_msg("Invalid shape function index i = " << i);
386 : }
387 : }
388 :
389 240800 : case C0POLYGON:
390 : {
391 : // C0Polygon requires using newer FE APIs
392 240800 : if (!elem)
393 0 : libmesh_error_msg("Code (see stack trace) used an outdated FE function overload.\n"
394 : "Shape functions on a C0Polygon are not defined by its ElemType alone.");
395 :
396 67025 : libmesh_assert(elem->type() == C0POLYGON);
397 :
398 67025 : const C0Polygon & poly = *cast_ptr<const C0Polygon *>(elem);
399 :
400 : // We can't use a small tolerance here, because in
401 : // inverse_map() Newton might hand us intermediate
402 : // iterates outside the polygon.
403 240800 : const auto [s, a, b] = poly.subtriangle_coordinates(p, 100);
404 240800 : if (s == invalid_uint)
405 0 : return 0;
406 67025 : libmesh_assert_less(s, poly.n_subtriangles());
407 :
408 240800 : const auto subtri = poly.subtriangle(s);
409 :
410 : // Avoid signed/unsigned comparison warnings
411 240800 : const int nodei = i;
412 240800 : if (nodei == subtri[0])
413 48160 : return 1-a-b;
414 192640 : if (nodei == subtri[1])
415 48160 : return a;
416 144480 : if (nodei == subtri[2])
417 48160 : return b;
418 :
419 : // Basis function i is not supported on p's subtriangle
420 26810 : return 0;
421 : }
422 :
423 0 : default:
424 0 : libmesh_error_msg("ERROR: Unsupported 2D element type: " << Utility::enum_to_string(type));
425 : }
426 : }
427 :
428 :
429 : // quadratic Lagrange shape functions
430 554683758 : case SECOND:
431 : {
432 407900867 : switch (type)
433 : {
434 9106464 : case QUAD8:
435 : case QUADSHELL8:
436 : {
437 9106464 : const Real xi = p(0);
438 9106464 : const Real eta = p(1);
439 :
440 2632048 : libmesh_assert_less (i, 8);
441 :
442 9106464 : switch (i)
443 : {
444 1138308 : case 0:
445 1138308 : return .25*(1. - xi)*(1. - eta)*(-1. - xi - eta);
446 :
447 1138308 : case 1:
448 1138308 : return .25*(1. + xi)*(1. - eta)*(-1. + xi - eta);
449 :
450 1138308 : case 2:
451 1138308 : return .25*(1. + xi)*(1. + eta)*(-1. + xi + eta);
452 :
453 1138308 : case 3:
454 1138308 : return .25*(1. - xi)*(1. + eta)*(-1. - xi + eta);
455 :
456 1138308 : case 4:
457 1138308 : return .5*(1. - xi*xi)*(1. - eta);
458 :
459 1138308 : case 5:
460 1138308 : return .5*(1. + xi)*(1. - eta*eta);
461 :
462 1138308 : case 6:
463 1138308 : return .5*(1. - xi*xi)*(1. + eta);
464 :
465 1138308 : case 7:
466 1138308 : return .5*(1. - xi)*(1. - eta*eta);
467 :
468 0 : default:
469 0 : libmesh_error_msg("Invalid shape function index i = " << i);
470 : }
471 : }
472 :
473 214380405 : case QUAD4:
474 4743 : libmesh_assert_msg(T == L2_LAGRANGE,
475 : "High order on first order elements only supported for L2 families");
476 : libmesh_fallthrough();
477 101863098 : case QUAD9:
478 : case QUADSHELL9:
479 : {
480 :
481 : // Compute quad shape functions as a tensor-product
482 316248246 : const Real xi = p(0);
483 316248246 : const Real eta = p(1);
484 :
485 101872584 : libmesh_assert_less (i, 9);
486 :
487 : // 0 1 2 3 4 5 6 7 8
488 : static const unsigned int i0[] = {0, 1, 1, 0, 2, 1, 2, 0, 2};
489 : static const unsigned int i1[] = {0, 0, 1, 1, 0, 2, 1, 2, 2};
490 :
491 316248246 : return (fe_lagrange_1D_quadratic_shape(i0[i], xi)*
492 398438469 : fe_lagrange_1D_quadratic_shape(i1[i], eta));
493 : }
494 :
495 157795866 : case TRI3:
496 3582 : libmesh_assert_msg(T == L2_LAGRANGE,
497 : "High order on first order elements only supported for L2 families");
498 : libmesh_fallthrough();
499 71288418 : case TRI6:
500 : case TRI7:
501 : {
502 229329048 : const Real zeta1 = p(0);
503 229329048 : const Real zeta2 = p(1);
504 229329048 : const Real zeta0 = 1. - zeta1 - zeta2;
505 :
506 71536764 : libmesh_assert_less (i, 6);
507 :
508 229329048 : switch(i)
509 : {
510 38221508 : case 0:
511 38221508 : return 2.*zeta0*(zeta0-0.5);
512 :
513 38221508 : case 1:
514 38221508 : return 2.*zeta1*(zeta1-0.5);
515 :
516 38221508 : case 2:
517 38221508 : return 2.*zeta2*(zeta2-0.5);
518 :
519 38221508 : case 3:
520 38221508 : return 4.*zeta0*zeta1;
521 :
522 38221508 : case 4:
523 38221508 : return 4.*zeta1*zeta2;
524 :
525 38221508 : case 5:
526 38221508 : return 4.*zeta2*zeta0;
527 :
528 0 : default:
529 0 : libmesh_error_msg("Invalid shape function index i = " << i);
530 : }
531 : }
532 :
533 0 : default:
534 0 : libmesh_error_msg("ERROR: Unsupported 2D element type: " << Utility::enum_to_string(type));
535 : }
536 : }
537 :
538 :
539 :
540 : // "cubic" (one cubic bubble) Lagrange shape functions on TRI7
541 307509370 : case THIRD:
542 : {
543 307509370 : switch (type)
544 : {
545 307509370 : case TRI7:
546 : {
547 307509370 : const Real zeta1 = p(0);
548 307509370 : const Real zeta2 = p(1);
549 307509370 : const Real zeta0 = 1. - zeta1 - zeta2;
550 307509370 : const Real bubble_27th = zeta0*zeta1*zeta2;
551 :
552 84860769 : libmesh_assert_less (i, 7);
553 :
554 307509370 : switch(i)
555 : {
556 41899595 : case 0:
557 41899595 : return 2.*zeta0*(zeta0-0.5) + 3.*bubble_27th;
558 :
559 41899595 : case 1:
560 41899595 : return 2.*zeta1*(zeta1-0.5) + 3.*bubble_27th;
561 :
562 41899595 : case 2:
563 41899595 : return 2.*zeta2*(zeta2-0.5) + 3.*bubble_27th;
564 :
565 41899595 : case 3:
566 41899595 : return 4.*zeta0*zeta1 - 12.*bubble_27th;
567 :
568 41899595 : case 4:
569 41899595 : return 4.*zeta1*zeta2 - 12.*bubble_27th;
570 :
571 41899595 : case 5:
572 41899595 : return 4.*zeta2*zeta0 - 12.*bubble_27th;
573 :
574 56111800 : case 6:
575 56111800 : return 27.*bubble_27th;
576 :
577 0 : default:
578 0 : libmesh_error_msg("Invalid shape function index i = " << i);
579 : }
580 : }
581 :
582 0 : default:
583 0 : libmesh_error_msg("ERROR: Unsupported 2D element type: " << Utility::enum_to_string(type));
584 : }
585 : } // end case THIRD
586 :
587 :
588 :
589 : // unsupported order
590 0 : default:
591 0 : libmesh_error_msg("ERROR: Unsupported 2D FE order: " << order);
592 : }
593 : #else // LIBMESH_DIM > 1
594 : libmesh_ignore(type, order, i, p);
595 : libmesh_not_implemented();
596 : #endif
597 : }
598 :
599 :
600 :
601 : template <FEFamily T>
602 5069517304 : Real fe_lagrange_2D_shape_deriv(const ElemType type,
603 : const Elem * elem,
604 : const Order order,
605 : const unsigned int i,
606 : const unsigned int j,
607 : const Point & p)
608 : {
609 : #if LIBMESH_DIM > 1
610 :
611 1308068594 : libmesh_assert_less (j, 2);
612 :
613 5069517304 : switch (order)
614 : {
615 : // linear Lagrange shape functions
616 147520014 : case FIRST:
617 : {
618 104764360 : switch (type)
619 : {
620 103526296 : case QUAD4:
621 : case QUADSHELL4:
622 : case QUAD8:
623 : case QUADSHELL8:
624 : case QUAD9:
625 : case QUADSHELL9:
626 : {
627 : // Compute quad shape functions as a tensor-product
628 103526296 : const Real xi = p(0);
629 103526296 : const Real eta = p(1);
630 :
631 31659976 : libmesh_assert_less (i, 4);
632 :
633 : // 0 1 2 3
634 : static const unsigned int i0[] = {0, 1, 1, 0};
635 : static const unsigned int i1[] = {0, 0, 1, 1};
636 :
637 103526296 : switch (j)
638 : {
639 : // d()/dxi
640 51763148 : case 0:
641 51763148 : return (fe_lagrange_1D_linear_shape_deriv(i0[i], 0, xi)*
642 87696308 : fe_lagrange_1D_linear_shape (i1[i], eta));
643 :
644 : // d()/deta
645 51763148 : case 1:
646 51763148 : return (fe_lagrange_1D_linear_shape (i0[i], xi)*
647 69729728 : fe_lagrange_1D_linear_shape_deriv(i1[i], 0, eta));
648 :
649 0 : default:
650 0 : libmesh_error_msg("ERROR: Invalid derivative index j = " << j);
651 : }
652 : }
653 :
654 43568118 : case TRI3:
655 : case TRISHELL3:
656 : case TRI6:
657 : case TRI7:
658 : {
659 11286084 : libmesh_assert_less (i, 3);
660 :
661 11286084 : const Real dzeta0dxi = -1.;
662 11286084 : const Real dzeta1dxi = 1.;
663 11286084 : const Real dzeta2dxi = 0.;
664 :
665 11286084 : const Real dzeta0deta = -1.;
666 11286084 : const Real dzeta1deta = 0.;
667 11286084 : const Real dzeta2deta = 1.;
668 :
669 43568118 : switch (j)
670 : {
671 : // d()/dxi
672 21784059 : case 0:
673 : {
674 16141011 : switch(i)
675 : {
676 1881014 : case 0:
677 1881014 : return dzeta0dxi;
678 :
679 1881014 : case 1:
680 1881014 : return dzeta1dxi;
681 :
682 1881014 : case 2:
683 1881014 : return dzeta2dxi;
684 :
685 0 : default:
686 0 : libmesh_error_msg("Invalid shape function index i = " << i);
687 : }
688 : }
689 : // d()/deta
690 21784059 : case 1:
691 : {
692 16141011 : switch(i)
693 : {
694 1881014 : case 0:
695 1881014 : return dzeta0deta;
696 :
697 1881014 : case 1:
698 1881014 : return dzeta1deta;
699 :
700 1881014 : case 2:
701 1881014 : return dzeta2deta;
702 :
703 0 : default:
704 0 : libmesh_error_msg("Invalid shape function index i = " << i);
705 : }
706 : }
707 0 : default:
708 0 : libmesh_error_msg("ERROR: Invalid derivative index j = " << j);
709 : }
710 : }
711 :
712 425600 : case C0POLYGON:
713 : {
714 : // C0Polygon requires using newer FE APIs
715 425600 : if (!elem)
716 0 : libmesh_error_msg("Code (see stack trace) used an outdated FE function overload.\n"
717 : "Shape functions on a C0Polygon are not defined by its ElemType alone.");
718 :
719 106510 : libmesh_assert(elem->type() == C0POLYGON);
720 :
721 106510 : const C0Polygon & poly = *cast_ptr<const C0Polygon *>(elem);
722 :
723 : // We can't use a small tolerance here, because in
724 : // inverse_map() Newton might hand us intermediate
725 : // iterates outside the polygon.
726 425600 : const auto [s, a, b] = poly.subtriangle_coordinates(p, 100);
727 425600 : if (s == invalid_uint)
728 0 : return 0;
729 106510 : libmesh_assert_less(s, poly.n_subtriangles());
730 :
731 425600 : const auto subtri = poly.subtriangle(s);
732 :
733 : // Find derivatives w.r.t. subtriangle barycentric
734 : // coordinates
735 106510 : Real du_da = 0, du_db = 0;
736 :
737 : // Avoid signed/unsigned comparison warnings
738 425600 : const int nodei = i;
739 425600 : if (nodei == subtri[0])
740 21302 : du_da = du_db = -1;
741 340480 : else if (nodei == subtri[1])
742 21302 : du_da = 1;
743 255360 : else if (nodei == subtri[2])
744 21302 : du_db = 1;
745 : else
746 : // Basis function i is not supported on p's subtriangle
747 42604 : return 0;
748 :
749 : // We want to return derivatives with respect to xi and
750 : // eta in master space for the polygon, but what we
751 : // calculated above are with respect to xi and eta
752 : // coordinates for a master *triangle*. We need to
753 : // convert from one to the other.
754 :
755 255360 : const auto master_points = poly.master_subtriangle(s);
756 :
757 255360 : const Real dxi_da = master_points[1](0) - master_points[0](0);
758 255360 : const Real dxi_db = master_points[2](0) - master_points[0](0);
759 255360 : const Real deta_da = master_points[1](1) - master_points[0](1);
760 255360 : const Real deta_db = master_points[2](1) - master_points[0](1);
761 255360 : const Real jac = dxi_da*deta_db - dxi_db*deta_da;
762 :
763 255360 : switch (j)
764 : {
765 : // d()/dxi
766 127680 : case 0:
767 : {
768 127680 : const Real da_dxi = deta_db / jac;
769 127680 : const Real db_dxi = -deta_da / jac;
770 127680 : return du_da*da_dxi + du_db*db_dxi;
771 : }
772 : // d()/deta
773 127680 : case 1:
774 : {
775 127680 : const Real da_deta = -dxi_db / jac;
776 127680 : const Real db_deta = dxi_da / jac;
777 127680 : return du_da*da_deta + du_db*db_deta;
778 : }
779 0 : default:
780 0 : libmesh_error_msg("ERROR: Invalid derivative index j = " << j);
781 : }
782 : }
783 :
784 0 : default:
785 0 : libmesh_error_msg("ERROR: Unsupported 2D element type: " << Utility::enum_to_string(type));
786 : }
787 : }
788 :
789 :
790 : // quadratic Lagrange shape functions
791 4447514318 : case SECOND:
792 : {
793 3309716598 : switch (type)
794 : {
795 17481008 : case QUAD8:
796 : case QUADSHELL8:
797 : {
798 17481008 : const Real xi = p(0);
799 17481008 : const Real eta = p(1);
800 :
801 4668912 : libmesh_assert_less (i, 8);
802 :
803 17481008 : switch (j)
804 : {
805 : // d/dxi
806 8740504 : case 0:
807 8740504 : switch (i)
808 : {
809 1092563 : case 0:
810 1676185 : return .25*(1. - eta)*((1. - xi)*(-1.) +
811 1092563 : (-1.)*(-1. - xi - eta));
812 :
813 1092563 : case 1:
814 1676185 : return .25*(1. - eta)*((1. + xi)*(1.) +
815 1092563 : (1.)*(-1. + xi - eta));
816 :
817 1092563 : case 2:
818 1676185 : return .25*(1. + eta)*((1. + xi)*(1.) +
819 1092563 : (1.)*(-1. + xi + eta));
820 :
821 1092563 : case 3:
822 1676185 : return .25*(1. + eta)*((1. - xi)*(-1.) +
823 1092563 : (-1.)*(-1. - xi + eta));
824 :
825 1092563 : case 4:
826 1092563 : return .5*(-2.*xi)*(1. - eta);
827 :
828 1092563 : case 5:
829 1092563 : return .5*(1.)*(1. - eta*eta);
830 :
831 1092563 : case 6:
832 1092563 : return .5*(-2.*xi)*(1. + eta);
833 :
834 1092563 : case 7:
835 1092563 : return .5*(-1.)*(1. - eta*eta);
836 :
837 0 : default:
838 0 : libmesh_error_msg("Invalid shape function index i = " << i);
839 : }
840 :
841 : // d/deta
842 8740504 : case 1:
843 8740504 : switch (i)
844 : {
845 1092563 : case 0:
846 1676185 : return .25*(1. - xi)*((1. - eta)*(-1.) +
847 1092563 : (-1.)*(-1. - xi - eta));
848 :
849 1092563 : case 1:
850 1676185 : return .25*(1. + xi)*((1. - eta)*(-1.) +
851 1092563 : (-1.)*(-1. + xi - eta));
852 :
853 1092563 : case 2:
854 1676185 : return .25*(1. + xi)*((1. + eta)*(1.) +
855 1092563 : (1.)*(-1. + xi + eta));
856 :
857 1092563 : case 3:
858 1676185 : return .25*(1. - xi)*((1. + eta)*(1.) +
859 1092563 : (1.)*(-1. - xi + eta));
860 :
861 1092563 : case 4:
862 1092563 : return .5*(1. - xi*xi)*(-1.);
863 :
864 1092563 : case 5:
865 1092563 : return .5*(1. + xi)*(-2.*eta);
866 :
867 1092563 : case 6:
868 1092563 : return .5*(1. - xi*xi)*(1.);
869 :
870 1092563 : case 7:
871 1092563 : return .5*(1. - xi)*(-2.*eta);
872 :
873 0 : default:
874 0 : libmesh_error_msg("Invalid shape function index i = " << i);
875 : }
876 :
877 0 : default:
878 0 : libmesh_error_msg("ERROR: Invalid derivative index j = " << j);
879 : }
880 : }
881 :
882 425703114 : case QUAD4:
883 7866 : libmesh_assert_msg(T == L2_LAGRANGE,
884 : "High order on first order elements only supported for L2 families");
885 : libmesh_fallthrough();
886 158957046 : case QUAD9:
887 : case QUADSHELL9:
888 : {
889 : // Compute quad shape functions as a tensor-product
890 584668026 : const Real xi = p(0);
891 584668026 : const Real eta = p(1);
892 :
893 158972778 : libmesh_assert_less (i, 9);
894 :
895 : // 0 1 2 3 4 5 6 7 8
896 : static const unsigned int i0[] = {0, 1, 1, 0, 2, 1, 2, 0, 2};
897 : static const unsigned int i1[] = {0, 0, 1, 1, 0, 2, 1, 2, 2};
898 :
899 584668026 : switch (j)
900 : {
901 : // d()/dxi
902 292334013 : case 0:
903 292334013 : return (fe_lagrange_1D_quadratic_shape_deriv(i0[i], 0, xi)*
904 370829421 : fe_lagrange_1D_quadratic_shape (i1[i], eta));
905 :
906 : // d()/deta
907 292334013 : case 1:
908 292334013 : return (fe_lagrange_1D_quadratic_shape (i0[i], xi)*
909 505181637 : fe_lagrange_1D_quadratic_shape_deriv(i1[i], 0, eta));
910 :
911 0 : default:
912 0 : libmesh_error_msg("ERROR: Invalid derivative index j = " << j);
913 : }
914 : }
915 :
916 2868932688 : case TRI3:
917 5484 : libmesh_assert_msg(T == L2_LAGRANGE,
918 : "High order on first order elements only supported for L2 families");
919 : libmesh_fallthrough();
920 976334076 : case TRI6:
921 : case TRI7:
922 : {
923 976438080 : libmesh_assert_less (i, 6);
924 :
925 3845365284 : const Real zeta1 = p(0);
926 3845365284 : const Real zeta2 = p(1);
927 3845365284 : const Real zeta0 = 1. - zeta1 - zeta2;
928 :
929 976438080 : const Real dzeta0dxi = -1.;
930 976438080 : const Real dzeta1dxi = 1.;
931 976438080 : const Real dzeta2dxi = 0.;
932 :
933 976438080 : const Real dzeta0deta = -1.;
934 976438080 : const Real dzeta1deta = 0.;
935 976438080 : const Real dzeta2deta = 1.;
936 :
937 3845365284 : switch(j)
938 : {
939 1922682642 : case 0:
940 : {
941 1922682642 : switch(i)
942 : {
943 320447107 : case 0:
944 320447107 : return (4.*zeta0-1.)*dzeta0dxi;
945 :
946 320447107 : case 1:
947 320447107 : return (4.*zeta1-1.)*dzeta1dxi;
948 :
949 320447107 : case 2:
950 320447107 : return (4.*zeta2-1.)*dzeta2dxi;
951 :
952 320447107 : case 3:
953 320447107 : return 4.*zeta1*dzeta0dxi + 4.*zeta0*dzeta1dxi;
954 :
955 320447107 : case 4:
956 320447107 : return 4.*zeta2*dzeta1dxi + 4.*zeta1*dzeta2dxi;
957 :
958 320447107 : case 5:
959 320447107 : return 4.*zeta2*dzeta0dxi + 4*zeta0*dzeta2dxi;
960 :
961 0 : default:
962 0 : libmesh_error_msg("Invalid shape function index i = " << i);
963 : }
964 : }
965 :
966 1922682642 : case 1:
967 : {
968 1922682642 : switch(i)
969 : {
970 320447107 : case 0:
971 320447107 : return (4.*zeta0-1.)*dzeta0deta;
972 :
973 320447107 : case 1:
974 320447107 : return (4.*zeta1-1.)*dzeta1deta;
975 :
976 320447107 : case 2:
977 320447107 : return (4.*zeta2-1.)*dzeta2deta;
978 :
979 320447107 : case 3:
980 320447107 : return 4.*zeta1*dzeta0deta + 4.*zeta0*dzeta1deta;
981 :
982 320447107 : case 4:
983 320447107 : return 4.*zeta2*dzeta1deta + 4.*zeta1*dzeta2deta;
984 :
985 320447107 : case 5:
986 320447107 : return 4.*zeta2*dzeta0deta + 4*zeta0*dzeta2deta;
987 :
988 0 : default:
989 0 : libmesh_error_msg("Invalid shape function index i = " << i);
990 : }
991 : }
992 0 : default:
993 0 : libmesh_error_msg("ERROR: Invalid derivative index j = " << j);
994 : }
995 : }
996 :
997 0 : default:
998 0 : libmesh_error_msg("ERROR: Unsupported 2D element type: " << Utility::enum_to_string(type));
999 : }
1000 : }
1001 :
1002 :
1003 :
1004 : // "cubic" (one cubic bubble) Lagrange shape functions on TRI7
1005 474482972 : case THIRD:
1006 : {
1007 474482972 : switch (type)
1008 : {
1009 474482972 : case TRI7:
1010 : {
1011 124936254 : libmesh_assert_less (i, 7);
1012 :
1013 474482972 : const Real zeta1 = p(0);
1014 474482972 : const Real zeta2 = p(1);
1015 474482972 : const Real zeta0 = 1. - zeta1 - zeta2;
1016 : // const Real bubble_27th = zeta0*zeta1*zeta2;
1017 :
1018 124936254 : const Real dzeta0dxi = -1.;
1019 124936254 : const Real dzeta1dxi = 1.;
1020 124936254 : const Real dzeta2dxi = 0.;
1021 474482972 : const Real dbubbledxi = zeta2 * (1. - 2.*zeta1 - zeta2);
1022 :
1023 124936254 : const Real dzeta0deta = -1.;
1024 124936254 : const Real dzeta1deta = 0.;
1025 124936254 : const Real dzeta2deta = 1.;
1026 474482972 : const Real dbubbledeta= zeta1 * (1. - zeta1 - 2.*zeta2);
1027 :
1028 474482972 : switch(j)
1029 : {
1030 237241486 : case 0:
1031 : {
1032 237241486 : switch(i)
1033 : {
1034 32754748 : case 0:
1035 32754748 : return (4.*zeta0-1.)*dzeta0dxi + 3.*dbubbledxi;
1036 :
1037 32754748 : case 1:
1038 32754748 : return (4.*zeta1-1.)*dzeta1dxi + 3.*dbubbledxi;
1039 :
1040 32754748 : case 2:
1041 32754748 : return (4.*zeta2-1.)*dzeta2dxi + 3.*dbubbledxi;
1042 :
1043 32754748 : case 3:
1044 32754748 : return 4.*zeta1*dzeta0dxi + 4.*zeta0*dzeta1dxi - 12.*dbubbledxi;
1045 :
1046 32754748 : case 4:
1047 32754748 : return 4.*zeta2*dzeta1dxi + 4.*zeta1*dzeta2dxi - 12.*dbubbledxi;
1048 :
1049 32754748 : case 5:
1050 32754748 : return 4.*zeta2*dzeta0dxi + 4*zeta0*dzeta2dxi - 12.*dbubbledxi;
1051 :
1052 40712998 : case 6:
1053 40712998 : return 27.*dbubbledxi;
1054 :
1055 0 : default:
1056 0 : libmesh_error_msg("Invalid shape function index i = " << i);
1057 : }
1058 : }
1059 :
1060 237241486 : case 1:
1061 : {
1062 237241486 : switch(i)
1063 : {
1064 32754748 : case 0:
1065 32754748 : return (4.*zeta0-1.)*dzeta0deta + 3.*dbubbledeta;
1066 :
1067 32754748 : case 1:
1068 32754748 : return (4.*zeta1-1.)*dzeta1deta + 3.*dbubbledeta;
1069 :
1070 32754748 : case 2:
1071 32754748 : return (4.*zeta2-1.)*dzeta2deta + 3.*dbubbledeta;
1072 :
1073 32754748 : case 3:
1074 32754748 : return 4.*zeta1*dzeta0deta + 4.*zeta0*dzeta1deta - 12.*dbubbledeta;
1075 :
1076 32754748 : case 4:
1077 32754748 : return 4.*zeta2*dzeta1deta + 4.*zeta1*dzeta2deta - 12.*dbubbledeta;
1078 :
1079 32754748 : case 5:
1080 32754748 : return 4.*zeta2*dzeta0deta + 4*zeta0*dzeta2deta - 12.*dbubbledeta;
1081 :
1082 40712998 : case 6:
1083 40712998 : return 27.*dbubbledeta;
1084 :
1085 0 : default:
1086 0 : libmesh_error_msg("Invalid shape function index i = " << i);
1087 : }
1088 : }
1089 0 : default:
1090 0 : libmesh_error_msg("ERROR: Invalid derivative index j = " << j);
1091 : }
1092 : }
1093 :
1094 0 : default:
1095 0 : libmesh_error_msg("ERROR: Unsupported 2D element type: " << Utility::enum_to_string(type));
1096 : }
1097 : } // end case THIRD
1098 :
1099 :
1100 :
1101 : // unsupported order
1102 0 : default:
1103 0 : libmesh_error_msg("ERROR: Unsupported 2D FE order: " << order);
1104 : }
1105 : #else // LIBMESH_DIM > 1
1106 : libmesh_ignore(type, order, i, j, p);
1107 : libmesh_not_implemented();
1108 : #endif
1109 : }
1110 :
1111 :
1112 :
1113 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1114 :
1115 : template <FEFamily T>
1116 71791671 : Real fe_lagrange_2D_shape_second_deriv(const ElemType type,
1117 : const Elem *,
1118 : const Order order,
1119 : const unsigned int i,
1120 : const unsigned int j,
1121 : const Point & p)
1122 : {
1123 : #if LIBMESH_DIM > 1
1124 :
1125 : // j = 0 ==> d^2 phi / dxi^2
1126 : // j = 1 ==> d^2 phi / dxi deta
1127 : // j = 2 ==> d^2 phi / deta^2
1128 25579572 : libmesh_assert_less (j, 3);
1129 :
1130 71791671 : switch (order)
1131 : {
1132 : // linear Lagrange shape functions
1133 732303 : case FIRST:
1134 : {
1135 732303 : switch (type)
1136 : {
1137 456588 : case QUAD4:
1138 : case QUADSHELL4:
1139 : case QUAD8:
1140 : case QUADSHELL8:
1141 : case QUAD9:
1142 : case QUADSHELL9:
1143 : {
1144 : // Compute quad shape functions as a tensor-product
1145 157200 : const Real xi = p(0);
1146 157200 : const Real eta = p(1);
1147 :
1148 157200 : libmesh_assert_less (i, 4);
1149 :
1150 : // 0 1 2 3
1151 : static const unsigned int i0[] = {0, 1, 1, 0};
1152 : static const unsigned int i1[] = {0, 0, 1, 1};
1153 :
1154 456588 : switch (j)
1155 : {
1156 : // d^2() / dxi^2
1157 52400 : case 0:
1158 52400 : return 0.;
1159 :
1160 : // d^2() / dxi deta
1161 152196 : case 1:
1162 152196 : return (fe_lagrange_1D_linear_shape_deriv(i0[i], 0, xi)*
1163 202094 : fe_lagrange_1D_linear_shape_deriv(i1[i], 0, eta));
1164 :
1165 : // d^2() / deta^2
1166 52400 : case 2:
1167 52400 : return 0.;
1168 :
1169 0 : default:
1170 0 : libmesh_error_msg("ERROR: Invalid derivative index j = " << j);
1171 : }
1172 : }
1173 :
1174 : // All second derivatives for linear triangles are zero.
1175 101913 : case TRI3:
1176 : case TRISHELL3:
1177 : case TRI6:
1178 : case TRI7:
1179 :
1180 : // All second derivatives for piecewise-linear polygons are
1181 : // zero or dirac-type distributions, but we can't put the
1182 : // latter in a Real, so beware when integrating...
1183 : case C0POLYGON:
1184 : {
1185 101913 : return 0.;
1186 : }
1187 :
1188 0 : default:
1189 0 : libmesh_error_msg("ERROR: Unsupported 2D element type: " << Utility::enum_to_string(type));
1190 :
1191 : } // end switch (type)
1192 : } // end case FIRST
1193 :
1194 :
1195 : // quadratic Lagrange shape functions
1196 31287258 : case SECOND:
1197 : {
1198 31287258 : switch (type)
1199 : {
1200 1318176 : case QUAD8:
1201 : case QUADSHELL8:
1202 : {
1203 1318176 : const Real xi = p(0);
1204 1318176 : const Real eta = p(1);
1205 :
1206 439488 : libmesh_assert_less (j, 3);
1207 :
1208 1318176 : switch (j)
1209 : {
1210 : // d^2() / dxi^2
1211 439392 : case 0:
1212 : {
1213 439392 : switch (i)
1214 : {
1215 109848 : case 0:
1216 : case 1:
1217 109848 : return 0.5*(1.-eta);
1218 :
1219 109848 : case 2:
1220 : case 3:
1221 109848 : return 0.5*(1.+eta);
1222 :
1223 54924 : case 4:
1224 54924 : return eta - 1.;
1225 :
1226 36624 : case 5:
1227 : case 7:
1228 36624 : return 0.0;
1229 :
1230 54924 : case 6:
1231 54924 : return -1. - eta;
1232 :
1233 0 : default:
1234 0 : libmesh_error_msg("Invalid shape function index i = " << i);
1235 : }
1236 : }
1237 :
1238 : // d^2() / dxi deta
1239 439392 : case 1:
1240 : {
1241 439392 : switch (i)
1242 : {
1243 54924 : case 0:
1244 54924 : return 0.25*( 1. - 2.*xi - 2.*eta);
1245 :
1246 54924 : case 1:
1247 54924 : return 0.25*(-1. - 2.*xi + 2.*eta);
1248 :
1249 54924 : case 2:
1250 54924 : return 0.25*( 1. + 2.*xi + 2.*eta);
1251 :
1252 54924 : case 3:
1253 54924 : return 0.25*(-1. + 2.*xi - 2.*eta);
1254 :
1255 18312 : case 4:
1256 18312 : return xi;
1257 :
1258 54924 : case 5:
1259 54924 : return -eta;
1260 :
1261 54924 : case 6:
1262 54924 : return -xi;
1263 :
1264 54924 : case 7:
1265 54924 : return eta;
1266 :
1267 0 : default:
1268 0 : libmesh_error_msg("Invalid shape function index i = " << i);
1269 : }
1270 : }
1271 :
1272 : // d^2() / deta^2
1273 439392 : case 2:
1274 : {
1275 439392 : switch (i)
1276 : {
1277 109848 : case 0:
1278 : case 3:
1279 109848 : return 0.5*(1.-xi);
1280 :
1281 109848 : case 1:
1282 : case 2:
1283 109848 : return 0.5*(1.+xi);
1284 :
1285 36624 : case 4:
1286 : case 6:
1287 36624 : return 0.0;
1288 :
1289 54924 : case 5:
1290 54924 : return -1.0 - xi;
1291 :
1292 54924 : case 7:
1293 54924 : return xi - 1.0;
1294 :
1295 0 : default:
1296 0 : libmesh_error_msg("Invalid shape function index i = " << i);
1297 : }
1298 : }
1299 :
1300 0 : default:
1301 0 : libmesh_error_msg("ERROR: Invalid derivative index j = " << j);
1302 : } // end switch (j)
1303 : } // end case QUAD8
1304 :
1305 4339197 : case QUAD4:
1306 11799 : libmesh_assert_msg(T == L2_LAGRANGE,
1307 : "High order on first order elements only supported for L2 families");
1308 : libmesh_fallthrough();
1309 2428218 : case QUAD9:
1310 : case QUADSHELL9:
1311 : {
1312 : // Compute QUAD9 second derivatives as tensor product
1313 6779214 : const Real xi = p(0);
1314 6779214 : const Real eta = p(1);
1315 :
1316 2451816 : libmesh_assert_less (i, 9);
1317 :
1318 : // 0 1 2 3 4 5 6 7 8
1319 : static const unsigned int i0[] = {0, 1, 1, 0, 2, 1, 2, 0, 2};
1320 : static const unsigned int i1[] = {0, 0, 1, 1, 0, 2, 1, 2, 2};
1321 :
1322 6779214 : switch (j)
1323 : {
1324 : // d^2() / dxi^2
1325 2259738 : case 0:
1326 2259738 : return (fe_lagrange_1D_quadratic_shape_second_deriv(i0[i], 0, xi)*
1327 2259738 : fe_lagrange_1D_quadratic_shape (i1[i], eta));
1328 :
1329 : // d^2() / dxi deta
1330 2259738 : case 1:
1331 2259738 : return (fe_lagrange_1D_quadratic_shape_deriv(i0[i], 0, xi)*
1332 3702204 : fe_lagrange_1D_quadratic_shape_deriv(i1[i], 0, eta));
1333 :
1334 : // d^2() / deta^2
1335 2259738 : case 2:
1336 2259738 : return (fe_lagrange_1D_quadratic_shape (i0[i], xi)*
1337 2740560 : fe_lagrange_1D_quadratic_shape_second_deriv(i1[i], 0, eta));
1338 :
1339 0 : default:
1340 0 : libmesh_error_msg("ERROR: Invalid derivative index j = " << j);
1341 : } // end switch (j)
1342 : } // end case QUAD9
1343 :
1344 14969826 : case TRI3:
1345 8226 : libmesh_assert_msg(T == L2_LAGRANGE,
1346 : "High order on first order elements only supported for L2 families");
1347 : libmesh_fallthrough();
1348 8201862 : case TRI6:
1349 : case TRI7:
1350 : {
1351 8228268 : const Real dzeta0dxi = -1.;
1352 8228268 : const Real dzeta1dxi = 1.;
1353 8228268 : const Real dzeta2dxi = 0.;
1354 :
1355 8228268 : const Real dzeta0deta = -1.;
1356 8228268 : const Real dzeta1deta = 0.;
1357 8228268 : const Real dzeta2deta = 1.;
1358 :
1359 8228268 : libmesh_assert_less (j, 3);
1360 :
1361 23189868 : switch (j)
1362 : {
1363 : // d^2() / dxi^2
1364 7729956 : case 0:
1365 : {
1366 7729956 : switch (i)
1367 : {
1368 457126 : case 0:
1369 457126 : return 4.*dzeta0dxi*dzeta0dxi;
1370 :
1371 457126 : case 1:
1372 457126 : return 4.*dzeta1dxi*dzeta1dxi;
1373 :
1374 457126 : case 2:
1375 457126 : return 4.*dzeta2dxi*dzeta2dxi;
1376 :
1377 457126 : case 3:
1378 457126 : return 8.*dzeta0dxi*dzeta1dxi;
1379 :
1380 457126 : case 4:
1381 457126 : return 8.*dzeta1dxi*dzeta2dxi;
1382 :
1383 457126 : case 5:
1384 457126 : return 8.*dzeta0dxi*dzeta2dxi;
1385 :
1386 0 : default:
1387 0 : libmesh_error_msg("Invalid shape function index i = " << i);
1388 : }
1389 : }
1390 :
1391 : // d^2() / dxi deta
1392 7729956 : case 1:
1393 : {
1394 7729956 : switch (i)
1395 : {
1396 457126 : case 0:
1397 457126 : return 4.*dzeta0dxi*dzeta0deta;
1398 :
1399 457126 : case 1:
1400 457126 : return 4.*dzeta1dxi*dzeta1deta;
1401 :
1402 457126 : case 2:
1403 457126 : return 4.*dzeta2dxi*dzeta2deta;
1404 :
1405 457126 : case 3:
1406 457126 : return 4.*dzeta1deta*dzeta0dxi + 4.*dzeta0deta*dzeta1dxi;
1407 :
1408 457126 : case 4:
1409 457126 : return 4.*dzeta2deta*dzeta1dxi + 4.*dzeta1deta*dzeta2dxi;
1410 :
1411 457126 : case 5:
1412 457126 : return 4.*dzeta2deta*dzeta0dxi + 4.*dzeta0deta*dzeta2dxi;
1413 :
1414 0 : default:
1415 0 : libmesh_error_msg("Invalid shape function index i = " << i);
1416 : }
1417 : }
1418 :
1419 : // d^2() / deta^2
1420 7729956 : case 2:
1421 : {
1422 7729956 : switch (i)
1423 : {
1424 457126 : case 0:
1425 457126 : return 4.*dzeta0deta*dzeta0deta;
1426 :
1427 457126 : case 1:
1428 457126 : return 4.*dzeta1deta*dzeta1deta;
1429 :
1430 457126 : case 2:
1431 457126 : return 4.*dzeta2deta*dzeta2deta;
1432 :
1433 457126 : case 3:
1434 457126 : return 8.*dzeta0deta*dzeta1deta;
1435 :
1436 457126 : case 4:
1437 457126 : return 8.*dzeta1deta*dzeta2deta;
1438 :
1439 457126 : case 5:
1440 457126 : return 8.*dzeta0deta*dzeta2deta;
1441 :
1442 0 : default:
1443 0 : libmesh_error_msg("Invalid shape function index i = " << i);
1444 : }
1445 : }
1446 :
1447 0 : default:
1448 0 : libmesh_error_msg("ERROR: Invalid derivative index j = " << j);
1449 : } // end switch (j)
1450 : } // end case TRI6+TRI7
1451 :
1452 0 : default:
1453 0 : libmesh_error_msg("ERROR: Unsupported 2D element type: " << Utility::enum_to_string(type));
1454 : }
1455 : } // end case SECOND
1456 :
1457 :
1458 :
1459 : // "cubic" (one cubic bubble) Lagrange shape functions on TRI7
1460 39772110 : case THIRD:
1461 : {
1462 39772110 : switch (type)
1463 : {
1464 25571223 : case TRI3:
1465 0 : libmesh_assert_msg(T == L2_LAGRANGE,
1466 : "High order on first order elements only supported for L2 families");
1467 : libmesh_fallthrough();
1468 14190660 : case TRI6:
1469 : case TRI7:
1470 : {
1471 39772110 : const Real zeta1 = p(0);
1472 39772110 : const Real zeta2 = p(1);
1473 : // const Real zeta0 = 1. - zeta1 - zeta2;
1474 :
1475 14200887 : const Real dzeta0dxi = -1.;
1476 14200887 : const Real dzeta1dxi = 1.;
1477 14200887 : const Real dzeta2dxi = 0.;
1478 : // const Real dbubbledxi = zeta2 * (1. - 2.*zeta1 - zeta2);
1479 39772110 : const Real d2bubbledxi2 = -2. * zeta2;
1480 :
1481 14200887 : const Real dzeta0deta = -1.;
1482 14200887 : const Real dzeta1deta = 0.;
1483 14200887 : const Real dzeta2deta = 1.;
1484 : // const Real dbubbledeta= zeta1 * (1. - zeta1 - 2.*zeta2);
1485 39772110 : const Real d2bubbledeta2 = -2. * zeta1;
1486 :
1487 39772110 : const Real d2bubbledxideta = (1. - 2.*zeta1 - 2.*zeta2);
1488 :
1489 14200887 : libmesh_assert_less (j, 3);
1490 :
1491 39772110 : switch (j)
1492 : {
1493 : // d^2() / dxi^2
1494 13257370 : case 0:
1495 : {
1496 13257370 : switch (i)
1497 : {
1498 1771285 : case 0:
1499 1771285 : return 4.*dzeta0dxi*dzeta0dxi + 3.*d2bubbledxi2;
1500 :
1501 1771285 : case 1:
1502 1771285 : return 4.*dzeta1dxi*dzeta1dxi + 3.*d2bubbledxi2;
1503 :
1504 1771285 : case 2:
1505 1771285 : return 4.*dzeta2dxi*dzeta2dxi + 3.*d2bubbledxi2;
1506 :
1507 1771285 : case 3:
1508 1771285 : return 8.*dzeta0dxi*dzeta1dxi - 12.*d2bubbledxi2;
1509 :
1510 1771285 : case 4:
1511 1771285 : return 8.*dzeta1dxi*dzeta2dxi - 12.*d2bubbledxi2;
1512 :
1513 1771285 : case 5:
1514 1771285 : return 8.*dzeta0dxi*dzeta2dxi - 12.*d2bubbledxi2;
1515 :
1516 2629660 : case 6:
1517 2629660 : return 27.*d2bubbledxi2;
1518 :
1519 0 : default:
1520 0 : libmesh_error_msg("Invalid shape function index i = " << i);
1521 : }
1522 : }
1523 :
1524 : // d^2() / dxi deta
1525 13257370 : case 1:
1526 : {
1527 13257370 : switch (i)
1528 : {
1529 1771285 : case 0:
1530 1771285 : return 4.*dzeta0dxi*dzeta0deta + 3.*d2bubbledxideta;
1531 :
1532 1771285 : case 1:
1533 1771285 : return 4.*dzeta1dxi*dzeta1deta + 3.*d2bubbledxideta;
1534 :
1535 1771285 : case 2:
1536 1771285 : return 4.*dzeta2dxi*dzeta2deta + 3.*d2bubbledxideta;
1537 :
1538 1771285 : case 3:
1539 1771285 : return 4.*dzeta1deta*dzeta0dxi + 4.*dzeta0deta*dzeta1dxi - 12.*d2bubbledxideta;
1540 :
1541 1771285 : case 4:
1542 1771285 : return 4.*dzeta2deta*dzeta1dxi + 4.*dzeta1deta*dzeta2dxi - 12.*d2bubbledxideta;
1543 :
1544 1771285 : case 5:
1545 1771285 : return 4.*dzeta2deta*dzeta0dxi + 4.*dzeta0deta*dzeta2dxi - 12.*d2bubbledxideta;
1546 :
1547 2629660 : case 6:
1548 2629660 : return 27.*d2bubbledxideta;
1549 :
1550 0 : default:
1551 0 : libmesh_error_msg("Invalid shape function index i = " << i);
1552 : }
1553 : }
1554 :
1555 : // d^2() / deta^2
1556 13257370 : case 2:
1557 : {
1558 13257370 : switch (i)
1559 : {
1560 1771285 : case 0:
1561 1771285 : return 4.*dzeta0deta*dzeta0deta + 3.*d2bubbledeta2;
1562 :
1563 1771285 : case 1:
1564 1771285 : return 4.*dzeta1deta*dzeta1deta + 3.*d2bubbledeta2;
1565 :
1566 1771285 : case 2:
1567 1771285 : return 4.*dzeta2deta*dzeta2deta + 3.*d2bubbledeta2;
1568 :
1569 1771285 : case 3:
1570 1771285 : return 8.*dzeta0deta*dzeta1deta - 12.*d2bubbledeta2;
1571 :
1572 1771285 : case 4:
1573 1771285 : return 8.*dzeta1deta*dzeta2deta - 12.*d2bubbledeta2;
1574 :
1575 1771285 : case 5:
1576 1771285 : return 8.*dzeta0deta*dzeta2deta - 12.*d2bubbledeta2;
1577 :
1578 2629660 : case 6:
1579 2629660 : return 27.*d2bubbledeta2;
1580 :
1581 0 : default:
1582 0 : libmesh_error_msg("Invalid shape function index i = " << i);
1583 : }
1584 : }
1585 :
1586 0 : default:
1587 0 : libmesh_error_msg("ERROR: Invalid derivative index j = " << j);
1588 : } // end switch (j)
1589 : } // end case TRI6+TRI7
1590 :
1591 0 : default:
1592 0 : libmesh_error_msg("ERROR: Unsupported 2D element type: " << Utility::enum_to_string(type));
1593 : }
1594 : } // end case THIRD
1595 :
1596 : // unsupported order
1597 0 : default:
1598 0 : libmesh_error_msg("ERROR: Unsupported 2D FE order: " << order);
1599 :
1600 : } // end switch (order)
1601 :
1602 : #else // LIBMESH_DIM > 1
1603 : libmesh_ignore(type, order, i, j, p);
1604 : libmesh_not_implemented();
1605 : #endif
1606 : }
1607 :
1608 : #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
1609 :
1610 : } // anonymous namespace
|
