19 #include "libmesh/libmesh_config.h"    21 #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES    24 #include "libmesh/fe.h"    25 #include "libmesh/elem.h"    26 #include "libmesh/number_lookups.h"    27 #include "libmesh/utility.h"    28 #include "libmesh/enum_to_string.h"    40 std::pair<unsigned int, unsigned int> quad_i0_i1 (
const unsigned int i,
    41                                                   const Order totalorder,
    44   libmesh_assert_less (i, (totalorder+1u)*(totalorder+1u));
    60   else if (i < totalorder + 3u)
    61     { i0 = i - 2; i1 = 0; }
    62   else if (i < 2u*totalorder + 2)
    63     { i0 = 1; i1 = i - totalorder - 1; }
    64   else if (i < 3u*totalorder + 1)
    65     { i0 = i - 2u*totalorder; i1 = 1; }
    66   else if (i < 4u*totalorder)
    67     { i0 = 0; i1 = i - 3u*totalorder + 1; }
    71       unsigned int basisnum = i - 4*totalorder;
    79   if      ((i>= 4                 && i<= 4+  totalorder-2u) && elem.
point(0) > elem.
point(1)) i0=totalorder+2-i0;
    80   else if ((i>= 4+  totalorder-1u && i<= 4+2*totalorder-3u) && elem.
point(1) > elem.
point(2)) i1=totalorder+2-i1;
    81   else if ((i>= 4+2*totalorder-2u && i<= 4+3*totalorder-4u) && elem.
point(3) > elem.
point(2)) i0=totalorder+2-i0;
    82   else if ((i>= 4+3*totalorder-3u && i<= 4+4*totalorder-5u) && elem.
point(0) > elem.
point(3)) i1=totalorder+2-i1;
    84   return std::make_pair(i0, i1);
   101                             const bool add_p_level)
   107   const Order totalorder =
   108     order + add_p_level*elem->
p_level();
   123         auto [i0, i1] = quad_i0_i1(i, totalorder, *elem);
   132         libmesh_assert_less (totalorder, 3);
   134         const Real xi  = p(0);
   135         const Real eta = p(1);
   137         libmesh_assert_less (i, 8);
   140         static const unsigned int i0[] = {0, 1, 1, 0, 2, 1, 2, 0, 2};
   141         static const unsigned int i1[] = {0, 0, 1, 1, 0, 2, 1, 2, 2};
   142         static const Real scal[] = {-0.25, -0.25, -0.25, -0.25, 0.5, 0.5, 0.5, 0.5};
   155       libmesh_assert_less (totalorder, 2);
   156       libmesh_fallthrough();
   167             libmesh_assert_less (i, 3);
   176                 libmesh_error_msg(
"Invalid shape function index i = " << i);
   185             libmesh_assert_less (i, 6);
   193               case 3: 
return 2.*x*r;
   194               case 4: 
return 2.*x*y;
   195               case 5: 
return 2.*r*y;
   198                 libmesh_error_msg(
"Invalid shape function index i = " << i);
   206             libmesh_assert_less (i, 10);
   208             unsigned int shape=i;
   217               case 0: 
return r*r*r;
   218               case 1: 
return x*x*x;
   219               case 2: 
return y*y*y;
   221               case 3: 
return 3.*x*r*r;
   222               case 4: 
return 3.*x*x*r;
   224               case 5: 
return 3.*y*x*x;
   225               case 6: 
return 3.*y*y*x;
   227               case 7: 
return 3.*y*r*r;
   228               case 8: 
return 3.*y*y*r;
   230               case 9: 
return 6.*x*y*r;
   233                 libmesh_error_msg(
"Invalid shape function index shape = " << shape);
   241             unsigned int shape=i;
   243             libmesh_assert_less (i, 15);
   253               case  0: 
return r*r*r*r;
   254               case  1: 
return x*x*x*x;
   255               case  2: 
return y*y*y*y;
   258               case  3: 
return 4.*x*r*r*r;
   259               case  4: 
return 6.*x*x*r*r;
   260               case  5: 
return 4.*x*x*x*r;
   262               case  6: 
