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/hdiv_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     4648324 : void HDivFETransformation<OutputShape>::init_map_phi(const FEGenericBase<OutputShape> & fe) const
      27             : {
      28      388054 :   fe.get_fe_map().get_dxidx();
      29     4648324 : }
      30             : 
      31             : 
      32             : 
      33             : template<typename OutputShape>
      34     1799364 : void HDivFETransformation<OutputShape>::init_map_dphi(const FEGenericBase<OutputShape> & fe) const
      35             : {
      36      149928 :   fe.get_fe_map().get_dxidx();
      37     1799364 : }
      38             : 
      39             : 
      40             : 
      41             : template<typename OutputShape>
      42           0 : void HDivFETransformation<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     1815188 : void HDivFETransformation<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     1815188 :   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      509270 :     case 2:
      67             :       {
      68       42574 :         const std::vector<RealGradient> & dxyz_dxi   = fe.get_fe_map().get_dxyzdxi();
      69       42574 :         const std::vector<RealGradient> & dxyz_deta  = fe.get_fe_map().get_dxyzdeta();
      70             : 
      71       42574 :         const std::vector<Real> & J = fe.get_fe_map().get_jacobian();
      72             : 
      73             :         // phi = J^{-1} * (dx/dxi) * \hat{phi}
      74     3472402 :         for (auto i : index_range(phi))
      75    28148692 :           for (auto p : index_range(phi[i]))
      76             :             {
      77    25185560 :               Real dx_dxi   = dxyz_dxi[p](0);
      78    25185560 :               Real dx_deta  = dxyz_deta[p](0);
      79             : 
      80    25185560 :               Real dy_dxi   = dxyz_dxi[p](1);
      81    25185560 :               Real dy_deta  = dxyz_deta[p](1);
      82             : 
      83    25185560 :               Real dz_dxi   = dxyz_dxi[p](2);
      84    25185560 :               Real dz_deta  = dxyz_deta[p](2);
      85             : 
      86             :               // Need to temporarily cache reference shape functions
      87             :               // We are computing mapping basis functions, so we explicitly ignore
      88             :               // any non-zero p_level() the Elem might have.
      89     2096882 :               OutputShape phi_ref;
      90    29379324 :               FEInterface::shape(fe.get_fe_type(), /*extra_order=*/0, elem, i, qp[p], phi_ref);
      91             : 
      92    29379324 :               phi[i][p](0) = (dx_dxi*phi_ref(0) + dx_deta*phi_ref(1))/J[p];
      93    25185560 :               phi[i][p](1) = (dy_dxi*phi_ref(0) + dy_deta*phi_ref(1))/J[p];
      94    25185560 :               phi[i][p](2) = (dz_dxi*phi_ref(0) + dz_deta*phi_ref(1))/J[p];
      95             :             }
      96             : 
      97       42574 :         break;
      98             :       }
      99     1305918 :     case 3:
     100             :       {
     101      108847 :         const std::vector<RealGradient> & dxyz_dxi   = fe.get_fe_map().get_dxyzdxi();
     102      108847 :         const std::vector<RealGradient> & dxyz_deta  = fe.get_fe_map().get_dxyzdeta();
     103      108847 :         const std::vector<RealGradient> & dxyz_dzeta = fe.get_fe_map().get_dxyzdzeta();
     104             : 
     105      108847 :         const std::vector<Real> & J = fe.get_fe_map().get_jacobian();
     106             : 
     107             :         // phi = J^{-1} * (dx/dxi) * \hat{phi}
     108     7369698 :         for (auto i : index_range(phi))
     109    88092162 :           for (auto p : index_range(phi[i]))
     110             :             {
     111    82028382 :               Real dx_dxi   = dxyz_dxi[p](0);
     112    82028382 :               Real dx_deta  = dxyz_deta[p](0);
     113    82028382 :               Real dx_dzeta = dxyz_dzeta[p](0);
     114             : 
     115    82028382 :               Real dy_dxi   = dxyz_dxi[p](1);
     116    82028382 :               Real dy_deta  = dxyz_deta[p](1);
     117    82028382 :               Real dy_dzeta = dxyz_dzeta[p](1);
     118             : 
     119    82028382 :               Real dz_dxi   = dxyz_dxi[p](2);
     120    82028382 :               Real dz_deta  = dxyz_deta[p](2);
     121    82028382 :               Real dz_dzeta = dxyz_dzeta[p](2);
     122             : 
     123             :               // Need to temporarily cache reference shape functions
     124             :               // We are computing mapping basis functions, so we explicitly ignore
     125             :               // any non-zero p_level() the Elem might have.
     126     6836142 :               OutputShape phi_ref;
     127    95700666 :               FEInterface::shape(fe.get_fe_type(), /*extra_order=*/0, elem, i, qp[p], phi_ref);
     128             : 
     129    95700666 :               phi[i][p](0) = (dx_dxi*phi_ref(0) + dx_deta*phi_ref(1) + dx_dzeta*phi_ref(2))/J[p];
     130    82028382 :               phi[i][p](1) = (dy_dxi*phi_ref(0) + dy_deta*phi_ref(1) + dy_dzeta*phi_ref(2))/J[p];
     131    82028382 :               phi[i][p](2) = (dz_dxi*phi_ref(0) + dz_deta*phi_ref(1) + dz_dzeta*phi_ref(2))/J[p];
     132             :             }
     133             : 
     134      108847 :         break;
     135             :       }
     136             : 
     137           0 :     default:
     138           0 :       libmesh_error_msg("Invalid dim = " << dim);
     139             :     } // switch(dim)
     140     1815188 : }
     141             : 
     142             : template<typename OutputShape>
     143      882516 : void HDivFETransformation<OutputShape>::map_div(const unsigned int dim,
     144             :                                                 const Elem * const,
     145             :                                                 const std::vector<Point> &,
     146             :                                                 const FEGenericBase<OutputShape> & fe,
     147             :                                                 std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputDivergence>> & div_phi) const
     148             : {
     149      882516 :   switch (dim)
     150             :     {
     151           0 :     case 0:
     152             :     case 1:
     153           0 :       libmesh_error_msg("These element transformations only make sense in 2D and 3D.");
     154             : 
     155       18180 :     case 2:
     156             :       {
     157       18180 :         const std::vector<std::vector<OutputShape>> & dphi_dxi = fe.get_dphidxi();
     158       18180 :         const std::vector<std::vector<OutputShape>> & dphi_deta = fe.get_dphideta();
     159             : 
     160       18180 :         const std::vector<Real> & J = fe.get_fe_map().get_jacobian();
     161             : 
     162             :         // div(phi) = J^{-1} * div(\hat{phi})
     163     1683222 :         for (auto i : index_range(div_phi))
     164    18605052 :           for (auto p : index_range(div_phi[i]))
     165             :             {
     166    25688016 :               div_phi[i][p] = (dphi_dxi[i][p](0) + dphi_deta[i][p](1))/J[p];
     167             :             }
     168             : 
     169       18180 :         break;
     170             :       }
     171       55344 :     case 3:
     172             :       {
     173       55344 :         const std::vector<std::vector<OutputShape>> & dphi_dxi = fe.get_dphidxi();
     174       55344 :         const std::vector<std::vector<OutputShape>> & dphi_deta = fe.get_dphideta();
     175       55344 :         const std::vector<std::vector<OutputShape>> & dphi_dzeta = fe.get_dphidzeta();
     176             : 
     177       55344 :         const std::vector<Real> & J = fe.get_fe_map().get_jacobian();
     178             : 
     179             :         // div(phi) = J^{-1} * div(\hat{phi})
     180     3792822 :         for (auto i : index_range(div_phi))
     181    45283410 :           for (auto p : index_range(div_phi[i]))
     182             :             {
     183    70247814 :               div_phi[i][p] = (dphi_dxi[i][p](0) + dphi_deta[i][p](1) + dphi_dzeta[i][p](2))/J[p];
     184             :             }
     185             : 
     186       55344 :         break;
     187             :       }
     188             : 
     189           0 :     default:
     190           0 :       libmesh_error_msg("Invalid dim = " << dim);
     191             :     } // switch(dim)
     192      882516 : }
     193             : 
     194             : template class LIBMESH_EXPORT HDivFETransformation<RealGradient>;
     195             : 
     196             : template<>
     197           0 : void HDivFETransformation<Real>::init_map_phi(const FEGenericBase<Real> & ) const
     198             : {
     199           0 :   libmesh_error_msg("HDiv transformations only make sense for vector-valued elements.");
     200             : }
     201             : 
     202             : template<>
     203           0 : void HDivFETransformation<Real>::init_map_dphi(const FEGenericBase<Real> & ) const
     204             : {
     205           0 :   libmesh_error_msg("HDiv transformations only make sense for vector-valued elements.");
     206             : }
     207             : 
     208             : template<>
     209           0 : void HDivFETransformation<Real>::init_map_d2phi(const FEGenericBase<Real> & ) const
     210             : {
     211           0 :   libmesh_error_msg("HDiv transformations only make sense for vector-valued elements.");
     212             : }
     213             : 
     214             : template<>
     215           0 : void HDivFETransformation<Real>::map_phi(const unsigned int,
     216             :                                          const Elem * const,
     217             :                                          const std::vector<Point> &,
     218             :                                          const FEGenericBase<Real> &,
     219             :                                          std::vector<std::vector<Real>> &,
     220             :                                          bool) const
     221             : {
     222           0 :   libmesh_error_msg("HDiv transformations only make sense for vector-valued elements.");
     223             : }
     224             : 
     225             : template<>
     226           0 : void HDivFETransformation<Real>::map_div(const unsigned int,
     227             :                                          const Elem * const,
     228             :                                          const std::vector<Point> &,
     229             :                                          const FEGenericBase<Real> &,
     230             :                                          std::vector<std::vector<Real>> &) const
     231             : {
     232           0 :   libmesh_error_msg("HDiv transformations only make sense for vector-valued elements.");
     233             : }
     234             : 
     235             : 
     236             : } // namespace libMesh