www.mooseframework.org
MooseLagrangeHelpers.h
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://www.mooseframework.org
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 
15 namespace Moose
16 {
17 using libMesh::Order;
18 
19 // Copy in libmesh's lagrange helper functions, but we template it
20 template <typename T>
21 T
22 fe_lagrange_1D_shape(const Order order, const unsigned int i, const T & xi)
23 {
24  switch (order)
25  {
26  // Lagrange linears
27  case 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 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 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 
94 template <typename T>
95 T
96 fe_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 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 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 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
169 template <typename T, template <typename> class VectorType>
170 T
171 fe_lagrange_2D_shape(const 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 FIRST:
180  {
181  switch (type)
182  {
183  case QUAD4:
184  case QUADSHELL4:
185  case QUAD8:
186  case QUADSHELL8:
187  case 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 TRI3:
203  case TRISHELL3:
204  case TRI6:
205  case 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 SECOND:
236  {
237  switch (type)
238  {
239  case 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 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(SECOND, i0[i], xi) *
282  fe_lagrange_1D_shape(SECOND, i1[i], eta));
283  }
284  case TRI6:
285  case 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 THIRD:
325  {
326  switch (type)
327  {
328  case 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 
376 template <typename T, template <typename> class VectorType>
377 T
378 fe_lagrange_2D_shape_deriv(const 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 FIRST:
390  {
391  switch (type)
392  {
393  case QUAD4:
394  case QUADSHELL4:
395  case QUAD8:
396  case QUADSHELL8:
397  case 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 TRI3:
427  case TRISHELL3:
428  case TRI6:
429  case 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 SECOND:
491  {
492  switch (type)
493  {
494  case QUAD8:
495  case 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 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(SECOND, i0[i], xi) *
590  fe_lagrange_1D_shape(SECOND, i1[i], eta));
591 
592  // d()/deta
593  case 1:
594  return (fe_lagrange_1D_shape(SECOND, i0[i], xi) *
595  fe_lagrange_1D_shape_deriv(SECOND, i1[i], eta));
596 
597  default:
598  mooseError("Invalid derivative index j = ", j);
599  }
600  }
601 
602  case TRI6:
603  case 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 THIRD:
686  {
687  switch (type)
688  {
689  case 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 }
ElemType
ElemType
QUAD8
QUAD8
mooseError
void mooseError(Args &&... args)
Emit an error message with the given stringified, concatenated args and terminate the application...
Definition: MooseError.h:299
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
FIRST
FIRST
SECOND
SECOND
TRI3
TRI3
QUAD4
QUAD4
Moose::fe_lagrange_2D_shape_deriv
T fe_lagrange_2D_shape_deriv(const ElemType type, const Order order, const unsigned int i, const unsigned int j, const VectorType< T > &p)
Definition: MooseLagrangeHelpers.h:378
TRI6
TRI6
QUADSHELL4
QUADSHELL4
QUADSHELL8
QUADSHELL8
TRISHELL3
TRISHELL3
TRI7
TRI7
Moose::fe_lagrange_1D_shape
T fe_lagrange_1D_shape(const Order order, const unsigned int i, const T &xi)
Definition: MooseLagrangeHelpers.h:22
QUAD9
QUAD9
Moose
MOOSE now contains C++17 code, so give a reasonable error message stating what the user can do to add...
Definition: SetAdaptivityOptionsAction.h:14
Moose::fe_lagrange_2D_shape
T fe_lagrange_2D_shape(const ElemType type, const Order order, const unsigned int i, const VectorType< T > &p)
Definition: MooseLagrangeHelpers.h:171
THIRD
THIRD
Generated on Wed Apr 17 2024 11:09:24 for www.mooseframework.org by   doxygen 1.8.14