Branch data     Line data    Source code

       1                 :            : // *****************************************************************************
       2                 :            : /*!
       3                 :            :   \file      src/PDE/Integrate/Mass.cpp
       4                 :            :   \copyright 2012-2015 J. Bakosi,
       5                 :            :              2016-2018 Los Alamos National Security, LLC.,
       6                 :            :              2019-2021 Triad National Security, LLC.
       7                 :            :              All rights reserved. See the LICENSE file for details.
       8                 :            :   \brief     Functions for computing the mass matrix for a system of PDEs
       9                 :            :   \details   This file contains functionality for computing the mass matrix for
      10                 :            :      a system of PDEs used in discontinuous and continuous Galerkin methods for
      11                 :            :      various orders of numerical representation.
      12                 :            : */
      13                 :            : // *****************************************************************************
      14                 :            : 
      15                 :            : #include "Mass.hpp"
      16                 :            : #include "Vector.hpp"
      17                 :            : 
      18                 :            : void
      19                 :        933 : tk::mass( ncomp_t ncomp,
      20                 :            :           ncomp_t ndof,
      21                 :            :           const Fields& geoElem,
      22                 :            :           Fields& l )
      23                 :            : // *****************************************************************************
      24                 :            : //  Compute the block-diagonal mass matrix for DG
      25                 :            : //! \param[in] ncomp Number of scalar components in this PDE system
      26                 :            : //! \param[in] ndof Number of solution degrees of freedom
      27                 :            : //! \param[in] geoElem Element geometry array
      28                 :            : //! \param[in,out] l Block diagonal mass matrix
      29                 :            : // *****************************************************************************
      30                 :            : {
      31                 :            :   Assert( geoElem.nunk() == l.nunk(), "Size mismatch" );
      32                 :        933 :   const auto nelem = geoElem.nunk();
      33                 :            :   Assert( l.nprop() == ndof*ncomp, "lhs mass matrix size incorrect" );
      34                 :            : 
      35                 :            :   // Compute LHS
      36         [ +  + ]:     628173 :   for (std::size_t e=0; e<nelem; ++e) {
      37                 :     627240 :     auto M = massMatrixDubiner(ndof, geoElem(e,0));
      38         [ +  + ]:    2074316 :     for (ncomp_t c=0; c<ncomp; ++c) {
      39                 :    1447076 :       const auto mark = c * ndof;
      40         [ +  + ]:    6582280 :       for (std::size_t i=0; i<ndof; ++i) {
      41                 :    5135204 :         l(e, mark+i) = M[i];
      42                 :            :       }
      43                 :            :     }
      44                 :            :   }
      45                 :        933 : }
      46                 :            : 
      47                 :            : tk::Fields
      48                 :          0 : tk::lump( ncomp_t ncomp,
      49                 :            :           const std::array< std::vector< tk::real >, 3 >& coord,
      50                 :            :           const std::vector< std::size_t >& inpoel )
      51                 :            : // *****************************************************************************
      52                 :            : //  Compute lumped mass matrix for CG
      53                 :            : //! \param[in] ncomp Total number of scalar components in the eq system
      54                 :            : //! \param[in] coord Mesh node coordinates
      55                 :            : //! \param[in] inpoel Mesh element connectivity
      56                 :            : //! \return Lumped mass matrix
      57                 :            : // *****************************************************************************
      58                 :            : {
      59                 :            :   const auto& x = coord[0];
      60                 :            :   const auto& y = coord[1];
      61                 :            :   const auto& z = coord[2];
      62                 :            : 
      63                 :          0 :   tk::Fields L( coord[0].size(), ncomp );
      64                 :            :   L.fill( 0.0 );
      65                 :            : 
      66         [ -  - ]:          0 :   for (std::size_t e=0; e<inpoel.size()/4; ++e) {
      67         [ -  - ]:          0 :     const std::array< std::size_t, 4 > N{{ inpoel[e*4+0], inpoel[e*4+1],
      68                 :          0 :                                            inpoel[e*4+2], inpoel[e*4+3] }};
      69                 :            :     // compute element Jacobi determinant
      70                 :            :     const std::array< tk::real, 3 >
      71                 :          0 :       ba{{ x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]] }},
      72                 :          0 :       ca{{ x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]] }},
      73         [ -  - ]:          0 :       da{{ x[N[3]]-x[N[0]], y[N[3]]-y[N[0]], z[N[3]]-z[N[0]] }};
      74                 :          0 :     const auto J = tk::triple( ba, ca, da ) * 5.0 / 120.0;
      75                 :            :     Assert( J > 0, "Element Jacobian non-positive" );
      76                 :            : 
      77                 :            :     // access pointer to lumped mass left hand side at element nodes
      78 [ -  - ][ -  - ]:          0 :     std::vector< const tk::real* > l( ncomp );
      79         [ -  - ]:          0 :     for (ncomp_t c=0; c<ncomp; ++c) l[c] = L.cptr( c );
      80                 :            : 
      81                 :            :     // scatter-add lumped mass element contributions to lhs nodes
      82         [ -  - ]:          0 :     for (ncomp_t c=0; c<ncomp; ++c)
      83         [ -  - ]:          0 :       for (std::size_t j=0; j<4; ++j)
      84                 :          0 :         L.var(l[c],N[j]) += J;
      85                 :            :   }
      86                 :            : 
      87                 :          0 :   return L;
      88                 :            : }
      89                 :            : 
      90                 :            : std::vector< tk::real >
      91                 :     672720 : tk::massMatrixDubiner( std::size_t dof,
      92                 :            :   tk::real vol )
      93                 :            : // *****************************************************************************
      94                 :            : //  Compute the diagonal mass matrix for DG with Dubiner basis functions
      95                 :            : //! \param[in] dof Number of degrees of freedom
      96                 :            : //! \param[in] vol Volume of element whose mass matrix is desired
      97                 :            : //! \return Diagonal of the mass matrix
      98                 :            : // *****************************************************************************
      99                 :            : {
     100                 :     672720 :   std::vector<tk::real> M(dof, 0.0);
     101                 :            : 
     102                 :            :   // Mass matrix for DG(P0)
     103                 :     672720 :   M[0] = vol;
     104                 :            : 
     105                 :            :   // Augment mass matrix for DG(P1)
     106         [ +  + ]:     672720 :   if(dof > 1) {
     107                 :     179016 :     M[1] = vol / 10.0;
     108                 :     179016 :     M[2] = vol * 3.0/10.0;
     109                 :     179016 :     M[3] = vol * 3.0/5.0;
     110                 :            :   }
     111                 :            : 
     112                 :            :   // Augment mass matrix for DG(P2)
     113         [ +  + ]:     672720 :   if(dof > 4) {
     114                 :      74944 :     M[4] = vol / 35.0;
     115                 :      74944 :     M[5] = vol / 21.0;
     116                 :      74944 :     M[6] = vol / 14.0;
     117                 :      74944 :     M[7] = vol / 7.0;
     118                 :      74944 :     M[8] = vol * 3.0/14.0;
     119                 :      74944 :     M[9] = vol * 3.0/7.0;
     120                 :            :   }
     121                 :            : 
     122                 :     672720 :   return M;
     123                 :            : }