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