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