18 #include "libmesh/fe_interface.h" 
   19 #include "libmesh/h1_fe_transformation.h" 
   20 #include "libmesh/tensor_value.h" 
   21 #include "libmesh/elem.h" 
   22 #include "libmesh/int_range.h" 
   26 template<
typename OutputShape>
 
   34 template<
typename OutputShape>
 
   42 template<
typename OutputShape>
 
   46 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 
   51 template<
typename OutputShape>
 
   53                                                const Elem * 
const elem,
 
   54                                                const std::vector<Point> & qp,
 
   56                                                std::vector<std::vector<OutputShape>> & phi )
 const 
   64             libmesh_assert_equal_to ( qp.size(), phi[i].size() );
 
   66               FEInterface::shape<OutputShape>(0, fe.
get_fe_type(), elem, i, qp[p], phi[i][p]);
 
   74             libmesh_assert_equal_to ( qp.size(), phi[i].size() );
 
   76               FEInterface::shape<OutputShape>(1, fe.
get_fe_type(), elem, i, qp[p], phi[i][p]);
 
   84             libmesh_assert_equal_to ( qp.size(), phi[i].size() );
 
   86               FEInterface::shape<OutputShape>(2, fe.
get_fe_type(), elem, i, qp[p], phi[i][p]);
 
   94             libmesh_assert_equal_to ( qp.size(), phi[i].size() );
 
   96               FEInterface::shape<OutputShape>(3, fe.
get_fe_type(), elem, i, qp[p], phi[i][p]);
 
  102       libmesh_error_msg(
"Invalid dim = " << 
dim);
 
  107 template<
typename OutputShape>
 
  110                                                 const std::vector<Point> &,
 
  113                                                 std::vector<std::vector<OutputShape>> & dphidx,
 
  114                                                 std::vector<std::vector<OutputShape>> & dphidy,
 
  115                                                 std::vector<std::vector<OutputShape>> & dphidz )
 const 
  129         const std::vector<std::vector<OutputShape>> & dphidxi = fe.
get_dphidxi();
 
  143               dphi[i][p].slice(0) = dphidx[i][p] = dphidxi[i][p]*dxidx_map[p];
 
  146               dphi[i][p].slice(1)  = dphidy[i][p] = dphidxi[i][p]*dxidy_map[p];
 
  149               dphi[i][p].slice(2) = dphidz[i][p] = dphidxi[i][p]*dxidz_map[p];
 
  158         const std::vector<std::vector<OutputShape>> & dphidxi = fe.
get_dphidxi();
 
  159         const std::vector<std::vector<OutputShape>> & dphideta = fe.
get_dphideta();
 
  177               dphi[i][p].slice(0) = dphidx[i][p] = (dphidxi[i][p]*dxidx_map[p] +
 
  178                                                     dphideta[i][p]*detadx_map[p]);
 
  181               dphi[i][p].slice(1) = dphidy[i][p] = (dphidxi[i][p]*dxidy_map[p] +
 
  182                                                     dphideta[i][p]*detady_map[p]);
 
  186               dphi[i][p].slice(2) = dphidz[i][p] = (dphidxi[i][p]*dxidz_map[p] +
 
  187                                                     dphideta[i][p]*detadz_map[p]);
 
  196         const std::vector<std::vector<OutputShape>> & dphidxi = fe.
get_dphidxi();
 
  197         const std::vector<std::vector<OutputShape>> & dphideta = fe.
get_dphideta();
 
  198         const std::vector<std::vector<OutputShape>> & dphidzeta = fe.
get_dphidzeta();
 
  216               dphi[i][p].slice(0) = dphidx[i][p] = (dphidxi[i][p]*dxidx_map[p] +
 
  217                                                     dphideta[i][p]*detadx_map[p] +
 
  218                                                     dphidzeta[i][p]*dzetadx_map[p]);
 
  221               dphi[i][p].slice(1) = dphidy[i][p] = (dphidxi[i][p]*dxidy_map[p] +
 
  222                                                     dphideta[i][p]*detady_map[p] +
 
  223                                                     dphidzeta[i][p]*dzetady_map[p]);
 
  226               dphi[i][p].slice(2) = dphidz[i][p] = (dphidxi[i][p]*dxidz_map[p] +
 
  227                                                     dphideta[i][p]*detadz_map[p] +
 
  228                                                     dphidzeta[i][p]*dzetadz_map[p]);
 
  234       libmesh_error_msg(
"Invalid dim = " << 
dim);
 
  238 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 
  239 template<
