20 #include "libmesh/libmesh_config.h" 
   21 #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES 
   23 #include "libmesh/fe.h" 
   24 #include "libmesh/elem.h" 
   25 #include "libmesh/number_lookups.h" 
   26 #include "libmesh/utility.h" 
   39   libmesh_error_msg(
"Bernstein polynomials require the element type \nbecause edge orientation is needed.");
 
   50                             const bool add_p_level)
 
   56   const Order totalorder =
 
   57     static_cast<Order>(order + add_p_level * elem->
p_level());
 
   72         const Real eta = p(1);
 
   74         libmesh_assert_less (i, (totalorder+1u)*(totalorder+1u));
 
   95         else if (i < totalorder + 3u)
 
   96           { i0 = i - 2; i1 = 0; }
 
   97         else if (i < 2u*totalorder + 2)
 
   98           { i0 = 1; i1 = i - totalorder - 1; }
 
   99         else if (i < 3u*totalorder + 1)
 
  100           { i0 = i - 2u*totalorder; i1 = 1; }
 
  101         else if (i < 4u*totalorder)
 
  102           { i0 = 0; i1 = i - 3u*totalorder + 1; }
 
  106             unsigned int basisnum = i - 4*totalorder;
 
  114         if     ((i>= 4                 && i<= 4+  totalorder-2u) && elem->
point(0) > elem->
point(1)) i0=totalorder+2-i0;
 
  115         else if ((i>= 4+  totalorder-1u && i<= 4+2*totalorder-3u) && elem->
point(1) > elem->
point(2)) i1=totalorder+2-i1;
 
  116         else if ((i>= 4+2*totalorder-2u && i<= 4+3*totalorder-4u) && elem->
point(3) > elem->
point(2)) i0=totalorder+2-i0;
 
  117         else if ((i>= 4+3*totalorder-3u && i<= 4+4*totalorder-5u) && elem->
point(0) > elem->
point(3)) i1=totalorder+2-i1;
 
  127         libmesh_assert_less (totalorder, 3);
 
  129         const Real xi  = p(0);
 
  130         const Real eta = p(1);
 
  132         libmesh_assert_less (i, 8);
 
  135         static const unsigned int i0[] = {0, 1, 1, 0, 2, 1, 2, 0, 2};
 
  136         static const unsigned int i1[] = {0, 0, 1, 1, 0, 2, 1, 2, 2};
 
  137         static const Real scal[] = {-0.25, -0.25, -0.25, -0.25, 0.5, 0.5, 0.5, 0.5};
 
  150       libmesh_assert_less (totalorder, 2);
 
  151       libmesh_fallthrough();
 
  161             libmesh_assert_less (i, 3);
 
  170                 libmesh_error_msg(
"Invalid shape function index i = " << i);
 
  179             libmesh_assert_less (i, 6);
 
  187               case 3: 
return 2.*x*r;
 
  188               case 4: 
return 2.*x*y;
 
  189               case 5: 
return 2.*r*y;
 
  192                 libmesh_error_msg(
"Invalid shape function index i = " << i);
 
  200             libmesh_assert_less (i, 10);
 
  202             unsigned int shape=i;
 
  205             if ((i==3||i==4) && elem->
point(0) > elem->
point(1)) shape=7-i;
 
  206             if ((i==5||i==6) && elem->
point(1) > elem->
point(2)) shape=11-i;
 
  207             if ((i==7||i==8) && elem->
point(0) > elem->
point(2)) shape=15-i;
 
  211               case 0: 
return r*r*r;
 
  212               case 1: 
return x*x*x;
 
  213               case 2: 
return y*y*y;
 
  215               case 3: 
return 3.*x*r*r;
 
  216               case 4: 
return 3.*x*x*r;
 
  218               case 5: 
return 3.*y*x*x;
 
  219               case 6: 
return 3.*y*y*x;
 
  221               case 7: 
return 3.*y*r*r;
 
  222               case 8: 
return 3.*y*y*r;
 
  224               case 9: 
return 6.*x*y*r;
 
  227                 libmesh_error_msg(
"Invalid shape function index shape = " << shape);
 
  235             unsigned int shape=i;
 
  237             libmesh_assert_less (i, 15);
 
  239             if ((i==3||i== 5) && elem->
point(0) > elem->
point(1)) shape=8-i;
 
  240             if ((i==6||i== 8) && elem->
point(1) > elem->
point(2)) shape=14-i;
 
  241             if ((i==9||i==11) && elem->
point(0) > elem->
point(2)) shape=20-i;
 
  247               case  0: 
return r*r*r*r;
 
  248               case  1: 
return x*x*x*x;
 
  249               case  2: 
return y*y*y*y;
 
  252               case  3: 
return 4.*x*r*r*r;
 
  253               case  4: 
