Line data    Source code

       1             : // The libMesh Finite Element Library.
       2             : // Copyright (C) 2002-2025 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   300656408 : void H1FETransformation<OutputShape>::init_map_phi(const FEGenericBase<OutputShape> & /* fe */ ) const
      28             : {
      29             :   // Nothing needed
      30   300656408 : }
      31             : 
      32             : 
      33             : 
      34             : template<typename OutputShape>
      35   217812516 : void H1FETransformation<OutputShape>::init_map_dphi(const FEGenericBase<OutputShape> & fe) const
      36             : {
      37    19354772 :   fe.get_fe_map().get_dxidx();
      38   217812516 : }
      39             : 
      40             : 
      41             : 
      42             : template<typename OutputShape>
      43    58107956 : void H1FETransformation<OutputShape>::init_map_d2phi(const FEGenericBase<OutputShape> & fe) const
      44             : {
      45     4101306 :   fe.get_fe_map().get_dxidx();
      46             : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
      47     4101306 :   fe.get_fe_map().get_d2xidxyz2();
      48             : #endif
      49    58107956 : }
      50             : 
      51             : template<typename OutputShape>
      52   103801535 : 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   113113247 :   FEInterface::all_shapes<OutputShape>(dim, fe.get_fe_type(), elem, qp, phi, add_p_level);
      60   103801535 : }
      61             : 
      62             : 
      63             : template<typename OutputShape>
      64    79160232 : 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    79160232 :   switch(dim)
      74             :     {
      75       18792 :     case 0: // No derivatives in 0D
      76             :       {
      77      391392 :         for (auto i : index_range(dphi))
      78      391392 :           for (auto p : index_range(dphi[i]))
      79       18792 :             dphi[i][p] = 0.;
      80       18792 :         break;
      81             :       }
      82             : 
      83       13500 :     case 1:
      84             :       {
      85       13500 :         const std::vector<std::vector<OutputShape>> & dphidxi = fe.get_dphidxi();
      86             : 
      87       13500 :         const std::vector<Real> & dxidx_map = fe.get_fe_map().get_dxidx();
      88             : #if LIBMESH_DIM>1
      89       13500 :         const std::vector<Real> & dxidy_map = fe.get_fe_map().get_dxidy();
      90             : #endif
      91             : #if LIBMESH_DIM>2
      92       13500 :         const std::vector<Real> & dxidz_map = fe.get_fe_map().get_dxidz();
      93             : #endif
      94             : 
      95     1854250 :         for (auto i : index_range(dphi))
      96     6548502 :           for (auto p : index_range(dphi[i]))
      97             :             {
      98             :               // dphi/dx    = (dphi/dxi)*(dxi/dx)
      99     5736892 :               dphi[i][p].slice(0) = dphidx[i][p] = dphidxi[i][p]*dxidx_map[p];
     100             : 
     101             : #if LIBMESH_DIM>1
     102     5335858 :               dphi[i][p].slice(1)  = dphidy[i][p] = dphidxi[i][p]*dxidy_map[p];
     103             : #endif
     104             : #if LIBMESH_DIM>2
     105     5469536 :               dphi[i][p].slice(2) = dphidz[i][p] = dphidxi[i][p]*dxidz_map[p];
     106             : #endif
     107             :             }
     108             : 
     109       13500 :         break;
     110             :       }
     111             : 
     112     6618760 :     case 2:
     113             :       {
     114     6618760 :         const std::vector<std::vector<OutputShape>> & dphidxi = fe.get_dphidxi();
     115     6618760 :         const std::vector<std::vector<OutputShape>> & dphideta = fe.get_dphideta();
     116             : 
     117     6618760 :         const std::vector<Real> & dxidx_map = fe.get_fe_map().get_dxidx();
     118     6618760 :         const std::vector<Real> & dxidy_map = fe.get_fe_map().get_dxidy();
     119             : #if LIBMESH_DIM > 2
     120     6618760 :         const std::vector<Real> & dxidz_map = fe.get_fe_map().get_dxidz();
     121             : #endif
     122             : 
     123     6618760 :         const std::vector<Real> & detadx_map = fe.get_fe_map().get_detadx();
     124     6618760 :         const std::vector<Real> & detady_map = fe.get_fe_map().get_detady();
     125             : #if LIBMESH_DIM > 2
     126     6618760 :         const std::vector<Real> & detadz_map = fe.get_fe_map().get_detadz();
     127             : #endif
     128             : 
     129   372031837 :         for (auto i : index_range(dphi))
     130  1735124509 :           for (auto p : index_range(dphi[i]))
     131             :             {
     132             :               // dphi/dx    = (dphi/dxi)*(dxi/dx) + (dphi/deta)*(deta/dx)
     133  2100867612 :               dphi[i][p].slice(0) = dphidx[i][p] = (dphidxi[i][p]*dxidx_map[p] +
     134  1824899800 :                                                     dphideta[i][p]*detadx_map[p]);
     135             : 
     136             :               // dphi/dy    = (dphi/dxi)*(dxi/dy) + (dphi/deta)*(deta/dy)
     137  1701723234 :               dphi[i][p].slice(1) = dphidy[i][p] = (dphidxi[i][p]*dxidy_map[p] +
     138  1558803548 :                                                     dphideta[i][p]*detady_map[p]);
     139             : 
     140             : #if LIBMESH_DIM > 2
     141             :               // dphi/dz    = (dphi/dxi)*(dxi/dz) + (dphi/deta)*(deta/dz)
     142  1834034202 :               dphi[i][p].slice(2) = dphidz[i][p] = (dphidxi[i][p]*dxidz_map[p] +
     143  1559540706 :                                                     dphideta[i][p]*detadz_map[p]);
     144             : #endif
     145             :             }
     146             : 
     147     6618760 :         break;
     148             :       }
     149             : 
     150      714182 :     case 3:
     151             :       {
     152      714182 :         const std::vector<std::vector<OutputShape>> & dphidxi = fe.get_dphidxi();
     153      714182 :         const std::vector<std::vector<OutputShape>> & dphideta = fe.get_dphideta();
     154      714182 :         const std::vector<std::vector<OutputShape>> & dphidzeta = fe.get_dphidzeta();
     155             : 
     156      714182 :         const std::vector<Real> & dxidx_map = fe.get_fe_map().get_dxidx();
     157      714182 :         const std::vector<Real> & dxidy_map = fe.get_fe_map().get_dxidy();
     158      714182 :         const std::vector<Real> & dxidz_map = fe.get_fe_map().get_dxidz();
     159             : 
     160      714182 :         const std::vector<Real> & detadx_map = fe.get_fe_map().get_detadx();
     161      714182 :         const std::vector<Real> & detady_map = fe.get_fe_map().get_detady();
     162      714182 :         const std::vector<Real> & detadz_map = fe.get_fe_map().get_detadz();
     163             : 
     164      714182 :         const std::vector<Real> & dzetadx_map = fe.get_fe_map().get_dzetadx();
     165      714182 :         const std::vector<Real> & dzetady_map = fe.get_fe_map().get_dzetady();
     166      714182 :         const std::vector<Real> & dzetadz_map = fe.get_fe_map().get_dzetadz();
     167             : 
     168   149623668 :         for (auto i : index_range(dphi))
     169   638602216 :           for (auto p : index_range(dphi[i]))
     170             :             {
     171             :               // dphi/dx    = (dphi/dxi)*(dxi/dx) + (dphi/deta)*(deta/dx) + (dphi/dzeta)*(dzeta/dx);
     172   719306597 :               dphi[i][p].slice(0) = dphidx[i][p] = (dphidxi[i][p]*dxidx_map[p] +
     173   569742746 :                                                     dphideta[i][p]*detadx_map[p] +
     174   609463731 :                                                     dphidzeta[i][p]*dzetadx_map[p]);
     175             : 
     176             :               // dphi/dy    = (dphi/dxi)*(dxi/dy) + (dphi/deta)*(deta/dy) + (dphi/dzeta)*(dzeta/dy);
     177   589413914 :               dphi[i][p].slice(1) = dphidy[i][p] = (dphidxi[i][p]*dxidy_map[p] +
     178   483147624 :                                                     dphideta[i][p]*detady_map[p] +
     179   522868609 :                                                     dphidzeta[i][p]*dzetady_map[p]);
     180             : 
     181             :               // dphi/dz    = (dphi/dxi)*(dxi/dz) + (dphi/deta)*(deta/dz) + (dphi/dzeta)*(dzeta/dz);
     182   630923187 :               dphi[i][p].slice(2) = dphidz[i][p] = (dphidxi[i][p]*dxidz_map[p] +
     183   483147624 :                                                     dphideta[i][p]*detadz_map[p] +
     184   524656897 :                                                     dphidzeta[i][p]*dzetadz_map[p]);
     185             :             }
     186      714182 :         break;
     187             :       }
     188             : 
     189           0 :     default:
     190           0 :       libmesh_error_msg("Invalid dim = " << dim);
     191             :     } // switch(dim)
     192    79160232 : }
     193             : 
     194             : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
     195             : template<typename OutputShape>
     196     6568893 : 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     6568893 :   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     1507539 :         for (auto i : index_range(d2phi))
     237     5899178 :           for (auto p : index_range(d2phi[i]))
     238             :             {
     239             :               // phi_{x x}
     240     4797676 :               d2phi[i][p].slice(0).slice(0) = d2phidx2[i][p] =
     241     4900600 :                 d2phidxi2[i][p]*dxidx_map[p]*dxidx_map[p] + // (xi_x)^2 * phi_{xi xi}
     242     5003524 :                 d2xidxyz2[p][0]*dphidxi[i][p];              // xi_{x x} * phi_{xi}
     243             : 
     244             : #if LIBMESH_DIM>1
     245             :               // phi_{x y}
     246     4694752 :               d2phi[i][p].slice(0).slice(1) = d2phi[i][p].slice(1).slice(0) = d2phidxdy[i][p] =
     247     4694752 :                 d2phidxi2[i][p]*dxidx_map[p]*dxidy_map[p] + // xi_x * xi_y * phi_{xi xi}
     248     4694752 :                 d2xidxyz2[p][1]*dphidxi[i][p];              // xi_{x y} * phi_{xi}
     249             : #endif
     250             : 
     251             : #if LIBMESH_DIM>2
     252             :               // phi_{x z}
     253     4694752 :               d2phi[i][p].slice(0).slice(2) = d2phi[i][p].slice(2).slice(0) = d2phidxdz[i][p] =
     254     4694752 :                 d2phidxi2[i][p]*dxidx_map[p]*dxidz_map[p] + // xi_x * xi_z * phi_{xi xi}
     255     4694752 :                 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     4694752 :               d2phi[i][p].slice(1).slice(1) = d2phidy2[i][p] =
     262     4797676 :                 d2phidxi2[i][p]*dxidy_map[p]*dxidy_map[p] + // (xi_y)^2 * phi_{xi xi}
     263     4694752 :                 d2xidxyz2[p][3]*dphidxi[i][p];              // xi_{y y} * phi_{xi}
     264             : #endif
     265             : 
     266             : #if LIBMESH_DIM>2
     267             :               // phi_{y z}
     268     4694752 :               d2phi[i][p].slice(1).slice(2) = d2phi[i][p].slice(2).slice(1) = d2phidydz[i][p] =
     269     4797676 :                 d2phidxi2[i][p]*dxidy_map[p]*dxidz_map[p] + // xi_y * xi_z * phi_{xi xi}
     270     4694752 :                 d2xidxyz2[p][4]*dphidxi[i][p];              // xi_{y z} * phi_{xi}
     271             : 
     272             :               // phi_{z z}
     273     4797676 :               d2phi[i][p].slice(2).slice(2) = d2phidz2[i][p] =
     274     4797676 :                 d2phidxi2[i][p]*dxidz_map[p]*dxidz_map[p] + // (xi_z)^2 * phi_{xi xi}
     275     4694752 :                 d2xidxyz2[p][5]*dphidxi[i][p];              // xi_{z z} * phi_{xi}
     276             : #endif
     277             :             }
     278        7660 :         break;
     279             :       }
     280             : 
     281       81763 :     case 2:
     282             :       {
     283       81763 :         const std::vector<std::vector<OutputShape>> & d2phidxi2 = fe.get_d2phidxi2();
     284       81763 :         const std::vector<std::vector<OutputShape>> & d2phidxideta = fe.get_d2phidxideta();
     285       81763 :         const std::vector<std::vector<OutputShape>> & d2phideta2 = fe.get_d2phideta2();
     286             : 
     287       81763 :         const std::vector<Real> & dxidx_map = fe.get_fe_map().get_dxidx();
     288       81763 :         const std::vector<Real> & dxidy_map = fe.get_fe_map().get_dxidy();
     289             : #if LIBMESH_DIM > 2
     290       81763 :         const std::vector<Real> & dxidz_map = fe.get_fe_map().get_dxidz();
     291             : #endif
     292             : 
     293       81763 :         const std::vector<Real> & detadx_map = fe.get_fe_map().get_detadx();
     294       81763 :         const std::vector<Real> & detady_map = fe.get_fe_map().get_detady();
     295             : #if LIBMESH_DIM > 2
     296       81763 :         const std::vector<Real> & detadz_map = fe.get_fe_map().get_detadz();
     297             : #endif
     298             : 
     299             :         // Shape function derivatives in reference space
     300       81763 :         const std::vector<std::vector<OutputShape>> & dphidxi = fe.get_dphidxi();
     301       81763 :         const std::vector<std::vector<OutputShape>> & dphideta = fe.get_dphideta();
     302             : 
     303             :         // Inverse map second derivatives
     304       81763 :         const std::vector<std::vector<Real>> & d2xidxyz2 = fe.get_fe_map().get_d2xidxyz2();
     305       81763 :         const std::vector<std::vector<Real>> & d2etadxyz2 = fe.get_fe_map().get_d2etadxyz2();
     306             : 
     307     9706021 :         for (auto i : index_range(d2phi))
     308    58685558 :           for (auto p : index_range(d2phi[i]))
     309             :             {
     310             :               // phi_{x x}
     311    53725785 :               d2phi[i][p].slice(0).slice(0) = d2phidx2[i][p] =
     312    57501816 :                 d2phidxi2[i][p]*dxidx_map[p]*dxidx_map[p] +       // (xi_x)^2 * phi_{xi xi}
     313    57501816 :                 d2phideta2[i][p]*detadx_map[p]*detadx_map[p] +    // (eta_x)^2 * phi_{eta eta}
     314    53725785 :                 2*d2phidxideta[i][p]*dxidx_map[p]*detadx_map[p] + // 2 * xi_x * eta_x * phi_{xi eta}
     315    57501816 :                 d2xidxyz2[p][0]*dphidxi[i][p] +                   // xi_{x x} * phi_{xi}
     316    61277847 :                 d2etadxyz2[p][0]*dphideta[i][p];                  // eta_{x x} * phi_{eta}
     317             : 
     318             :               // phi_{x y}
     319    49949754 :               d2phi[i][p].slice(0).slice(1) = d2phi[i][p].slice(1).slice(0) = d2phidxdy[i][p] =
     320    49949754 :                 d2phidxi2[i][p]*dxidx_map[p]*dxidy_map[p] +                                    // xi_x * xi_y * phi_{xi xi}
     321    49949754 :                 d2phideta2[i][p]*detadx_map[p]*detady_map[p] +                                 // eta_x * eta_y * phi_{eta eta}
     322    53725785 :                 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    53725785 :                 d2xidxyz2[p][1]*dphidxi[i][p] +                                                // xi_{x y} * phi_{xi}
     324    49949754 :                 d2etadxyz2[p][1]*dphideta[i][p];                                               // eta_{x y} * phi_{eta}
     325             : 
     326             : #if LIBMESH_DIM > 2
     327             :               // phi_{x z}
     328    49949754 :               d2phi[i][p].slice(0).slice(2) = d2phi[i][p].slice(2).slice(0) = d2phidxdz[i][p] =
     329    49949754 :                 d2phidxi2[i][p]*dxidx_map[p]*dxidz_map[p] +                                    // xi_x * xi_z * phi_{xi xi}
     330    49949754 :                 d2phideta2[i][p]*detadx_map[p]*detadz_map[p] +                                 // eta_x * eta_z * phi_{eta eta}
     331    53725785 :                 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    53725785 :                 d2xidxyz2[p][2]*dphidxi[i][p] +                                                // xi_{x z} * phi_{xi}
     333    49949754 :                 d2etadxyz2[p][2]*dphideta[i][p];                                               // eta_{x z} * phi_{eta}
     334             : #endif
     335             : 
     336             :               // phi_{y y}
     337    49949754 :               d2phi[i][p].slice(1).slice(1) = d2phidy2[i][p] =
     338    53725785 :                 d2phidxi2[i][p]*dxidy_map[p]*dxidy_map[p] +       // (xi_y)^2 * phi_{xi xi}
     339    53725785 :                 d2phideta2[i][p]*detady_map[p]*detady_map[p] +    // (eta_y)^2 * phi_{eta eta}
     340    53725785 :                 2*d2phidxideta[i][p]*dxidy_map[p]*detady_map[p] + // 2 * xi_y * eta_y * phi_{xi eta}
     341    53725785 :                 d2xidxyz2[p][3]*dphidxi[i][p] +                   // xi_{y y} * phi_{xi}
     342    49949754 :                 d2etadxyz2[p][3]*dphideta[i][p];                  // eta_{y y} * phi_{eta}
     343             : 
     344             : #if LIBMESH_DIM > 2
     345             :               // phi_{y z}
     346    49949754 :               d2phi[i][p].slice(1).slice(2) = d2phi[i][p].slice(2).slice(1) = d2phidydz[i][p] =
     347    53725785 :                 d2phidxi2[i][p]*dxidy_map[p]*dxidz_map[p] +                                    // xi_y * xi_z * phi_{xi xi}
     348    53725785 :                 d2phideta2[i][p]*detady_map[p]*detadz_map[p] +                                 // eta_y * eta_z * phi_{eta eta}
     349    53725785 :                 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    53725785 :                 d2xidxyz2[p][4]*dphidxi[i][p] +                                                // xi_{y z} * phi_{xi}
     351    49949754 :                 d2etadxyz2[p][4]*dphideta[i][p];                                               // eta_{y z} * phi_{eta}
     352             : 
     353             :               // phi_{z z}
     354    53725785 :               d2phi[i][p].slice(2).slice(2) = d2phidz2[i][p] =
     355    53725785 :                 d2phidxi2[i][p]*dxidz_map[p]*dxidz_map[p] +       // (xi_z)^2 * phi_{xi xi}
     356    53725785 :                 d2phideta2[i][p]*detadz_map[p]*detadz_map[p] +    // (eta_z)^2 * phi_{eta eta}
     357    53725785 :                 2*d2phidxideta[i][p]*dxidz_map[p]*detadz_map[p] + // 2 * xi_z * eta_z * phi_{xi eta}
     358    53725785 :                 d2xidxyz2[p][5]*dphidxi[i][p] +                   // xi_{z z} * phi_{xi}
     359    49949754 :                 d2etadxyz2[p][5]*dphideta[i][p];                  // eta_{z z} * phi_{eta}
     360             : #endif
     361             :             }
     362             : 
     363       81763 :         break;
     364             :       }
     365             : 
     366      432120 :     case 3:
     367             :       {
     368      432120 :         const std::vector<std::vector<OutputShape>> & d2phidxi2 = fe.get_d2phidxi2();
     369      432120 :         const std::vector<std::vector<OutputShape>> & d2phidxideta = fe.get_d2phidxideta();
     370      432120 :         const std::vector<std::vector<OutputShape>> & d2phideta2 = fe.get_d2phideta2();
     371      432120 :         const std::vector<std::vector<OutputShape>> & d2phidxidzeta = fe.get_d2phidxidzeta();
     372      432120 :         const std::vector<std::vector<OutputShape>> & d2phidetadzeta = fe.get_d2phidetadzeta();
     373      432120 :         const std::vector<std::vector<OutputShape>> & d2phidzeta2 = fe.get_d2phidzeta2();
     374             : 
     375      432120 :         const std::vector<Real> & dxidx_map = fe.get_fe_map().get_dxidx();
     376      432120 :         const std::vector<Real> & dxidy_map = fe.get_fe_map().get_dxidy();
     377      432120 :         const std::vector<Real> & dxidz_map = fe.get_fe_map().get_dxidz();
     378             : 
     379      432120 :         const std::vector<Real> & detadx_map = fe.get_fe_map().get_detadx();
     380      432120 :         const std::vector<Real> & detady_map = fe.get_fe_map().get_detady();
     381      432120 :         const std::vector<Real> & detadz_map = fe.get_fe_map().get_detadz();
     382             : 
     383      432120 :         const std::vector<Real> & dzetadx_map = fe.get_fe_map().get_dzetadx();
     384      432120 :         const std::vector<Real> & dzetady_map = fe.get_fe_map().get_dzetady();
     385      432120 :         const std::vector<Real> & dzetadz_map = fe.get_fe_map().get_dzetadz();
     386             : 
     387             :         // Shape function derivatives in reference space
     388      432120 :         const std::vector<std::vector<OutputShape>> & dphidxi = fe.get_dphidxi();
     389      432120 :         const std::vector<std::vector<OutputShape>> & dphideta = fe.get_dphideta();
     390      432120 :         const std::vector<std::vector<OutputShape>> & dphidzeta = fe.get_dphidzeta();
     391             : 
     392             :         // Inverse map second derivatives
     393      432120 :         const std::vector<std::vector<Real>> & d2xidxyz2 = fe.get_fe_map().get_d2xidxyz2();
     394      432120 :         const std::vector<std::vector<Real>> & d2etadxyz2 = fe.get_fe_map().get_d2etadxyz2();
     395      432120 :         const std::vector<std::vector<Real>> & d2zetadxyz2 = fe.get_fe_map().get_d2zetadxyz2();
     396             : 
     397   125449572 :         for (auto i : index_range(d2phi))
     398   429680475 :           for (auto p : index_range(d2phi[i]))
     399             :             {
     400             :               // phi_{x x}
     401   335089672 :               d2phi[i][p].slice(0).slice(0) = d2phidx2[i][p] =
     402   360652878 :                 d2phidxi2[i][p]*dxidx_map[p]*dxidx_map[p] +           // (xi_x)^2 * phi_{xi xi}
     403   360652878 :                 d2phideta2[i][p]*detadx_map[p]*detadx_map[p] +        // (eta_x)^2 * phi_{eta eta}
     404   360652878 :                 d2phidzeta2[i][p]*dzetadx_map[p]*dzetadx_map[p] +     // (zeta_x)^2 * phi_{zeta zeta}
     405   335089672 :                 2*d2phidxideta[i][p]*dxidx_map[p]*detadx_map[p] +     // 2 * xi_x * eta_x * phi_{xi eta}
     406   335089672 :                 2*d2phidxidzeta[i][p]*dxidx_map[p]*dzetadx_map[p] +   // 2 * xi_x * zeta_x * phi_{xi zeta}
     407   335089672 :                 2*d2phidetadzeta[i][p]*detadx_map[p]*dzetadx_map[p] + // 2 * eta_x * zeta_x * phi_{eta zeta}
     408   360652878 :                 d2xidxyz2[p][0]*dphidxi[i][p] +                       // xi_{x x} * phi_{xi}
     409   360652878 :                 d2etadxyz2[p][0]*dphideta[i][p] +                     // eta_{x x} * phi_{eta}
     410   386216084 :                 d2zetadxyz2[p][0]*dphidzeta[i][p];                    // zeta_{x x} * phi_{zeta}
     411             : 
     412             :               // phi_{x y}
     413   309526466 :               d2phi[i][p].slice(0).slice(1) = d2phi[i][p].slice(1).slice(0) = d2phidxdy[i][p] =
     414   309526466 :                 d2phidxi2[i][p]*dxidx_map[p]*dxidy_map[p] +                                          // xi_x * xi_y * phi_{xi xi}
     415   309526466 :                 d2phideta2[i][p]*detadx_map[p]*detady_map[p] +                                       // eta_x * eta_y * phi_{eta eta}
     416   309526466 :                 d2phidzeta2[i][p]*dzetadx_map[p]*dzetady_map[p] +                                    // zeta_x * zeta_y * phi_{zeta zeta}
     417   335089672 :                 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   335089672 :                 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   335089672 :                 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   335089672 :                 d2xidxyz2[p][1]*dphidxi[i][p] +                                                      // xi_{x y} * phi_{xi}
     421   335089672 :                 d2etadxyz2[p][1]*dphideta[i][p] +                                                    // eta_{x y} * phi_{eta}
     422   309526466 :                 d2zetadxyz2[p][1]*dphidzeta[i][p];                                                   // zeta_{x y} * phi_{zeta}
     423             : 
     424             :               // phi_{x z}
     425   309526466 :               d2phi[i][p].slice(0).slice(2) = d2phi[i][p].slice(2).slice(0) = d2phidxdz[i][p] =
     426   309526466 :                 d2phidxi2[i][p]*dxidx_map[p]*dxidz_map[p] +                                          // xi_x * xi_z * phi_{xi xi}
     427   309526466 :                 d2phideta2[i][p]*detadx_map[p]*detadz_map[p] +                                       // eta_x * eta_z * phi_{eta eta}
     428   309526466 :                 d2phidzeta2[i][p]*dzetadx_map[p]*dzetadz_map[p] +                                    // zeta_x * zeta_z * phi_{zeta zeta}
     429   335089672 :                 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   335089672 :                 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   335089672 :                 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   335089672 :                 d2xidxyz2[p][2]*dphidxi[i][p] +                                                      // xi_{x z} * phi_{xi}
     433   335089672 :                 d2etadxyz2[p][2]*dphideta[i][p] +                                                    // eta_{x z} * phi_{eta}
     434   309526466 :                 d2zetadxyz2[p][2]*dphidzeta[i][p];                                                   // zeta_{x z} * phi_{zeta}
     435             : 
     436             :               // phi_{y y}
     437   309526466 :               d2phi[i][p].slice(1).slice(1) = d2phidy2[i][p] =
     438   335089672 :                 d2phidxi2[i][p]*dxidy_map[p]*dxidy_map[p] +           // (xi_y)^2 * phi_{xi xi}
     439   335089672 :                 d2phideta2[i][p]*detady_map[p]*detady_map[p] +        // (eta_y)^2 * phi_{eta eta}
     440   335089672 :                 d2phidzeta2[i][p]*dzetady_map[p]*dzetady_map[p] +     // (zeta_y)^2 * phi_{zeta zeta}
     441   335089672 :                 2*d2phidxideta[i][p]*dxidy_map[p]*detady_map[p] +     // 2 * xi_y * eta_y * phi_{xi eta}
     442   335089672 :                 2*d2phidxidzeta[i][p]*dxidy_map[p]*dzetady_map[p] +   // 2 * xi_y * zeta_y * phi_{xi zeta}
     443   335089672 :                 2*d2phidetadzeta[i][p]*detady_map[p]*dzetady_map[p] + // 2 * eta_y * zeta_y * phi_{eta zeta}
     444   335089672 :                 d2xidxyz2[p][3]*dphidxi[i][p] +                       // xi_{y y} * phi_{xi}
     445   335089672 :                 d2etadxyz2[p][3]*dphideta[i][p] +                     // eta_{y y} * phi_{eta}
     446   309526466 :                 d2zetadxyz2[p][3]*dphidzeta[i][p];                    // zeta_{y y} * phi_{zeta}
     447             : 
     448             :               // phi_{y z}
     449   309526466 :               d2phi[i][p].slice(1).slice(2) = d2phi[i][p].slice(2).slice(1) = d2phidydz[i][p] =
     450   335089672 :                 d2phidxi2[i][p]*dxidy_map[p]*dxidz_map[p] +                                          // xi_y * xi_z * phi_{xi xi}
     451   335089672 :                 d2phideta2[i][p]*detady_map[p]*detadz_map[p] +                                       // eta_y * eta_z * phi_{eta eta}
     452   335089672 :                 d2phidzeta2[i][p]*dzetady_map[p]*dzetadz_map[p] +                                    // zeta_y * zeta_z * phi_{zeta zeta}
     453   335089672 :                 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   335089672 :                 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   335089672 :                 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   335089672 :                 d2xidxyz2[p][4]*dphidxi[i][p] +                                                      // xi_{y z} * phi_{xi}
     457   335089672 :                 d2etadxyz2[p][4]*dphideta[i][p] +                                                    // eta_{y z} * phi_{eta}
     458   309526466 :                 d2zetadxyz2[p][4]*dphidzeta[i][p];                                                   // zeta_{y z} * phi_{zeta}
     459             : 
     460             :               // phi_{z z}
     461   335089672 :               d2phi[i][p].slice(2).slice(2) = d2phidz2[i][p] =
     462   335089672 :                 d2phidxi2[i][p]*dxidz_map[p]*dxidz_map[p] +           // (xi_z)^2 * phi_{xi xi}
     463   335089672 :                 d2phideta2[i][p]*detadz_map[p]*detadz_map[p] +        // (eta_z)^2 * phi_{eta eta}
     464   335089672 :                 d2phidzeta2[i][p]*dzetadz_map[p]*dzetadz_map[p] +     // (zeta_z)^2 * phi_{zeta zeta}
     465   335089672 :                 2*d2phidxideta[i][p]*dxidz_map[p]*detadz_map[p] +     // 2 * xi_z * eta_z * phi_{xi eta}
     466   335089672 :                 2*d2phidxidzeta[i][p]*dxidz_map[p]*dzetadz_map[p] +   // 2 * xi_z * zeta_z * phi_{xi zeta}
     467   335089672 :                 2*d2phidetadzeta[i][p]*detadz_map[p]*dzetadz_map[p] + // 2 * eta_z * zeta_z * phi_{eta zeta}
     468   335089672 :                 d2xidxyz2[p][5]*dphidxi[i][p] +                       // xi_{z z} * phi_{xi}
     469   335089672 :                 d2etadxyz2[p][5]*dphideta[i][p] +                     // eta_{z z} * phi_{eta}
     470   309526466 :                 d2zetadxyz2[p][5]*dphidzeta[i][p];                    // zeta_{z z} * phi_{zeta}
     471             :             }
     472             : 
     473      432120 :         break;
     474             :       }
     475             : 
     476           0 :     default:
     477           0 :       libmesh_error_msg("Invalid dim = " << dim);
     478             :     } // switch(dim)
     479     6568893 : }
     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      186876 : 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      186876 :   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      735408 :         for (auto i : index_range(curl_phi))
     525     7506480 :           for (auto p : index_range(curl_phi[i]))
     526             :             {
     527             : 
     528     7990048 :               Real dphiy_dx = (dphidxi[i][p].slice(1))*dxidx_map[p]
     529     7990048 :                 + (dphideta[i][p].slice(1))*detadx_map[p];
     530             : 
     531     6845856 :               Real dphix_dy = (dphidxi[i][p].slice(0))*dxidy_map[p]
     532     6845856 :                 + (dphideta[i][p].slice(0))*detady_map[p];
     533             : 
     534     6845856 :               curl_phi[i][p].slice(2) = dphiy_dx - dphix_dy;
     535             : 
     536             : #if LIBMESH_DIM > 2
     537     6845856 :               Real dphiy_dz = (dphidxi[i][p].slice(1))*dxidz_map[p]
     538     6845856 :                 + (dphideta[i][p].slice(1))*detadz_map[p];
     539             : 
     540     6845856 :               Real dphix_dz = (dphidxi[i][p].slice(0))*dxidz_map[p]
     541     6845856 :                 + (dphideta[i][p].slice(0))*detadz_map[p];
     542             : 
     543     6845856 :               curl_phi[i][p].slice(0) = -dphiy_dz;
     544     6845856 :               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     1475196 :         for (auto i : index_range(curl_phi))
     570    22246560 :           for (auto p : index_range(curl_phi[i]))
     571             :             {
     572    24364032 :               Real dphiz_dy = (dphidxi[i][p].slice(2))*dxidy_map[p]
     573    24364032 :                 + (dphideta[i][p].slice(2))*detady_map[p]
     574    24364032 :                 + (dphidzeta[i][p].slice(2))*dzetady_map[p];
     575             : 
     576    20883456 :               Real dphiy_dz = (dphidxi[i][p].slice(1))*dxidz_map[p]
     577    20883456 :                 + (dphideta[i][p].slice(1))*detadz_map[p]
     578    20883456 :                 + (dphidzeta[i][p].slice(1))*dzetadz_map[p];
     579             : 
     580    20883456 :               Real dphix_dz = (dphidxi[i][p].slice(0))*dxidz_map[p]
     581    20883456 :                 + (dphideta[i][p].slice(0))*detadz_map[p]
     582    20883456 :                 + (dphidzeta[i][p].slice(0))*dzetadz_map[p];
     583             : 
     584    20883456 :               Real dphiz_dx = (dphidxi[i][p].slice(2))*dxidx_map[p]
     585    20883456 :                 + (dphideta[i][p].slice(2))*detadx_map[p]
     586    20883456 :                 + (dphidzeta[i][p].slice(2))*dzetadx_map[p];
     587             : 
     588    20883456 :               Real dphiy_dx = (dphidxi[i][p].slice(1))*dxidx_map[p]
     589    20883456 :                 + (dphideta[i][p].slice(1))*detadx_map[p]
     590    20883456 :                 + (dphidzeta[i][p].slice(1))*dzetadx_map[p];
     591             : 
     592    20883456 :               Real dphix_dy = (dphidxi[i][p].slice(0))*dxidy_map[p]
     593    20883456 :                 + (dphideta[i][p].slice(0))*detady_map[p]
     594    20883456 :                 + (dphidzeta[i][p].slice(0))*dzetady_map[p];
     595             : 
     596    20883456 :               curl_phi[i][p].slice(0) = dphiz_dy - dphiy_dz;
     597             : 
     598    20883456 :               curl_phi[i][p].slice(1) = dphix_dz - dphiz_dx;
     599             : 
     600    20883456 :               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      186876 : }
     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      570380 : 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      570380 :   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     2295488 :         for (auto i : index_range(div_phi))
     648    18878448 :           for (auto p : index_range(div_phi[i]))
     649             :             {
     650    19638592 :               Real dphix_dx = (dphidxi[i][p].slice(0))*dxidx_map[p]
     651    19638592 :                 + (dphideta[i][p].slice(0))*detadx_map[p];
     652             : 
     653    16820064 :               Real dphiy_dy = (dphidxi[i][p].slice(1))*dxidy_map[p]
     654    16820064 :                 + (dphideta[i][p].slice(1))*detady_map[p];
     655             : 
     656    18229328 :               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     4350588 :         for (auto i : index_range(div_phi))
     680    38171808 :           for (auto p : index_range(div_phi[i]))
     681             :             {
     682    39846912 :               Real dphix_dx = (dphidxi[i][p].slice(0))*dxidx_map[p]
     683    39846912 :                 + (dphideta[i][p].slice(0))*detadx_map[p]
     684    39846912 :                 + (dphidzeta[i][p].slice(0))*dzetadx_map[p];
     685             : 
     686    34154496 :               Real dphiy_dy = (dphidxi[i][p].slice(1))*dxidy_map[p]
     687    34154496 :                 + (dphideta[i][p].slice(1))*detady_map[p]
     688    34154496 :                 + (dphidzeta[i][p].slice(1))*dzetady_map[p];
     689             : 
     690    34154496 :               Real dphiz_dz = (dphidxi[i][p].slice(2))*dxidz_map[p]
     691    34154496 :                 + (dphideta[i][p].slice(2))*detadz_map[p]
     692    34154496 :                 + (dphidzeta[i][p].slice(2))*dzetadz_map[p];
     693             : 
     694    37000704 :               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      570380 :   return;
     705             : }
     706             : 
     707             : // Explicit Instantiations
     708             : template class LIBMESH_EXPORT H1FETransformation<Real>;
     709             : template class LIBMESH_EXPORT H1FETransformation<RealGradient>;
     710             : 
     711             : } //namespace libMesh