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 : 985 : 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 [ - + ][ - - ]: 985 : Assert( geoElem.nunk() == l.nunk(), "Size mismatch" );
[ - - ][ - - ]
32 : 985 : const auto nelem = geoElem.nunk();
33 [ - + ][ - - ]: 985 : Assert( l.nprop() == ndof*ncomp, "lhs mass matrix size incorrect" );
[ - - ][ - - ]
34 : :
35 : : // Compute LHS
36 [ + + ]: 650701 : for (std::size_t e=0; e<nelem; ++e) {
37 [ + - ][ + - ]: 1299432 : auto M = massMatrixDubiner(ndof, geoElem(e,0));
38 [ + + ]: 2202032 : for (ncomp_t c=0; c<ncomp; ++c) {
39 : 1552316 : const auto mark = c * ndof;
40 [ + + ]: 6792760 : for (std::size_t i=0; i<ndof; ++i) {
41 [ + - ]: 5240444 : l(e, mark+i) = M[i];
42 : : }
43 : : }
44 : : }
45 : 985 : }
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 : 0 : const auto& x = coord[0];
60 : 0 : const auto& y = coord[1];
61 : 0 : const auto& z = coord[2];
62 : :
63 : 0 : tk::Fields L( coord[0].size(), ncomp );
64 [ - - ]: 0 : 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 [ - - ][ - - ]: 0 : 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 : 767556 : 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 [ + - ]: 767556 : std::vector<tk::real> M(dof, 0.0);
101 : :
102 : : // Mass matrix for DG(P0)
103 : 767556 : M[0] = vol;
104 : :
105 : : // Augment mass matrix for DG(P1)
106 [ + + ]: 767556 : if(dof > 1) {
107 : 251376 : M[1] = vol / 10.0;
108 : 251376 : M[2] = vol * 3.0/10.0;
109 : 251376 : M[3] = vol * 3.0/5.0;
110 : : }
111 : :
112 : : // Augment mass matrix for DG(P2)
113 [ + + ]: 767556 : 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 : 767556 : return M;
123 : : }
|