typename OutputShape>
 
  241                                                 const std::vector<Point> &,
 
  244                                                 std::vector<std::vector<OutputShape>> & d2phidx2,
 
  245                                                 std::vector<std::vector<OutputShape>> & d2phidxdy,
 
  246                                                 std::vector<std::vector<OutputShape>> & d2phidxdz,
 
  247                                                 std::vector<std::vector<OutputShape>> & d2phidy2,
 
  248                                                 std::vector<std::vector<OutputShape>> & d2phidydz,
 
  249                                                 std::vector<std::vector<OutputShape>> & d2phidz2)
 const 
  264         const std::vector<std::vector<OutputShape>> & d2phidxi2 = fe.
get_d2phidxi2();
 
  275         const std::vector<std::vector<OutputShape>> & dphidxi = fe.
get_dphidxi();
 
  284               d2phi[i][p].slice(0).slice(0) = d2phidx2[i][p] =
 
  285                 d2phidxi2[i][p]*dxidx_map[p]*dxidx_map[p] + 
 
  286                 d2xidxyz2[p][0]*dphidxi[i][p];              
 
  290               d2phi[i][p].slice(0).slice(1) = d2phi[i][p].slice(1).slice(0) = d2phidxdy[i][p] =
 
  291                 d2phidxi2[i][p]*dxidx_map[p]*dxidy_map[p] + 
 
  292                 d2xidxyz2[p][1]*dphidxi[i][p];              
 
  297               d2phi[i][p].slice(0).slice(2) = d2phi[i][p].slice(2).slice(0) = d2phidxdz[i][p] =
 
  298                 d2phidxi2[i][p]*dxidx_map[p]*dxidz_map[p] + 
 
  299                 d2xidxyz2[p][2]*dphidxi[i][p];              
 
  305               d2phi[i][p].slice(1).slice(1) = d2phidy2[i][p] =
 
  306                 d2phidxi2[i][p]*dxidy_map[p]*dxidy_map[p] + 
 
  307                 d2xidxyz2[p][3]*dphidxi[i][p];              
 
  312               d2phi[i][p].slice(1).slice(2) = d2phi[i][p].slice(2).slice(1) = d2phidydz[i][p] =
 
  313                 d2phidxi2[i][p]*dxidy_map[p]*dxidz_map[p] + 
 
  314                 d2xidxyz2[p][4]*dphidxi[i][p];              
 
  317               d2phi[i][p].slice(2).slice(2) = d2phidz2[i][p] =
 
  318                 d2phidxi2[i][p]*dxidz_map[p]*dxidz_map[p] + 
 
  319                 d2xidxyz2[p][5]*dphidxi[i][p];              
 
  327         const std::vector<std::vector<OutputShape>> & d2phidxi2 = fe.
get_d2phidxi2();
 
  328         const std::vector<std::vector<OutputShape>> & d2phidxideta = fe.
get_d2phidxideta();
 
  329         const std::vector<std::vector<OutputShape>> & d2phideta2 = fe.
get_d2phideta2();
 
  344         const std::vector<std::vector<OutputShape>> & dphidxi = fe.
get_dphidxi();
 
  345         const std::vector<std::vector<OutputShape>> & dphideta = fe.
