https://mooseframework.inl.gov
Loading...
Searching...
No Matches
MooseLagrangeHelpers.h
Go to the documentation of this file.
1//* This file is part of the MOOSE framework
2//* https://mooseframework.inl.gov
3//*
4//* All rights reserved, see COPYRIGHT for full restrictions
5//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6//*
7//* Licensed under LGPL 2.1, please see LICENSE for details
8//* https://www.gnu.org/licenses/lgpl-2.1.html
9
10#pragma once
11
12#include "MooseError.h"
13#include "libmesh/fe_type.h"
14
15namespace Moose
16{
17// Copy in libmesh's lagrange helper functions, but we template it
18template <typename T>
19T
20fe_lagrange_1D_shape(const Order order, const unsigned int i, const T & xi)
21{
22 switch (order)
23 {
24 // Lagrange linears
25 case libMesh::FIRST:
26 {
27 libmesh_assert_less(i, 2);
28
29 switch (i)
30 {
31 case 0:
32 return .5 * (1. - xi);
33
34 case 1:
35 return .5 * (1. + xi);
36
37 default:
38 mooseError("Invalid shape function index i = ", i);
39 }
40 }
41
42 // Lagrange quadratics
43 case libMesh::SECOND:
44 {
45 libmesh_assert_less(i, 3);
46
47 switch (i)
48 {
49 case 0:
50 return .5 * xi * (xi - 1.);
51
52 case 1:
53 return .5 * xi * (xi + 1);
54
55 case 2:
56 return (1. - xi * xi);
57
58 default:
59 mooseError("Invalid shape function index i = ", i);
60 }
61 }
62
63 // Lagrange cubics
64 case libMesh::THIRD:
65 {
66 libmesh_assert_less(i, 4);
67
68 switch (i)
69 {
70 case 0:
71 return 9. / 16. * (1. / 9. - xi * xi) * (xi - 1.);
72
73 case 1:
74 return -9. / 16. * (1. / 9. - xi * xi) * (xi + 1.);
75
76 case 2:
77 return 27. / 16. * (1. - xi * xi) * (1. / 3. - xi);
78
79 case 3:
80 return 27. / 16. * (1. - xi * xi) * (1. / 3. + xi);
81
82 default:
83 mooseError("Invalid shape function index i = ", i);
84 }
85 }
86
87 default:
88 mooseError("Unsupported order");
89 }
90}
91
92template <typename T>
93T
94fe_lagrange_1D_shape_deriv(const Order order, const unsigned int i, const T & xi)
95{
96 switch (order)
97 {
98 // Lagrange linear shape function derivatives
99 case libMesh::FIRST:
100 {
101 libmesh_assert_less(i, 2);
102
103 switch (i)
104 {
105 case 0:
106 return -.5;
107
108 case 1:
109 return .5;
110
111 default:
112 mooseError("Invalid shape function index i = ", i);
113 }
114 }
115
116 // Lagrange quadratic shape function derivatives
117 case libMesh::SECOND:
118 {
119 libmesh_assert_less(i, 3);
120
121 switch (i)
122 {
123 case 0:
124 return xi - .5;
125
126 case 1:
127 return xi + .5;
128
129 case 2:
130 return -2. * xi;
131
132 default:
133 mooseError("Invalid shape function index i = ", i);
134 }
135 }
136
137 // Lagrange cubic shape function derivatives
138 case libMesh::THIRD:
139 {
140 libmesh_assert_less(i, 4);
141
142 switch (i)
143 {
144 case 0:
145 return -9. / 16. * (3. * xi * xi - 2. * xi - 1. / 9.);
146
147 case 1:
148 return -9. / 16. * (-3. * xi * xi - 2. * xi + 1. / 9.);
149
150 case 2:
151 return 27. / 16. * (3. * xi * xi - 2. / 3. * xi - 1.);
152
153 case 3:
154 return 27. / 16. * (-3. * xi * xi - 2. / 3. * xi + 1.);
155
156 default:
157 mooseError("Invalid shape function index i = ", i);
158 }
159 }
160
161 default:
162 mooseError("Unsupported order");
163 }
164}
165
166// Copy of libMesh function but templated to enable calling with DualNumber vectors
167template <typename T, template <typename> class VectorType>
168T
169fe_lagrange_2D_shape(const libMesh::ElemType type,
170 const Order order,
171 const unsigned int i,
172 const VectorType<T> & p)
173{
174 switch (order)
175 {
176 // linear Lagrange shape functions
177 case libMesh::FIRST:
178 {
179 switch (type)
180 {
181 case libMesh::QUAD4:
182 case libMesh::QUADSHELL4:
183 case libMesh::QUAD8:
184 case libMesh::QUADSHELL8:
185 case libMesh::QUAD9:
186 {
187 // Compute quad shape functions as a tensor-product
188 const T xi = p(0);
189 const T eta = p(1);
190
191 libmesh_assert_less(i, 4);
192
193 // 0 1 2 3
194 static const unsigned int i0[] = {0, 1, 1, 0};
195 static const unsigned int i1[] = {0, 0, 1, 1};
196
197 return (fe_lagrange_1D_shape(FIRST, i0[i], xi) * fe_lagrange_1D_shape(FIRST, i1[i], eta));
198 }
199
200 case libMesh::TRI3:
201 case libMesh::TRISHELL3:
202 case libMesh::TRI6:
203 case libMesh::TRI7:
204 {
205 const T zeta1 = p(0);
206 const T zeta2 = p(1);
207 const T zeta0 = 1. - zeta1 - zeta2;
208
209 libmesh_assert_less(i, 3);
210
211 switch (i)
212 {
213 case 0:
214 return zeta0;
215
216 case 1:
217 return zeta1;
218
219 case 2:
220 return zeta2;
221
222 default:
223 mooseError("Invalid shape function index i = ", i);
224 }
225 }
226
227 default:
228 mooseError("Unsupported element type:", type);
229 }
230 }
231
232 // quadratic Lagrange shape functions
233 case libMesh::SECOND:
234 {
235 switch (type)
236 {
237 case libMesh::QUAD8:
238 {
239 // Compute quad shape functions as a tensor-product
240 const T xi = p(0);
241 const T eta = p(1);
242
243 libmesh_assert_less(i, 8);
244
245 switch (i)
246 {
247 case 0:
248 return .25 * (1. - xi) * (1. - eta) * (-1. - xi - eta);
249 case 1:
250 return .25 * (1. + xi) * (1. - eta) * (-1. + xi - eta);
251 case 2:
252 return .25 * (1. + xi) * (eta + 1.) * (-1. + xi + eta);
253 case 3:
254 return .25 * (1. - xi) * (eta + 1.) * (-1. - xi + eta);
255 case 4:
256 return .5 * (1. - xi * xi) * (1. - eta);
257 case 5:
258 return .5 * (1. + xi) * (1. - eta * eta);
259 case 6:
260 return .5 * (1. - xi * xi) * (1. + eta);
261 case 7:
262 return .5 * (1. - xi) * (1. - eta * eta);
263 default:
264 mooseError("Invalid shape function index i = ", i);
265 }
266 }
267 case libMesh::QUAD9:
268 {
269 // Compute quad shape functions as a tensor-product
270 const T xi = p(0);
271 const T eta = p(1);
272
273 libmesh_assert_less(i, 9);
274
275 // 0 1 2 3 4 5 6 7 8
276 static const unsigned int i0[] = {0, 1, 1, 0, 2, 1, 2, 0, 2};
277 static const unsigned int i1[] = {0, 0, 1, 1, 0, 2, 1, 2, 2};
278
279 return (fe_lagrange_1D_shape(libMesh::SECOND, i0[i], xi) *
280 fe_lagrange_1D_shape(libMesh::SECOND, i1[i], eta));
281 }
282 case libMesh::TRI6:
283 case libMesh::TRI7:
284 {
285 const T zeta1 = p(0);
286 const T zeta2 = p(1);
287 const T zeta0 = 1. - zeta1 - zeta2;
288
289 libmesh_assert_less(i, 6);
290
291 switch (i)
292 {
293 case 0:
294 return 2. * zeta0 * (zeta0 - 0.5);
295
296 case 1:
297 return 2. * zeta1 * (zeta1 - 0.5);
298
299 case 2:
300 return 2. * zeta2 * (zeta2 - 0.5);
301
302 case 3:
303 return 4. * zeta0 * zeta1;
304
305 case 4:
306 return 4. * zeta1 * zeta2;
307
308 case 5:
309 return 4. * zeta2 * zeta0;
310
311 default:
312 mooseError("Invalid shape function index i = ", i);
313 }
314 }
315
316 default:
317 mooseError("Unsupported 2D element type");
318 }
319 }
320
321 // "cubic" (one cubic bubble) Lagrange shape functions
322 case libMesh::THIRD:
323 {
324 switch (type)
325 {
326 case libMesh::TRI7:
327 {
328 const T zeta1 = p(0);
329 const T zeta2 = p(1);
330 const T zeta0 = 1. - zeta1 - zeta2;
331 const T bubble_27th = zeta0 * zeta1 * zeta2;
332
333 libmesh_assert_less(i, 7);
334
335 switch (i)
336 {
337 case 0:
338 return 2. * zeta0 * (zeta0 - 0.5) + 3. * bubble_27th;
339
340 case 1:
341 return 2. * zeta1 * (zeta1 - 0.5) + 3. * bubble_27th;
342
343 case 2:
344 return 2. * zeta2 * (zeta2 - 0.5) + 3. * bubble_27th;
345
346 case 3:
347 return 4. * zeta0 * zeta1 - 12. * bubble_27th;
348
349 case 4:
350 return 4. * zeta1 * zeta2 - 12. * bubble_27th;
351
352 case 5:
353 return 4. * zeta2 * zeta0 - 12. * bubble_27th;
354
355 case 6:
356 return 27. * bubble_27th;
357
358 default:
359 mooseError("Invalid shape function index i = ", i);
360 }
361 }
362
363 default:
364 mooseError("Unsupported 2D element type");
365 }
366 }
367
368 // unsupported order
369 default:
370 mooseError("Unsupported order");
371 }
372}
373
374template <typename T, template <typename> class VectorType>
375T
376fe_lagrange_2D_shape_deriv(const libMesh::ElemType type,
377 const Order order,
378 const unsigned int i,
379 const unsigned int j,
380 const VectorType<T> & p)
381{
382 libmesh_assert_less(j, 2);
383
384 switch (order)
385 {
386 // linear Lagrange shape functions
387 case libMesh::FIRST:
388 {
389 switch (type)
390 {
391 case libMesh::QUAD4:
392 case libMesh::QUADSHELL4:
393 case libMesh::QUAD8:
394 case libMesh::QUADSHELL8:
395 case libMesh::QUAD9:
396 {
397 // Compute quad shape functions as a tensor-product
398 const T xi = p(0);
399 const T eta = p(1);
400
401 libmesh_assert_less(i, 4);
402
403 // 0 1 2 3
404 static const unsigned int i0[] = {0, 1, 1, 0};
405 static const unsigned int i1[] = {0, 0, 1, 1};
406
407 switch (j)
408 {
409 // d()/dxi
410 case 0:
411 return (fe_lagrange_1D_shape_deriv(FIRST, i0[i], xi) *
412 fe_lagrange_1D_shape(FIRST, i1[i], eta));
413
414 // d()/deta
415 case 1:
416 return (fe_lagrange_1D_shape(FIRST, i0[i], xi) *
417 fe_lagrange_1D_shape_deriv(FIRST, i1[i], eta));
418
419 default:
420 mooseError("Invalid derivative index j = ", j);
421 }
422 }
423
424 case libMesh::TRI3:
425 case libMesh::TRISHELL3:
426 case libMesh::TRI6:
427 case libMesh::TRI7:
428 {
429 libmesh_assert_less(i, 3);
430
431 const T dzeta0dxi = -1.;
432 const T dzeta1dxi = 1.;
433 const T dzeta2dxi = 0.;
434
435 const T dzeta0deta = -1.;
436 const T dzeta1deta = 0.;
437 const T dzeta2deta = 1.;
438
439 switch (j)
440 {
441 // d()/dxi
442 case 0:
443 {
444 switch (i)
445 {
446 case 0:
447 return dzeta0dxi;
448
449 case 1:
450 return dzeta1dxi;
451
452 case 2:
453 return dzeta2dxi;
454
455 default:
456 mooseError("Invalid shape function index i = ", i);
457 }
458 }
459 // d()/deta
460 case 1:
461 {
462 switch (i)
463 {
464 case 0:
465 return dzeta0deta;
466
467 case 1:
468 return dzeta1deta;
469
470 case 2:
471 return dzeta2deta;
472
473 default:
474 mooseError("Invalid shape function index i = ", i);
475 }
476 }
477 default:
478 mooseError("Invalid derivative index j = ", j);
479 }
480 }
481
482 default:
483 mooseError("Unsupported 2D element type");
484 }
485 }
486
487 // quadratic Lagrange shape functions
488 case libMesh::SECOND:
489 {
490 switch (type)
491 {
492 case libMesh::QUAD8:
493 case libMesh::QUADSHELL8:
494 {
495 const T xi = p(0);
496 const T eta = p(1);
497
498 libmesh_assert_less(i, 8);
499
500 switch (j)
501 {
502 // d/dxi
503 case 0:
504 switch (i)
505 {
506 case 0:
507 return .25 * (1. - eta) * ((1. - xi) * (-1.) + (-1.) * (-1. - xi - eta));
508
509 case 1:
510 return .25 * (1. - eta) * ((1. + xi) * (1.) + (1.) * (-1. + xi - eta));
511
512 case 2:
513 return .25 * (1. + eta) * ((1. + xi) * (1.) + (1.) * (-1. + xi + eta));
514
515 case 3:
516 return .25 * (1. + eta) * ((1. - xi) * (-1.) + (-1.) * (-1. - xi + eta));
517
518 case 4:
519 return .5 * (-2. * xi) * (1. - eta);
520
521 case 5:
522 return .5 * (1.) * (1. - eta * eta);
523
524 case 6:
525 return .5 * (-2. * xi) * (1. + eta);
526
527 case 7:
528 return .5 * (-1.) * (1. - eta * eta);
529
530 default:
531 mooseError("Invalid shape function index i = ", i);
532 }
533
534 // d/deta
535 case 1:
536 switch (i)
537 {
538 case 0:
539 return .25 * (1. - xi) * ((1. - eta) * (-1.) + (-1.) * (-1. - xi - eta));
540
541 case 1:
542 return .25 * (1. + xi) * ((1. - eta) * (-1.) + (-1.) * (-1. + xi - eta));
543
544 case 2:
545 return .25 * (1. + xi) * ((1. + eta) * (1.) + (1.) * (-1. + xi + eta));
546
547 case 3:
548 return .25 * (1. - xi) * ((1. + eta) * (1.) + (1.) * (-1. - xi + eta));
549
550 case 4:
551 return .5 * (1. - xi * xi) * (-1.);
552
553 case 5:
554 return .5 * (1. + xi) * (-2. * eta);
555
556 case 6:
557 return .5 * (1. - xi * xi) * (1.);
558
559 case 7:
560 return .5 * (1. - xi) * (-2. * eta);
561
562 default:
563 mooseError("Invalid shape function index i = ", i);
564 }
565
566 default:
567 mooseError("ERROR: Invalid derivative index j = ", j);
568 }
569 }
570
571 case libMesh::QUAD9:
572 {
573 // Compute quad shape functions as a tensor-product
574 const T xi = p(0);
575 const T eta = p(1);
576
577 libmesh_assert_less(i, 9);
578
579 // 0 1 2 3 4 5 6 7 8
580 static const unsigned int i0[] = {0, 1, 1, 0, 2, 1, 2, 0, 2};
581 static const unsigned int i1[] = {0, 0, 1, 1, 0, 2, 1, 2, 2};
582
583 switch (j)
584 {
585 // d()/dxi
586 case 0:
587 return (fe_lagrange_1D_shape_deriv(libMesh::SECOND, i0[i], xi) *
588 fe_lagrange_1D_shape(libMesh::SECOND, i1[i], eta));
589
590 // d()/deta
591 case 1:
592 return (fe_lagrange_1D_shape(libMesh::SECOND, i0[i], xi) *
593 fe_lagrange_1D_shape_deriv(libMesh::SECOND, i1[i], eta));
594
595 default:
596 mooseError("Invalid derivative index j = ", j);
597 }
598 }
599
600 case libMesh::TRI6:
601 case libMesh::TRI7:
602 {
603 libmesh_assert_less(i, 6);
604
605 const T zeta1 = p(0);
606 const T zeta2 = p(1);
607 const T zeta0 = 1. - zeta1 - zeta2;
608
609 const T dzeta0dxi = -1.;
610 const T dzeta1dxi = 1.;
611 const T dzeta2dxi = 0.;
612
613 const T dzeta0deta = -1.;
614 const T dzeta1deta = 0.;
615 const T dzeta2deta = 1.;
616
617 switch (j)
618 {
619 case 0:
620 {
621 switch (i)
622 {
623 case 0:
624 return (4. * zeta0 - 1.) * dzeta0dxi;
625
626 case 1:
627 return (4. * zeta1 - 1.) * dzeta1dxi;
628
629 case 2:
630 return (4. * zeta2 - 1.) * dzeta2dxi;
631
632 case 3:
633 return 4. * zeta1 * dzeta0dxi + 4. * zeta0 * dzeta1dxi;
634
635 case 4:
636 return 4. * zeta2 * dzeta1dxi + 4. * zeta1 * dzeta2dxi;
637
638 case 5:
639 return 4. * zeta2 * dzeta0dxi + 4 * zeta0 * dzeta2dxi;
640
641 default:
642 mooseError("Invalid shape function index i = ", i);
643 }
644 }
645
646 case 1:
647 {
648 switch (i)
649 {
650 case 0:
651 return (4. * zeta0 - 1.) * dzeta0deta;
652
653 case 1:
654 return (4. * zeta1 - 1.) * dzeta1deta;
655
656 case 2:
657 return (4. * zeta2 - 1.) * dzeta2deta;
658
659 case 3:
660 return 4. * zeta1 * dzeta0deta + 4. * zeta0 * dzeta1deta;
661
662 case 4:
663 return 4. * zeta2 * dzeta1deta + 4. * zeta1 * dzeta2deta;
664
665 case 5:
666 return 4. * zeta2 * dzeta0deta + 4 * zeta0 * dzeta2deta;
667
668 default:
669 mooseError("Invalid shape function index i = ", i);
670 }
671 }
672 default:
673 mooseError("ERROR: Invalid derivative index j = ", j);
674 }
675 }
676
677 default:
678 mooseError("ERROR: Unsupported 2D element type");
679 }
680 }
681
682 // "cubic" (one cubic bubble) Lagrange shape functions
683 case libMesh::THIRD:
684 {
685 switch (type)
686 {
687 case libMesh::TRI7:
688 {
689 libmesh_assert_less(i, 7);
690
691 const T zeta1 = p(0);
692 const T zeta2 = p(1);
693 const T zeta0 = 1. - zeta1 - zeta2;
694
695 const T dzeta0dxi = -1.;
696 const T dzeta1dxi = 1.;
697 const T dzeta2dxi = 0.;
698 const T dbubbledxi = zeta2 * (1. - 2. * zeta1 - zeta2);
699
700 const T dzeta0deta = -1.;
701 const T dzeta1deta = 0.;
702 const T dzeta2deta = 1.;
703 const T dbubbledeta = zeta1 * (1. - zeta1 - 2. * zeta2);
704
705 switch (j)
706 {
707 case 0:
708 {
709 switch (i)
710 {
711 case 0:
712 return (4. * zeta0 - 1.) * dzeta0dxi + 3. * dbubbledxi;
713
714 case 1:
715 return (4. * zeta1 - 1.) * dzeta1dxi + 3. * dbubbledxi;
716
717 case 2:
718 return (4. * zeta2 - 1.) * dzeta2dxi + 3. * dbubbledxi;
719
720 case 3:
721 return 4. * zeta1 * dzeta0dxi + 4. * zeta0 * dzeta1dxi - 12. * dbubbledxi;
722
723 case 4:
724 return 4. * zeta2 * dzeta1dxi + 4. * zeta1 * dzeta2dxi - 12. * dbubbledxi;
725
726 case 5:
727 return 4. * zeta2 * dzeta0dxi + 4 * zeta0 * dzeta2dxi - 12. * dbubbledxi;
728
729 case 6:
730 return 27. * dbubbledxi;
731
732 default:
733 mooseError("Invalid shape function index i = ", i);
734 }
735 }
736
737 case 1:
738 {
739 switch (i)
740 {
741 case 0:
742 return (4. * zeta0 - 1.) * dzeta0deta + 3. * dbubbledeta;
743
744 case 1:
745 return (4. * zeta1 - 1.) * dzeta1deta + 3. * dbubbledeta;
746
747 case 2:
748 return (4. * zeta2 - 1.) * dzeta2deta + 3. * dbubbledeta;
749
750 case 3:
751 return 4. * zeta1 * dzeta0deta + 4. * zeta0 * dzeta1deta - 12. * dbubbledeta;
752
753 case 4:
754 return 4. * zeta2 * dzeta1deta + 4. * zeta1 * dzeta2deta - 12. * dbubbledeta;
755
756 case 5:
757 return 4. * zeta2 * dzeta0deta + 4 * zeta0 * dzeta2deta - 12. * dbubbledeta;
758
759 case 6:
760 return 27. * dbubbledeta;
761
762 default:
763 mooseError("Invalid shape function index i = ", i);
764 }
765 }
766 default:
767 mooseError("ERROR: Invalid derivative index j = ", j);
768 }
769 }
770
771 default:
772 mooseError("ERROR: Unsupported 2D element type");
773 }
774 }
775
776 // unsupported order
777 default:
778 mooseError("Unsupported order");
779 }
780}
781}
mooseError
void mooseError(Args &&... args)
Emit an error message with the given stringified, concatenated args and terminate the application.
Definition MooseError.h:311
eta
Point eta
Definition MortarUtils.C:60
xi
Point xi
Definition MortarUtils.C:59
Moose
MOOSE now contains C++17 code, so give a reasonable error message stating what the user can do to add...
Definition SetAdaptivityOptionsAction.h:15
Moose::fe_lagrange_2D_shape
T fe_lagrange_2D_shape(const libMesh::ElemType type, const Order order, const unsigned int i, const VectorType< T > &p)
Definition MooseLagrangeHelpers.h:169
Moose::fe_lagrange_1D_shape
T fe_lagrange_1D_shape(const Order order, const unsigned int i, const T &xi)
Definition MooseLagrangeHelpers.h:20
Moose::fe_lagrange_2D_shape_deriv
T fe_lagrange_2D_shape_deriv(const libMesh::ElemType type, const Order order, const unsigned int i, const unsigned int j, const VectorType< T > &p)
Definition MooseLagrangeHelpers.h:376
Moose::fe_lagrange_1D_shape_deriv
T fe_lagrange_1D_shape_deriv(const Order order, const unsigned int i, const T &xi)
Definition MooseLagrangeHelpers.h:94
Generated on Mon Sep 14 2026 20:29:10 for https://mooseframework.inl.gov by doxygen 1.9.8