Line data Source code
1 0 : // Distributed under the MIT License. 2 : // See LICENSE.txt for details. 3 : 4 : #pragma once 5 : 6 : #include <cstddef> 7 : 8 : #include "DataStructures/Variables.hpp" 9 : #include "NumericalAlgorithms/Spectral/Mesh.hpp" 10 : #include "Utilities/ErrorHandling/Assert.hpp" 11 : #include "Utilities/Gsl.hpp" 12 : 13 1 : namespace dg { 14 : 15 : namespace detail { 16 : template <size_t Dim> 17 : void apply_mass_matrix_impl(gsl::not_null<double*> data, const Mesh<Dim>& mesh); 18 : template <size_t Dim> 19 : void apply_inverse_mass_matrix_impl(gsl::not_null<double*> data, 20 : const Mesh<Dim>& mesh); 21 : } // namespace detail 22 : 23 : /*! 24 : * \brief Apply the DG mass matrix to the data, in the diagonal mass-matrix 25 : * approximation ("mass-lumping") 26 : * 27 : * The DG mass matrix is: 28 : * 29 : * \f{equation} 30 : * M_{pq} = \int_{\Omega_k} \psi_p(\xi) \psi_q(\xi) \mathrm{d}V 31 : * \f} 32 : * 33 : * where \f$\psi_p(\xi)\f$ are the basis functions on the element 34 : * \f$\Omega_k\f$. In the diagonal mass-matrix approximation ("mass-lumping") we 35 : * evaluate the integral directly on the collocation points, i.e. with a Gauss 36 : * or Gauss-Lobatto quadrature determined by the element mesh. Then it reduces 37 : * to: 38 : * 39 : * \f{equation} 40 : * M_{pq} \approx \delta_{pq} \prod_{i=1}^d w_{p_i} 41 : * \f} 42 : * 43 : * where \f$d\f$ is the spatial dimension and \f$w_{p_i}\f$ are the quadrature 44 : * weights in dimension \f$i\f$. To apply the mass matrix in coordinates 45 : * different than logical, or to account for a curved background metric, the 46 : * data can be pre-multiplied with the Jacobian determinant and/or a metric 47 : * determinant. 48 : * 49 : * \note The mass-lumping is exact on Legendre-Gauss meshes, but omits a 50 : * correction term on Legendre-Gauss-Lobatto meshes. 51 : */ 52 : /// @{ 53 : template <size_t Dim> 54 1 : void apply_mass_matrix(const gsl::not_null<DataVector*> data, 55 : const Mesh<Dim>& mesh) { 56 : ASSERT(data->size() == mesh.number_of_grid_points(), 57 : "The DataVector has size " << data->size() << ", but expected size " 58 : << mesh.number_of_grid_points() 59 : << " on the given mesh."); 60 : detail::apply_mass_matrix_impl(data->data(), mesh); 61 : } 62 : 63 : template <size_t Dim, typename TagsList> 64 1 : void apply_mass_matrix(const gsl::not_null<Variables<TagsList>*> data, 65 : const Mesh<Dim>& mesh) { 66 : const size_t num_points = data->number_of_grid_points(); 67 : ASSERT(num_points == mesh.number_of_grid_points(), 68 : "The Variables data has " 69 : << num_points << " grid points, but expected " 70 : << mesh.number_of_grid_points() << " on the given mesh."); 71 : constexpr size_t num_comps = 72 : Variables<TagsList>::number_of_independent_components; 73 : for (size_t i = 0; i < num_comps; ++i) { 74 : detail::apply_mass_matrix_impl(data->data() + i * num_points, mesh); 75 : } 76 : } 77 : /// @} 78 : 79 : /*! 80 : * \brief Apply the inverse DG mass matrix to the data 81 : * 82 : * \see dg::apply_mass_matrix 83 : */ 84 : /// @{ 85 : template <size_t Dim> 86 1 : void apply_inverse_mass_matrix(const gsl::not_null<DataVector*> data, 87 : const Mesh<Dim>& mesh) { 88 : ASSERT(data->size() == mesh.number_of_grid_points(), 89 : "The DataVector has size " << data->size() << ", but expected size " 90 : << mesh.number_of_grid_points() 91 : << " on the given mesh."); 92 : detail::apply_inverse_mass_matrix_impl(data->data(), mesh); 93 : } 94 : 95 : template <size_t Dim, typename TagsList> 96 1 : void apply_inverse_mass_matrix(const gsl::not_null<Variables<TagsList>*> data, 97 : const Mesh<Dim>& mesh) { 98 : const size_t num_points = data->number_of_grid_points(); 99 : ASSERT(num_points == mesh.number_of_grid_points(), 100 : "The Variables data has " 101 : << num_points << " grid points, but expected " 102 : << mesh.number_of_grid_points() << " on the given mesh."); 103 : constexpr size_t num_comps = 104 : Variables<TagsList>::number_of_independent_components; 105 : for (size_t i = 0; i < num_comps; ++i) { 106 : detail::apply_inverse_mass_matrix_impl(data->data() + i * num_points, mesh); 107 : } 108 : } 109 : /// @} 110 : 111 : } // namespace dg