return 6.*x*x*r*r;
 
  254               case  5: 
return 4.*x*x*x*r;
 
  256               case  6: 
return 4.*y*x*x*x;
 
  257               case  7: 
return 6.*y*y*x*x;
 
  258               case  8: 
return 4.*y*y*y*x;
 
  260               case  9: 
return 4.*y*r*r*r;
 
  261               case 10: 
return 6.*y*y*r*r;
 
  262               case 11: 
return 4.*y*y*y*r;
 
  265               case 12: 
return 12.*x*y*r*r;
 
  266               case 13: 
return 12.*x*x*y*r;
 
  267               case 14: 
return 12.*x*y*y*r;
 
  270                 libmesh_error_msg(
"Invalid shape function index shape = " << shape);
 
  278             unsigned int shape=i;
 
  280             libmesh_assert_less (i, 21);
 
  282             if ((i>= 3&&i<= 6) && elem->
point(0) > elem->
point(1)) shape=9-i;
 
  283             if ((i>= 7&&i<=10) && elem->
point(1) > elem->
point(2)) shape=17-i;
 
  284             if ((i>=11&&i<=14) && elem->
point(0) > elem->
point(2)) shape=25-i;
 
  289               case  0: 
return pow<5>(r);
 
  290               case  1: 
return pow<5>(x);
 
  291               case  2: 
return pow<5>(y);
 
  294               case  3: 
return  5.*x        *pow<4>(r);
 
  295               case  4: 
return 10.*pow<2>(x)*pow<3>(r);
 
  296               case  5: 
return 10.*pow<3>(x)*pow<2>(r);
 
  297               case  6: 
return  5.*pow<4>(x)*r;
 
  299               case  7: 
return  5.*y   *pow<4>(x);
 
  300               case  8: 
return 10.*pow<2>(y)*pow<3>(x);
 
  301               case  9: 
return 10.*pow<3>(y)*pow<2>(x);
 
  302               case 10: 
return  5.*pow<4>(y)*x;
 
  304               case 11: 
return  5.*y   *pow<4>(r);
 
  305               case 12: 
return 10.*pow<2>(y)*pow<3>(r);
 
  306               case 13: 
return 10.*pow<3>(y)*pow<2>(r);
 
  307               case 14: 
return  5.*pow<4>(y)*r;
 
  310               case 15: 
return 20.*x*y*pow<3>(r);
 
  311               case 16: 
return 30.*x*pow<2>(y)*pow<2>(r);
 
  312               case 17: 
return 30.*pow<2>(x)*y*pow<2>(r);
 
  313               case 18: 
return 20.*x*pow<3>(y)*r;
 
  314               case 19: 
return 20.*pow<3>(x)*y*r;
 
  315               case 20: 
return 30.*pow<2>(x)*pow<2>(y)*r;
 
  318                 libmesh_error_msg(
"Invalid shape function index shape = " << shape);
 
  326             unsigned int shape=i;
 
  328             libmesh_assert_less (i, 28);
 
  330             if ((i>= 3&&i<= 7) && elem->
point(0) > elem->
point(1)) shape=10-i;
 
  331             if ((i>= 8&&i<=12) && elem->
point(1) > elem->
point(2)) shape=20-i;
 
  332             if ((i>=13&&i<=17) && elem->
point(0) > elem->
point(2)) shape=30-i;
 
  337               case  0: 
return pow<6>(r);
 
  338               case  1: 
return pow<6>(x);
 
  339               case  2: 
return pow<6>(y);
 
  342               case  3: 
return  6.*x        *pow<5>(r);
 
  343               case  4: 
return 15.*pow<2>(x)*pow<4>(r);
 
  344               case  5: 
return 20.*pow<3>(x)*pow<3>(r);
 
  345               case  6: 
return 15.*pow<4>(x)*pow<2>(r);
 
  346               case  7: 
return  6.*pow<5>(x)*r;
 
  348               case  8: 
return  6.*y        *pow<5>(x);
 
  349               case  9: 
return 15.*pow<2>(y)*pow<4>(x);
 
  350               case 10: 
return 20.*pow<3>(y)*pow<3>(x);
 
  351               case 11: 
return 15.*pow<4>(y)*pow<2>(x);
 
  352               case 12: 
return  6.*pow<5>(y)*x;
 
  354               case 13: 
return  6.*y        *pow<5>(r);
 
  355               case 14: 
return 15.*pow<2>(y)*pow<4>(r);
 
  356               case 15: 
return 20.*pow<3>(y)*pow<3>(r);
 
  357               case 16: 
return 15.*pow<4>(y)*pow<2>(r);
 
  358               case 17: 
return  6.*pow<5>(y)*r;
 
  361               case 18: 