return 4.*y*x*x*x;
   263               case  7: 
return 6.*y*y*x*x;
   264               case  8: 
return 4.*y*y*y*x;
   266               case  9: 
return 4.*y*r*r*r;
   267               case 10: 
return 6.*y*y*r*r;
   268               case 11: 
return 4.*y*y*y*r;
   271               case 12: 
return 12.*x*y*r*r;
   272               case 13: 
return 12.*x*x*y*r;
   273               case 14: 
return 12.*x*y*y*r;
   276                 libmesh_error_msg(
"Invalid shape function index shape = " << shape);
   284             unsigned int shape=i;
   286             libmesh_assert_less (i, 21);
   295               case  0: 
return pow<5>(r);
   296               case  1: 
return pow<5>(x);
   297               case  2: 
return pow<5>(y);
   300               case  3: 
return  5.*x        *pow<4>(r);
   301               case  4: 
return 10.*pow<2>(x)*pow<3>(r);
   302               case  5: 
return 10.*pow<3>(x)*pow<2>(r);
   303               case  6: 
return  5.*pow<4>(x)*r;
   305               case  7: 
return  5.*y   *pow<4>(x);
   306               case  8: 
return 10.*pow<2>(y)*pow<3>(x);
   307               case  9: 
return 10.*pow<3>(y)*pow<2>(x);
   308               case 10: 
return  5.*pow<4>(y)*x;
   310               case 11: 
return  5.*y   *pow<4>(r);
   311               case 12: 
return 10.*pow<2>(y)*pow<3>(r);
   312               case 13: 
return 10.*pow<3>(y)*pow<2>(r);
   313               case 14: 
return  5.*pow<4>(y)*r;
   316               case 15: 
return 20.*x*y*pow<3>(r);
   317               case 16: 
return 30.*x*pow<2>(y)*pow<2>(r);
   318               case 17: 
return 30.*pow<2>(x)*y*pow<2>(r);
   319               case 18: 
return 20.*x*pow<3>(y)*r;
   320               case 19: 
return 20.*pow<3>(x)*y*r;
   321               case 20: 
return 30.*pow<2>(x)*pow<2>(y)*r;
   324                 libmesh_error_msg(
"Invalid shape function index shape = " << shape);
   332             unsigned int shape=i;
   334             libmesh_assert_less (i, 28);
   343               case  0: 
return pow<6>(r);
   344               case  1: 
return pow<6>(x);
   345               case  2: 
return pow<6>(y);
   348               case  3: 
return  6.*x        *pow<5>(r);
   349               case  4: 
return 15.*pow<2>(x)*pow<4>(r);
   350               case  5: 
return 20.*pow<3>(x)*pow<3>(r);
   351               case  6: 
return 15.*pow<4>(x)*pow<2>(r);
   352               case  7: 
return  6.*pow<5>(x)*r;
   354               case  8: 
return  6.*y        *pow<5>(x);
   355               case  9: 
return 15.*pow<2>(y)*pow<4>(x);
   356               case 10: 
return 20.*pow<3>(y)*pow<3>(x);
   357               case 11: 
return 15.*pow<4>(y)*pow<2>(x);
   358               case 12: 
return  6.*pow<5>(y)*x;
   360               case 13: 
return  6.*y        *pow<5>(r);
   361               case 14: 
return 15.*pow<2>(y)*pow<4>(r);
   362               case 15: 
return 20.*pow<3>(y)*pow<3>(r);
   363               case 16: 
return 15.*pow<4>(y)*pow<2>(r);
   364               case 17: 
return  6.*pow<5>(y)*r;
   367               case 18: 
return 30.*x*y*pow<4>(r);
   368               case 19: 
return 60.*x*pow<2>(y)*pow<3>(r);
   369               case 20: 
return 60.*  pow<2>(x)*y*pow<3>(r);
   370               case 21: 
return 60.*x*pow<3>(y)*pow<2>(r);
   371               case 22: 
return 60.*pow<3>(x)*y*pow<2>(r);
   372               case 23: 
