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/hcurl_fe_transformation.h"
      19             : #include "libmesh/fe_interface.h"
      20             : #include "libmesh/int_range.h"
      21             : 
      22             : namespace libMesh
      23             : {
      24             : 
      25             : template<typename OutputShape>
      26      879482 : void HCurlFETransformation<OutputShape>::init_map_phi(const FEGenericBase<OutputShape> & fe) const
      27             : {
      28       74160 :   fe.get_fe_map().get_dxidx();
      29      879482 : }
      30             : 
      31             : 
      32             : 
      33             : template<typename OutputShape>
      34      539406 : void HCurlFETransformation<OutputShape>::init_map_dphi(const FEGenericBase<OutputShape> & fe) const
      35             : {
      36       45444 :   fe.get_fe_map().get_dxidx();
      37      539406 : }
      38             : 
      39             : 
      40             : 
      41             : template<typename OutputShape>
      42           0 : void HCurlFETransformation<OutputShape>::init_map_d2phi(const FEGenericBase<OutputShape> & fe) const
      43             : {
      44           0 :   fe.get_fe_map().get_dxidx();
      45             : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
      46           0 :   fe.get_fe_map().get_d2xidxyz2();
      47             : #endif
      48           0 : }
      49             : 
      50             : 
      51             : 
      52             : template<typename OutputShape>
      53      381956 : void HCurlFETransformation<OutputShape>::map_phi(const unsigned int dim,
      54             :                                                  const Elem * const elem,
      55             :                                                  const std::vector<Point> & qp,
      56             :                                                  const FEGenericBase<OutputShape> & fe,
      57             :                                                  std::vector<std::vector<OutputShape>> & phi,
      58             :                                                  const bool /*add_p_level*/) const
      59             : {
      60      381956 :   switch (dim)
      61             :     {
      62           0 :     case 0:
      63             :     case 1:
      64           0 :       libmesh_error_msg("These element transformations only make sense in 2D and 3D.");
      65             : 
      66      245698 :     case 2:
      67             :       {
      68       20470 :         const std::vector<Real> & dxidx_map = fe.get_fe_map().get_dxidx();
      69       20470 :         const std::vector<Real> & dxidy_map = fe.get_fe_map().get_dxidy();
      70       20470 :         const std::vector<Real> & dxidz_map = fe.get_fe_map().get_dxidz();
      71             : 
      72       20470 :         const std::vector<Real> & detadx_map = fe.get_fe_map().get_detadx();
      73       20470 :         const std::vector<Real> & detady_map = fe.get_fe_map().get_detady();
      74       20470 :         const std::vector<Real> & detadz_map = fe.get_fe_map().get_detadz();
      75             : 
      76             :         // phi = (dx/dxi)^-T * \hat{phi}
      77             :         // In 2D (embedded in 3D):
      78             :         // (dx/dxi)^{-1} = [ dxi/dx  dxi/dy  dxi/dz
      79             :         //                   deta/dx deta/dy deta/dz ]
      80             :         //
      81             :         // so: dxi/dx^{-T} * \hat{phi} = [ dxi/dx deta/dx   [ \hat{phi}_xi
      82             :         //                                 dxi/dy deta/dy     \hat{phi}_eta ]
      83             :         //                                 dxi/dz deta/dz ]
      84             :         //
      85             :         // or in indicial notation:  phi_j = xi_{i,j}*\hat{phi}_i
      86     1751402 :         for (auto i : index_range(phi))
      87    13123912 :           for (auto p : index_range(phi[i]))
      88             :             {
      89             :               // Need to temporarily cache reference shape functions
      90             :               // TODO: PB: Might be worth trying to build phi_ref separately to see
      91             :               //           if we can get vectorization
      92             :               // We are computing mapping basis functions, so we explicitly ignore
      93             :               // any non-zero p_level() the Elem might have.
      94      968734 :               OutputShape phi_ref;
      95    13555676 :               FEInterface::shape(fe.get_fe_type(), /*extra_order=*/0, elem, i, qp[p], phi_ref);
      96             : 
      97    14524410 :               phi[i][p](0) = dxidx_map[p]*phi_ref.slice(0) + detadx_map[p]*phi_ref.slice(1);
      98             : 
      99    12586942 :               phi[i][p](1) = dxidy_map[p]*phi_ref.slice(0) + detady_map[p]*phi_ref.slice(1);
     100             : 
     101    13555676 :               phi[i][p](2) = dxidz_map[p]*phi_ref.slice(0) + detadz_map[p]*phi_ref.slice(1);
     102             :             }
     103             : 
     104       20470 :         break;
     105             :       }
     106             : 
     107      136258 :     case 3:
     108             :       {
     109       11668 :         const std::vector<Real> & dxidx_map = fe.get_fe_map().get_dxidx();
     110       11668 :         const std::vector<Real> & dxidy_map = fe.get_fe_map().get_dxidy();
     111       11668 :         const std::vector<Real> & dxidz_map = fe.get_fe_map().get_dxidz();
     112             : 
     113       11668 :         const std::vector<Real> & detadx_map = fe.get_fe_map().get_detadx();
     114       11668 :         const std::vector<Real> & detady_map = fe.get_fe_map().get_detady();
     115       11668 :         const std::vector<Real> & detadz_map = fe.get_fe_map().get_detadz();
     116             : 
     117       11668 :         const std::vector<Real> & dzetadx_map = fe.get_fe_map().get_dzetadx();
     118       11668 :         const std::vector<Real> & dzetady_map = fe.get_fe_map().get_dzetady();
     119       11668 :         const std::vector<Real> & dzetadz_map = fe.get_fe_map().get_dzetadz();
     120             : 
     121             :         // phi = (dx/dxi)^-T * \hat{phi}
     122             :         // In 3D:
     123             :         // dx/dxi^-1 = [  dxi/dx    dxi/dy    dxi/dz
     124             :         // deta/dx   deta/dy   deta/dz
     125             :         // dzeta/dx  dzeta/dy  dzeta/dz]
     126             :         //
     127             :         // so: dxi/dx^-T * \hat{phi} = [ dxi/dx  deta/dx  dzeta/dx   [ \hat{phi}_xi
     128             :         // dxi/dy  deta/dy  dzeta/dy     \hat{phi}_eta
     129             :         // dxi/dz  deta/dz  dzeta/dz ]   \hat{phi}_zeta ]
     130             :         //
     131             :         // or in indicial notation:  phi_j = xi_{i,j}*\hat{phi}_i
     132     1020250 :         for (auto i : index_range(phi))
     133    11830944 :           for (auto p : index_range(phi[i]))
     134             :             {
     135             :               // Need to temporarily cache reference shape functions
     136             :               // TODO: PB: Might be worth trying to build phi_ref separately to see
     137             :               //           if we can get vectorization
     138             :               // We are computing mapping basis functions, so we explicitly ignore
     139             :               // any non-zero p_level() the Elem might have.
     140      927312 :               OutputShape phi_ref;
     141    12801576 :               FEInterface::shape(fe.get_fe_type(), /*extra_order=*/0, elem, i, qp[p], phi_ref);
     142             : 
     143    13728888 :               phi[i][p].slice(0) = dxidx_map[p]*phi_ref.slice(0) + detadx_map[p]*phi_ref.slice(1)
     144    11874264 :                 + dzetadx_map[p]*phi_ref.slice(2);
     145             : 
     146    12801576 :               phi[i][p].slice(1) = dxidy_map[p]*phi_ref.slice(0) + detady_map[p]*phi_ref.slice(1)
     147    10946952 :                 + dzetady_map[p]*phi_ref.slice(2);
     148             : 
     149    12801576 :               phi[i][p].slice(2) = dxidz_map[p]*phi_ref.slice(0) + detadz_map[p]*phi_ref.slice(1)
     150    11874264 :                 + dzetadz_map[p]*phi_ref.slice(2);
     151             :             }
     152       11668 :         break;
     153             :       }
     154             : 
     155           0 :     default:
     156           0 :       libmesh_error_msg("Invalid dim = " << dim);
     157             :     } // switch(dim)
     158      381956 : }
     159             : 
     160             : template<typename OutputShape>
     161      255694 : void HCurlFETransformation<OutputShape>::map_curl(const unsigned int dim,
     162             :                                                   const Elem * const,
     163             :                                                   const std::vector<Point> &,
     164             :                                                   const FEGenericBase<OutputShape> & fe,
     165             :                                                   std::vector<std::vector<OutputShape>> & curl_phi) const
     166             : {
     167      255694 :   switch (dim)
     168             :     {
     169           0 :     case 0:
     170             :     case 1:
     171           0 :       libmesh_error_msg("These element transformations only make sense in 2D and 3D.");
     172             : 
     173      165900 :     case 2:
     174             :       {
     175       13824 :         const std::vector<Real> & dxidx_map = fe.get_fe_map().get_dxidx();
     176       13824 :         const std::vector<Real> & dxidy_map = fe.get_fe_map().get_dxidy();
     177       13824 :         const std::vector<Real> & dxidz_map = fe.get_fe_map().get_dxidz();
     178             : 
     179       13824 :         const std::vector<Real> & detadx_map = fe.get_fe_map().get_detadx();
     180       13824 :         const std::vector<Real> & detady_map = fe.get_fe_map().get_detady();
     181       13824 :         const std::vector<Real> & detadz_map = fe.get_fe_map().get_detadz();
     182             : 
     183       13824 :         const std::vector<std::vector<OutputShape>> & dphi_dxi = fe.get_dphidxi();
     184       13824 :         const std::vector<std::vector<OutputShape>> & dphi_deta = fe.get_dphideta();
     185             : 
     186             :         // In 2D (embedded in 3D): curl(phi)_j = C(j) * curl(\hat{phi}) where C(j)
     187             :         // is the jth column cofactor of (dx/dxi)^{-1} = [ dxi/dx  dxi/dy  dxi/dz
     188             :         //                                                 deta/dx deta/dy deta/dz ]
     189     1175820 :         for (auto i : index_range(curl_phi))
     190     9476264 :           for (auto p : index_range(curl_phi[i]))
     191             :             {
     192    10584626 :               const Real curl_ref = dphi_dxi[i][p].slice(1) - dphi_deta[i][p].slice(0);
     193    11290720 :               curl_phi[i][p].slice(0) = curl_ref *  (dxidy_map[p]*detadz_map[p] - dxidz_map[p]*detady_map[p]);
     194     9172438 :               curl_phi[i][p].slice(1) = curl_ref * -(dxidx_map[p]*detadz_map[p] - dxidz_map[p]*detadx_map[p]);
     195     8466344 :               curl_phi[i][p].slice(2) = curl_ref *  (dxidx_map[p]*detady_map[p] - dxidy_map[p]*detadx_map[p]);
     196             :             }
     197             : 
     198       13824 :         break;
     199             :       }
     200             : 
     201        7700 :     case 3:
     202             :       {
     203        7700 :         const std::vector<std::vector<OutputShape>> & dphi_dxi = fe.get_dphidxi();
     204        7700 :         const std::vector<std::vector<OutputShape>> & dphi_deta = fe.get_dphideta();
     205        7700 :         const std::vector<std::vector<OutputShape>> & dphi_dzeta = fe.get_dphidzeta();
     206             : 
     207        7700 :         const std::vector<RealGradient> & dxyz_dxi   = fe.get_fe_map().get_dxyzdxi();
     208        7700 :         const std::vector<RealGradient> & dxyz_deta  = fe.get_fe_map().get_dxyzdeta();
     209        7700 :         const std::vector<RealGradient> & dxyz_dzeta = fe.get_fe_map().get_dxyzdzeta();
     210             : 
     211        7700 :         const std::vector<Real> & J = fe.get_fe_map().get_jacobian();
     212             : 
     213      632794 :         for (auto i : index_range(curl_phi))
     214     8255136 :           for (auto p : index_range(curl_phi[i]))
     215             :             {
     216     7712136 :               Real dx_dxi   = dxyz_dxi[p](0);
     217     7712136 :               Real dx_deta  = dxyz_deta[p](0);
     218     7712136 :               Real dx_dzeta = dxyz_dzeta[p](0);
     219             : 
     220     7712136 :               Real dy_dxi   = dxyz_dxi[p](1);
     221     7712136 :               Real dy_deta  = dxyz_deta[p](1);
     222     7712136 :               Real dy_dzeta = dxyz_dzeta[p](1);
     223             : 
     224     7712136 :               Real dz_dxi   = dxyz_dxi[p](2);
     225     7712136 :               Real dz_deta  = dxyz_deta[p](2);
     226     7712136 :               Real dz_dzeta = dxyz_dzeta[p](2);
     227             : 
     228     7712136 :               const Real inv_jac = 1.0/J[p];
     229             : 
     230             :               /* In 3D: curl(phi) = J^{-1} dx/dxi * curl(\hat{phi})
     231             : 
     232             :                  dx/dxi = [  dx/dxi  dx/deta  dx/dzeta
     233             :                  dy/dxi  dy/deta  dy/dzeta
     234             :                  dz/dxi  dz/deta  dz/dzeta ]
     235             : 
     236             :                  curl(u) = [ du_z/deta  - du_y/dzeta
     237             :                  du_x/dzeta - du_z/dxi
     238             :                  du_y/dxi   - du_x/deta ]
     239             :               */
     240     9666936 :               curl_phi[i][p].slice(0) = inv_jac*( dx_dxi*( dphi_deta[i][p].slice(2)  -
     241     8363736 :                                                            dphi_dzeta[i][p].slice(1)   ) +
     242     9015336 :                                                   dx_deta*( dphi_dzeta[i][p].slice(0) -
     243     8363736 :                                                             dphi_dxi[i][p].slice(2)     ) +
     244     8363736 :                                                   dx_dzeta*( dphi_dxi[i][p].slice(1) -
     245      651600 :                                                              dphi_deta[i][p].slice(0)    ) );
     246             : 
     247     9666936 :               curl_phi[i][p].slice(1) = inv_jac*( dy_dxi*( dphi_deta[i][p].slice(2) -
     248     7712136 :                                                            dphi_dzeta[i][p].slice(1)  ) +
     249     8363736 :                                                   dy_deta*( dphi_dzeta[i][p].slice(0)-
     250     7712136 :                                                             dphi_dxi[i][p].slice(2)    ) +
     251     8363736 :                                                   dy_dzeta*( dphi_dxi[i][p].slice(1) -
     252      651600 :                                                              dphi_deta[i][p].slice(0)   ) );
     253             : 
     254     9666936 :               curl_phi[i][p].slice(2) = inv_jac*( dz_dxi*( dphi_deta[i][p].slice(2) -
     255     7712136 :                                                            dphi_dzeta[i][p].slice(1)   ) +
     256     8363736 :                                                   dz_deta*( dphi_dzeta[i][p].slice(0) -
     257     7712136 :                                                             dphi_dxi[i][p].slice(2)     ) +
     258     8363736 :                                                   dz_dzeta*( dphi_dxi[i][p].slice(1) -
     259      651600 :                                                              dphi_deta[i][p].slice(0)    ) );
     260             :             }
     261             : 
     262        7700 :         break;
     263             :       }
     264             : 
     265           0 :     default:
     266           0 :       libmesh_error_msg("Invalid dim = " << dim);
     267             :     } // switch(dim)
     268      255694 : }
     269             : 
     270             : template class LIBMESH_EXPORT HCurlFETransformation<RealGradient>;
     271             : 
     272             : template<>
     273           0 : void HCurlFETransformation<Real>::init_map_phi(const FEGenericBase<Real> & ) const
     274             : {
     275           0 :   libmesh_error_msg("HCurl transformations only make sense for vector-valued elements.");
     276             : }
     277             : 
     278             : template<>
     279           0 : void HCurlFETransformation<Real>::init_map_dphi(const FEGenericBase<Real> & ) const
     280             : {
     281           0 :   libmesh_error_msg("HCurl transformations only make sense for vector-valued elements.");
     282             : }
     283             : 
     284             : template<>
     285           0 : void HCurlFETransformation<Real>::init_map_d2phi(const FEGenericBase<Real> & ) const
     286             : {
     287           0 :   libmesh_error_msg("HCurl transformations only make sense for vector-valued elements.");
     288             : }
     289             : 
     290             : template<>
     291           0 : void HCurlFETransformation<Real>::map_phi(const unsigned int,
     292             :                                           const Elem * const,
     293             :                                           const std::vector<Point> &,
     294             :                                           const FEGenericBase<Real> &,
     295             :                                           std::vector<std::vector<Real>> &,
     296             :                                           bool) const
     297             : {
     298           0 :   libmesh_error_msg("HCurl transformations only make sense for vector-valued elements.");
     299             : }
     300             : 
     301             : template<>
     302           0 : void HCurlFETransformation<Real>::map_curl(const unsigned int,
     303             :                                            const Elem * const,
     304             :                                            const std::vector<Point> &,
     305             :                                            const FEGenericBase<Real> &,
     306             :                                            std::vector<std::vector<Real>> &) const
     307             : {
     308           0 :   libmesh_error_msg("HCurl transformations only make sense for vector-valued elements.");
     309             : }
     310             : 
     311             : 
     312             : } // namespace libMesh