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                 :            : tk::Fields
      19                 :          0 : tk::lump( ncomp_t ncomp,
      20                 :            :           const std::array< std::vector< tk::real >, 3 >& coord,
      21                 :            :           const std::vector< std::size_t >& inpoel )
      22                 :            : // *****************************************************************************
      23                 :            : //  Compute lumped mass matrix for CG
      24                 :            : //! \param[in] ncomp Total number of scalar components in the eq system
      25                 :            : //! \param[in] coord Mesh node coordinates
      26                 :            : //! \param[in] inpoel Mesh element connectivity
      27                 :            : //! \return Lumped mass matrix
      28                 :            : // *****************************************************************************
      29                 :            : {
      30                 :          0 :   const auto& x = coord[0];
      31                 :          0 :   const auto& y = coord[1];
      32                 :          0 :   const auto& z = coord[2];
      33                 :            : 
      34                 :          0 :   tk::Fields L( coord[0].size(), ncomp );
      35         [ -  - ]:          0 :   L.fill( 0.0 );
      36                 :            : 
      37         [ -  - ]:          0 :   for (std::size_t e=0; e<inpoel.size()/4; ++e) {
      38                 :          0 :     const std::array< std::size_t, 4 > N{{ inpoel[e*4+0], inpoel[e*4+1],
      39                 :          0 :                                            inpoel[e*4+2], inpoel[e*4+3] }};
      40                 :            :     // compute element Jacobi determinant
      41                 :            :     const std::array< tk::real, 3 >
      42                 :          0 :       ba{{ x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]] }},
      43                 :          0 :       ca{{ x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]] }},
      44                 :          0 :       da{{ x[N[3]]-x[N[0]], y[N[3]]-y[N[0]], z[N[3]]-z[N[0]] }};
      45                 :          0 :     const auto J = tk::triple( ba, ca, da ) * 5.0 / 120.0;
      46 [ -  - ][ -  - ]:          0 :     Assert( J > 0, "Element Jacobian non-positive" );
         [ -  - ][ -  - ]
      47                 :            : 
      48                 :            :     // access pointer to lumped mass left hand side at element nodes
      49         [ -  - ]:          0 :     std::vector< const tk::real* > l( ncomp );
      50 [ -  - ][ -  - ]:          0 :     for (ncomp_t c=0; c<ncomp; ++c) l[c] = L.cptr( c );
      51                 :            : 
      52                 :            :     // scatter-add lumped mass element contributions to lhs nodes
      53         [ -  - ]:          0 :     for (ncomp_t c=0; c<ncomp; ++c)
      54         [ -  - ]:          0 :       for (std::size_t j=0; j<4; ++j)
      55         [ -  - ]:          0 :         L.var(l[c],N[j]) += J;
      56                 :            :   }
      57                 :            : 
      58                 :          0 :   return L;
      59                 :            : }