return 90.*pow<2>(x)*pow<2>(y)*pow<2>(r);
   373               case 24: 
return 30.*x*pow<4>(y)*r;
   374               case 25: 
return 60.*pow<2>(x)*pow<3>(y)*r;
   375               case 26: 
return 60.*pow<3>(x)*pow<2>(y)*r;
   376               case 27: 
return 30.*pow<4>(x)*y*r;
   379                 libmesh_error_msg(
"Invalid shape function index shape = " << shape);
   383           libmesh_error_msg(
"Invalid totalorder = " << totalorder);
   398   libmesh_error_msg(
"Bernstein polynomials require the element type \nbecause edge orientation is needed.");
   406                             const unsigned int i,
   408                             const bool add_p_level)
   418                                   const unsigned int i,
   419                                   const unsigned int j,
   421                                   const bool add_p_level)
   427   const Order totalorder =
   428     order + add_p_level*elem->
p_level();
   438         auto [i0, i1] = quad_i0_i1(i, totalorder, *elem);
   453             libmesh_error_msg(
"Invalid shape function derivative j = " << j);
   462         libmesh_assert_less (totalorder, 3);
   464         const Real xi  = p(0);
   465         const Real eta = p(1);
   467         libmesh_assert_less (i, 8);
   470         static const unsigned int i0[] = {0, 1, 1, 0, 2, 1, 2, 0, 2};
   471         static const unsigned int i1[] = {0, 0, 1, 1, 0, 2, 1, 2, 2};
   472         static const Real scal[] = {-0.25, -0.25, -0.25, -0.25, 0.5, 0.5, 0.5, 0.5};
   492             libmesh_error_msg(
"Invalid shape function derivative j = " << j);
   498       libmesh_assert_less (totalorder, 2);
   499       libmesh_fallthrough();
   521   libmesh_error_msg(
"Bernstein polynomials require the element type \nbecause edge orientation is needed.");
   528                                   const unsigned int i,
   529                                   const unsigned int j,
   531                                   const bool add_p_level)
   539 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES   546                                          const unsigned int i,
   547                                          const unsigned int j,
   549                                          const bool add_p_level)
   555   const Order totalorder =
   556     order + add_p_level*elem->
p_level();
   566         auto [i0, i1] = quad_i0_i1(i, totalorder, *elem);
   586             libmesh_error_msg(
"Invalid shape function derivative j = " << j);
   593       libmesh_assert_less (totalorder, 2);
   594       libmesh_fallthrough();
   617   libmesh_error_msg(
"Bernstein polynomials require the element type \nbecause edge orientation is needed.");
   625                                          const unsigned int i,
   626                                          const unsigned int j,
   628                                          const bool add_p_level)
   641 #endif// LIBMESH_ENABLE_HIGHER_ORDER_SHAPES class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
ElemType
Defines an enum for geometric element types. 
Order
defines an enum for polynomial orders. 
static OutputShape shape(const ElemType t, const Order o, const unsigned int i, const Point &p)
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. 
OutputShape fe_fdm_deriv(const Elem *elem, const Order order, const unsigned int i, const unsigned int j, const Point &p, const bool add_p_level, OutputShape(*shape_func)(const Elem *, const Order, const unsigned int, const Point &, const bool))
Helper functions for finite differenced derivatives in cases where analytical calculations haven't be...
const unsigned char square_number_column[]
LIBMESH_DEFAULT_VECTORIZED_FE(template<>Real FE< 0, BERNSTEIN)
A specific instantiation of the FEBase class. 
bool positive_edge_orientation(const unsigned int i) const
OutputShape fe_fdm_second_deriv(const Elem *elem, const Order order, const unsigned int i, const unsigned int j, const Point &p, const bool add_p_level, OutputShape(*deriv_func)(const Elem *, const Order, const unsigned int, const unsigned int, const Point &, const bool))
std::string enum_to_string(const T e)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const unsigned char square_number_row[]
virtual ElemType type() const =0
A Point defines a location in LIBMESH_DIM dimensional Real space. 
const Point & point(const unsigned int i) const
static OutputShape shape_second_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)