Line data    Source code

       1             : // The libMesh Finite Element Library.
       2             : // Copyright (C) 2002-2026 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             : #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"
      23             : 
      24             : namespace libMesh
      25             : {
      26             : template<typename OutputShape>
      27    95469641 : void H1FETransformation<OutputShape>::init_map_phi(const FEGenericBase<OutputShape> & /* fe */ ) const
      28             : {
      29             :   // Nothing needed
      30    95469641 : }
      31             : 
      32             : 
      33             : 
      34             : template<typename OutputShape>
      35    66281041 : void H1FETransformation<OutputShape>::init_map_dphi(const FEGenericBase<OutputShape> & fe) const
      36             : {
      37    19808122 :   fe.get_fe_map().get_dxidx();
      38    66281041 : }
      39             : 
      40             : 
      41             : 
      42             : template<typename OutputShape>
      43    14764532 : void H1FETransformation<OutputShape>::init_map_d2phi(const FEGenericBase<OutputShape> & fe) const
      44             : {
      45     4120951 :   fe.get_fe_map().get_dxidx();
      46             : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
      47     4120951 :   fe.get_fe_map().get_d2xidxyz2();
      48             : #endif
      49    14764532 : }
      50             : 
      51             : template<typename OutputShape>
      52    33527243 : void H1FETransformation<OutputShape>::map_phi( const unsigned int dim,
      53             :                                                const Elem * const elem,
      54             :                                                const std::vector<Point> & qp,
      55             :                                                const FEGenericBase<OutputShape> & fe,
      56             :                                                std::vector<std::vector<OutputShape>> & phi,
      57             :                                                const bool add_p_level) const
      58             : {
      59    43301838 :   FEInterface::all_shapes<OutputShape>(dim, fe.get_fe_type(), elem, qp, phi, add_p_level);
      60    33527243 : }
      61             : 
      62             : 
      63             : template<typename OutputShape>
      64    24584158 : void H1FETransformation<OutputShape>::map_dphi( const unsigned int dim,
      65             :                                                 const Elem * const,
      66             :                                                 const std::vector<Point> &,
      67             :                                                 const FEGenericBase<OutputShape> & fe,
      68             :                                                 std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputGradient>> & dphi,
      69             :                                                 std::vector<std::vector<OutputShape>> & dphidx,
      70             :                                                 std::vector<std::vector<OutputShape>> & dphidy,
      71             :                                                 std::vector<std::vector<OutputShape>> & dphidz ) const
      72             : {
      73    24584158 :   switch(dim)
      74             :     {
      75       18792 :     case 0: // No derivatives in 0D
      76             :       {
      77      134784 :         for (auto i : index_range(dphi))
      78      134784 :           for (auto p : index_range(dphi[i]))
      79       18792 :             dphi[i][p] = 0.;
      80       18792 :         break;
      81             :       }
      82             : 
      83       13458 :     case 1:
      84             :       {
      85       13458 :         const std::vector<std::vector<OutputShape>> & dphidxi = fe.get_dphidxi();
      86             : 
      87       13458 :         const std::vector<Real> & dxidx_map = fe.get_fe_map().get_dxidx();
      88             : #if LIBMESH_DIM>1
      89       13458 :         const std::vector<Real> & dxidy_map = fe.get_fe_map().get_dxidy();
      90             : #endif
      91             : #if LIBMESH_DIM>2
      92       13458 :         const std::vector<Real> & dxidz_map = fe.get_fe_map().get_dxidz();
      93             : #endif
      94             : 
      95      247611 :         for (auto i : index_range(dphi))
      96      700406 :           for (auto p : index_range(dphi[i]))
      97             :             {
      98             :               // dphi/dx    = (dphi/dxi)*(dxi/dx)
      99     1162572 :               dphi[i][p].slice(0) = dphidx[i][p] = dphidxi[i][p]*dxidx_map[p];
     100             : 
     101             : #if LIBMESH_DIM>1
     102      767130 :               dphi[i][p].slice(1)  = dphidy[i][p] = dphidxi[i][p]*dxidy_map[p];
     103             : #endif
     104             : #if LIBMESH_DIM>2
     105      898944 :               dphi[i][p].slice(2) = dphidz[i][p] = dphidxi[i][p]*dxidz_map[p];
     106             : #endif
     107             :             }
     108             : 
     109       13458 :         break;
     110             :       }
     111             : 
     112     6734999 :     case 2:
     113             :       {
     114     6734999 :         const std::vector<std::vector<OutputShape>> & dphidxi = fe.get_dphidxi();
     115     6734999 :         const std::vector<std::vector<OutputShape>> & dphideta = fe.get_dphideta();
     116             : 
     117     6734999 :         const std::vector<Real> & dxidx_map = fe.get_fe_map().get_dxidx();
     118     6734999 :         const std::vector<Real> & dxidy_map = fe.get_fe_map().get_dxidy();
     119             : #if LIBMESH_DIM > 2
     120     6734999 :         const std::vector<Real> & dxidz_map = fe.get_fe_map().get_dxidz();
     121             : #endif
     122             : 
     123     6734999 :         const std::vector<Real> & detadx_map = fe.get_fe_map().get_detadx();
     124     6734999 :         const std::vector<Real> & detady_map = fe.get_fe_map().get_detady();
     125             : #if LIBMESH_DIM > 2
     126     6734999 :         const std::vector<Real> & detadz_map = fe.get_fe_map().get_detadz();
     127             : #endif
     128             : 
     129   119803697 :         for (auto i : index_range(dphi))
     130   556665683 :           for (auto p : index_range(dphi[i]))
     131             :             {
     132             :               // dphi/dx    = (dphi/dxi)*(dxi/dx) + (dphi/deta)*(deta/dx)
     133  1148745635 :               dphi[i][p].slice(0) = dphidx[i][p] = (dphidxi[i][p]*dxidx_map[p] +
     134   869907707 :                                                     dphideta[i][p]*detadx_map[p]);
     135             : 
     136             :               // dphi/dy    = (dphi/dxi)*(dxi/dy) + (dphi/deta)*(deta/dy)
     137   736131539 :               dphi[i][p].slice(1) = dphidy[i][p] = (dphidxi[i][p]*dxidy_map[p] +
     138   594831643 :                                                     dphideta[i][p]*detady_map[p]);
     139             : 
     140             : #if LIBMESH_DIM > 2
     141             :               // dphi/dz    = (dphi/dxi)*(dxi/dz) + (dphi/deta)*(deta/dz)
     142   872932413 :               dphi[i][p].slice(2) = dphidz[i][p] = (dphidxi[i][p]*dxidz_map[p] +
     143   595568801 :                                                     dphideta[i][p]*detadz_map[p]);
     144             : #endif
     145             :             }
     146             : 
     147     6734999 :         break;
     148             :       }
     149             : 
     150      717455 :     case 3:
     151             :       {
     152      717455 :         const std::vector<std::vector<OutputShape>> & dphidxi = fe.get_dphidxi();
     153      717455 :         const std::vector<std::vector<OutputShape>> & dphideta = fe.get_dphideta();
     154      717455 :         const std::vector<std::vector<OutputShape>> & dphidzeta = fe.get_dphidzeta();
     155             : 
     156      717455 :         const std::vector<Real> & dxidx_map = fe.get_fe_map().get_dxidx();
     157      717455 :         const std::vector<Real> & dxidy_map = fe.get_fe_map().get_dxidy();
     158      717455 :         const std::vector<Real> & dxidz_map = fe.get_fe_map().get_dxidz();
     159             : 
     160      717455 :         const std::vector<Real> & detadx_map = fe.get_fe_map().get_detadx();
     161      717455 :         const std::vector<Real> & detady_map = fe.get_fe_map().get_detady();
     162      717455 :         const std::vector<Real> & detadz_map = fe.get_fe_map().get_detadz();
     163             : 
     164      717455 :         const std::vector<Real> & dzetadx_map = fe.get_fe_map().get_dzetadx();
     165      717455 :         const std::vector<Real> & dzetady_map = fe.get_fe_map().get_dzetady();
     166      717455 :         const std::vector<Real> & dzetadz_map = fe.get_fe_map().get_dzetadz();
     167             : 
     168    48055427 :         for (auto i : index_range(dphi))
     169   210382428 :           for (auto p : index_range(dphi[i]))
     170             :             {
     171             :               // dphi/dx    = (dphi/dxi)*(dxi/dx) + (dphi/deta)*(deta/dx) + (dphi/dzeta)*(dzeta/dx);
     172   386918180 :               dphi[i][p].slice(0) = dphidx[i][p] = (dphidxi[i][p]*dxidx_map[p] +
     173   251619233 :                                                     dphideta[i][p]*detadx_map[p] +
     174   291354018 :                                                     dphidzeta[i][p]*dzetadx_map[p]);
     175             : 
     176             :               // dphi/dy    = (dphi/dxi)*(dxi/dy) + (dphi/deta)*(deta/dy) + (dphi/dzeta)*(dzeta/dy);
     177   256984097 :               dphi[i][p].slice(1) = dphidy[i][p] = (dphidxi[i][p]*dxidy_map[p] +
     178   164996511 :                                                     dphideta[i][p]*detady_map[p] +
     179   204731296 :                                                     dphidzeta[i][p]*dzetady_map[p]);
     180             : 
     181             :               // dphi/dz    = (dphi/dxi)*(dxi/dz) + (dphi/deta)*(deta/dz) + (dphi/dzeta)*(dzeta/dz);
     182   298507170 :               dphi[i][p].slice(2) = dphidz[i][p] = (dphidxi[i][p]*dxidz_map[p] +
     183   164996511 :                                                     dphideta[i][p]*detadz_map[p] +
     184   206519584 :                                                     dphidzeta[i][p]*dzetadz_map[p]);
     185             :             }
     186      717455 :         break;
     187             :       }
     188             : 
     189           0 :     default:
     190           0 :       libmesh_error_msg("Invalid dim = " << dim);
     191             :     } // switch(dim)
     192    24584158 : }
     193             : 
     194             : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
     195             : template<typename OutputShape>
     196     2014772 : void H1FETransformation<OutputShape>::map_d2phi(const unsigned int dim,
     197             :                                                 const std::vector<Point> &,
     198             :                                                 const FEGenericBase<OutputShape> & fe,
     199             :                                                 std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputTensor>> & d2phi,
     200             :                                                 std::vector<std::vector<OutputShape>> & d2phidx2,
     201             :                                                 std::vector<std::vector<OutputShape>> & d2phidxdy,
     202             :                                                 std::vector<std::vector<OutputShape>> & d2phidxdz,
     203             :                                                 std::vector<std::vector<OutputShape>> & d2phidy2,
     204             :                                                 std::vector<std::vector<OutputShape>> & d2phidydz,
     205             :                                                 std::vector<std::vector<OutputShape>> & d2phidz2) const
     206             : {
     207     2014772 :   switch (dim)
     208             :     {
     209           0 :     case 0: // No derivatives in 0D
     210             :       {
     211           0 :         for (auto i : index_range(d2phi))
     212           0 :           for (auto p : index_range(d2phi[i]))
     213           0 :             d2phi[i][p] = 0.;
     214             : 
     215           0 :         break;
     216             :       }
     217             : 
     218        7660 :     case 1:
     219             :       {
     220        7660 :         const std::vector<std::vector<OutputShape>> & d2phidxi2 = fe.get_d2phidxi2();
     221             : 
     222        7660 :         const std::vector<Real> & dxidx_map = fe.get_fe_map().get_dxidx();
     223             : #if LIBMESH_DIM>1
     224        7660 :         const std::vector<Real> & dxidy_map = fe.get_fe_map().get_dxidy();
     225             : #endif
     226             : #if LIBMESH_DIM>2
     227        7660 :         const std::vector<Real> & dxidz_map = fe.get_fe_map().get_dxidz();
     228             : #endif
     229             : 
     230             :         // Shape function derivatives in reference space
     231        7660 :         const std::vector<std::vector<OutputShape>> & dphidxi = fe.get_dphidxi();
     232             : 
     233             :         // Inverse map second derivatives
     234        7660 :         const std::vector<std::vector<Real>> & d2xidxyz2 = fe.get_fe_map().get_d2xidxyz2();
     235             : 
     236      135055 :         for (auto i : index_range(d2phi))
     237      495047 :           for (auto p : index_range(d2phi[i]))
     238             :             {
     239             :               // phi_{x x}
     240      490340 :               d2phi[i][p].slice(0).slice(0) = d2phidx2[i][p] =
     241      593120 :                 d2phidxi2[i][p]*dxidx_map[p]*dxidx_map[p] + // (xi_x)^2 * phi_{xi xi}
     242      695900 :                 d2xidxyz2[p][0]*dphidxi[i][p];              // xi_{x x} * phi_{xi}
     243             : 
     244             : #if LIBMESH_DIM>1
     245             :               // phi_{x y}
     246      387560 :               d2phi[i][p].slice(0).slice(1) = d2phi[i][p].slice(1).slice(0) = d2phidxdy[i][p] =
     247      387560 :                 d2phidxi2[i][p]*dxidx_map[p]*dxidy_map[p] + // xi_x * xi_y * phi_{xi xi}
     248      387560 :                 d2xidxyz2[p][1]*dphidxi[i][p];              // xi_{x y} * phi_{xi}
     249             : #endif
     250             : 
     251             : #if LIBMESH_DIM>2
     252             :               // phi_{x z}
     253      387560 :               d2phi[i][p].slice(0).slice(2) = d2phi[i][p].slice(2).slice(0) = d2phidxdz[i][p] =
     254      387560 :                 d2phidxi2[i][p]*dxidx_map[p]*dxidz_map[p] + // xi_x * xi_z * phi_{xi xi}
     255      387560 :                 d2xidxyz2[p][2]*dphidxi[i][p];              // xi_{x z} * phi_{xi}
     256             : #endif
     257             : 
     258             : 
     259             : #if LIBMESH_DIM>1
     260             :               // phi_{y y}
     261      387560 :               d2phi[i][p].slice(1).slice(1) = d2phidy2[i][p] =
     262      490340 :                 d2phidxi2[i][p]*dxidy_map[p]*dxidy_map[p] + // (xi_y)^2 * phi_{xi xi}
     263      387560 :                 d2xidxyz2[p][3]*dphidxi[i][p];              // xi_{y y} * phi_{xi}
     264             : #endif
     265             : 
     266             : #if LIBMESH_DIM>2
     267             :               // phi_{y z}
     268      387560 :               d2phi[i][p].slice(1).slice(2) = d2phi[i][p].slice(2).slice(1) = d2phidydz[i][p] =
     269      490340 :                 d2phidxi2[i][p]*dxidy_map[p]*dxidz_map[p] + // xi_y * xi_z * phi_{xi xi}
     270      387560 :                 d2xidxyz2[p][4]*dphidxi[i][p];              // xi_{y z} * phi_{xi}
     271             : 
     272             :               // phi_{z z}
     273      490340 :               d2phi[i][p].slice(2).slice(2) = d2phidz2[i][p] =
     274      490340 :                 d2phidxi2[i][p]*dxidz_map[p]*dxidz_map[p] + // (xi_z)^2 * phi_{xi xi}
     275      387560 :                 d2xidxyz2[p][5]*dphidxi[i][p];              // xi_{z z} * phi_{xi}
     276             : #endif
     277             :             }
     278        7660 :         break;
     279             :       }
     280             : 
     281       81767 :     case 2:
     282             :       {
     283       81767 :         const std::vector<std::vector<OutputShape>> & d2phidxi2 = fe.get_d2phidxi2();
     284       81767 :         const std::vector<std::vector<OutputShape>> & d2phidxideta = fe.get_d2phidxideta();
     285       81767 :         const std::vector<std::vector<OutputShape>> & d2phideta2 = fe.get_d2phideta2();
     286             : 
     287       81767 :         const std::vector<Real> & dxidx_map = fe.get_fe_map().get_dxidx();
     288       81767 :         const std::vector<Real> & dxidy_map = fe.get_fe_map().get_dxidy();
     289             : #if LIBMESH_DIM > 2
     290       81767 :         const std::vector<Real> & dxidz_map = fe.get_fe_map().get_dxidz();
     291             : #endif
     292             : 
     293       81767 :         const std::vector<Real> & detadx_map = fe.get_fe_map().get_detadx();
     294       81767 :         const std::vector<Real> & detady_map = fe.get_fe_map().get_detady();
     295             : #if LIBMESH_DIM > 2
     296       81767 :         const std::vector<Real> & detadz_map = fe.get_fe_map().get_detadz();
     297             : #endif
     298             : 
     299             :         // Shape function derivatives in reference space
     300       81767 :         const std::vector<std::vector<OutputShape>> & dphidxi = fe.get_dphidxi();
     301       81767 :         const std::vector<std::vector<OutputShape>> & dphideta = fe.get_dphideta();
     302             : 
     303             :         // Inverse map second derivatives
     304       81767 :         const std::vector<std::vector<Real>> & d2xidxyz2 = fe.get_fe_map().get_d2xidxyz2();
     305       81767 :         const std::vector<std::vector<Real>> & d2etadxyz2 = fe.get_fe_map().get_d2etadxyz2();
     306             : 
     307     3077039 :         for (auto i : index_range(d2phi))
     308    17626277 :           for (auto p : index_range(d2phi[i]))
     309             :             {
     310             :               // phi_{x x}
     311    18633236 :               d2phi[i][p].slice(0).slice(0) = d2phidx2[i][p] =
     312    22409287 :                 d2phidxi2[i][p]*dxidx_map[p]*dxidx_map[p] +       // (xi_x)^2 * phi_{xi xi}
     313    22409287 :                 d2phideta2[i][p]*detadx_map[p]*detadx_map[p] +    // (eta_x)^2 * phi_{eta eta}
     314    18633236 :                 2*d2phidxideta[i][p]*dxidx_map[p]*detadx_map[p] + // 2 * xi_x * eta_x * phi_{xi eta}
     315    22409287 :                 d2xidxyz2[p][0]*dphidxi[i][p] +                   // xi_{x x} * phi_{xi}
     316    26185338 :                 d2etadxyz2[p][0]*dphideta[i][p];                  // eta_{x x} * phi_{eta}
     317             : 
     318             :               // phi_{x y}
     319    14857185 :               d2phi[i][p].slice(0).slice(1) = d2phi[i][p].slice(1).slice(0) = d2phidxdy[i][p] =
     320    14857185 :                 d2phidxi2[i][p]*dxidx_map[p]*dxidy_map[p] +                                    // xi_x * xi_y * phi_{xi xi}
     321    14857185 :                 d2phideta2[i][p]*detadx_map[p]*detady_map[p] +                                 // eta_x * eta_y * phi_{eta eta}
     322    18633236 :                 d2phidxideta[i][p]*(dxidx_map[p]*detady_map[p] + detadx_map[p]*dxidy_map[p]) + // (xi_x*eta_y + eta_x*xi_y) * phi_{xi eta}
     323    18633236 :                 d2xidxyz2[p][1]*dphidxi[i][p] +                                                // xi_{x y} * phi_{xi}
     324    14857185 :                 d2etadxyz2[p][1]*dphideta[i][p];                                               // eta_{x y} * phi_{eta}
     325             : 
     326             : #if LIBMESH_DIM > 2
     327             :               // phi_{x z}
     328    14857185 :               d2phi[i][p].slice(0).slice(2) = d2phi[i][p].slice(2).slice(0) = d2phidxdz[i][p] =
     329    14857185 :                 d2phidxi2[i][p]*dxidx_map[p]*dxidz_map[p] +                                    // xi_x * xi_z * phi_{xi xi}
     330    14857185 :                 d2phideta2[i][p]*detadx_map[p]*detadz_map[p] +                                 // eta_x * eta_z * phi_{eta eta}
     331    18633236 :                 d2phidxideta[i][p]*(dxidx_map[p]*detadz_map[p] + detadx_map[p]*dxidz_map[p]) + // (xi_x*eta_z + eta_x*xi_z) * phi_{xi eta}
     332    18633236 :                 d2xidxyz2[p][2]*dphidxi[i][p] +                                                // xi_{x z} * phi_{xi}
     333    14857185 :                 d2etadxyz2[p][2]*dphideta[i][p];                                               // eta_{x z} * phi_{eta}
     334             : #endif
     335             : 
     336             :               // phi_{y y}
     337    14857185 :               d2phi[i][p].slice(1).slice(1) = d2phidy2[i][p] =
     338    18633236 :                 d2phidxi2[i][p]*dxidy_map[p]*dxidy_map[p] +       // (xi_y)^2 * phi_{xi xi}
     339    18633236 :                 d2phideta2[i][p]*detady_map[p]*detady_map[p] +    // (eta_y)^2 * phi_{eta eta}
     340    18633236 :                 2*d2phidxideta[i][p]*dxidy_map[p]*detady_map[p] + // 2 * xi_y * eta_y * phi_{xi eta}
     341    18633236 :                 d2xidxyz2[p][3]*dphidxi[i][p] +                   // xi_{y y} * phi_{xi}
     342    14857185 :                 d2etadxyz2[p][3]*dphideta[i][p];                  // eta_{y y} * phi_{eta}
     343             : 
     344             : #if LIBMESH_DIM > 2
     345             :               // phi_{y z}
     346    14857185 :               d2phi[i][p].slice(1).slice(2) = d2phi[i][p].slice(2).slice(1) = d2phidydz[i][p] =
     347    18633236 :                 d2phidxi2[i][p]*dxidy_map[p]*dxidz_map[p] +                                    // xi_y * xi_z * phi_{xi xi}
     348    18633236 :                 d2phideta2[i][p]*detady_map[p]*detadz_map[p] +                                 // eta_y * eta_z * phi_{eta eta}
     349    18633236 :                 d2phidxideta[i][p]*(dxidy_map[p]*detadz_map[p] + detady_map[p]*dxidz_map[p]) + // (xi_y*eta_z + eta_y*xi_z) * phi_{xi eta}
     350    18633236 :                 d2xidxyz2[p][4]*dphidxi[i][p] +                                                // xi_{y z} * phi_{xi}
     351    14857185 :                 d2etadxyz2[p][4]*dphideta[i][p];                                               // eta_{y z} * phi_{eta}
     352             : 
     353             :               // phi_{z z}
     354    18633236 :               d2phi[i][p].slice(2).slice(2) = d2phidz2[i][p] =
     355    18633236 :                 d2phidxi2[i][p]*dxidz_map[p]*dxidz_map[p] +       // (xi_z)^2 * phi_{xi xi}
     356    18633236 :                 d2phideta2[i][p]*detadz_map[p]*detadz_map[p] +    // (eta_z)^2 * phi_{eta eta}
     357    18633236 :                 2*d2phidxideta[i][p]*dxidz_map[p]*detadz_map[p] + // 2 * xi_z * eta_z * phi_{xi eta}
     358    18633236 :                 d2xidxyz2[p][5]*dphidxi[i][p] +                   // xi_{z z} * phi_{xi}
     359    14857185 :                 d2etadxyz2[p][5]*dphideta[i][p];                  // eta_{z z} * phi_{eta}
     360             : #endif
     361             :             }
     362             : 
     363       81767 :         break;
     364             :       }
     365             : 
     366      435393 :     case 3:
     367             :       {
     368      435393 :         const std::vector<std::vector<OutputShape>> & d2phidxi2 = fe.get_d2phidxi2();
     369      435393 :         const std::vector<std::vector<OutputShape>> & d2phidxideta = fe.get_d2phidxideta();
     370      435393 :         const std::vector<std::vector<OutputShape>> & d2phideta2 = fe.get_d2phideta2();
     371      435393 :         const std::vector<std::vector<OutputShape>> & d2phidxidzeta = fe.get_d2phidxidzeta();
     372      435393 :         const std::vector<std::vector<OutputShape>> & d2phidetadzeta = fe.get_d2phidetadzeta();
     373      435393 :         const std::vector<std::vector<OutputShape>> & d2phidzeta2 = fe.get_d2phidzeta2();
     374             : 
     375      435393 :         const std::vector<Real> & dxidx_map = fe.get_fe_map().get_dxidx();
     376      435393 :         const std::vector<Real> & dxidy_map = fe.get_fe_map().get_dxidy();
     377      435393 :         const std::vector<Real> & dxidz_map = fe.get_fe_map().get_dxidz();
     378             : 
     379      435393 :         const std::vector<Real> & detadx_map = fe.get_fe_map().get_detadx();
     380      435393 :         const std::vector<Real> & detady_map = fe.get_fe_map().get_detady();
     381      435393 :         const std::vector<Real> & detadz_map = fe.get_fe_map().get_detadz();
     382             : 
     383      435393 :         const std::vector<Real> & dzetadx_map = fe.get_fe_map().get_dzetadx();
     384      435393 :         const std::vector<Real> & dzetady_map = fe.get_fe_map().get_dzetady();
     385      435393 :         const std::vector<Real> & dzetadz_map = fe.get_fe_map().get_dzetadz();
     386             : 
     387             :         // Shape function derivatives in reference space
     388      435393 :         const std::vector<std::vector<OutputShape>> & dphidxi = fe.get_dphidxi();
     389      435393 :         const std::vector<std::vector<OutputShape>> & dphideta = fe.get_dphideta();
     390      435393 :         const std::vector<std::vector<OutputShape>> & dphidzeta = fe.get_dphidzeta();
     391             : 
     392             :         // Inverse map second derivatives
     393      435393 :         const std::vector<std::vector<Real>> & d2xidxyz2 = fe.get_fe_map().get_d2xidxyz2();
     394      435393 :         const std::vector<std::vector<Real>> & d2etadxyz2 = fe.get_fe_map().get_d2etadxyz2();
     395      435393 :         const std::vector<std::vector<Real>> & d2zetadxyz2 = fe.get_fe_map().get_d2zetadxyz2();
     396             : 
     397    40151510 :         for (auto i : index_range(d2phi))
     398   138761472 :           for (auto p : index_range(d2phi[i]))
     399             :             {
     400             :               // phi_{x x}
     401   125866225 :               d2phi[i][p].slice(0).slice(0) = d2phidx2[i][p] =
     402   151443231 :                 d2phidxi2[i][p]*dxidx_map[p]*dxidx_map[p] +           // (xi_x)^2 * phi_{xi xi}
     403   151443231 :                 d2phideta2[i][p]*detadx_map[p]*detadx_map[p] +        // (eta_x)^2 * phi_{eta eta}
     404   151443231 :                 d2phidzeta2[i][p]*dzetadx_map[p]*dzetadx_map[p] +     // (zeta_x)^2 * phi_{zeta zeta}
     405   125866225 :                 2*d2phidxideta[i][p]*dxidx_map[p]*detadx_map[p] +     // 2 * xi_x * eta_x * phi_{xi eta}
     406   125866225 :                 2*d2phidxidzeta[i][p]*dxidx_map[p]*dzetadx_map[p] +   // 2 * xi_x * zeta_x * phi_{xi zeta}
     407   125866225 :                 2*d2phidetadzeta[i][p]*detadx_map[p]*dzetadx_map[p] + // 2 * eta_x * zeta_x * phi_{eta zeta}
     408   151443231 :                 d2xidxyz2[p][0]*dphidxi[i][p] +                       // xi_{x x} * phi_{xi}
     409   151443231 :                 d2etadxyz2[p][0]*dphideta[i][p] +                     // eta_{x x} * phi_{eta}
     410   177020237 :                 d2zetadxyz2[p][0]*dphidzeta[i][p];                    // zeta_{x x} * phi_{zeta}
     411             : 
     412             :               // phi_{x y}
     413   100289219 :               d2phi[i][p].slice(0).slice(1) = d2phi[i][p].slice(1).slice(0) = d2phidxdy[i][p] =
     414   100289219 :                 d2phidxi2[i][p]*dxidx_map[p]*dxidy_map[p] +                                          // xi_x * xi_y * phi_{xi xi}
     415   100289219 :                 d2phideta2[i][p]*detadx_map[p]*detady_map[p] +                                       // eta_x * eta_y * phi_{eta eta}
     416   100289219 :                 d2phidzeta2[i][p]*dzetadx_map[p]*dzetady_map[p] +                                    // zeta_x * zeta_y * phi_{zeta zeta}
     417   125866225 :                 d2phidxideta[i][p]*(dxidx_map[p]*detady_map[p] + detadx_map[p]*dxidy_map[p]) +       // (xi_x*eta_y + eta_x*xi_y) * phi_{xi eta}
     418   125866225 :                 d2phidxidzeta[i][p]*(dxidx_map[p]*dzetady_map[p] + dzetadx_map[p]*dxidy_map[p]) +    // (zeta_x*xi_y + xi_x*zeta_y) * phi_{xi zeta}
     419   125866225 :                 d2phidetadzeta[i][p]*(detadx_map[p]*dzetady_map[p] + dzetadx_map[p]*detady_map[p]) + // (zeta_x*eta_y + eta_x*zeta_y) * phi_{eta zeta}
     420   125866225 :                 d2xidxyz2[p][1]*dphidxi[i][p] +                                                      // xi_{x y} * phi_{xi}
     421   125866225 :                 d2etadxyz2[p][1]*dphideta[i][p] +                                                    // eta_{x y} * phi_{eta}
     422   100289219 :                 d2zetadxyz2[p][1]*dphidzeta[i][p];                                                   // zeta_{x y} * phi_{zeta}
     423             : 
     424             :               // phi_{x z}
     425   100289219 :               d2phi[i][p].slice(0).slice(2) = d2phi[i][p].slice(2).slice(0) = d2phidxdz[i][p] =
     426   100289219 :                 d2phidxi2[i][p]*dxidx_map[p]*dxidz_map[p] +                                          // xi_x * xi_z * phi_{xi xi}
     427   100289219 :                 d2phideta2[i][p]*detadx_map[p]*detadz_map[p] +                                       // eta_x * eta_z * phi_{eta eta}
     428   100289219 :                 d2phidzeta2[i][p]*dzetadx_map[p]*dzetadz_map[p] +                                    // zeta_x * zeta_z * phi_{zeta zeta}
     429   125866225 :                 d2phidxideta[i][p]*(dxidx_map[p]*detadz_map[p] + detadx_map[p]*dxidz_map[p]) +       // (xi_x*eta_z + eta_x*xi_z) * phi_{xi eta}
     430   125866225 :                 d2phidxidzeta[i][p]*(dxidx_map[p]*dzetadz_map[p] + dzetadx_map[p]*dxidz_map[p]) +    // (zeta_x*xi_z + xi_x*zeta_z) * phi_{xi zeta}
     431   125866225 :                 d2phidetadzeta[i][p]*(detadx_map[p]*dzetadz_map[p] + dzetadx_map[p]*detadz_map[p]) + // (zeta_x*eta_z + eta_x*zeta_z) * phi_{eta zeta}
     432   125866225 :                 d2xidxyz2[p][2]*dphidxi[i][p] +                                                      // xi_{x z} * phi_{xi}
     433   125866225 :                 d2etadxyz2[p][2]*dphideta[i][p] +                                                    // eta_{x z} * phi_{eta}
     434   100289219 :                 d2zetadxyz2[p][2]*dphidzeta[i][p];                                                   // zeta_{x z} * phi_{zeta}
     435             : 
     436             :               // phi_{y y}
     437   100289219 :               d2phi[i][p].slice(1).slice(1) = d2phidy2[i][p] =
     438   125866225 :                 d2phidxi2[i][p]*dxidy_map[p]*dxidy_map[p] +           // (xi_y)^2 * phi_{xi xi}
     439   125866225 :                 d2phideta2[i][p]*detady_map[p]*detady_map[p] +        // (eta_y)^2 * phi_{eta eta}
     440   125866225 :                 d2phidzeta2[i][p]*dzetady_map[p]*dzetady_map[p] +     // (zeta_y)^2 * phi_{zeta zeta}
     441   125866225 :                 2*d2phidxideta[i][p]*dxidy_map[p]*detady_map[p] +     // 2 * xi_y * eta_y * phi_{xi eta}
     442   125866225 :                 2*d2phidxidzeta[i][p]*dxidy_map[p]*dzetady_map[p] +   // 2 * xi_y * zeta_y * phi_{xi zeta}
     443   125866225 :                 2*d2phidetadzeta[i][p]*detady_map[p]*dzetady_map[p] + // 2 * eta_y * zeta_y * phi_{eta zeta}
     444   125866225 :                 d2xidxyz2[p][3]*dphidxi[i][p] +                       // xi_{y y} * phi_{xi}
     445   125866225 :                 d2etadxyz2[p][3]*dphideta[i][p] +                     // eta_{y y} * phi_{eta}
     446   100289219 :                 d2zetadxyz2[p][3]*dphidzeta[i][p];                    // zeta_{y y} * phi_{zeta}
     447             : 
     448             :               // phi_{y z}
     449   100289219 :               d2phi[i][p].slice(1).slice(2) = d2phi[i][p].slice(2).slice(1) = d2phidydz[i][p] =
     450   125866225 :                 d2phidxi2[i][p]*dxidy_map[p]*dxidz_map[p] +                                          // xi_y * xi_z * phi_{xi xi}
     451   125866225 :                 d2phideta2[i][p]*detady_map[p]*detadz_map[p] +                                       // eta_y * eta_z * phi_{eta eta}
     452   125866225 :                 d2phidzeta2[i][p]*dzetady_map[p]*dzetadz_map[p] +                                    // zeta_y * zeta_z * phi_{zeta zeta}
     453   125866225 :                 d2phidxideta[i][p]*(dxidy_map[p]*detadz_map[p] + detady_map[p]*dxidz_map[p]) +       // (xi_y*eta_z + eta_y*xi_z) * phi_{xi eta}
     454   125866225 :                 d2phidxidzeta[i][p]*(dxidy_map[p]*dzetadz_map[p] + dzetady_map[p]*dxidz_map[p]) +    // (zeta_y*xi_z + xi_y*zeta_z) * phi_{xi zeta}
     455   125866225 :                 d2phidetadzeta[i][p]*(detady_map[p]*dzetadz_map[p] + dzetady_map[p]*detadz_map[p]) + // (zeta_y*eta_z + eta_y*zeta_z) * phi_{eta zeta}
     456   125866225 :                 d2xidxyz2[p][4]*dphidxi[i][p] +                                                      // xi_{y z} * phi_{xi}
     457   125866225 :                 d2etadxyz2[p][4]*dphideta[i][p] +                                                    // eta_{y z} * phi_{eta}
     458   100289219 :                 d2zetadxyz2[p][4]*dphidzeta[i][p];                                                   // zeta_{y z} * phi_{zeta}
     459             : 
     460             :               // phi_{z z}
     461   125866225 :               d2phi[i][p].slice(2).slice(2) = d2phidz2[i][p] =
     462   125866225 :                 d2phidxi2[i][p]*dxidz_map[p]*dxidz_map[p] +           // (xi_z)^2 * phi_{xi xi}
     463   125866225 :                 d2phideta2[i][p]*detadz_map[p]*detadz_map[p] +        // (eta_z)^2 * phi_{eta eta}
     464   125866225 :                 d2phidzeta2[i][p]*dzetadz_map[p]*dzetadz_map[p] +     // (zeta_z)^2 * phi_{zeta zeta}
     465   125866225 :                 2*d2phidxideta[i][p]*dxidz_map[p]*detadz_map[p] +     // 2 * xi_z * eta_z * phi_{xi eta}
     466   125866225 :                 2*d2phidxidzeta[i][p]*dxidz_map[p]*dzetadz_map[p] +   // 2 * xi_z * zeta_z * phi_{xi zeta}
     467   125866225 :                 2*d2phidetadzeta[i][p]*detadz_map[p]*dzetadz_map[p] + // 2 * eta_z * zeta_z * phi_{eta zeta}
     468   125866225 :                 d2xidxyz2[p][5]*dphidxi[i][p] +                       // xi_{z z} * phi_{xi}
     469   125866225 :                 d2etadxyz2[p][5]*dphideta[i][p] +                     // eta_{z z} * phi_{eta}
     470   100289219 :                 d2zetadxyz2[p][5]*dphidzeta[i][p];                    // zeta_{z z} * phi_{zeta}
     471             :             }
     472             : 
     473      435393 :         break;
     474             :       }
     475             : 
     476           0 :     default:
     477           0 :       libmesh_error_msg("Invalid dim = " << dim);
     478             :     } // switch(dim)
     479     2014772 : }
     480             : #endif //LIBMESH_ENABLE_SECOND_DERIVATIVES
     481             : 
     482             : template<>
     483           0 : void H1FETransformation<Real>::map_curl(const unsigned int,
     484             :                                         const Elem * const,
     485             :                                         const std::vector<Point> &,
     486             :                                         const FEGenericBase<Real> &,
     487             :                                         std::vector<std::vector<Real>> &) const
     488             : {
     489           0 :   libmesh_error_msg("Computing the curl of a shape function only \nmakes sense for vector-valued elements.");
     490             : }
     491             : 
     492             : template<>
     493       62396 : void H1FETransformation<RealGradient>::map_curl(const unsigned int dim,
     494             :                                                 const Elem * const,
     495             :                                                 const std::vector<Point> &,
     496             :                                                 const FEGenericBase<RealGradient> & fe,
     497             :                                                 std::vector<std::vector<RealGradient>> & curl_phi ) const
     498             : {
     499       62396 :   switch (dim)
     500             :     {
     501           0 :     case 0:
     502             :     case 1:
     503           0 :       libmesh_error_msg("The curl only make sense in 2D and 3D");
     504             : 
     505        6244 :     case 2:
     506             :       {
     507        6244 :         const std::vector<std::vector<RealGradient>> & dphidxi = fe.get_dphidxi();
     508        6244 :         const std::vector<std::vector<RealGradient>> & dphideta = fe.get_dphideta();
     509             : 
     510        6244 :         const std::vector<Real> & dxidx_map = fe.get_fe_map().get_dxidx();
     511        6244 :         const std::vector<Real> & dxidy_map = fe.get_fe_map().get_dxidy();
     512             : #if LIBMESH_DIM > 2
     513        6244 :         const std::vector<Real> & dxidz_map = fe.get_fe_map().get_dxidz();
     514             : #endif
     515             : 
     516        6244 :         const std::vector<Real> & detadx_map = fe.get_fe_map().get_detadx();
     517        6244 :         const std::vector<Real> & detady_map = fe.get_fe_map().get_detady();
     518             : #if LIBMESH_DIM > 2
     519        6244 :         const std::vector<Real> & detadz_map = fe.get_fe_map().get_detadz();
     520             : #endif
     521             : 
     522             :         // For 2D elements in 3D space:
     523             :         // curl = ( -dphi_y/dz, dphi_x/dz, dphi_y/dx - dphi_x/dy )
     524      246384 :         for (auto i : index_range(curl_phi))
     525     2517240 :           for (auto p : index_range(curl_phi[i]))
     526             :             {
     527             : 
     528     3440080 :               Real dphiy_dx = (dphidxi[i][p].slice(1))*dxidx_map[p]
     529     3440080 :                 + (dphideta[i][p].slice(1))*detadx_map[p];
     530             : 
     531     2295888 :               Real dphix_dy = (dphidxi[i][p].slice(0))*dxidy_map[p]
     532     2295888 :                 + (dphideta[i][p].slice(0))*detady_map[p];
     533             : 
     534     2295888 :               curl_phi[i][p].slice(2) = dphiy_dx - dphix_dy;
     535             : 
     536             : #if LIBMESH_DIM > 2
     537     2295888 :               Real dphiy_dz = (dphidxi[i][p].slice(1))*dxidz_map[p]
     538     2295888 :                 + (dphideta[i][p].slice(1))*detadz_map[p];
     539             : 
     540     2295888 :               Real dphix_dz = (dphidxi[i][p].slice(0))*dxidz_map[p]
     541     2295888 :                 + (dphideta[i][p].slice(0))*detadz_map[p];
     542             : 
     543     2295888 :               curl_phi[i][p].slice(0) = -dphiy_dz;
     544     2295888 :               curl_phi[i][p].slice(1) = dphix_dz;
     545             : #endif
     546             :             }
     547             : 
     548        6244 :         break;
     549             :       }
     550        9341 :     case 3:
     551             :       {
     552        9341 :         const std::vector<std::vector<RealGradient>> & dphidxi = fe.get_dphidxi();
     553        9341 :         const std::vector<std::vector<RealGradient>> & dphideta = fe.get_dphideta();
     554        9341 :         const std::vector<std::vector<RealGradient>> & dphidzeta = fe.get_dphidzeta();
     555             : 
     556        9341 :         const std::vector<Real> & dxidx_map = fe.get_fe_map().get_dxidx();
     557        9341 :         const std::vector<Real> & dxidy_map = fe.get_fe_map().get_dxidy();
     558        9341 :         const std::vector<Real> & dxidz_map = fe.get_fe_map().get_dxidz();
     559             : 
     560        9341 :         const std::vector<Real> & detadx_map = fe.get_fe_map().get_detadx();
     561        9341 :         const std::vector<Real> & detady_map = fe.get_fe_map().get_detady();
     562        9341 :         const std::vector<Real> & detadz_map = fe.get_fe_map().get_detadz();
     563             : 
     564        9341 :         const std::vector<Real> & dzetadx_map = fe.get_fe_map().get_dzetadx();
     565        9341 :         const std::vector<Real> & dzetady_map = fe.get_fe_map().get_dzetady();
     566        9341 :         const std::vector<Real> & dzetadz_map = fe.get_fe_map().get_dzetadz();
     567             : 
     568             :         // In 3D: curl = ( dphi_z/dy - dphi_y/dz, dphi_x/dz - dphi_z/dx, dphi_y/dx - dphi_x/dy )
     569      491732 :         for (auto i : index_range(curl_phi))
     570     7415520 :           for (auto p : index_range(curl_phi[i]))
     571             :             {
     572    10441728 :               Real dphiz_dy = (dphidxi[i][p].slice(2))*dxidy_map[p]
     573    10441728 :                 + (dphideta[i][p].slice(2))*detady_map[p]
     574    10441728 :                 + (dphidzeta[i][p].slice(2))*dzetady_map[p];
     575             : 
     576     6961152 :               Real dphiy_dz = (dphidxi[i][p].slice(1))*dxidz_map[p]
     577     6961152 :                 + (dphideta[i][p].slice(1))*detadz_map[p]
     578     6961152 :                 + (dphidzeta[i][p].slice(1))*dzetadz_map[p];
     579             : 
     580     6961152 :               Real dphix_dz = (dphidxi[i][p].slice(0))*dxidz_map[p]
     581     6961152 :                 + (dphideta[i][p].slice(0))*detadz_map[p]
     582     6961152 :                 + (dphidzeta[i][p].slice(0))*dzetadz_map[p];
     583             : 
     584     6961152 :               Real dphiz_dx = (dphidxi[i][p].slice(2))*dxidx_map[p]
     585     6961152 :                 + (dphideta[i][p].slice(2))*detadx_map[p]
     586     6961152 :                 + (dphidzeta[i][p].slice(2))*dzetadx_map[p];
     587             : 
     588     6961152 :               Real dphiy_dx = (dphidxi[i][p].slice(1))*dxidx_map[p]
     589     6961152 :                 + (dphideta[i][p].slice(1))*detadx_map[p]
     590     6961152 :                 + (dphidzeta[i][p].slice(1))*dzetadx_map[p];
     591             : 
     592     6961152 :               Real dphix_dy = (dphidxi[i][p].slice(0))*dxidy_map[p]
     593     6961152 :                 + (dphideta[i][p].slice(0))*detady_map[p]
     594     6961152 :                 + (dphidzeta[i][p].slice(0))*dzetady_map[p];
     595             : 
     596     6961152 :               curl_phi[i][p].slice(0) = dphiz_dy - dphiy_dz;
     597             : 
     598     6961152 :               curl_phi[i][p].slice(1) = dphix_dz - dphiz_dx;
     599             : 
     600     6961152 :               curl_phi[i][p].slice(2) = dphiy_dx - dphix_dy;
     601             :             }
     602             : 
     603        9341 :         break;
     604             :       }
     605           0 :     default:
     606           0 :       libmesh_error_msg("Invalid dim = " << dim);
     607             :     }
     608       62396 : }
     609             : 
     610             : 
     611             : template<>
     612           0 : void H1FETransformation<Real>::map_div(const unsigned int,
     613             :                                        const Elem * const,
     614             :                                        const std::vector<Point> &,
     615             :                                        const FEGenericBase<Real> &,
     616             :                                        std::vector<std::vector<FEGenericBase<Real>::OutputDivergence>> & ) const
     617             : {
     618           0 :   libmesh_error_msg("Computing the divergence of a shape function only \nmakes sense for vector-valued elements.");
     619             : }
     620             : 
     621             : 
     622             : template<>
     623      189644 : void H1FETransformation<RealGradient>::map_div(const unsigned int dim,
     624             :                                                const Elem * const,
     625             :                                                const std::vector<Point> &,
     626             :                                                const FEGenericBase<RealGradient> & fe,
     627             :                                                std::vector<std::vector<FEGenericBase<RealGradient>::OutputDivergence>> & div_phi) const
     628             : {
     629      189644 :   switch (dim)
     630             :     {
     631           0 :     case 0:
     632             :     case 1:
     633           0 :       libmesh_error_msg("The divergence only make sense in 2D and 3D");
     634             : 
     635       19884 :     case 2:
     636             :       {
     637       19884 :         const std::vector<std::vector<RealGradient>> & dphidxi = fe.get_dphidxi();
     638       19884 :         const std::vector<std::vector<RealGradient>> & dphideta = fe.get_dphideta();
     639             : 
     640       19884 :         const std::vector<Real> & dxidx_map = fe.get_fe_map().get_dxidx();
     641       19884 :         const std::vector<Real> & dxidy_map = fe.get_fe_map().get_dxidy();
     642             : 
     643       19884 :         const std::vector<Real> & detadx_map = fe.get_fe_map().get_detadx();
     644       19884 :         const std::vector<Real> & detady_map = fe.get_fe_map().get_detady();
     645             : 
     646             :         // In 2D: div = dphi_x/dx + dphi_y/dy
     647      763344 :         for (auto i : index_range(div_phi))
     648     6291384 :           for (auto p : index_range(div_phi[i]))
     649             :             {
     650     8425120 :               Real dphix_dx = (dphidxi[i][p].slice(0))*dxidx_map[p]
     651     8425120 :                 + (dphideta[i][p].slice(0))*detadx_map[p];
     652             : 
     653     5606592 :               Real dphiy_dy = (dphidxi[i][p].slice(1))*dxidy_map[p]
     654     5606592 :                 + (dphideta[i][p].slice(1))*detady_map[p];
     655             : 
     656     7015856 :               div_phi[i][p] = dphix_dx + dphiy_dy;
     657             :             }
     658       19884 :         break;
     659             :       }
     660       27773 :     case 3:
     661             :       {
     662       27773 :         const std::vector<std::vector<RealGradient>> & dphidxi = fe.get_dphidxi();
     663       27773 :         const std::vector<std::vector<RealGradient>> & dphideta = fe.get_dphideta();
     664       27773 :         const std::vector<std::vector<RealGradient>> & dphidzeta = fe.get_dphidzeta();
     665             : 
     666       27773 :         const std::vector<Real> & dxidx_map = fe.get_fe_map().get_dxidx();
     667       27773 :         const std::vector<Real> & dxidy_map = fe.get_fe_map().get_dxidy();
     668       27773 :         const std::vector<Real> & dxidz_map = fe.get_fe_map().get_dxidz();
     669             : 
     670       27773 :         const std::vector<Real> & detadx_map = fe.get_fe_map().get_detadx();
     671       27773 :         const std::vector<Real> & detady_map = fe.get_fe_map().get_detady();
     672       27773 :         const std::vector<Real> & detadz_map = fe.get_fe_map().get_detadz();
     673             : 
     674       27773 :         const std::vector<Real> & dzetadx_map = fe.get_fe_map().get_dzetadx();
     675       27773 :         const std::vector<Real> & dzetady_map = fe.get_fe_map().get_dzetady();
     676       27773 :         const std::vector<Real> & dzetadz_map = fe.get_fe_map().get_dzetadz();
     677             : 
     678             :         // In 3D: div = dphi_x/dx + dphi_y/dy + dphi_z/dz
     679     1450196 :         for (auto i : index_range(div_phi))
     680    12723936 :           for (auto p : index_range(div_phi[i]))
     681             :             {
     682    17077248 :               Real dphix_dx = (dphidxi[i][p].slice(0))*dxidx_map[p]
     683    17077248 :                 + (dphideta[i][p].slice(0))*detadx_map[p]
     684    17077248 :                 + (dphidzeta[i][p].slice(0))*dzetadx_map[p];
     685             : 
     686    11384832 :               Real dphiy_dy = (dphidxi[i][p].slice(1))*dxidy_map[p]
     687    11384832 :                 + (dphideta[i][p].slice(1))*detady_map[p]
     688    11384832 :                 + (dphidzeta[i][p].slice(1))*dzetady_map[p];
     689             : 
     690    11384832 :               Real dphiz_dz = (dphidxi[i][p].slice(2))*dxidz_map[p]
     691    11384832 :                 + (dphideta[i][p].slice(2))*detadz_map[p]
     692    11384832 :                 + (dphidzeta[i][p].slice(2))*dzetadz_map[p];
     693             : 
     694    14231040 :               div_phi[i][p] = dphix_dx + dphiy_dy + dphiz_dz;
     695             :             }
     696             : 
     697       27773 :         break;
     698             :       }
     699             : 
     700           0 :     default:
     701           0 :       libmesh_error_msg("Invalid dim = " << dim);
     702             :     } // switch(dim)
     703             : 
     704      189644 :   return;
     705             : }
     706             : 
     707             : // Explicit Instantiations
     708             : template class LIBMESH_EXPORT H1FETransformation<Real>;
     709             : template class LIBMESH_EXPORT H1FETransformation<RealGradient>;
     710             : 
     711             : } //namespace libMesh