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/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 692933346 : LIBMESH_DEFAULT_VECTORIZED_FE(2,LAGRANGE)
69 4072422 : LIBMESH_DEFAULT_VECTORIZED_FE(2,L2_LAGRANGE)
70 :
71 :
72 : template <>
73 801607763 : Real FE<2,LAGRANGE>::shape(const ElemType type,
74 : const Order order,
75 : const unsigned int i,
76 : const Point & p)
77 : {
78 801607763 : 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 1704317561 : 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 166266064 : libmesh_assert(elem);
101 :
102 : // call the orientation-independent shape functions
103 2013986183 : 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 7215116 : 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 8428924 : 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 1133853900 : 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 105345897 : libmesh_assert(elem);
130 1310318300 : 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 851765168 : 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 851765168 : 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 679542334 : 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 62143074 : libmesh_assert(elem);
180 :
181 : // call the orientation-independent shape functions
182 803804958 : 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 1719616 : 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 1999696 : 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 15335067878 : 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 1130308006 : libmesh_assert(elem);
211 17591582478 : 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 150865560 : 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 150865560 : 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 3916107 : 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 336354 : libmesh_assert(elem);
265 :
266 : // call the orientation-independent shape functions
267 4588815 : 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 974544 : 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 1132296 : 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 151963722 : 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 12044028 : libmesh_assert(elem);
296 176051778 : 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 3646994340 : 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 3646994340 : switch (order)
335 : {
336 : // linear Lagrange shape functions
337 1236011793 : case FIRST:
338 : {
339 1121034684 : switch (type)
340 : {
341 1069182392 : 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 1069182392 : const Real xi = p(0);
350 1069182392 : const Real eta = p(1);
351 :
352 104683048 : 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 1069182392 : return (fe_lagrange_1D_linear_shape(i0[i], xi)*
359 1069182392 : fe_lagrange_1D_linear_shape(i1[i], eta));
360 : }
361 :
362 166125201 : case TRI3:
363 : case TRISHELL3:
364 : case TRI6:
365 : case TRI7:
366 : {
367 166125201 : const Real zeta1 = p(0);
368 166125201 : const Real zeta2 = p(1);
369 166125201 : const Real zeta0 = 1. - zeta1 - zeta2;
370 :
371 14899806 : libmesh_assert_less (i, 3);
372 :
373 166125201 : switch(i)
374 : {
375 4966602 : case 0:
376 4966602 : return zeta0;
377 :
378 55375067 : case 1:
379 55375067 : return zeta1;
380 :
381 55375067 : case 2:
382 55375067 : return zeta2;
383 :
384 0 : default:
385 0 : libmesh_error_msg("Invalid shape function index i = " << i);
386 : }
387 : }
388 :
389 704200 : case C0POLYGON:
390 : {
391 : // C0Polygon requires using newer FE APIs
392 704200 : 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 704200 : const auto [s, a, b] = poly.subtriangle_coordinates(p, 100);
404 704200 : if (s == invalid_uint)
405 0 : return 0;
406 67025 : libmesh_assert_less(s, poly.n_subtriangles());
407 :
408 704200 : const auto subtri = poly.subtriangle(s);
409 :
410 : // Avoid signed/unsigned comparison warnings
411 704200 : const int nodei = i;
412 704200 : if (nodei == subtri[0])
413 140840 : return 1-a-b;
414 563360 : if (nodei == subtri[1])
415 140840 : return a;
416 422520 : if (nodei == subtri[2])
417 140840 : 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 1418870298 : case SECOND:
431 : {
432 1302423609 : switch (type)
433 : {
434 31236744 : case QUAD8:
435 : case QUADSHELL8:
436 : {
437 31236744 : const Real xi = p(0);
438 31236744 : const Real eta = p(1);
439 :
440 2687088 : libmesh_assert_less (i, 8);
441 :
442 31236744 : switch (i)
443 : {
444 3904593 : case 0:
445 3904593 : return .25*(1. - xi)*(1. - eta)*(-1. - xi - eta);
446 :
447 3904593 : case 1:
448 3904593 : return .25*(1. + xi)*(1. - eta)*(-1. + xi - eta);
449 :
450 3904593 : case 2:
451 3904593 : return .25*(1. + xi)*(1. + eta)*(-1. + xi + eta);
452 :
453 3904593 : case 3:
454 3904593 : return .25*(1. - xi)*(1. + eta)*(-1. - xi + eta);
455 :
456 3904593 : case 4:
457 3904593 : return .5*(1. - xi*xi)*(1. - eta);
458 :
459 3904593 : case 5:
460 3904593 : return .5*(1. + xi)*(1. - eta*eta);
461 :
462 3904593 : case 6:
463 3904593 : return .5*(1. - xi*xi)*(1. + eta);
464 :
465 3904593 : case 7:
466 3904593 : 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 748143342 : 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 80335647 : case QUAD9:
478 : case QUADSHELL9:
479 : {
480 :
481 : // Compute quad shape functions as a tensor-product
482 828483732 : const Real xi = p(0);
483 828483732 : const Real eta = p(1);
484 :
485 80345133 : 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 828483732 : return (fe_lagrange_1D_quadratic_shape(i0[i], xi)*
492 828483732 : fe_lagrange_1D_quadratic_shape(i1[i], eta));
493 : }
494 :
495 502500960 : 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 56404098 : case TRI6:
500 : case TRI7:
501 : {
502 559149822 : const Real zeta1 = p(0);
503 559149822 : const Real zeta2 = p(1);
504 559149822 : const Real zeta0 = 1. - zeta1 - zeta2;
505 :
506 56652444 : libmesh_assert_less (i, 6);
507 :
508 559149822 : switch(i)
509 : {
510 93191637 : case 0:
511 93191637 : return 2.*zeta0*(zeta0-0.5);
512 :
513 93191637 : case 1:
514 93191637 : return 2.*zeta1*(zeta1-0.5);
515 :
516 93191637 : case 2:
517 93191637 : return 2.*zeta2*(zeta2-0.5);
518 :
519 93191637 : case 3:
520 93191637 : return 4.*zeta0*zeta1;
521 :
522 93191637 : case 4:
523 93191637 : return 4.*zeta1*zeta2;
524 :
525 93191637 : case 5:
526 93191637 : 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 992112249 : case THIRD:
542 : {
543 992112249 : switch (type)
544 : {
545 992112249 : case TRI7:
546 : {
547 992112249 : const Real zeta1 = p(0);
548 992112249 : const Real zeta2 = p(1);
549 992112249 : const Real zeta0 = 1. - zeta1 - zeta2;
550 992112249 : const Real bubble_27th = zeta0*zeta1*zeta2;
551 :
552 84155654 : libmesh_assert_less (i, 7);
553 :
554 992112249 : switch(i)
555 : {
556 135393747 : case 0:
557 135393747 : return 2.*zeta0*(zeta0-0.5) + 3.*bubble_27th;
558 :
559 135393747 : case 1:
560 135393747 : return 2.*zeta1*(zeta1-0.5) + 3.*bubble_27th;
561 :
562 135393747 : case 2:
563 135393747 : return 2.*zeta2*(zeta2-0.5) + 3.*bubble_27th;
564 :
565 135393747 : case 3:
566 135393747 : return 4.*zeta0*zeta1 - 12.*bubble_27th;
567 :
568 135393747 : case 4:
569 135393747 : return 4.*zeta1*zeta2 - 12.*bubble_27th;
570 :
571 135393747 : case 5:
572 135393747 : return 4.*zeta2*zeta0 - 12.*bubble_27th;
573 :
574 179749767 : case 6:
575 179749767 : 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 16868094996 : 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 1265045440 : libmesh_assert_less (j, 2);
612 :
613 16868094996 : switch (order)
614 : {
615 : // linear Lagrange shape functions
616 461591924 : case FIRST:
617 : {
618 411842878 : switch (type)
619 : {
620 326901376 : 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 326901376 : const Real xi = p(0);
629 326901376 : const Real eta = p(1);
630 :
631 38651416 : 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 326901376 : switch (j)
638 : {
639 : // d()/dxi
640 163450688 : case 0:
641 163450688 : return (fe_lagrange_1D_linear_shape_deriv(i0[i], 0, xi)*
642 163450688 : fe_lagrange_1D_linear_shape (i1[i], eta));
643 :
644 : // d()/deta
645 163450688 : case 1:
646 163450688 : return (fe_lagrange_1D_linear_shape (i0[i], xi)*
647 163450688 : 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 133391748 : case TRI3:
655 : case TRISHELL3:
656 : case TRI6:
657 : case TRI7:
658 : {
659 11264616 : libmesh_assert_less (i, 3);
660 :
661 11264616 : const Real dzeta0dxi = -1.;
662 11264616 : const Real dzeta1dxi = 1.;
663 11264616 : const Real dzeta2dxi = 0.;
664 :
665 11264616 : const Real dzeta0deta = -1.;
666 11264616 : const Real dzeta1deta = 0.;
667 11264616 : const Real dzeta2deta = 1.;
668 :
669 133391748 : switch (j)
670 : {
671 : // d()/dxi
672 66695874 : case 0:
673 : {
674 61063638 : switch(i)
675 : {
676 1877436 : case 0:
677 1877436 : return dzeta0dxi;
678 :
679 1877436 : case 1:
680 1877436 : return dzeta1dxi;
681 :
682 1877436 : case 2:
683 1877436 : return dzeta2dxi;
684 :
685 0 : default:
686 0 : libmesh_error_msg("Invalid shape function index i = " << i);
687 : }
688 : }
689 : // d()/deta
690 66695874 : case 1:
691 : {
692 61063638 : switch(i)
693 : {
694 1877436 : case 0:
695 1877436 : return dzeta0deta;
696 :
697 1877436 : case 1:
698 1877436 : return dzeta1deta;
699 :
700 1877436 : case 2:
701 1877436 : 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 1298800 : case C0POLYGON:
713 : {
714 : // C0Polygon requires using newer FE APIs
715 1298800 : 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 1298800 : const auto [s, a, b] = poly.subtriangle_coordinates(p, 100);
727 1298800 : if (s == invalid_uint)
728 0 : return 0;
729 106510 : libmesh_assert_less(s, poly.n_subtriangles());
730 :
731 1298800 : 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 1298800 : const int nodei = i;
739 1298800 : if (nodei == subtri[0])
740 21302 : du_da = du_db = -1;
741 1039040 : else if (nodei == subtri[1])
742 21302 : du_da = 1;
743 779280 : 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 779280 : const auto master_points = poly.master_subtriangle(s);
756 :
757 779280 : const Real dxi_da = master_points[1](0) - master_points[0](0);
758 779280 : const Real dxi_db = master_points[2](0) - master_points[0](0);
759 779280 : const Real deta_da = master_points[1](1) - master_points[0](1);
760 779280 : const Real deta_db = master_points[2](1) - master_points[0](1);
761 779280 : const Real jac = dxi_da*deta_db - dxi_db*deta_da;
762 :
763 779280 : switch (j)
764 : {
765 : // d()/dxi
766 389640 : case 0:
767 : {
768 389640 : const Real da_dxi = deta_db / jac;
769 389640 : const Real db_dxi = -deta_da / jac;
770 389640 : return du_da*da_dxi + du_db*db_dxi;
771 : }
772 : // d()/deta
773 389640 : case 1:
774 : {
775 389640 : const Real da_deta = -dxi_db / jac;
776 389640 : const Real db_deta = dxi_da / jac;
777 389640 : 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 14898884448 : case SECOND:
792 : {
793 13809998144 : switch (type)
794 : {
795 59970336 : case QUAD8:
796 : case QUADSHELL8:
797 : {
798 59970336 : const Real xi = p(0);
799 59970336 : const Real eta = p(1);
800 :
801 4863104 : libmesh_assert_less (i, 8);
802 :
803 59970336 : switch (j)
804 : {
805 : // d/dxi
806 29985168 : case 0:
807 29985168 : switch (i)
808 : {
809 3748146 : case 0:
810 4356030 : return .25*(1. - eta)*((1. - xi)*(-1.) +
811 3748146 : (-1.)*(-1. - xi - eta));
812 :
813 3748146 : case 1:
814 4356030 : return .25*(1. - eta)*((1. + xi)*(1.) +
815 3748146 : (1.)*(-1. + xi - eta));
816 :
817 3748146 : case 2:
818 4356030 : return .25*(1. + eta)*((1. + xi)*(1.) +
819 3748146 : (1.)*(-1. + xi + eta));
820 :
821 3748146 : case 3:
822 4356030 : return .25*(1. + eta)*((1. - xi)*(-1.) +
823 3748146 : (-1.)*(-1. - xi + eta));
824 :
825 3748146 : case 4:
826 3748146 : return .5*(-2.*xi)*(1. - eta);
827 :
828 3748146 : case 5:
829 3748146 : return .5*(1.)*(1. - eta*eta);
830 :
831 3748146 : case 6:
832 3748146 : return .5*(-2.*xi)*(1. + eta);
833 :
834 3748146 : case 7:
835 3748146 : 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 29985168 : case 1:
843 29985168 : switch (i)
844 : {
845 3748146 : case 0:
846 4356030 : return .25*(1. - xi)*((1. - eta)*(-1.) +
847 3748146 : (-1.)*(-1. - xi - eta));
848 :
849 3748146 : case 1:
850 4356030 : return .25*(1. + xi)*((1. - eta)*(-1.) +
851 3748146 : (-1.)*(-1. + xi - eta));
852 :
853 3748146 : case 2:
854 4356030 : return .25*(1. + xi)*((1. + eta)*(1.) +
855 3748146 : (1.)*(-1. + xi + eta));
856 :
857 3748146 : case 3:
858 4356030 : return .25*(1. - xi)*((1. + eta)*(1.) +
859 3748146 : (1.)*(-1. - xi + eta));
860 :
861 3748146 : case 4:
862 3748146 : return .5*(1. - xi*xi)*(-1.);
863 :
864 3748146 : case 5:
865 3748146 : return .5*(1. + xi)*(-2.*eta);
866 :
867 3748146 : case 6:
868 3748146 : return .5*(1. - xi*xi)*(1.);
869 :
870 3748146 : case 7:
871 3748146 : 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 1850488434 : 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 129656016 : case QUAD9:
887 : case QUADSHELL9:
888 : {
889 : // Compute quad shape functions as a tensor-product
890 1980152316 : const Real xi = p(0);
891 1980152316 : const Real eta = p(1);
892 :
893 129671748 : 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 1980152316 : switch (j)
900 : {
901 : // d()/dxi
902 990076158 : case 0:
903 990076158 : return (fe_lagrange_1D_quadratic_shape_deriv(i0[i], 0, xi)*
904 990076158 : fe_lagrange_1D_quadratic_shape (i1[i], eta));
905 :
906 : // d()/deta
907 990076158 : case 1:
908 990076158 : return (fe_lagrange_1D_quadratic_shape (i0[i], xi)*
909 990076158 : 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 11902635432 : 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 956027844 : case TRI6:
921 : case TRI7:
922 : {
923 956131848 : libmesh_assert_less (i, 6);
924 :
925 12858761796 : const Real zeta1 = p(0);
926 12858761796 : const Real zeta2 = p(1);
927 12858761796 : const Real zeta0 = 1. - zeta1 - zeta2;
928 :
929 956131848 : const Real dzeta0dxi = -1.;
930 956131848 : const Real dzeta1dxi = 1.;
931 956131848 : const Real dzeta2dxi = 0.;
932 :
933 956131848 : const Real dzeta0deta = -1.;
934 956131848 : const Real dzeta1deta = 0.;
935 956131848 : const Real dzeta2deta = 1.;
936 :
937 12858761796 : switch(j)
938 : {
939 6429380898 : case 0:
940 : {
941 6429380898 : switch(i)
942 : {
943 1071563483 : case 0:
944 1071563483 : return (4.*zeta0-1.)*dzeta0dxi;
945 :
946 1071563483 : case 1:
947 1071563483 : return (4.*zeta1-1.)*dzeta1dxi;
948 :
949 1071563483 : case 2:
950 1071563483 : return (4.*zeta2-1.)*dzeta2dxi;
951 :
952 1071563483 : case 3:
953 1071563483 : return 4.*zeta1*dzeta0dxi + 4.*zeta0*dzeta1dxi;
954 :
955 1071563483 : case 4:
956 1071563483 : return 4.*zeta2*dzeta1dxi + 4.*zeta1*dzeta2dxi;
957 :
958 1071563483 : case 5:
959 1071563483 : 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 6429380898 : case 1:
967 : {
968 6429380898 : switch(i)
969 : {
970 1071563483 : case 0:
971 1071563483 : return (4.*zeta0-1.)*dzeta0deta;
972 :
973 1071563483 : case 1:
974 1071563483 : return (4.*zeta1-1.)*dzeta1deta;
975 :
976 1071563483 : case 2:
977 1071563483 : return (4.*zeta2-1.)*dzeta2deta;
978 :
979 1071563483 : case 3:
980 1071563483 : return 4.*zeta1*dzeta0deta + 4.*zeta0*dzeta1deta;
981 :
982 1071563483 : case 4:
983 1071563483 : return 4.*zeta2*dzeta1deta + 4.*zeta1*dzeta2deta;
984 :
985 1071563483 : case 5:
986 1071563483 : 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 1507618624 : case THIRD:
1006 : {
1007 1507618624 : switch (type)
1008 : {
1009 1507618624 : case TRI7:
1010 : {
1011 124356198 : libmesh_assert_less (i, 7);
1012 :
1013 1507618624 : const Real zeta1 = p(0);
1014 1507618624 : const Real zeta2 = p(1);
1015 1507618624 : const Real zeta0 = 1. - zeta1 - zeta2;
1016 : // const Real bubble_27th = zeta0*zeta1*zeta2;
1017 :
1018 124356198 : const Real dzeta0dxi = -1.;
1019 124356198 : const Real dzeta1dxi = 1.;
1020 124356198 : const Real dzeta2dxi = 0.;
1021 1507618624 : const Real dbubbledxi = zeta2 * (1. - 2.*zeta1 - zeta2);
1022 :
1023 124356198 : const Real dzeta0deta = -1.;
1024 124356198 : const Real dzeta1deta = 0.;
1025 124356198 : const Real dzeta2deta = 1.;
1026 1507618624 : const Real dbubbledeta= zeta1 * (1. - zeta1 - 2.*zeta2);
1027 :
1028 1507618624 : switch(j)
1029 : {
1030 753809312 : case 0:
1031 : {
1032 753809312 : switch(i)
1033 : {
1034 104126846 : case 0:
1035 104126846 : return (4.*zeta0-1.)*dzeta0dxi + 3.*dbubbledxi;
1036 :
1037 104126846 : case 1:
1038 104126846 : return (4.*zeta1-1.)*dzeta1dxi + 3.*dbubbledxi;
1039 :
1040 104126846 : case 2:
1041 104126846 : return (4.*zeta2-1.)*dzeta2dxi + 3.*dbubbledxi;
1042 :
1043 104126846 : case 3:
1044 104126846 : return 4.*zeta1*dzeta0dxi + 4.*zeta0*dzeta1dxi - 12.*dbubbledxi;
1045 :
1046 104126846 : case 4:
1047 104126846 : return 4.*zeta2*dzeta1dxi + 4.*zeta1*dzeta2dxi - 12.*dbubbledxi;
1048 :
1049 104126846 : case 5:
1050 104126846 : return 4.*zeta2*dzeta0dxi + 4*zeta0*dzeta2dxi - 12.*dbubbledxi;
1051 :
1052 129048236 : case 6:
1053 129048236 : return 27.*dbubbledxi;
1054 :
1055 0 : default:
1056 0 : libmesh_error_msg("Invalid shape function index i = " << i);
1057 : }
1058 : }
1059 :
1060 753809312 : case 1:
1061 : {
1062 753809312 : switch(i)
1063 : {
1064 104126846 : case 0:
1065 104126846 : return (4.*zeta0-1.)*dzeta0deta + 3.*dbubbledeta;
1066 :
1067 104126846 : case 1:
1068 104126846 : return (4.*zeta1-1.)*dzeta1deta + 3.*dbubbledeta;
1069 :
1070 104126846 : case 2:
1071 104126846 : return (4.*zeta2-1.)*dzeta2deta + 3.*dbubbledeta;
1072 :
1073 104126846 : case 3:
1074 104126846 : return 4.*zeta1*dzeta0deta + 4.*zeta0*dzeta1deta - 12.*dbubbledeta;
1075 :
1076 104126846 : case 4:
1077 104126846 : return 4.*zeta2*dzeta1deta + 4.*zeta1*dzeta2deta - 12.*dbubbledeta;
1078 :
1079 104126846 : case 5:
1080 104126846 : return 4.*zeta2*dzeta0deta + 4*zeta0*dzeta2deta - 12.*dbubbledeta;
1081 :
1082 129048236 : case 6:
1083 129048236 : 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 282140523 : 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 25579410 : libmesh_assert_less (j, 3);
1129 :
1130 282140523 : switch (order)
1131 : {
1132 : // linear Lagrange shape functions
1133 2833725 : case FIRST:
1134 : {
1135 2833725 : switch (type)
1136 : {
1137 1739052 : 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 1739052 : const Real xi = p(0);
1146 1739052 : 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 1739052 : switch (j)
1155 : {
1156 : // d^2() / dxi^2
1157 52400 : case 0:
1158 52400 : return 0.;
1159 :
1160 : // d^2() / dxi deta
1161 579684 : case 1:
1162 579684 : return (fe_lagrange_1D_linear_shape_deriv(i0[i], 0, xi)*
1163 579684 : 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 127588464 : case SECOND:
1197 : {
1198 127588464 : switch (type)
1199 : {
1200 4866912 : case QUAD8:
1201 : case QUADSHELL8:
1202 : {
1203 4866912 : const Real xi = p(0);
1204 4866912 : const Real eta = p(1);
1205 :
1206 439488 : libmesh_assert_less (j, 3);
1207 :
1208 4866912 : switch (j)
1209 : {
1210 : // d^2() / dxi^2
1211 1622304 : case 0:
1212 : {
1213 1622304 : switch (i)
1214 : {
1215 405576 : case 0:
1216 : case 1:
1217 405576 : return 0.5*(1.-eta);
1218 :
1219 405576 : case 2:
1220 : case 3:
1221 405576 : return 0.5*(1.+eta);
1222 :
1223 202788 : case 4:
1224 202788 : return eta - 1.;
1225 :
1226 36624 : case 5:
1227 : case 7:
1228 36624 : return 0.0;
1229 :
1230 202788 : case 6:
1231 202788 : 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 1622304 : case 1:
1240 : {
1241 1622304 : switch (i)
1242 : {
1243 202788 : case 0:
1244 202788 : return 0.25*( 1. - 2.*xi - 2.*eta);
1245 :
1246 202788 : case 1:
1247 202788 : return 0.25*(-1. - 2.*xi + 2.*eta);
1248 :
1249 202788 : case 2:
1250 202788 : return 0.25*( 1. + 2.*xi + 2.*eta);
1251 :
1252 202788 : case 3:
1253 202788 : return 0.25*(-1. + 2.*xi - 2.*eta);
1254 :
1255 18312 : case 4:
1256 18312 : return xi;
1257 :
1258 202788 : case 5:
1259 202788 : return -eta;
1260 :
1261 202788 : case 6:
1262 202788 : return -xi;
1263 :
1264 202788 : case 7:
1265 202788 : 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 1622304 : case 2:
1274 : {
1275 1622304 : switch (i)
1276 : {
1277 405576 : case 0:
1278 : case 3:
1279 405576 : return 0.5*(1.-xi);
1280 :
1281 405576 : case 1:
1282 : case 2:
1283 405576 : return 0.5*(1.+xi);
1284 :
1285 36624 : case 4:
1286 : case 6:
1287 36624 : return 0.0;
1288 :
1289 202788 : case 5:
1290 202788 : return -1.0 - xi;
1291 :
1292 202788 : case 7:
1293 202788 : 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 27317979 : 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 2428056 : case QUAD9:
1310 : case QUADSHELL9:
1311 : {
1312 : // Compute QUAD9 second derivatives as tensor product
1313 29757834 : const Real xi = p(0);
1314 29757834 : const Real eta = p(1);
1315 :
1316 2451654 : 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 29757834 : switch (j)
1323 : {
1324 : // d^2() / dxi^2
1325 9919278 : case 0:
1326 9919278 : return (fe_lagrange_1D_quadratic_shape_second_deriv(i0[i], 0, xi)*
1327 9919278 : fe_lagrange_1D_quadratic_shape (i1[i], eta));
1328 :
1329 : // d^2() / dxi deta
1330 9919278 : case 1:
1331 9919278 : return (fe_lagrange_1D_quadratic_shape_deriv(i0[i], 0, xi)*
1332 9919278 : fe_lagrange_1D_quadratic_shape_deriv(i1[i], 0, eta));
1333 :
1334 : // d^2() / deta^2
1335 9919278 : case 2:
1336 9919278 : return (fe_lagrange_1D_quadratic_shape (i0[i], xi)*
1337 9919278 : 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 84743676 : 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 92963718 : switch (j)
1362 : {
1363 : // d^2() / dxi^2
1364 30987906 : case 0:
1365 : {
1366 30987906 : 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 30987906 : case 1:
1393 : {
1394 30987906 : 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 30987906 : case 2:
1421 : {
1422 30987906 : 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 151718334 : case THIRD:
1461 : {
1462 151718334 : switch (type)
1463 : {
1464 137517447 : 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 151718334 : const Real zeta1 = p(0);
1472 151718334 : 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 151718334 : 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 151718334 : const Real d2bubbledeta2 = -2. * zeta1;
1486 :
1487 151718334 : const Real d2bubbledxideta = (1. - 2.*zeta1 - 2.*zeta2);
1488 :
1489 14200887 : libmesh_assert_less (j, 3);
1490 :
1491 151718334 : switch (j)
1492 : {
1493 : // d^2() / dxi^2
1494 50572778 : case 0:
1495 : {
1496 50572778 : switch (i)
1497 : {
1498 6754954 : case 0:
1499 6754954 : return 4.*dzeta0dxi*dzeta0dxi + 3.*d2bubbledxi2;
1500 :
1501 6754954 : case 1:
1502 6754954 : return 4.*dzeta1dxi*dzeta1dxi + 3.*d2bubbledxi2;
1503 :
1504 6754954 : case 2:
1505 6754954 : return 4.*dzeta2dxi*dzeta2dxi + 3.*d2bubbledxi2;
1506 :
1507 6754954 : case 3:
1508 6754954 : return 8.*dzeta0dxi*dzeta1dxi - 12.*d2bubbledxi2;
1509 :
1510 6754954 : case 4:
1511 6754954 : return 8.*dzeta1dxi*dzeta2dxi - 12.*d2bubbledxi2;
1512 :
1513 6754954 : case 5:
1514 6754954 : return 8.*dzeta0dxi*dzeta2dxi - 12.*d2bubbledxi2;
1515 :
1516 10043054 : case 6:
1517 10043054 : 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 50572778 : case 1:
1526 : {
1527 50572778 : switch (i)
1528 : {
1529 6754954 : case 0:
1530 6754954 : return 4.*dzeta0dxi*dzeta0deta + 3.*d2bubbledxideta;
1531 :
1532 6754954 : case 1:
1533 6754954 : return 4.*dzeta1dxi*dzeta1deta + 3.*d2bubbledxideta;
1534 :
1535 6754954 : case 2:
1536 6754954 : return 4.*dzeta2dxi*dzeta2deta + 3.*d2bubbledxideta;
1537 :
1538 6754954 : case 3:
1539 6754954 : return 4.*dzeta1deta*dzeta0dxi + 4.*dzeta0deta*dzeta1dxi - 12.*d2bubbledxideta;
1540 :
1541 6754954 : case 4:
1542 6754954 : return 4.*dzeta2deta*dzeta1dxi + 4.*dzeta1deta*dzeta2dxi - 12.*d2bubbledxideta;
1543 :
1544 6754954 : case 5:
1545 6754954 : return 4.*dzeta2deta*dzeta0dxi + 4.*dzeta0deta*dzeta2dxi - 12.*d2bubbledxideta;
1546 :
1547 10043054 : case 6:
1548 10043054 : 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 50572778 : case 2:
1557 : {
1558 50572778 : switch (i)
1559 : {
1560 6754954 : case 0:
1561 6754954 : return 4.*dzeta0deta*dzeta0deta + 3.*d2bubbledeta2;
1562 :
1563 6754954 : case 1:
1564 6754954 : return 4.*dzeta1deta*dzeta1deta + 3.*d2bubbledeta2;
1565 :
1566 6754954 : case 2:
1567 6754954 : return 4.*dzeta2deta*dzeta2deta + 3.*d2bubbledeta2;
1568 :
1569 6754954 : case 3:
1570 6754954 : return 8.*dzeta0deta*dzeta1deta - 12.*d2bubbledeta2;
1571 :
1572 6754954 : case 4:
1573 6754954 : return 8.*dzeta1deta*dzeta2deta - 12.*d2bubbledeta2;
1574 :
1575 6754954 : case 5:
1576 6754954 : return 8.*dzeta0deta*dzeta2deta - 12.*d2bubbledeta2;
1577 :
1578 10043054 : case 6:
1579 10043054 : 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
|
