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 : 1158 : tk::mass( ncomp_t ncomp,
20 : : ncomp_t offset,
21 : : ncomp_t ndof,
22 : : const Fields& geoElem,
23 : : Fields& l )
24 : : // *****************************************************************************
25 : : // Compute the block-diagonal mass matrix for DG
26 : : //! \param[in] ncomp Number of scalar components in this PDE system
27 : : //! \param[in] offset Offset this PDE system operates from
28 : : //! \param[in] ndof Number of solution degrees of freedom
29 : : //! \param[in] geoElem Element geometry array
30 : : //! \param[in,out] l Block diagonal mass matrix
31 : : // *****************************************************************************
32 : : {
33 [ - + ][ - - ]: 1158 : Assert( geoElem.nunk() == l.nunk(), "Size mismatch" );
[ - - ][ - - ]
34 : 1158 : const auto nelem = geoElem.nunk();
35 [ - + ][ - - ]: 1158 : Assert( l.nprop() == ndof*ncomp, "lhs mass matrix size incorrect" );
[ - - ][ - - ]
36 : :
37 : : // Compute LHS for DG(P0)
38 [ + + ]: 1622610 : for (std::size_t e=0; e<nelem; ++e)
39 [ + + ]: 4100724 : for (ncomp_t c=0; c<ncomp; ++c)
40 : 2479272 : l(e, c*ndof, offset) = geoElem(e,0,0);
41 : :
42 : : // Augment LHS for DG(P1)
43 [ + + ]: 1158 : if (ndof > 1) {
44 [ + + ]: 140344 : for (std::size_t e=0; e<nelem; ++e) {
45 [ + + ]: 710032 : for (ncomp_t c=0; c<ncomp; ++c) {
46 : 570032 : const auto mark = c * ndof;
47 : 570032 : l(e, mark+1, offset) = geoElem(e,0,0) / 10.0;
48 : 570032 : l(e, mark+2, offset) = geoElem(e,0,0) * 3.0/10.0;
49 : 570032 : l(e, mark+3, offset) = geoElem(e,0,0) * 3.0/5.0;
50 : : }
51 : : }
52 : : }
53 : :
54 : : // Augment LHS for DG(P2)
55 [ + + ]: 1158 : if (ndof > 4){
56 [ + + ]: 83083 : for (std::size_t e=0; e<nelem; ++e) {
57 [ + + ]: 465624 : for (ncomp_t c=0; c<ncomp; ++c) {
58 : 382700 : const auto mark = c * ndof;
59 : 382700 : l(e, mark+4, offset) = geoElem(e,0,0) / 35.0;
60 : 382700 : l(e, mark+5, offset) = geoElem(e,0,0) / 21.0;
61 : 382700 : l(e, mark+6, offset) = geoElem(e,0,0) / 14.0;
62 : 382700 : l(e, mark+7, offset) = geoElem(e,0,0) / 7.0;
63 : 382700 : l(e, mark+8, offset) = geoElem(e,0,0) * 3.0/14.0;
64 : 382700 : l(e, mark+9, offset) = geoElem(e,0,0) * 3.0/7.0;
65 : : }
66 : : }
67 : : }
68 : 1158 : }
69 : :
70 : : tk::Fields
71 : 730 : tk::lump( ncomp_t ncomp,
72 : : const std::array< std::vector< tk::real >, 3 >& coord,
73 : : const std::vector< std::size_t >& inpoel )
74 : : // *****************************************************************************
75 : : // Compute lumped mass matrix for CG
76 : : //! \param[in] ncomp Total number of scalar components in the eq system
77 : : //! \param[in] coord Mesh node coordinates
78 : : //! \param[in] inpoel Mesh element connectivity
79 : : //! \return Lumped mass matrix
80 : : // *****************************************************************************
81 : : {
82 : 730 : const auto& x = coord[0];
83 : 730 : const auto& y = coord[1];
84 : 730 : const auto& z = coord[2];
85 : :
86 : 730 : tk::Fields L( coord[0].size(), ncomp );
87 [ + - ]: 730 : L.fill( 0.0 );
88 : :
89 [ + + ]: 1159348 : for (std::size_t e=0; e<inpoel.size()/4; ++e) {
90 : 1158618 : const std::array< std::size_t, 4 > N{{ inpoel[e*4+0], inpoel[e*4+1],
91 : 1158618 : inpoel[e*4+2], inpoel[e*4+3] }};
92 : : // compute element Jacobi determinant
93 : : const std::array< tk::real, 3 >
94 : 1158618 : ba{{ x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]] }},
95 : 1158618 : ca{{ x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]] }},
96 : 1158618 : da{{ x[N[3]]-x[N[0]], y[N[3]]-y[N[0]], z[N[3]]-z[N[0]] }};
97 : 1158618 : const auto J = tk::triple( ba, ca, da ) * 5.0 / 120.0;
98 [ - + ][ - - ]: 1158618 : Assert( J > 0, "Element Jacobian non-positive" );
[ - - ][ - - ]
99 : :
100 : : // access pointer to lumped mass left hand side at element nodes
101 [ + - ]: 2317236 : std::vector< const tk::real* > l( ncomp );
102 [ + + ][ + - ]: 3265038 : for (ncomp_t c=0; c<ncomp; ++c) l[c] = L.cptr( c, 0 );
103 : :
104 : : // scatter-add lumped mass element contributions to lhs nodes
105 [ + + ]: 3265038 : for (ncomp_t c=0; c<ncomp; ++c)
106 [ + + ]: 10532100 : for (std::size_t j=0; j<4; ++j)
107 [ + - ]: 8425680 : L.var(l[c],N[j]) += J;
108 : : }
109 : :
110 : 730 : return L;
111 : : }
|