get_dphideta();
 
  355               d2phi[i][p].slice(0).slice(0) = d2phidx2[i][p] =
 
  356                 d2phidxi2[i][p]*dxidx_map[p]*dxidx_map[p] +       
 
  357                 d2phideta2[i][p]*detadx_map[p]*detadx_map[p] +    
 
  358                 2*d2phidxideta[i][p]*dxidx_map[p]*detadx_map[p] + 
 
  359                 d2xidxyz2[p][0]*dphidxi[i][p] +                   
 
  360                 d2etadxyz2[p][0]*dphideta[i][p];                  
 
  363               d2phi[i][p].slice(0).slice(1) = d2phi[i][p].slice(1).slice(0) = d2phidxdy[i][p] =
 
  364                 d2phidxi2[i][p]*dxidx_map[p]*dxidy_map[p] +                                    
 
  365                 d2phideta2[i][p]*detadx_map[p]*detady_map[p] +                                 
 
  366                 d2phidxideta[i][p]*(dxidx_map[p]*detady_map[p] + detadx_map[p]*dxidy_map[p]) + 
 
  367                 d2xidxyz2[p][1]*dphidxi[i][p] +                                                
 
  368                 d2etadxyz2[p][1]*dphideta[i][p];                                               
 
  372               d2phi[i][p].slice(0).slice(2) = d2phi[i][p].slice(2).slice(0) = d2phidxdz[i][p] =
 
  373                 d2phidxi2[i][p]*dxidx_map[p]*dxidz_map[p] +                                    
 
  374                 d2phideta2[i][p]*detadx_map[p]*detadz_map[p] +                                 
 
  375                 d2phidxideta[i][p]*(dxidx_map[p]*detadz_map[p] + detadx_map[p]*dxidz_map[p]) + 
 
  376                 d2xidxyz2[p][2]*dphidxi[i][p] +                                                
 
  377                 d2etadxyz2[p][2]*dphideta[i][p];                                               
 
  381               d2phi[i][p].slice(1).slice(1) = d2phidy2[i][p] =
 
  382                 d2phidxi2[i][p]*dxidy_map[p]*dxidy_map[p] +       
 
  383                 d2phideta2[i][p]*detady_map[p]*detady_map[p] +    
 
  384                 2*d2phidxideta[i][p]*dxidy_map[p]*detady_map[p] + 
 
  385                 d2xidxyz2[p][3]*dphidxi[i][p] +                   
 
  386                 d2etadxyz2[p][3]*dphideta[i][p];                  
 
  390               d2phi[i][p].slice(1).slice(2) = d2phi[i][p].slice(2).slice(1) = d2phidydz[i][p] =
 
  391                 d2phidxi2[i][p]*dxidy_map[p]*dxidz_map[p] +                                    
 
  392                 d2phideta2[i][p]*detady_map[p]*detadz_map[p] +                                 
 
  393                 d2phidxideta[i][p]*(dxidy_map[p]*detadz_map[p] + detady_map[p]*dxidz_map[p]) + 
 
  394                 d2xidxyz2[p][4]*dphidxi[i][p] +                                                
 
  395                 d2etadxyz2[p][4]*dphideta[i][p];                                               
 
  398               d2phi[i][p].slice(2).slice(2) = d2phidz2[i][p] =
 
  399                 d2phidxi2[i][p]*dxidz_map[p]*dxidz_map[p] +       
 
  400                 d2phideta2[i][p]*detadz_map[p]*detadz_map[p] +    
 
  401                 2*d2phidxideta[i][p]*dxidz_map[p]*detadz_map[p] + 
 
  402                 d2xidxyz2[p][5]*dphidxi[i][p] +                   
 
  403                 d2etadxyz2[p][5]*dphideta[i][p];                  
 
  412         const std::vector<std::vector<OutputShape>> & d2phidxi2 = fe.
get_d2phidxi2();
 
  413         const std::vector<std::vector<OutputShape>> & d2phidxideta = fe.
get_d2phidxideta();
 
  414         const std::vector<std::vector<OutputShape>> & d2phideta2 = fe.
get_d2phideta2();
 
  415         const std::vector<std::vector<OutputShape>> & d2phidxidzeta = fe.
get_d2phidxidzeta();
 
  416         const std::vector<std::vector<OutputShape>> & d2phidetadzeta = fe.
get_d2phidetadzeta();
 
  417         const std::vector<std::vector<OutputShape>> & d2phidzeta2 = fe.
get_d2phidzeta2();
 
  432         const std::vector<std::vector<OutputShape>> & dphidxi = fe.
get_dphidxi();
 
  433         const std::vector<std::vector<OutputShape>> & dphideta = fe.
get_dphideta();
 
  434         const std::vector<std::vector<OutputShape>> & dphidzeta = fe.
