Line data    Source code

       1            1 : // Distributed under the MIT License.
       2              : // See LICENSE.txt for details.
       3              : 
       4              : /// \file
       5              : /// Defines functions mean_value and mean_value_on_boundary.
       6              : 
       7              : #pragma once
       8              : 
       9              : #include <cstddef>
      10              : #include <utility>
      11              : 
      12              : #include "Domain/Structure/Side.hpp"
      13              : #include "NumericalAlgorithms/LinearOperators/DefiniteIntegral.hpp"
      14              : #include "Utilities/ConstantExpressions.hpp"
      15              : #include "Utilities/Gsl.hpp"
      16              : 
      17              : /// \cond
      18              : class DataVector;
      19              : template <size_t>
      20              : class Mesh;
      21              : /// \endcond
      22              : 
      23              : /*!
      24              :  * \ingroup NumericalAlgorithmsGroup
      25              :  * \brief Compute the mean value of a function over a manifold.
      26              :  *
      27              :  * Given a function \f$f\f$, compute its mean value \f$\bar{f}\f$ with respect
      28              :  * to the logical coordinates \f$\boldsymbol{\xi} = (\xi, \eta, \zeta)\f$. E.g.,
      29              :  * in 1 dimension, \f$\bar{f} = \int_{-1}^1 f d\xi \Big/ \int_{-1}^1 d\xi\f$.
      30              :  *
      31              :  * \note
      32              :  * The mean w.r.t. a different set of coordinates
      33              :  * \f$\boldsymbol{x} = \boldsymbol{x}(\boldsymbol{\xi})\f$ can't be directly
      34              :  * computed using this function. Before calling `mean_value`, \f$f\f$ must be
      35              :  * pre-multiplied by the Jacobian determinant
      36              :  * \f$J = \det d\boldsymbol{x}/d\boldsymbol{\xi}\f$ of the mapping
      37              :  * \f$\boldsymbol{x}(\boldsymbol{\xi})\f$. Additionally, the output of
      38              :  * `mean_value` must be multiplied by a factor
      39              :  * \f$2^{\text{d}} / \int J d^{\text{d}}\xi\f$ (in \f$d\f$ dimensions), to
      40              :  * account for the different volume of the manifold in the \f$\boldsymbol{x}\f$
      41              :  * coordinates.
      42              :  *
      43              :  * \param f the function to average.
      44              :  * \param mesh the Mesh defining the grid points on the manifold.
      45              :  */
      46              : template <size_t Dim>
      47            1 : double mean_value(const DataVector& f, const Mesh<Dim>& mesh) {
      48              :   return definite_integral(f, mesh) / two_to_the(Dim);
      49              : }
      50              : 
      51              : /// @{
      52              : /*!
      53              :  * \ingroup NumericalAlgorithmsGroup
      54              :  * \brief Compute the mean value of a function over a boundary of a manifold.
      55              :  *
      56              :  * Given a function \f$f\f$, compute its mean value \f$\bar{f}\f$, over a
      57              :  * boundary, with respect to the logical coordinates
      58              :  * \f$\boldsymbol{\xi} = (\xi, \eta, \zeta)\f$.
      59              :  *
      60              :  * \see `mean_value` for notes about means w.r.t. other coordinates.
      61              :  *
      62              :  * - `f` the function to average.
      63              :  * - `mesh` the Mesh defining the grid points on the manifold.
      64              :  * - `d` the dimension which is sliced away to get the boundary.
      65              :  * - `side` whether it is the lower or upper boundary in the d-th dimension.
      66              :  * - `boundary_buffer` is a pointer to a DataVector of size
      67              :  *   `mesh.slice_away(d).number_of_grid_points()` used as a temporary buffer
      68              :  *   when slicing the data to the boundary.
      69              :  * - `volume_and_slice_indices` a pair of `(volume_index_for_point,
      70              :  *   slice_index_for_point)` computed using the `SliceIterator`. Because
      71              :  *   `SliceIterator` is somewhat expensive, if computing the mean value on the
      72              :  *   same boundary for many different tensor components, prefer computing the
      73              :  *   slice indices once.
      74              :  */
      75              : template <size_t Dim>
      76            1 : double mean_value_on_boundary(const DataVector& f, const Mesh<Dim>& mesh,
      77              :                               size_t d, Side side);
      78              : 
      79              : template <size_t Dim>
      80            1 : double mean_value_on_boundary(gsl::not_null<DataVector*> boundary_buffer,
      81              :                               const DataVector& f, const Mesh<Dim>& mesh,
      82              :                               size_t d, Side side);
      83              : 
      84            1 : double mean_value_on_boundary(gsl::not_null<DataVector*> /*boundary_buffer*/,
      85              :                               const DataVector& f, const Mesh<1>& mesh,
      86              :                               size_t d, Side side);
      87              : 
      88              : template <size_t Dim>
      89            1 : double mean_value_on_boundary(
      90              :     gsl::not_null<DataVector*> boundary_buffer,
      91              :     gsl::span<std::pair<size_t, size_t>> volume_and_slice_indices,
      92              :     const DataVector& f, const Mesh<Dim>& mesh, size_t d, Side /*side*/);
      93              : 
      94            1 : double mean_value_on_boundary(
      95              :     gsl::not_null<DataVector*> /*boundary_buffer*/,
      96              :     gsl::span<std::pair<size_t, size_t>> /*volume_and_slice_indices*/,
      97              :     const DataVector& f, const Mesh<1>& mesh, size_t d, Side side);
      98              : /// @}