ApplyMassMatrix.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <cstddef>
7 
11 #include "Utilities/Gsl.hpp"
12 
13 namespace dg {
14 
15 namespace detail {
16 template <size_t Dim>
17 void apply_mass_matrix_impl(gsl::not_null<double*> data,
18  const Mesh<Dim>& mesh) noexcept;
19 } // namespace detail
20 
21 /*!
22  * \brief Apply the DG mass matrix to the data, in the diagonal mass-matrix
23  * approximation ("mass-lumping")
24  *
25  * The DG mass matrix is:
26  *
27  * \f{equation}
28  * M_{pq} = \int_{\Omega_k} \psi_p(\xi) \psi_q(\xi) \mathrm{d}V
29  * \f}
30  *
31  * where \f$\psi_p(\xi)\f$ are the basis functions on the element
32  * \f$\Omega_k\f$. In the diagonal mass-matrix approximation ("mass-lumping") we
33  * evaluate the integral directly on the collocation points, i.e. with a Gauss
34  * or Gauss-Lobatto quadrature determined by the element mesh. Then it reduces
35  * to:
36  *
37  * \f{equation}
38  * M_{pq} \approx \delta_{pq} \prod_{i=1}^d w_{p_i}
39  * \f}
40  *
41  * where \f$d\f$ is the spatial dimension and \f$w_{p_i}\f$ are the quadrature
42  * weights in dimension \f$i\f$. To apply the mass matrix in coordinates
43  * different than logical, or to account for a curved background metric, the
44  * data can be pre-multiplied with the Jacobian determinant and/or a metric
45  * determinant.
46  *
47  * \note The mass-lumping is exact on Legendre-Gauss meshes, but omits a
48  * correction term on Legendre-Gauss-Lobatto meshes.
49  */
50 /// @{
51 template <size_t Dim>
53  const Mesh<Dim>& mesh) noexcept {
54  ASSERT(data->size() == mesh.number_of_grid_points(),
55  "The DataVector has size " << data->size() << ", but expected size "
56  << mesh.number_of_grid_points()
57  << " on the given mesh.");
58  detail::apply_mass_matrix_impl(data->data(), mesh);
59 }
60 
61 template <size_t Dim, typename TagsList>
62 void apply_mass_matrix(const gsl::not_null<Variables<TagsList>*> data,
63  const Mesh<Dim>& mesh) noexcept {
64  const size_t num_points = data->number_of_grid_points();
65  ASSERT(num_points == mesh.number_of_grid_points(),
66  "The Variables data has "
67  << num_points << " grid points, but expected "
68  << mesh.number_of_grid_points() << " on the given mesh.");
69  constexpr size_t num_comps =
70  Variables<TagsList>::number_of_independent_components;
71  for (size_t i = 0; i < num_comps; ++i) {
72  detail::apply_mass_matrix_impl(data->data() + i * num_points, mesh);
73  }
74 }
75 /// @}
76 
77 } // namespace dg
dg::apply_mass_matrix
void apply_mass_matrix(const gsl::not_null< DataVector * > data, const Mesh< Dim > &mesh) noexcept
Apply the DG mass matrix to the data, in the diagonal mass-matrix approximation ("mass-lumping")
Definition: ApplyMassMatrix.hpp:52
dg
Functionality related to discontinuous Galerkin schemes.
Definition: ApplyMassMatrix.hpp:13
cstddef
Assert.hpp
ASSERT
#define ASSERT(a, m)
Assert that an expression should be true.
Definition: Assert.hpp:49
Variables.hpp
Mesh
Holds the number of grid points, basis, and quadrature in each direction of the computational grid.
Definition: Mesh.hpp:49
Gsl.hpp
Mesh.hpp
gsl::not_null
Require a pointer to not be a nullptr
Definition: ReadSpecPiecewisePolynomial.hpp:13