get_dphidzeta();
 
  445               d2phi[i][p].slice(0).slice(0) = d2phidx2[i][p] =
 
  446                 d2phidxi2[i][p]*dxidx_map[p]*dxidx_map[p] +           
 
  447                 d2phideta2[i][p]*detadx_map[p]*detadx_map[p] +        
 
  448                 d2phidzeta2[i][p]*dzetadx_map[p]*dzetadx_map[p] +     
 
  449                 2*d2phidxideta[i][p]*dxidx_map[p]*detadx_map[p] +     
 
  450                 2*d2phidxidzeta[i][p]*dxidx_map[p]*dzetadx_map[p] +   
 
  451                 2*d2phidetadzeta[i][p]*detadx_map[p]*dzetadx_map[p] + 
 
  452                 d2xidxyz2[p][0]*dphidxi[i][p] +                       
 
  453                 d2etadxyz2[p][0]*dphideta[i][p] +                     
 
  454                 d2zetadxyz2[p][0]*dphidzeta[i][p];                    
 
  457               d2phi[i][p].slice(0).slice(1) = d2phi[i][p].slice(1).slice(0) = d2phidxdy[i][p] =
 
  458                 d2phidxi2[i][p]*dxidx_map[p]*dxidy_map[p] +                                          
 
  459                 d2phideta2[i][p]*detadx_map[p]*detady_map[p] +                                       
 
  460                 d2phidzeta2[i][p]*dzetadx_map[p]*dzetady_map[p] +                                    
 
  461                 d2phidxideta[i][p]*(dxidx_map[p]*detady_map[p] + detadx_map[p]*dxidy_map[p]) +       
 
  462                 d2phidxidzeta[i][p]*(dxidx_map[p]*dzetady_map[p] + dzetadx_map[p]*dxidy_map[p]) +    
 
  463                 d2phidetadzeta[i][p]*(detadx_map[p]*dzetady_map[p] + dzetadx_map[p]*detady_map[p]) + 
 
  464                 d2xidxyz2[p][1]*dphidxi[i][p] +                                                      
 
  465                 d2etadxyz2[p][1]*dphideta[i][p] +                                                    
 
  466                 d2zetadxyz2[p][1]*dphidzeta[i][p];                                                   
 
  469               d2phi[i][p].slice(0).slice(2) = d2phi[i][p].slice(2).slice(0) = d2phidy2[i][p] =
 
  470                 d2phidxi2[i][p]*dxidx_map[p]*dxidz_map[p] +                                          
 
  471                 d2phideta2[i][p]*detadx_map[p]*detadz_map[p] +                                       
 
  472                 d2phidzeta2[i][p]*dzetadx_map[p]*dzetadz_map[p] +                                    
 
  473                 d2phidxideta[i][p]*(dxidx_map[p]*detadz_map[p] + detadx_map[p]*dxidz_map[p]) +       
 
  474                 d2phidxidzeta[i][p]*(dxidx_map[p]*dzetadz_map[p] + dzetadx_map[p]*dxidz_map[p]) +    
 
  475                 d2phidetadzeta[i][p]*(detadx_map[p]*dzetadz_map[p] + dzetadx_map[p]*detadz_map[p]) + 
 
  476                 d2xidxyz2[p][2]*dphidxi[i][p] +                                                      
 
  477                 d2etadxyz2[p][2]*dphideta[i][p] +                                                    
 
  478                 d2zetadxyz2[p][2]*dphidzeta[i][p];                                                   
 
  481               d2phi[i][p].slice(1).slice(1) = d2phidxdz[i][p] =
 
  482                 d2phidxi2[i][p]*dxidy_map[p]*dxidy_map[p] +           
 
  483                 d2phideta2[i][p]*detady_map[p]*detady_map[p] +        
 
  484                 d2phidzeta2[i][p]*dzetady_map[p]*dzetady_map[p] +     
 
  485                 2*d2phidxideta[i][p]*dxidy_map[p]*detady_map[p] +     
 
  486                 2*d2phidxidzeta[i][p]*dxidy_map[p]*dzetady_map[p] +   
 
  487                 2*d2phidetadzeta[i][p]*detady_map[p]*dzetady_map[p] + 
 
  488                 d2xidxyz2[p][3]*dphidxi[i][p] +                       
 
  489                 d2etadxyz2[p][3]*dphideta[i][p] +                     
 
  490                 d2zetadxyz2[p][3]*dphidzeta[i][p];                    
 
  493               d2phi[i][p].slice(1).slice(2) = d2phi[i][p].slice(2).slice(1) = d2phidydz[i][p] =
 
  494                 d2phidxi2[i][p]*dxidy_map[p]*dxidz_map[p] +                                          
 
  495                 d2phideta2[i][p]*detady_map[p]*detadz_map[p] +                                       
 
  496                 d2phidzeta2[i][p]*dzetady_map[p]*dzetadz_map[p] +                                    
 
  497                 d2phidxideta[i][p]*(dxidy_map[p]*detadz_map[p] + detady_map[p]*dxidz_map[p]) +       
 
  498                 d2phidxidzeta[i][p]*(dxidy_map[p]*dzetadz_map[p] + dzetady_map[p]*dxidz_map[p]) +    
 
  499                 d2phidetadzeta[i][p]*(detady_map[p]*dzetadz_map[p] + dzetady_map[p]*detadz_map[p]) + 
 
  500                 d2xidxyz2[p][4]*dphidxi[i][p] +                                                      
 
  501                 d2etadxyz2[p][4]*dphideta[i][p] +                                                    
 
  502                 d2zetadxyz2[p][4]*dphidzeta[i][p];                                                   
 
  505               d2phi[i][p].slice(2).slice(2) = d2phidz2[i][p] =
 
  506                 d2phidxi2[i][p]*dxidz_map[p]*dxidz_map[p] +           
 
  507                 d2phideta2[i][p]*detadz_map[p]*detadz_map[p] +        
 
  508                 d2phidzeta2[i][p]*dzetadz_map[p]*dzetadz_map[p] +     
 
  509                 2*d2phidxideta[i][p]*dxidz_map[p]*detadz_map[p] +     
 
  510                 2*d2phidxidzeta[i][p]*dxidz_map[p]*dzetadz_map[p] +   
 
  511                 2*d2phidetadzeta[i][p]*detadz_map[p]*dzetadz_map[p] + 
 
  512                 d2xidxyz2[p][5]*dphidxi[i][p] +                       
 
  513                 d2etadxyz2[p][5]*dphideta[i][p] +                     
 
  514                 d2zetadxyz2[p][5]*dphidzeta[i][p];                    
 
  521       libmesh_error_msg(
"Invalid dim = " << 
dim);
 
  524 #endif //LIBMESH_ENABLE_SECOND_DERIVATIVES 
  529                                         const std::vector<Point> &,
 
  531                                         std::vector<std::vector<Real>> &)
 const 
  533   libmesh_error_msg(
"Computing the curl of a shape function only \nmakes sense for vector-valued elements.");
 
  539                                                 const std::vector<Point> &,
 
  541                                                 std::vector<std::vector<RealGradient>> & curl_phi )
 const 
  547       libmesh_error_msg(
"The curl only make sense in 2D and 3D");
 
  551         const std::vector<std::vector<RealGradient>> & dphidxi = fe.
