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