return 30.*x*y*pow<4>(r);
 
  362               case 19: 
return 60.*x*pow<2>(y)*pow<3>(r);
 
  363               case 20: 
return 60.*  pow<2>(x)*y*pow<3>(r);
 
  364               case 21: 
return 60.*x*pow<3>(y)*pow<2>(r);
 
  365               case 22: 
return 60.*pow<3>(x)*y*pow<2>(r);
 
  366               case 23: 
return 90.*pow<2>(x)*pow<2>(y)*pow<2>(r);
 
  367               case 24: 
return 30.*x*pow<4>(y)*r;
 
  368               case 25: 
return 60.*pow<2>(x)*pow<3>(y)*r;
 
  369               case 26: 
return 60.*pow<3>(x)*pow<2>(y)*r;
 
  370               case 27: 
return 30.*pow<4>(x)*y*r;
 
  373                 libmesh_error_msg(
"Invalid shape function index shape = " << shape);
 
  377           libmesh_error_msg(
"Invalid totalorder = " << totalorder);
 
  381       libmesh_error_msg(
"ERROR: Unsupported element type = " << type);
 
  394   libmesh_error_msg(
"Bernstein polynomials require the element type \nbecause edge orientation is needed.");
 
  403                                   const unsigned int i,
 
  404                                   const unsigned int j,
 
  406                                   const bool add_p_level)
 
  412   const Order totalorder =
 
  413     static_cast<Order>(order + add_p_level * elem->
p_level());
 
  422         const Real xi  = p(0);
 
  423         const Real eta = p(1);
 
  425         libmesh_assert_less (i, (totalorder+1u)*(totalorder+1u));
 
  441         else if (i < totalorder + 3u)
 
  442           { i0 = i - 2; i1 = 0; }
 
  443         else if (i < 2u*totalorder + 2)
 
  444           { i0 = 1; i1 = i - totalorder - 1; }
 
  445         else if (i < 3u*totalorder + 1)
 
  446           { i0 = i - 2u*totalorder; i1 = 1; }
 
  447         else if (i < 4u*totalorder)
 
  448           { i0 = 0; i1 = i - 3u*totalorder + 1; }
 
  452             unsigned int basisnum = i - 4*totalorder;
 
  460         if      ((i>= 4                 && i<= 4+  totalorder-2u) && elem->
point(0) > elem->
point(1)) i0=totalorder+2-i0;
 
  461         else if ((i>= 4+  totalorder-1u && i<= 4+2*totalorder-3u) && elem->
point(1) > elem->
point(2)) i1=totalorder+2-i1;
 
  462         else if ((i>= 4+2*totalorder-2u && i<= 4+3*totalorder-4u) && elem->
point(3) > elem->
point(2)) i0=totalorder+2-i0;
 
  463         else if ((i>= 4+3*totalorder-3u && i<= 4+4*totalorder-5u) && elem->
point(0) > elem->
point(3)) i1=totalorder+2-i1;
 
  478             libmesh_error_msg(
"Invalid shape function derivative j = " << j);
 
  487         libmesh_assert_less (totalorder, 3);
 
  489         const Real xi  = p(0);
 
  490         const Real eta = p(1);
 
  492         libmesh_assert_less (i, 8);
 
  495         static const unsigned int i0[] = {0, 1, 1, 0, 2, 1, 2, 0, 2};
 
  496         static const unsigned int i1[] = {0, 0, 1, 1, 0, 2, 1, 2, 2};
 
  497         static const Real scal[] = {-0.25, -0.25, -0.25, -0.25, 0.5, 0.5, 0.5, 0.5};
 
  517             libmesh_error_msg(
"Invalid shape function derivative j = " << j);
 
  523       libmesh_assert_less (totalorder, 2);
 
  524       libmesh_fallthrough();
 
  529         const Real eps = 1.e-6;
 
  536               const Point pp(p(0)+eps, p(1));
 
  537               const Point pm(p(0)-eps, p(1));
 
  546               const Point pp(p(0), p(1)+eps);
 
  547               const Point pm(p(0), p(1)-eps);
 
  555             libmesh_error_msg(
"Invalid shape function derivative j = " << j);
 
  560       libmesh_error_msg(
"ERROR: Unsupported element type = " << type);
 
  566 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 
  574   static bool warning_given = 
false;
 
  577     libMesh::err << 
"Second derivatives for Bernstein elements " 
  578                  << 
"are not yet implemented!" 
  581   warning_given = 
true;
 
  595   static bool warning_given = 
false;
 
  598     libMesh::err << 
"Second derivatives for Bernstein elements " 
  599                  << 
"are not yet implemented!" 
  602   warning_given = 
true;
 
  611 #endif// LIBMESH_ENABLE_HIGHER_ORDER_SHAPES