get_dphidxi();
 
  552         const std::vector<std::vector<RealGradient>> & dphideta = fe.
get_dphideta();
 
  572               Real dphiy_dx = (dphidxi[i][p].slice(1))*dxidx_map[p]
 
  573                 + (dphideta[i][p].slice(1))*detadx_map[p];
 
  575               Real dphix_dy = (dphidxi[i][p].slice(0))*dxidy_map[p]
 
  576                 + (dphideta[i][p].slice(0))*detady_map[p];
 
  578               curl_phi[i][p].slice(2) = dphiy_dx - dphix_dy;
 
  581               Real dphiy_dz = (dphidxi[i][p].slice(1))*dxidz_map[p]
 
  582                 + (dphideta[i][p].slice(1))*detadz_map[p];
 
  584               Real dphix_dz = (dphidxi[i][p].slice(0))*dxidz_map[p]
 
  585                 + (dphideta[i][p].slice(0))*detadz_map[p];
 
  587               curl_phi[i][p].slice(0) = -dphiy_dz;
 
  588               curl_phi[i][p].slice(1) = dphix_dz;
 
  596         const std::vector<std::vector<RealGradient>> & dphidxi = fe.
get_dphidxi();
 
  597         const std::vector<std::vector<RealGradient>> & dphideta = fe.
get_dphideta();
 
  598         const std::vector<std::vector<RealGradient>> & dphidzeta = fe.
