libMesh
fe_lagrange_shape_2D.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2024 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 
25 // Anonymous namespace for functions shared by LAGRANGE and
26 // L2_LAGRANGE implementations. Implementations appear at the bottom
27 // of this file.
28 namespace
29 {
30 using namespace libMesh;
31 
32 template <FEFamily T>
33 Real fe_lagrange_2D_shape(const ElemType,
34  const Order order,
35  const unsigned int i,
36  const Point & p);
37 
38 template <FEFamily T>
39 Real fe_lagrange_2D_shape_deriv(const ElemType type,
40  const Order order,
41  const unsigned int i,
42  const unsigned int j,
43  const Point & p);
44 
45 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
46 
47 template <FEFamily T>
48 Real fe_lagrange_2D_shape_second_deriv(const ElemType type,
49  const Order order,
50  const unsigned int i,
51  const unsigned int j,
52  const Point & p);
53 
54 #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
55 
56 } // anonymous namespace
57 
58 
59 
60 namespace libMesh
61 {
62 
63 
66 
67 
68 template <>
70  const Order order,
71  const unsigned int i,
72  const Point & p)
73 {
74  return fe_lagrange_2D_shape<LAGRANGE>(type, order, i, p);
75 }
76 
77 
78 
79 template <>
81  const Order order,
82  const unsigned int i,
83  const Point & p)
84 {
85  return fe_lagrange_2D_shape<L2_LAGRANGE>(type, order, i, p);
86 }
87 
88 
89 template <>
91  const Order order,
92  const unsigned int i,
93  const Point & p,
94  const bool add_p_level)
95 {
96  libmesh_assert(elem);
97 
98  // call the orientation-independent shape functions
99  return fe_lagrange_2D_shape<LAGRANGE>(elem->type(), static_cast<Order>(order + add_p_level * elem->p_level()), i, p);
100 }
101 
102 
103 
104 template <>
106  const Order order,
107  const unsigned int i,
108  const Point & p,
109  const bool add_p_level)
110 {
111  libmesh_assert(elem);
112 
113  // call the orientation-independent shape functions
114  return fe_lagrange_2D_shape<L2_LAGRANGE>(elem->type(), static_cast<Order>(order + add_p_level * elem->p_level()), i, p);
115 }
116 
117 
118 template <>
120  const Elem * elem,
121  const unsigned int i,
122  const Point & p,
123  const bool add_p_level)
124 {
125  libmesh_assert(elem);
126  return fe_lagrange_2D_shape<LAGRANGE>(elem->type(), static_cast<Order>(fet.order + add_p_level * elem->p_level()), i, p);
127 }
128 
129 
130 
131 template <>
133  const Elem * elem,
134  const unsigned int i,
135  const Point & p,
136  const bool add_p_level)
137 {
138  libmesh_assert(elem);
139  return fe_lagrange_2D_shape<L2_LAGRANGE>(elem->type(), static_cast<Order>(fet.order + add_p_level * elem->p_level()), i, p);
140 }
141 
142 
143 template <>
145  const Order order,
146  const unsigned int i,
147  const unsigned int j,
148  const Point & p)
149 {
150  return fe_lagrange_2D_shape_deriv<LAGRANGE>(type, order, i, j, p);
151 }
152 
153 
154 
155 template <>
157  const Order order,
158  const unsigned int i,
159  const unsigned int j,
160  const Point & p)
161 {
162  return fe_lagrange_2D_shape_deriv<L2_LAGRANGE>(type, order, i, j, p);
163 }
164 
165 
166 
167 template <>
169  const Order order,
170  const unsigned int i,
171  const unsigned int j,
172  const Point & p,
173  const bool add_p_level)
174 {
175  libmesh_assert(elem);
176 
177  // call the orientation-independent shape functions
178  return fe_lagrange_2D_shape_deriv<LAGRANGE>(elem->type(), static_cast<Order>(order + add_p_level * elem->p_level()), i, j, p);
179 }
180 
181 
182 
183 template <>
185  const Order order,
186  const unsigned int i,
187  const unsigned int j,
188  const Point & p,
189  const bool add_p_level)
190 {
191  libmesh_assert(elem);
192 
193 
194  // call the orientation-independent shape functions
195  return fe_lagrange_2D_shape_deriv<L2_LAGRANGE>(elem->type(), static_cast<Order>(order + add_p_level * elem->p_level()), i, j, p);
196 }
197 
198 
199 template <>
201  const Elem * elem,
202  const unsigned int i,
203  const unsigned int j,
204  const Point & p,
205  const bool add_p_level)
206 {
207  libmesh_assert(elem);
208  return fe_lagrange_2D_shape_deriv<LAGRANGE>(elem->type(), static_cast<Order>(fet.order + add_p_level * elem->p_level()), i, j, p);
209 }
210 
211 
212 
213 template <>
215  const Elem * elem,
216  const unsigned int i,
217  const unsigned int j,
218  const Point & p,
219  const bool add_p_level)
220 {
221  libmesh_assert(elem);
222  return fe_lagrange_2D_shape_deriv<L2_LAGRANGE>(elem->type(), static_cast<Order>(fet.order + add_p_level * elem->p_level()), i, j, p);
223 }
224 
225 
226 
227 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
228 
229 template <>
231  const Order order,
232  const unsigned int i,
233  const unsigned int j,
234  const Point & p)
235 {
236  return fe_lagrange_2D_shape_second_deriv<LAGRANGE>(type, order, i, j, p);
237 }
238 
239 
240 
241 template <>
243  const Order order,
244  const unsigned int i,
245  const unsigned int j,
246  const Point & p)
247 {
248  return fe_lagrange_2D_shape_second_deriv<L2_LAGRANGE>(type, order, i, j, p);
249 }
250 
251 
252 
253 template <>
255  const Order order,
256  const unsigned int i,
257  const unsigned int j,
258  const Point & p,
259  const bool add_p_level)
260 {
261  libmesh_assert(elem);
262 
263  // call the orientation-independent shape functions
264  return fe_lagrange_2D_shape_second_deriv<LAGRANGE>(elem->type(), static_cast<Order>(order + add_p_level * elem->p_level()), i, j, p);
265 }
266 
267 
268 
269 template <>
271  const Order order,
272  const unsigned int i,
273  const unsigned int j,
274  const Point & p,
275  const bool add_p_level)
276 {
277  libmesh_assert(elem);
278 
279  // call the orientation-independent shape functions
280  return fe_lagrange_2D_shape_second_deriv<L2_LAGRANGE>(elem->type(), static_cast<Order>(order + add_p_level * elem->p_level()), i, j, p);
281 }
282 
283 
284 template <>
286  const Elem * elem,
287  const unsigned int i,
288  const unsigned int j,
289  const Point & p,
290  const bool add_p_level)
291 {
292  libmesh_assert(elem);
293  return fe_lagrange_2D_shape_second_deriv<LAGRANGE>(elem->type(), static_cast<Order>(fet.order + add_p_level * elem->p_level()), i, j, p);
294 }
295 
296 
297 
298 template <>
300  const Elem * elem,
301  const unsigned int i,
302  const unsigned int j,
303  const Point & p,
304  const bool add_p_level)
305 {
306  libmesh_assert(elem);
307  return fe_lagrange_2D_shape_second_deriv<L2_LAGRANGE>(elem->type(), static_cast<Order>(fet.order + add_p_level * elem->p_level()), i, j, p);
308 }
309 
310 #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
311 
312 } // namespace libMesh
313 
314 
315 
316 
317 // Anonymous namespace function definitions
318 namespace
319 {
320 using namespace libMesh;
321 
322 template <FEFamily T>
323 Real fe_lagrange_2D_shape(const ElemType type,
324  const Order order,
325  const unsigned int i,
326  const Point & p)
327 {
328 #if LIBMESH_DIM > 1
329 
330  switch (order)
331  {
332  // linear Lagrange shape functions
333  case FIRST:
334  {
335  switch (type)
336  {
337  case QUAD4:
338  case QUADSHELL4:
339  case QUAD8:
340  case QUADSHELL8:
341  case QUAD9:
342  {
343  // Compute quad shape functions as a tensor-product
344  const Real xi = p(0);
345  const Real eta = p(1);
346 
347  libmesh_assert_less (i, 4);
348 
349  // 0 1 2 3
350  static const unsigned int i0[] = {0, 1, 1, 0};
351  static const unsigned int i1[] = {0, 0, 1, 1};
352 
353  return (fe_lagrange_1D_linear_shape(i0[i], xi)*
354  fe_lagrange_1D_linear_shape(i1[i], eta));
355  }
356 
357  case TRI3:
358  case TRISHELL3:
359  case TRI6:
360  case TRI7:
361  {
362  const Real zeta1 = p(0);
363  const Real zeta2 = p(1);
364  const Real zeta0 = 1. - zeta1 - zeta2;
365 
366  libmesh_assert_less (i, 3);
367 
368  switch(i)
369  {
370  case 0:
371  return zeta0;
372 
373  case 1:
374  return zeta1;
375 
376  case 2:
377  return zeta2;
378 
379  default:
380  libmesh_error_msg("Invalid shape function index i = " << i);
381  }
382  }
383 
384  default:
385  libmesh_error_msg("ERROR: Unsupported 2D element type: " << Utility::enum_to_string(type));
386  }
387  }
388 
389 
390  // quadratic Lagrange shape functions
391  case SECOND:
392  {
393  switch (type)
394  {
395  case QUAD8:
396  case QUADSHELL8:
397  {
398  const Real xi = p(0);
399  const Real eta = p(1);
400 
401  libmesh_assert_less (i, 8);
402 
403  switch (i)
404  {
405  case 0:
406  return .25*(1. - xi)*(1. - eta)*(-1. - xi - eta);
407 
408  case 1:
409  return .25*(1. + xi)*(1. - eta)*(-1. + xi - eta);
410 
411  case 2:
412  return .25*(1. + xi)*(1. + eta)*(-1. + xi + eta);
413 
414  case 3:
415  return .25*(1. - xi)*(1. + eta)*(-1. - xi + eta);
416 
417  case 4:
418  return .5*(1. - xi*xi)*(1. - eta);
419 
420  case 5:
421  return .5*(1. + xi)*(1. - eta*eta);
422 
423  case 6:
424  return .5*(1. - xi*xi)*(1. + eta);
425 
426  case 7:
427  return .5*(1. - xi)*(1. - eta*eta);
428 
429  default:
430  libmesh_error_msg("Invalid shape function index i = " << i);
431  }
432  }
433 
434  case QUAD4:
435  libmesh_assert_msg(T == L2_LAGRANGE,
436  "High order on first order elements only supported for L2 families");
437  libmesh_fallthrough();
438  case QUAD9:
439  {
440 
441  // Compute quad shape functions as a tensor-product
442  const Real xi = p(0);
443  const Real eta = p(1);
444 
445  libmesh_assert_less (i, 9);
446 
447  // 0 1 2 3 4 5 6 7 8
448  static const unsigned int i0[] = {0, 1, 1, 0, 2, 1, 2, 0, 2};
449  static const unsigned int i1[] = {0, 0, 1, 1, 0, 2, 1, 2, 2};
450 
451  return (fe_lagrange_1D_quadratic_shape(i0[i], xi)*
452  fe_lagrange_1D_quadratic_shape(i1[i], eta));
453  }
454 
455  case TRI3:
456  libmesh_assert_msg(T == L2_LAGRANGE,
457  "High order on first order elements only supported for L2 families");
458  libmesh_fallthrough();
459  case TRI6:
460  case TRI7:
461  {
462  const Real zeta1 = p(0);
463  const Real zeta2 = p(1);
464  const Real zeta0 = 1. - zeta1 - zeta2;
465 
466  libmesh_assert_less (i, 6);
467 
468  switch(i)
469  {
470  case 0:
471  return 2.*zeta0*(zeta0-0.5);
472 
473  case 1:
474  return 2.*zeta1*(zeta1-0.5);
475 
476  case 2:
477  return 2.*zeta2*(zeta2-0.5);
478 
479  case 3:
480  return 4.*zeta0*zeta1;
481 
482  case 4:
483  return 4.*zeta1*zeta2;
484 
485  case 5:
486  return 4.*zeta2*zeta0;
487 
488  default:
489  libmesh_error_msg("Invalid shape function index i = " << i);
490  }
491  }
492 
493  default:
494  libmesh_error_msg("ERROR: Unsupported 2D element type: " << Utility::enum_to_string(type));
495  }
496  }
497 
498 
499 
500  // "cubic" (one cubic bubble) Lagrange shape functions on TRI7
501  case THIRD:
502  {
503  switch (type)
504  {
505  case TRI7:
506  {
507  const Real zeta1 = p(0);
508  const Real zeta2 = p(1);
509  const Real zeta0 = 1. - zeta1 - zeta2;
510  const Real bubble_27th = zeta0*zeta1*zeta2;
511 
512  libmesh_assert_less (i, 7);
513 
514  switch(i)
515  {
516  case 0:
517  return 2.*zeta0*(zeta0-0.5) + 3.*bubble_27th;
518 
519  case 1:
520  return 2.*zeta1*(zeta1-0.5) + 3.*bubble_27th;
521 
522  case 2:
523  return 2.*zeta2*(zeta2-0.5) + 3.*bubble_27th;
524 
525  case 3:
526  return 4.*zeta0*zeta1 - 12.*bubble_27th;
527 
528  case 4:
529  return 4.*zeta1*zeta2 - 12.*bubble_27th;
530 
531  case 5:
532  return 4.*zeta2*zeta0 - 12.*bubble_27th;
533 
534  case 6:
535  return 27.*bubble_27th;
536 
537  default:
538  libmesh_error_msg("Invalid shape function index i = " << i);
539  }
540  }
541 
542  default:
543  libmesh_error_msg("ERROR: Unsupported 2D element type: " << Utility::enum_to_string(type));
544  }
545  } // end case THIRD
546 
547 
548 
549  // unsupported order
550  default:
551  libmesh_error_msg("ERROR: Unsupported 2D FE order: " << order);
552  }
553 #else // LIBMESH_DIM > 1
554  libmesh_ignore(type, order, i, p);
555  libmesh_not_implemented();
556 #endif
557 }
558 
559 
560 
561 template <FEFamily T>
562 Real fe_lagrange_2D_shape_deriv(const ElemType type,
563  const Order order,
564  const unsigned int i,
565  const unsigned int j,
566  const Point & p)
567 {
568 #if LIBMESH_DIM > 1
569 
570  libmesh_assert_less (j, 2);
571 
572  switch (order)
573  {
574  // linear Lagrange shape functions
575  case FIRST:
576  {
577  switch (type)
578  {
579  case QUAD4:
580  case QUADSHELL4:
581  case QUAD8:
582  case QUADSHELL8:
583  case QUAD9:
584  {
585  // Compute quad shape functions as a tensor-product
586  const Real xi = p(0);
587  const Real eta = p(1);
588 
589  libmesh_assert_less (i, 4);
590 
591  // 0 1 2 3
592  static const unsigned int i0[] = {0, 1, 1, 0};
593  static const unsigned int i1[] = {0, 0, 1, 1};
594 
595  switch (j)
596  {
597  // d()/dxi
598  case 0:
599  return (fe_lagrange_1D_linear_shape_deriv(i0[i], 0, xi)*
600  fe_lagrange_1D_linear_shape (i1[i], eta));
601 
602  // d()/deta
603  case 1:
604  return (fe_lagrange_1D_linear_shape (i0[i], xi)*
605  fe_lagrange_1D_linear_shape_deriv(i1[i], 0, eta));
606 
607  default:
608  libmesh_error_msg("ERROR: Invalid derivative index j = " << j);
609  }
610  }
611 
612  case TRI3:
613  case TRISHELL3:
614  case TRI6:
615  case TRI7:
616  {
617  libmesh_assert_less (i, 3);
618 
619  const Real dzeta0dxi = -1.;
620  const Real dzeta1dxi = 1.;
621  const Real dzeta2dxi = 0.;
622 
623  const Real dzeta0deta = -1.;
624  const Real dzeta1deta = 0.;
625  const Real dzeta2deta = 1.;
626 
627  switch (j)
628  {
629  // d()/dxi
630  case 0:
631  {
632  switch(i)
633  {
634  case 0:
635  return dzeta0dxi;
636 
637  case 1:
638  return dzeta1dxi;
639 
640  case 2:
641  return dzeta2dxi;
642 
643  default:
644  libmesh_error_msg("Invalid shape function index i = " << i);
645  }
646  }
647  // d()/deta
648  case 1:
649  {
650  switch(i)
651  {
652  case 0:
653  return dzeta0deta;
654 
655  case 1:
656  return dzeta1deta;
657 
658  case 2:
659  return dzeta2deta;
660 
661  default:
662  libmesh_error_msg("Invalid shape function index i = " << i);
663  }
664  }
665  default:
666  libmesh_error_msg("ERROR: Invalid derivative index j = " << j);
667  }
668  }
669 
670  default:
671  libmesh_error_msg("ERROR: Unsupported 2D element type: " << Utility::enum_to_string(type));
672  }
673  }
674 
675 
676  // quadratic Lagrange shape functions
677  case SECOND:
678  {
679  switch (type)
680  {
681  case QUAD8:
682  case QUADSHELL8:
683  {
684  const Real xi = p(0);
685  const Real eta = p(1);
686 
687  libmesh_assert_less (i, 8);
688 
689  switch (j)
690  {
691  // d/dxi
692  case 0:
693  switch (i)
694  {
695  case 0:
696  return .25*(1. - eta)*((1. - xi)*(-1.) +
697  (-1.)*(-1. - xi - eta));
698 
699  case 1:
700  return .25*(1. - eta)*((1. + xi)*(1.) +
701  (1.)*(-1. + xi - eta));
702 
703  case 2:
704  return .25*(1. + eta)*((1. + xi)*(1.) +
705  (1.)*(-1. + xi + eta));
706 
707  case 3:
708  return .25*(1. + eta)*((1. - xi)*(-1.) +
709  (-1.)*(-1. - xi + eta));
710 
711  case 4:
712  return .5*(-2.*xi)*(1. - eta);
713 
714  case 5:
715  return .5*(1.)*(1. - eta*eta);
716 
717  case 6:
718  return .5*(-2.*xi)*(1. + eta);
719 
720  case 7:
721  return .5*(-1.)*(1. - eta*eta);
722 
723  default:
724  libmesh_error_msg("Invalid shape function index i = " << i);
725  }
726 
727  // d/deta
728  case 1:
729  switch (i)
730  {
731  case 0:
732  return .25*(1. - xi)*((1. - eta)*(-1.) +
733  (-1.)*(-1. - xi - eta));
734 
735  case 1:
736  return .25*(1. + xi)*((1. - eta)*(-1.) +
737  (-1.)*(-1. + xi - eta));
738 
739  case 2:
740  return .25*(1. + xi)*((1. + eta)*(1.) +
741  (1.)*(-1. + xi + eta));
742 
743  case 3:
744  return .25*(1. - xi)*((1. + eta)*(1.) +
745  (1.)*(-1. - xi + eta));
746 
747  case 4:
748  return .5*(1. - xi*xi)*(-1.);
749 
750  case 5:
751  return .5*(1. + xi)*(-2.*eta);
752 
753  case 6:
754  return .5*(1. - xi*xi)*(1.);
755 
756  case 7:
757  return .5*(1. - xi)*(-2.*eta);
758 
759  default:
760  libmesh_error_msg("Invalid shape function index i = " << i);
761  }
762 
763  default:
764  libmesh_error_msg("ERROR: Invalid derivative index j = " << j);
765  }
766  }
767 
768  case QUAD4:
769  libmesh_assert_msg(T == L2_LAGRANGE,
770  "High order on first order elements only supported for L2 families");
771  libmesh_fallthrough();
772  case QUAD9:
773  {
774  // Compute quad shape functions as a tensor-product
775  const Real xi = p(0);
776  const Real eta = p(1);
777 
778  libmesh_assert_less (i, 9);
779 
780  // 0 1 2 3 4 5 6 7 8
781  static const unsigned int i0[] = {0, 1, 1, 0, 2, 1, 2, 0, 2};
782  static const unsigned int i1[] = {0, 0, 1, 1, 0, 2, 1, 2, 2};
783 
784  switch (j)
785  {
786  // d()/dxi
787  case 0:
788  return (fe_lagrange_1D_quadratic_shape_deriv(i0[i], 0, xi)*
789  fe_lagrange_1D_quadratic_shape (i1[i], eta));
790 
791  // d()/deta
792  case 1:
793  return (fe_lagrange_1D_quadratic_shape (i0[i], xi)*
794  fe_lagrange_1D_quadratic_shape_deriv(i1[i], 0, eta));
795 
796  default:
797  libmesh_error_msg("ERROR: Invalid derivative index j = " << j);
798  }
799  }
800 
801  case TRI3:
802  libmesh_assert_msg(T == L2_LAGRANGE,
803  "High order on first order elements only supported for L2 families");
804  libmesh_fallthrough();
805  case TRI6:
806  case TRI7:
807  {
808  libmesh_assert_less (i, 6);
809 
810  const Real zeta1 = p(0);
811  const Real zeta2 = p(1);
812  const Real zeta0 = 1. - zeta1 - zeta2;
813 
814  const Real dzeta0dxi = -1.;
815  const Real dzeta1dxi = 1.;
816  const Real dzeta2dxi = 0.;
817 
818  const Real dzeta0deta = -1.;
819  const Real dzeta1deta = 0.;
820  const Real dzeta2deta = 1.;
821 
822  switch(j)
823  {
824  case 0:
825  {
826  switch(i)
827  {
828  case 0:
829  return (4.*zeta0-1.)*dzeta0dxi;
830 
831  case 1:
832  return (4.*zeta1-1.)*dzeta1dxi;
833 
834  case 2:
835  return (4.*zeta2-1.)*dzeta2dxi;
836 
837  case 3:
838  return 4.*zeta1*dzeta0dxi + 4.*zeta0*dzeta1dxi;
839 
840  case 4:
841  return 4.*zeta2*dzeta1dxi + 4.*zeta1*dzeta2dxi;
842 
843  case 5:
844  return 4.*zeta2*dzeta0dxi + 4*zeta0*dzeta2dxi;
845 
846  default:
847  libmesh_error_msg("Invalid shape function index i = " << i);
848  }
849  }
850 
851  case 1:
852  {
853  switch(i)
854  {
855  case 0:
856  return (4.*zeta0-1.)*dzeta0deta;
857 
858  case 1:
859  return (4.*zeta1-1.)*dzeta1deta;
860 
861  case 2:
862  return (4.*zeta2-1.)*dzeta2deta;
863 
864  case 3:
865  return 4.*zeta1*dzeta0deta + 4.*zeta0*dzeta1deta;
866 
867  case 4:
868  return 4.*zeta2*dzeta1deta + 4.*zeta1*dzeta2deta;
869 
870  case 5:
871  return 4.*zeta2*dzeta0deta + 4*zeta0*dzeta2deta;
872 
873  default:
874  libmesh_error_msg("Invalid shape function index i = " << i);
875  }
876  }
877  default:
878  libmesh_error_msg("ERROR: Invalid derivative index j = " << j);
879  }
880  }
881 
882  default:
883  libmesh_error_msg("ERROR: Unsupported 2D element type: " << Utility::enum_to_string(type));
884  }
885  }
886 
887 
888 
889  // "cubic" (one cubic bubble) Lagrange shape functions on TRI7
890  case THIRD:
891  {
892  switch (type)
893  {
894  case TRI7:
895  {
896  libmesh_assert_less (i, 7);
897 
898  const Real zeta1 = p(0);
899  const Real zeta2 = p(1);
900  const Real zeta0 = 1. - zeta1 - zeta2;
901  // const Real bubble_27th = zeta0*zeta1*zeta2;
902 
903  const Real dzeta0dxi = -1.;
904  const Real dzeta1dxi = 1.;
905  const Real dzeta2dxi = 0.;
906  const Real dbubbledxi = zeta2 * (1. - 2.*zeta1 - zeta2);
907 
908  const Real dzeta0deta = -1.;
909  const Real dzeta1deta = 0.;
910  const Real dzeta2deta = 1.;
911  const Real dbubbledeta= zeta1 * (1. - zeta1 - 2.*zeta2);
912 
913  switch(j)
914  {
915  case 0:
916  {
917  switch(i)
918  {
919  case 0:
920  return (4.*zeta0-1.)*dzeta0dxi + 3.*dbubbledxi;
921 
922  case 1:
923  return (4.*zeta1-1.)*dzeta1dxi + 3.*dbubbledxi;
924 
925  case 2:
926  return (4.*zeta2-1.)*dzeta2dxi + 3.*dbubbledxi;
927 
928  case 3:
929  return 4.*zeta1*dzeta0dxi + 4.*zeta0*dzeta1dxi - 12.*dbubbledxi;
930 
931  case 4:
932  return 4.*zeta2*dzeta1dxi + 4.*zeta1*dzeta2dxi - 12.*dbubbledxi;
933 
934  case 5:
935  return 4.*zeta2*dzeta0dxi + 4*zeta0*dzeta2dxi - 12.*dbubbledxi;
936 
937  case 6:
938  return 27.*dbubbledxi;
939 
940  default:
941  libmesh_error_msg("Invalid shape function index i = " << i);
942  }
943  }
944 
945  case 1:
946  {
947  switch(i)
948  {
949  case 0:
950  return (4.*zeta0-1.)*dzeta0deta + 3.*dbubbledeta;
951 
952  case 1:
953  return (4.*zeta1-1.)*dzeta1deta + 3.*dbubbledeta;
954 
955  case 2:
956  return (4.*zeta2-1.)*dzeta2deta + 3.*dbubbledeta;
957 
958  case 3:
959  return 4.*zeta1*dzeta0deta + 4.*zeta0*dzeta1deta - 12.*dbubbledeta;
960 
961  case 4:
962  return 4.*zeta2*dzeta1deta + 4.*zeta1*dzeta2deta - 12.*dbubbledeta;
963 
964  case 5:
965  return 4.*zeta2*dzeta0deta + 4*zeta0*dzeta2deta - 12.*dbubbledeta;
966 
967  case 6:
968  return 27.*dbubbledeta;
969 
970  default:
971  libmesh_error_msg("Invalid shape function index i = " << i);
972  }
973  }
974  default:
975  libmesh_error_msg("ERROR: Invalid derivative index j = " << j);
976  }
977  }
978 
979  default:
980  libmesh_error_msg("ERROR: Unsupported 2D element type: " << Utility::enum_to_string(type));
981  }
982  } // end case THIRD
983 
984 
985 
986  // unsupported order
987  default:
988  libmesh_error_msg("ERROR: Unsupported 2D FE order: " << order);
989  }
990 #else // LIBMESH_DIM > 1
991  libmesh_ignore(type, order, i, j, p);
992  libmesh_not_implemented();
993 #endif
994 }
995 
996 
997 
998 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
999 
1000 template <FEFamily T>
1001 Real fe_lagrange_2D_shape_second_deriv(const ElemType type,
1002  const Order order,
1003  const unsigned int i,
1004  const unsigned int j,
1005  const Point & p)
1006 {
1007 #if LIBMESH_DIM > 1
1008 
1009  // j = 0 ==> d^2 phi / dxi^2
1010  // j = 1 ==> d^2 phi / dxi deta
1011  // j = 2 ==> d^2 phi / deta^2
1012  libmesh_assert_less (j, 3);
1013 
1014  switch (order)
1015  {
1016  // linear Lagrange shape functions
1017  case FIRST:
1018  {
1019  switch (type)
1020  {
1021  case QUAD4:
1022  case QUADSHELL4:
1023  case QUAD8:
1024  case QUADSHELL8:
1025  case QUAD9:
1026  {
1027  // Compute quad shape functions as a tensor-product
1028  const Real xi = p(0);
1029  const Real eta = p(1);
1030 
1031  libmesh_assert_less (i, 4);
1032 
1033  // 0 1 2 3
1034  static const unsigned int i0[] = {0, 1, 1, 0};
1035  static const unsigned int i1[] = {0, 0, 1, 1};
1036 
1037  switch (j)
1038  {
1039  // d^2() / dxi^2
1040  case 0:
1041  return 0.;
1042 
1043  // d^2() / dxi deta
1044  case 1:
1045  return (fe_lagrange_1D_linear_shape_deriv(i0[i], 0, xi)*
1046  fe_lagrange_1D_linear_shape_deriv(i1[i], 0, eta));
1047 
1048  // d^2() / deta^2
1049  case 2:
1050  return 0.;
1051 
1052  default:
1053  libmesh_error_msg("ERROR: Invalid derivative index j = " << j);
1054  }
1055  }
1056 
1057  case TRI3:
1058  case TRISHELL3:
1059  case TRI6:
1060  case TRI7:
1061  {
1062  // All second derivatives for linear triangles are zero.
1063  return 0.;
1064  }
1065 
1066  default:
1067  libmesh_error_msg("ERROR: Unsupported 2D element type: " << Utility::enum_to_string(type));
1068 
1069  } // end switch (type)
1070  } // end case FIRST
1071 
1072 
1073  // quadratic Lagrange shape functions
1074  case SECOND:
1075  {
1076  switch (type)
1077  {
1078  case QUAD8:
1079  case QUADSHELL8:
1080  {
1081  const Real xi = p(0);
1082  const Real eta = p(1);
1083 
1084  libmesh_assert_less (j, 3);
1085 
1086  switch (j)
1087  {
1088  // d^2() / dxi^2
1089  case 0:
1090  {
1091  switch (i)
1092  {
1093  case 0:
1094  case 1:
1095  return 0.5*(1.-eta);
1096 
1097  case 2:
1098  case 3:
1099  return 0.5*(1.+eta);
1100 
1101  case 4:
1102  return eta - 1.;
1103 
1104  case 5:
1105  case 7:
1106  return 0.0;
1107 
1108  case 6:
1109  return -1. - eta;
1110 
1111  default:
1112  libmesh_error_msg("Invalid shape function index i = " << i);
1113  }
1114  }
1115 
1116  // d^2() / dxi deta
1117  case 1:
1118  {
1119  switch (i)
1120  {
1121  case 0:
1122  return 0.25*( 1. - 2.*xi - 2.*eta);
1123 
1124  case 1:
1125  return 0.25*(-1. - 2.*xi + 2.*eta);
1126 
1127  case 2:
1128  return 0.25*( 1. + 2.*xi + 2.*eta);
1129 
1130  case 3:
1131  return 0.25*(-1. + 2.*xi - 2.*eta);
1132 
1133  case 4:
1134  return xi;
1135 
1136  case 5:
1137  return -eta;
1138 
1139  case 6:
1140  return -xi;
1141 
1142  case 7:
1143  return eta;
1144 
1145  default:
1146  libmesh_error_msg("Invalid shape function index i = " << i);
1147  }
1148  }
1149 
1150  // d^2() / deta^2
1151  case 2:
1152  {
1153  switch (i)
1154  {
1155  case 0:
1156  case 3:
1157  return 0.5*(1.-xi);
1158 
1159  case 1:
1160  case 2:
1161  return 0.5*(1.+xi);
1162 
1163  case 4:
1164  case 6:
1165  return 0.0;
1166 
1167  case 5:
1168  return -1.0 - xi;
1169 
1170  case 7:
1171  return xi - 1.0;
1172 
1173  default:
1174  libmesh_error_msg("Invalid shape function index i = " << i);
1175  }
1176  }
1177 
1178  default:
1179  libmesh_error_msg("ERROR: Invalid derivative index j = " << j);
1180  } // end switch (j)
1181  } // end case QUAD8
1182 
1183  case QUAD4:
1184  libmesh_assert_msg(T == L2_LAGRANGE,
1185  "High order on first order elements only supported for L2 families");
1186  libmesh_fallthrough();
1187  case QUAD9:
1188  {
1189  // Compute QUAD9 second derivatives as tensor product
1190  const Real xi = p(0);
1191  const Real eta = p(1);
1192 
1193  libmesh_assert_less (i, 9);
1194 
1195  // 0 1 2 3 4 5 6 7 8
1196  static const unsigned int i0[] = {0, 1, 1, 0, 2, 1, 2, 0, 2};
1197  static const unsigned int i1[] = {0, 0, 1, 1, 0, 2, 1, 2, 2};
1198 
1199  switch (j)
1200  {
1201  // d^2() / dxi^2
1202  case 0:
1203  return (fe_lagrange_1D_quadratic_shape_second_deriv(i0[i], 0, xi)*
1204  fe_lagrange_1D_quadratic_shape (i1[i], eta));
1205 
1206  // d^2() / dxi deta
1207  case 1:
1208  return (fe_lagrange_1D_quadratic_shape_deriv(i0[i], 0, xi)*
1209  fe_lagrange_1D_quadratic_shape_deriv(i1[i], 0, eta));
1210 
1211  // d^2() / deta^2
1212  case 2:
1213  return (fe_lagrange_1D_quadratic_shape (i0[i], xi)*
1215 
1216  default:
1217  libmesh_error_msg("ERROR: Invalid derivative index j = " << j);
1218  } // end switch (j)
1219  } // end case QUAD9
1220 
1221  case TRI3:
1222  libmesh_assert_msg(T == L2_LAGRANGE,
1223  "High order on first order elements only supported for L2 families");
1224  libmesh_fallthrough();
1225  case TRI6:
1226  case TRI7:
1227  {
1228  const Real dzeta0dxi = -1.;
1229  const Real dzeta1dxi = 1.;
1230  const Real dzeta2dxi = 0.;
1231 
1232  const Real dzeta0deta = -1.;
1233  const Real dzeta1deta = 0.;
1234  const Real dzeta2deta = 1.;
1235 
1236  libmesh_assert_less (j, 3);
1237 
1238  switch (j)
1239  {
1240  // d^2() / dxi^2
1241  case 0:
1242  {
1243  switch (i)
1244  {
1245  case 0:
1246  return 4.*dzeta0dxi*dzeta0dxi;
1247 
1248  case 1:
1249  return 4.*dzeta1dxi*dzeta1dxi;
1250 
1251  case 2:
1252  return 4.*dzeta2dxi*dzeta2dxi;
1253 
1254  case 3:
1255  return 8.*dzeta0dxi*dzeta1dxi;
1256 
1257  case 4:
1258  return 8.*dzeta1dxi*dzeta2dxi;
1259 
1260  case 5:
1261  return 8.*dzeta0dxi*dzeta2dxi;
1262 
1263  default:
1264  libmesh_error_msg("Invalid shape function index i = " << i);
1265  }
1266  }
1267 
1268  // d^2() / dxi deta
1269  case 1:
1270  {
1271  switch (i)
1272  {
1273  case 0:
1274  return 4.*dzeta0dxi*dzeta0deta;
1275 
1276  case 1:
1277  return 4.*dzeta1dxi*dzeta1deta;
1278 
1279  case 2:
1280  return 4.*dzeta2dxi*dzeta2deta;
1281 
1282  case 3:
1283  return 4.*dzeta1deta*dzeta0dxi + 4.*dzeta0deta*dzeta1dxi;
1284 
1285  case 4:
1286  return 4.*dzeta2deta*dzeta1dxi + 4.*dzeta1deta*dzeta2dxi;
1287 
1288  case 5:
1289  return 4.*dzeta2deta*dzeta0dxi + 4.*dzeta0deta*dzeta2dxi;
1290 
1291  default:
1292  libmesh_error_msg("Invalid shape function index i = " << i);
1293  }
1294  }
1295 
1296  // d^2() / deta^2
1297  case 2:
1298  {
1299  switch (i)
1300  {
1301  case 0:
1302  return 4.*dzeta0deta*dzeta0deta;
1303 
1304  case 1:
1305  return 4.*dzeta1deta*dzeta1deta;
1306 
1307  case 2:
1308  return 4.*dzeta2deta*dzeta2deta;
1309 
1310  case 3:
1311  return 8.*dzeta0deta*dzeta1deta;
1312 
1313  case 4:
1314  return 8.*dzeta1deta*dzeta2deta;
1315 
1316  case 5:
1317  return 8.*dzeta0deta*dzeta2deta;
1318 
1319  default:
1320  libmesh_error_msg("Invalid shape function index i = " << i);
1321  }
1322  }
1323 
1324  default:
1325  libmesh_error_msg("ERROR: Invalid derivative index j = " << j);
1326  } // end switch (j)
1327  } // end case TRI6+TRI7
1328 
1329  default:
1330  libmesh_error_msg("ERROR: Unsupported 2D element type: " << Utility::enum_to_string(type));
1331  }
1332  } // end case SECOND
1333 
1334 
1335 
1336  // "cubic" (one cubic bubble) Lagrange shape functions on TRI7
1337  case THIRD:
1338  {
1339  switch (type)
1340  {
1341  case TRI3:
1342  libmesh_assert_msg(T == L2_LAGRANGE,
1343  "High order on first order elements only supported for L2 families");
1344  libmesh_fallthrough();
1345  case TRI6:
1346  case TRI7:
1347  {
1348  const Real zeta1 = p(0);
1349  const Real zeta2 = p(1);
1350  // const Real zeta0 = 1. - zeta1 - zeta2;
1351 
1352  const Real dzeta0dxi = -1.;
1353  const Real dzeta1dxi = 1.;
1354  const Real dzeta2dxi = 0.;
1355  // const Real dbubbledxi = zeta2 * (1. - 2.*zeta1 - zeta2);
1356  const Real d2bubbledxi2 = -2. * zeta2;
1357 
1358  const Real dzeta0deta = -1.;
1359  const Real dzeta1deta = 0.;
1360  const Real dzeta2deta = 1.;
1361  // const Real dbubbledeta= zeta1 * (1. - zeta1 - 2.*zeta2);
1362  const Real d2bubbledeta2 = -2. * zeta1;
1363 
1364  const Real d2bubbledxideta = (1. - 2.*zeta1 - 2.*zeta2);
1365 
1366  libmesh_assert_less (j, 3);
1367 
1368  switch (j)
1369  {
1370  // d^2() / dxi^2
1371  case 0:
1372  {
1373  switch (i)
1374  {
1375  case 0:
1376  return 4.*dzeta0dxi*dzeta0dxi + 3.*d2bubbledxi2;
1377 
1378  case 1:
1379  return 4.*dzeta1dxi*dzeta1dxi + 3.*d2bubbledxi2;
1380 
1381  case 2:
1382  return 4.*dzeta2dxi*dzeta2dxi + 3.*d2bubbledxi2;
1383 
1384  case 3:
1385  return 8.*dzeta0dxi*dzeta1dxi - 12.*d2bubbledxi2;
1386 
1387  case 4:
1388  return 8.*dzeta1dxi*dzeta2dxi - 12.*d2bubbledxi2;
1389 
1390  case 5:
1391  return 8.*dzeta0dxi*dzeta2dxi - 12.*d2bubbledxi2;
1392 
1393  case 6:
1394  return 27.*d2bubbledxi2;
1395 
1396  default:
1397  libmesh_error_msg("Invalid shape function index i = " << i);
1398  }
1399  }
1400 
1401  // d^2() / dxi deta
1402  case 1:
1403  {
1404  switch (i)
1405  {
1406  case 0:
1407  return 4.*dzeta0dxi*dzeta0deta + 3.*d2bubbledxideta;
1408 
1409  case 1:
1410  return 4.*dzeta1dxi*dzeta1deta + 3.*d2bubbledxideta;
1411 
1412  case 2:
1413  return 4.*dzeta2dxi*dzeta2deta + 3.*d2bubbledxideta;
1414 
1415  case 3:
1416  return 4.*dzeta1deta*dzeta0dxi + 4.*dzeta0deta*dzeta1dxi - 12.*d2bubbledxideta;
1417 
1418  case 4:
1419  return 4.*dzeta2deta*dzeta1dxi + 4.*dzeta1deta*dzeta2dxi - 12.*d2bubbledxideta;
1420 
1421  case 5:
1422  return 4.*dzeta2deta*dzeta0dxi + 4.*dzeta0deta*dzeta2dxi - 12.*d2bubbledxideta;
1423 
1424  case 6:
1425  return 27.*d2bubbledxideta;
1426 
1427  default:
1428  libmesh_error_msg("Invalid shape function index i = " << i);
1429  }
1430  }
1431 
1432  // d^2() / deta^2
1433  case 2:
1434  {
1435  switch (i)
1436  {
1437  case 0:
1438  return 4.*dzeta0deta*dzeta0deta + 3.*d2bubbledeta2;
1439 
1440  case 1:
1441  return 4.*dzeta1deta*dzeta1deta + 3.*d2bubbledeta2;
1442 
1443  case 2:
1444  return 4.*dzeta2deta*dzeta2deta + 3.*d2bubbledeta2;
1445 
1446  case 3:
1447  return 8.*dzeta0deta*dzeta1deta - 12.*d2bubbledeta2;
1448 
1449  case 4:
1450  return 8.*dzeta1deta*dzeta2deta - 12.*d2bubbledeta2;
1451 
1452  case 5:
1453  return 8.*dzeta0deta*dzeta2deta - 12.*d2bubbledeta2;
1454 
1455  case 6:
1456  return 27.*d2bubbledeta2;
1457 
1458  default:
1459  libmesh_error_msg("Invalid shape function index i = " << i);
1460  }
1461  }
1462 
1463  default:
1464  libmesh_error_msg("ERROR: Invalid derivative index j = " << j);
1465  } // end switch (j)
1466  } // end case TRI6+TRI7
1467 
1468  default:
1469  libmesh_error_msg("ERROR: Unsupported 2D element type: " << Utility::enum_to_string(type));
1470  }
1471  } // end case THIRD
1472 
1473  // unsupported order
1474  default:
1475  libmesh_error_msg("ERROR: Unsupported 2D FE order: " << order);
1476 
1477  } // end switch (order)
1478 
1479 #else // LIBMESH_DIM > 1
1480  libmesh_ignore(type, order, i, j, p);
1481  libmesh_not_implemented();
1482 #endif
1483 }
1484 
1485 #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
1486 
1487 } // anonymous namespace
Real fe_lagrange_1D_quadratic_shape_second_deriv(const unsigned int i, const unsigned int libmesh_dbg_var(j), const Real)
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:182
Real fe_lagrange_1D_linear_shape_deriv(const unsigned int i, const unsigned int libmesh_dbg_var(j), const Real)
ElemType
Defines an enum for geometric element types.
Order
defines an enum for polynomial orders.
Definition: enum_order.h:40
Real fe_lagrange_1D_quadratic_shape_deriv(const unsigned int i, const unsigned int libmesh_dbg_var(j), const Real xi)
static OutputShape shape(const ElemType t, const Order o, const unsigned int i, const Point &p)
Real fe_lagrange_1D_quadratic_shape(const unsigned int i, const Real xi)
This is the base class from which all geometric element types are derived.
Definition: elem.h:94
static OutputShape shape_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)
unsigned int p_level() const
Definition: elem.h:2945
OrderWrapper order
The approximation order of the element.
Definition: fe_type.h:201
The libMesh namespace provides an interface to certain functionality in the library.
LIBMESH_DEFAULT_VECTORIZED_FE(template<>Real FE< 0, BERNSTEIN)
void libmesh_ignore(const Args &...)
libmesh_assert(ctx)
std::string enum_to_string(const T e)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual ElemType type() const =0
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
static OutputShape shape_second_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)
Real fe_lagrange_1D_linear_shape(const unsigned int i, const Real xi)