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