get_dphidzeta();
 
  616               Real dphiz_dy = (dphidxi[i][p].slice(2))*dxidy_map[p]
 
  617                 + (dphideta[i][p].slice(2))*detady_map[p]
 
  618                 + (dphidzeta[i][p].slice(2))*dzetady_map[p];
 
  620               Real dphiy_dz = (dphidxi[i][p].slice(1))*dxidz_map[p]
 
  621                 + (dphideta[i][p].slice(1))*detadz_map[p]
 
  622                 + (dphidzeta[i][p].slice(1))*dzetadz_map[p];
 
  624               Real dphix_dz = (dphidxi[i][p].slice(0))*dxidz_map[p]
 
  625                 + (dphideta[i][p].slice(0))*detadz_map[p]
 
  626                 + (dphidzeta[i][p].slice(0))*dzetadz_map[p];
 
  628               Real dphiz_dx = (dphidxi[i][p].slice(2))*dxidx_map[p]
 
  629                 + (dphideta[i][p].slice(2))*detadx_map[p]
 
  630                 + (dphidzeta[i][p].slice(2))*dzetadx_map[p];
 
  632               Real dphiy_dx = (dphidxi[i][p].slice(1))*dxidx_map[p]
 
  633                 + (dphideta[i][p].slice(1))*detadx_map[p]
 
  634                 + (dphidzeta[i][p].slice(1))*dzetadx_map[p];
 
  636               Real dphix_dy = (dphidxi[i][p].slice(0))*dxidy_map[p]
 
  637                 + (dphideta[i][p].slice(0))*detady_map[p]
 
  638                 + (dphidzeta[i][p].slice(0))*dzetady_map[p];
 
  640               curl_phi[i][p].slice(0) = dphiz_dy - dphiy_dz;
 
  642               curl_phi[i][p].slice(1) = dphix_dz - dphiz_dx;
 
  644               curl_phi[i][p].slice(2) = dphiy_dx - dphix_dy;
 
  650       libmesh_error_msg(
"Invalid dim = " << 
dim);
 
  658                                        const std::vector<Point> &,
 
  662   libmesh_error_msg(
"Computing the divergence of a shape function only \nmakes sense for vector-valued elements.");
 
  669                                                const std::vector<Point> &,
 
  677       libmesh_error_msg(
"The divergence only make sense in 2D and 3D");
 
  681         const std::vector<std::vector<RealGradient>> & dphidxi = fe.
get_dphidxi();
 
  682         const std::vector<std::vector<RealGradient>> & dphideta = fe.
get_dphideta();
 
  694               Real dphix_dx = (dphidxi[i][p].slice(0))*dxidx_map[p]
 
  695                 + (dphideta[i][p].slice(0))*detadx_map[p];
 
  697               Real dphiy_dy = (dphidxi[i][p].slice(1))*dxidy_map[p]
 
  698                 + (dphideta[i][p].slice(1))*detady_map[p];
 
  700               div_phi[i][p] = dphix_dx + dphiy_dy;
 
  706         const std::vector<std::vector<RealGradient>> & dphidxi = fe.
get_dphidxi();
 
  707         const std::vector<std::vector<RealGradient>> & dphideta = fe.
get_dphideta();
 
  708         const std::vector<std::vector<RealGradient>> & dphidzeta = fe.
get_dphidzeta();
 
  726               Real dphix_dx = (dphidxi[i][p].slice(0))*dxidx_map[p]
 
  727                 + (dphideta[i][p].slice(0))*detadx_map[p]
 
  728                 + (dphidzeta[i][p].slice(0))*dzetadx_map[p];
 
  730               Real dphiy_dy = (dphidxi[i][p].slice(1))*dxidy_map[p]
 
  731                 + (dphideta[i][p].slice(1))*detady_map[p]
 
  732                 + (dphidzeta[i][p].slice(1))*dzetady_map[p];
 
  734               Real dphiz_dz = (dphidxi[i][p].slice(2))*dxidz_map[p]
 
  735                 + (dphideta[i][p].slice(2))*detadz_map[p]
 
  736                 + (dphidzeta[i][p].slice(2))*dzetadz_map[p];
 
  738               div_phi[i][p] = dphix_dx + dphiy_dy + dphiz_dz;
 
  745       libmesh_error_msg(
"Invalid dim = " << 
dim);