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" 45 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 48 Real fe_lagrange_2D_shape_second_deriv(
const ElemType type,
54 #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES 74 return fe_lagrange_2D_shape<LAGRANGE>(type, order, i, p);
85 return fe_lagrange_2D_shape<L2_LAGRANGE>(type, order, i, p);
94 const bool add_p_level)
99 return fe_lagrange_2D_shape<LAGRANGE>(elem->
type(),
static_cast<Order>(order + add_p_level * elem->
p_level()), i, p);
107 const unsigned int i,
109 const bool add_p_level)
114 return fe_lagrange_2D_shape<L2_LAGRANGE>(elem->
type(),
static_cast<Order>(order + add_p_level * elem->
p_level()), i, p);
121 const unsigned int i,
123 const bool add_p_level)
126 return fe_lagrange_2D_shape<LAGRANGE>(elem->
type(),
static_cast<Order>(fet.
order + add_p_level * elem->
p_level()), i, p);
134 const unsigned int i,
136 const bool add_p_level)
139 return fe_lagrange_2D_shape<L2_LAGRANGE>(elem->
type(),
static_cast<Order>(fet.
order + add_p_level * elem->
p_level()), i, p);
146 const unsigned int i,
147 const unsigned int j,
150 return fe_lagrange_2D_shape_deriv<LAGRANGE>(type, order, i, j, p);
158 const unsigned int i,
159 const unsigned int j,
162 return fe_lagrange_2D_shape_deriv<L2_LAGRANGE>(type, order, i, j, p);
170 const unsigned int i,
171 const unsigned int j,
173 const bool add_p_level)
178 return fe_lagrange_2D_shape_deriv<LAGRANGE>(elem->
type(),
static_cast<Order>(order + add_p_level * elem->
p_level()), i, j, p);
186 const unsigned int i,
187 const unsigned int j,
189 const bool add_p_level)
195 return fe_lagrange_2D_shape_deriv<L2_LAGRANGE>(elem->
type(),
static_cast<Order>(order + add_p_level * elem->
p_level()), i, j, p);
202 const unsigned int i,
203 const unsigned int j,
205 const bool add_p_level)
208 return fe_lagrange_2D_shape_deriv<LAGRANGE>(elem->
type(),
static_cast<Order>(fet.
order + add_p_level * elem->
p_level()), i, j, p);
216 const unsigned int i,
217 const unsigned int j,
219 const bool add_p_level)
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);
227 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 232 const unsigned int i,
233 const unsigned int j,
236 return fe_lagrange_2D_shape_second_deriv<LAGRANGE>(type, order, i, j, p);
244 const unsigned int i,
245 const unsigned int j,
248 return fe_lagrange_2D_shape_second_deriv<L2_LAGRANGE>(type, order, i, j, p);
256 const unsigned int i,
257 const unsigned int j,
259 const bool add_p_level)
264 return fe_lagrange_2D_shape_second_deriv<LAGRANGE>(elem->
type(),
static_cast<Order>(order + add_p_level * elem->
p_level()), i, j, p);
272 const unsigned int i,
273 const unsigned int j,
275 const bool add_p_level)
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);
287 const unsigned int i,
288 const unsigned int j,
290 const bool add_p_level)
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);
301 const unsigned int i,
302 const unsigned int j,
304 const bool add_p_level)
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);
310 #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES 322 template <FEFamily T>
325 const unsigned int i,
344 const Real xi = p(0);
345 const Real eta = p(1);
347 libmesh_assert_less (i, 4);
350 static const unsigned int i0[] = {0, 1, 1, 0};
351 static const unsigned int i1[] = {0, 0, 1, 1};
362 const Real zeta1 = p(0);
363 const Real zeta2 = p(1);
364 const Real zeta0 = 1. - zeta1 - zeta2;
366 libmesh_assert_less (i, 3);
380 libmesh_error_msg(
"Invalid shape function index i = " << i);
398 const Real xi = p(0);
399 const Real eta = p(1);
401 libmesh_assert_less (i, 8);
406 return .25*(1. - xi)*(1. - eta)*(-1. - xi - eta);
409 return .25*(1. + xi)*(1. - eta)*(-1. + xi - eta);
412 return .25*(1. + xi)*(1. + eta)*(-1. + xi + eta);
415 return .25*(1. - xi)*(1. + eta)*(-1. - xi + eta);
418 return .5*(1. - xi*xi)*(1. - eta);
421 return .5*(1. + xi)*(1. - eta*eta);
424 return .5*(1. - xi*xi)*(1. + eta);
427 return .5*(1. - xi)*(1. - eta*eta);
430 libmesh_error_msg(
"Invalid shape function index i = " << i);
436 "High order on first order elements only supported for L2 families");
437 libmesh_fallthrough();
442 const Real xi = p(0);
443 const Real eta = p(1);
445 libmesh_assert_less (i, 9);
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};
457 "High order on first order elements only supported for L2 families");
458 libmesh_fallthrough();
462 const Real zeta1 = p(0);
463 const Real zeta2 = p(1);
464 const Real zeta0 = 1. - zeta1 - zeta2;
466 libmesh_assert_less (i, 6);
471 return 2.*zeta0*(zeta0-0.5);
474 return 2.*zeta1*(zeta1-0.5);
477 return 2.*zeta2*(zeta2-0.5);
480 return 4.*zeta0*zeta1;
483 return 4.*zeta1*zeta2;
486 return 4.*zeta2*zeta0;
489 libmesh_error_msg(
"Invalid shape function index i = " << i);
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;
512 libmesh_assert_less (i, 7);
517 return 2.*zeta0*(zeta0-0.5) + 3.*bubble_27th;
520 return 2.*zeta1*(zeta1-0.5) + 3.*bubble_27th;
523 return 2.*zeta2*(zeta2-0.5) + 3.*bubble_27th;
526 return 4.*zeta0*zeta1 - 12.*bubble_27th;
529 return 4.*zeta1*zeta2 - 12.*bubble_27th;
532 return 4.*zeta2*zeta0 - 12.*bubble_27th;
535 return 27.*bubble_27th;
538 libmesh_error_msg(
"Invalid shape function index i = " << i);
551 libmesh_error_msg(
"ERROR: Unsupported 2D FE order: " << order);
553 #else // LIBMESH_DIM > 1 555 libmesh_not_implemented();
561 template <FEFamily T>
564 const unsigned int i,
565 const unsigned int j,
570 libmesh_assert_less (j, 2);
586 const Real xi = p(0);
587 const Real eta = p(1);
589 libmesh_assert_less (i, 4);
592 static const unsigned int i0[] = {0, 1, 1, 0};
593 static const unsigned int i1[] = {0, 0, 1, 1};
608 libmesh_error_msg(
"ERROR: Invalid derivative index j = " << j);
617 libmesh_assert_less (i, 3);
619 const Real dzeta0dxi = -1.;
620 const Real dzeta1dxi = 1.;
621 const Real dzeta2dxi = 0.;
623 const Real dzeta0deta = -1.;
624 const Real dzeta1deta = 0.;
625 const Real dzeta2deta = 1.;
644 libmesh_error_msg(
"Invalid shape function index i = " << i);
662 libmesh_error_msg(
"Invalid shape function index i = " << i);
666 libmesh_error_msg(
"ERROR: Invalid derivative index j = " << j);
684 const Real xi = p(0);
685 const Real eta = p(1);
687 libmesh_assert_less (i, 8);
696 return .25*(1. - eta)*((1. - xi)*(-1.) +
697 (-1.)*(-1. - xi - eta));
700 return .25*(1. - eta)*((1. + xi)*(1.) +
701 (1.)*(-1. + xi - eta));
704 return .25*(1. + eta)*((1. + xi)*(1.) +
705 (1.)*(-1. + xi + eta));
708 return .25*(1. + eta)*((1. - xi)*(-1.) +
709 (-1.)*(-1. - xi + eta));
712 return .5*(-2.*xi)*(1. - eta);
715 return .5*(1.)*(1. - eta*eta);
718 return .5*(-2.*xi)*(1. + eta);
721 return .5*(-1.)*(1. - eta*eta);
724 libmesh_error_msg(
"Invalid shape function index i = " << i);
732 return .25*(1. - xi)*((1. - eta)*(-1.) +
733 (-1.)*(-1. - xi - eta));
736 return .25*(1. + xi)*((1. - eta)*(-1.) +
737 (-1.)*(-1. + xi - eta));
740 return .25*(1. + xi)*((1. + eta)*(1.) +
741 (1.)*(-1. + xi + eta));
744 return .25*(1. - xi)*((1. + eta)*(1.) +
745 (1.)*(-1. - xi + eta));
748 return .5*(1. - xi*xi)*(-1.);
751 return .5*(1. + xi)*(-2.*eta);
754 return .5*(1. - xi*xi)*(1.);
757 return .5*(1. - xi)*(-2.*eta);
760 libmesh_error_msg(
"Invalid shape function index i = " << i);
764 libmesh_error_msg(
"ERROR: Invalid derivative index j = " << j);
770 "High order on first order elements only supported for L2 families");
771 libmesh_fallthrough();
775 const Real xi = p(0);
776 const Real eta = p(1);
778 libmesh_assert_less (i, 9);
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};
797 libmesh_error_msg(
"ERROR: Invalid derivative index j = " << j);
803 "High order on first order elements only supported for L2 families");
804 libmesh_fallthrough();
808 libmesh_assert_less (i, 6);
810 const Real zeta1 = p(0);
811 const Real zeta2 = p(1);
812 const Real zeta0 = 1. - zeta1 - zeta2;
814 const Real dzeta0dxi = -1.;
815 const Real dzeta1dxi = 1.;
816 const Real dzeta2dxi = 0.;
818 const Real dzeta0deta = -1.;
819 const Real dzeta1deta = 0.;
820 const Real dzeta2deta = 1.;
829 return (4.*zeta0-1.)*dzeta0dxi;
832 return (4.*zeta1-1.)*dzeta1dxi;
835 return (4.*zeta2-1.)*dzeta2dxi;
838 return 4.*zeta1*dzeta0dxi + 4.*zeta0*dzeta1dxi;
841 return 4.*zeta2*dzeta1dxi + 4.*zeta1*dzeta2dxi;
844 return 4.*zeta2*dzeta0dxi + 4*zeta0*dzeta2dxi;
847 libmesh_error_msg(
"Invalid shape function index i = " << i);
856 return (4.*zeta0-1.)*dzeta0deta;
859 return (4.*zeta1-1.)*dzeta1deta;
862 return (4.*zeta2-1.)*dzeta2deta;
865 return 4.*zeta1*dzeta0deta + 4.*zeta0*dzeta1deta;
868 return 4.*zeta2*dzeta1deta + 4.*zeta1*dzeta2deta;
871 return 4.*zeta2*dzeta0deta + 4*zeta0*dzeta2deta;
874 libmesh_error_msg(
"Invalid shape function index i = " << i);
878 libmesh_error_msg(
"ERROR: Invalid derivative index j = " << j);
896 libmesh_assert_less (i, 7);
898 const Real zeta1 = p(0);
899 const Real zeta2 = p(1);
900 const Real zeta0 = 1. - zeta1 - zeta2;
903 const Real dzeta0dxi = -1.;
904 const Real dzeta1dxi = 1.;
905 const Real dzeta2dxi = 0.;
906 const Real dbubbledxi = zeta2 * (1. - 2.*zeta1 - zeta2);
908 const Real dzeta0deta = -1.;
909 const Real dzeta1deta = 0.;
910 const Real dzeta2deta = 1.;
911 const Real dbubbledeta= zeta1 * (1. - zeta1 - 2.*zeta2);
920 return (4.*zeta0-1.)*dzeta0dxi + 3.*dbubbledxi;
923 return (4.*zeta1-1.)*dzeta1dxi + 3.*dbubbledxi;
926 return (4.*zeta2-1.)*dzeta2dxi + 3.*dbubbledxi;
929 return 4.*zeta1*dzeta0dxi + 4.*zeta0*dzeta1dxi - 12.*dbubbledxi;
932 return 4.*zeta2*dzeta1dxi + 4.*zeta1*dzeta2dxi - 12.*dbubbledxi;
935 return 4.*zeta2*dzeta0dxi + 4*zeta0*dzeta2dxi - 12.*dbubbledxi;
938 return 27.*dbubbledxi;
941 libmesh_error_msg(
"Invalid shape function index i = " << i);
950 return (4.*zeta0-1.)*dzeta0deta + 3.*dbubbledeta;
953 return (4.*zeta1-1.)*dzeta1deta + 3.*dbubbledeta;
956 return (4.*zeta2-1.)*dzeta2deta + 3.*dbubbledeta;
959 return 4.*zeta1*dzeta0deta + 4.*zeta0*dzeta1deta - 12.*dbubbledeta;
962 return 4.*zeta2*dzeta1deta + 4.*zeta1*dzeta2deta - 12.*dbubbledeta;
965 return 4.*zeta2*dzeta0deta + 4*zeta0*dzeta2deta - 12.*dbubbledeta;
968 return 27.*dbubbledeta;
971 libmesh_error_msg(
"Invalid shape function index i = " << i);
975 libmesh_error_msg(
"ERROR: Invalid derivative index j = " << j);
988 libmesh_error_msg(
"ERROR: Unsupported 2D FE order: " << order);
990 #else // LIBMESH_DIM > 1 992 libmesh_not_implemented();
998 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 1000 template <FEFamily T>
1001 Real fe_lagrange_2D_shape_second_deriv(
const ElemType type,
1003 const unsigned int i,
1004 const unsigned int j,
1012 libmesh_assert_less (j, 3);
1028 const Real xi = p(0);
1029 const Real eta = p(1);
1031 libmesh_assert_less (i, 4);
1034 static const unsigned int i0[] = {0, 1, 1, 0};
1035 static const unsigned int i1[] = {0, 0, 1, 1};
1053 libmesh_error_msg(
"ERROR: Invalid derivative index j = " << j);
1081 const Real xi = p(0);
1082 const Real eta = p(1);
1084 libmesh_assert_less (j, 3);
1095 return 0.5*(1.-eta);
1099 return 0.5*(1.+eta);
1112 libmesh_error_msg(
"Invalid shape function index i = " << i);
1122 return 0.25*( 1. - 2.*xi - 2.*eta);
1125 return 0.25*(-1. - 2.*xi + 2.*eta);
1128 return 0.25*( 1. + 2.*xi + 2.*eta);
1131 return 0.25*(-1. + 2.*xi - 2.*eta);
1146 libmesh_error_msg(
"Invalid shape function index i = " << i);
1174 libmesh_error_msg(
"Invalid shape function index i = " << i);
1179 libmesh_error_msg(
"ERROR: Invalid derivative index j = " << j);
1185 "High order on first order elements only supported for L2 families");
1186 libmesh_fallthrough();
1190 const Real xi = p(0);
1191 const Real eta = p(1);
1193 libmesh_assert_less (i, 9);
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};
1217 libmesh_error_msg(
"ERROR: Invalid derivative index j = " << j);
1223 "High order on first order elements only supported for L2 families");
1224 libmesh_fallthrough();
1228 const Real dzeta0dxi = -1.;
1229 const Real dzeta1dxi = 1.;
1230 const Real dzeta2dxi = 0.;
1232 const Real dzeta0deta = -1.;
1233 const Real dzeta1deta = 0.;
1234 const Real dzeta2deta = 1.;
1236 libmesh_assert_less (j, 3);
1246 return 4.*dzeta0dxi*dzeta0dxi;
1249 return 4.*dzeta1dxi*dzeta1dxi;
1252 return 4.*dzeta2dxi*dzeta2dxi;
1255 return 8.*dzeta0dxi*dzeta1dxi;
1258 return 8.*dzeta1dxi*dzeta2dxi;
1261 return 8.*dzeta0dxi*dzeta2dxi;
1264 libmesh_error_msg(
"Invalid shape function index i = " << i);
1274 return 4.*dzeta0dxi*dzeta0deta;
1277 return 4.*dzeta1dxi*dzeta1deta;
1280 return 4.*dzeta2dxi*dzeta2deta;
1283 return 4.*dzeta1deta*dzeta0dxi + 4.*dzeta0deta*dzeta1dxi;
1286 return 4.*dzeta2deta*dzeta1dxi + 4.*dzeta1deta*dzeta2dxi;
1289 return 4.*dzeta2deta*dzeta0dxi + 4.*dzeta0deta*dzeta2dxi;
1292 libmesh_error_msg(
"Invalid shape function index i = " << i);
1302 return 4.*dzeta0deta*dzeta0deta;
1305 return 4.*dzeta1deta*dzeta1deta;
1308 return 4.*dzeta2deta*dzeta2deta;
1311 return 8.*dzeta0deta*dzeta1deta;
1314 return 8.*dzeta1deta*dzeta2deta;
1317 return 8.*dzeta0deta*dzeta2deta;
1320 libmesh_error_msg(
"Invalid shape function index i = " << i);
1325 libmesh_error_msg(
"ERROR: Invalid derivative index j = " << j);
1343 "High order on first order elements only supported for L2 families");
1344 libmesh_fallthrough();
1348 const Real zeta1 = p(0);
1349 const Real zeta2 = p(1);
1352 const Real dzeta0dxi = -1.;
1353 const Real dzeta1dxi = 1.;
1354 const Real dzeta2dxi = 0.;
1356 const Real d2bubbledxi2 = -2. * zeta2;
1358 const Real dzeta0deta = -1.;
1359 const Real dzeta1deta = 0.;
1360 const Real dzeta2deta = 1.;
1362 const Real d2bubbledeta2 = -2. * zeta1;
1364 const Real d2bubbledxideta = (1. - 2.*zeta1 - 2.*zeta2);
1366 libmesh_assert_less (j, 3);
1376 return 4.*dzeta0dxi*dzeta0dxi + 3.*d2bubbledxi2;
1379 return 4.*dzeta1dxi*dzeta1dxi + 3.*d2bubbledxi2;
1382 return 4.*dzeta2dxi*dzeta2dxi + 3.*d2bubbledxi2;
1385 return 8.*dzeta0dxi*dzeta1dxi - 12.*d2bubbledxi2;
1388 return 8.*dzeta1dxi*dzeta2dxi - 12.*d2bubbledxi2;
1391 return 8.*dzeta0dxi*dzeta2dxi - 12.*d2bubbledxi2;
1394 return 27.*d2bubbledxi2;
1397 libmesh_error_msg(
"Invalid shape function index i = " << i);
1407 return 4.*dzeta0dxi*dzeta0deta + 3.*d2bubbledxideta;
1410 return 4.*dzeta1dxi*dzeta1deta + 3.*d2bubbledxideta;
1413 return 4.*dzeta2dxi*dzeta2deta + 3.*d2bubbledxideta;
1416 return 4.*dzeta1deta*dzeta0dxi + 4.*dzeta0deta*dzeta1dxi - 12.*d2bubbledxideta;
1419 return 4.*dzeta2deta*dzeta1dxi + 4.*dzeta1deta*dzeta2dxi - 12.*d2bubbledxideta;
1422 return 4.*dzeta2deta*dzeta0dxi + 4.*dzeta0deta*dzeta2dxi - 12.*d2bubbledxideta;
1425 return 27.*d2bubbledxideta;
1428 libmesh_error_msg(
"Invalid shape function index i = " << i);
1438 return 4.*dzeta0deta*dzeta0deta + 3.*d2bubbledeta2;
1441 return 4.*dzeta1deta*dzeta1deta + 3.*d2bubbledeta2;
1444 return 4.*dzeta2deta*dzeta2deta + 3.*d2bubbledeta2;
1447 return 8.*dzeta0deta*dzeta1deta - 12.*d2bubbledeta2;
1450 return 8.*dzeta1deta*dzeta2deta - 12.*d2bubbledeta2;
1453 return 8.*dzeta0deta*dzeta2deta - 12.*d2bubbledeta2;
1456 return 27.*d2bubbledeta2;
1459 libmesh_error_msg(
"Invalid shape function index i = " << i);
1464 libmesh_error_msg(
"ERROR: Invalid derivative index j = " << j);
1475 libmesh_error_msg(
"ERROR: Unsupported 2D FE order: " << order);
1479 #else // LIBMESH_DIM > 1 1481 libmesh_not_implemented();
1485 #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES 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...
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.
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.
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
OrderWrapper order
The approximation order of the element.
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 &...)
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.
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)