Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <array>
       7             : #include <cstddef>
       8             : #include <functional>
       9             : #include <ostream>
      10             : 
      11             : #include "NumericalAlgorithms/Spectral/SegmentSize.hpp"
      12             : #include "Utilities/ConstantExpressions.hpp"
      13             : 
      14             : /// \cond
      15             : class Matrix;
      16             : template <size_t Dim>
      17             : class Mesh;
      18             : /// \endcond
      19             : 
      20             : namespace Spectral {
      21             : /// Determine whether data needs to be projected between a child mesh and its
      22             : /// parent mesh. If no projection is necessary the data may be used as-is.
      23             : /// Projection is necessary if the child is either p-refined or h-refined
      24             : /// relative to its parent, or both. This operation is symmetric, i.e. it is
      25             : /// irrelevant in which order the child and the parent mesh are passed in.
      26             : template <size_t Dim>
      27           1 : bool needs_projection(const Mesh<Dim>& mesh1, const Mesh<Dim>& mesh2,
      28             :                       const std::array<SegmentSize, Dim>& child_sizes);
      29             : 
      30             : /*!
      31             :  * \brief The projection matrix from a child mesh to its parent.
      32             :  *
      33             :  * The projection matrices returned by this function (and by
      34             :  * projection_matrix_parent_to_child()) define orthogonal projection operators
      35             :  * between the spaces of functions on a parent mesh and its children. These
      36             :  * projections are usually the correct way to transfer data between meshes in
      37             :  * a mesh-refinement hierarchy, as well as between an element face and its
      38             :  * adjacent mortars.
      39             :  *
      40             :  * These functions assume that the `child_mesh` is at least as fine as the
      41             :  * `parent_mesh`, i.e. functions on the `parent_mesh` can be represented exactly
      42             :  * on the `child_mesh`. In practice this means that functions can be projected
      43             :  * to a mortar (the `child_mesh`) from both adjacent element faces (the
      44             :  * `parent_mesh`) without losing accuracy. Similarly, functions in a
      45             :  * mesh-refinement hierarchy don't lose accuracy when an element is split
      46             :  * (h-refined). For this reason, the `projection_matrix_child_to_parent` is
      47             :  * sometimes referred to as a "restriction operator" and the
      48             :  * `projection_matrix_parent_to_child` as a "prolongation operator".
      49             :  *
      50             :  * \par Massive quantities
      51             :  * If the quantity that should be projected is not a function over the
      52             :  * computational grid but a "massive" residual, i.e. a quantity
      53             :  * \f$\int_{\Omega_k} f(x) \psi_p(x) \mathrm{d}V\f$ where \f$\psi_p\f$ are the
      54             :  * basis functions on the mesh, then pass `true` for the parameter
      55             :  * `operand_is_massive` (default is `false`). The restriction operator for this
      56             :  * case is just the transpose of the prolongation operator, i.e. just an
      57             :  * interpolation matrix transpose. Note that the "massive" residual already
      58             :  * takes the difference in element size between parent and children into account
      59             :  * by including a Jacobian in the volume element of the integral.
      60             :  *
      61             :  * \par Implementation details
      62             :  * The half-interval projections are based on an equation derived by
      63             :  * Saul.  This shows that the projection from the spectral basis for
      64             :  * the entire interval to the spectral basis for the upper half
      65             :  * interval is
      66             :  * \f{equation*}
      67             :  * T_{jk} = \frac{2 j + 1}{2} 2^j \sum_{n=0}^{j-k} \binom{j}{k+n}
      68             :  * \binom{(j + k + n - 1)/2}{j} \frac{(k + n)!^2}{(2 k + n + 1)! n!}
      69             :  * \f}
      70             :  */
      71           1 : const Matrix& projection_matrix_child_to_parent(
      72             :     const Mesh<1>& child_mesh, const Mesh<1>& parent_mesh, SegmentSize size,
      73             :     bool operand_is_massive = false);
      74             : 
      75             : /// The projection matrix from a child mesh to its parent, in `Dim` dimensions.
      76             : template <size_t Dim>
      77             : std::array<std::reference_wrapper<const Matrix>, Dim>
      78           1 : projection_matrix_child_to_parent(
      79             :     const Mesh<Dim>& child_mesh, const Mesh<Dim>& parent_mesh,
      80             :     const std::array<SegmentSize, Dim>& child_sizes,
      81             :     bool operand_is_massive = false);
      82             : 
      83             : /// The projection matrix from a parent mesh to one of its children.
      84             : ///
      85             : /// \see projection_matrix_child_to_parent()
      86           1 : const Matrix& projection_matrix_parent_to_child(const Mesh<1>& parent_mesh,
      87             :                                                 const Mesh<1>& child_mesh,
      88             :                                                 SegmentSize size);
      89             : 
      90             : /// The projection matrix from a parent mesh to one of its children, in `Dim`
      91             : /// dimensions
      92             : template <size_t Dim>
      93             : std::array<std::reference_wrapper<const Matrix>, Dim>
      94           1 : projection_matrix_parent_to_child(
      95             :     const Mesh<Dim>& parent_mesh, const Mesh<Dim>& child_mesh,
      96             :     const std::array<SegmentSize, Dim>& child_sizes);
      97             : 
      98             : /// The projection matrices from a source mesh to a target mesh
      99             : /// covering given portions of an element
     100             : template <size_t Dim>
     101           1 : std::array<std::reference_wrapper<const Matrix>, Dim> projection_matrices(
     102             :     const Mesh<Dim>& source_mesh, const Mesh<Dim>& target_mesh,
     103             :     const std::array<SegmentSize, Dim>& source_sizes,
     104             :     const std::array<SegmentSize, Dim>& target_sizes);
     105             : 
     106             : /// The projection matrices from a source mesh to a target mesh where the
     107             : /// meshes cover the same physical volume
     108             : template <size_t Dim>
     109           1 : std::array<std::reference_wrapper<const Matrix>, Dim> p_projection_matrices(
     110             :     const Mesh<Dim>& source_mesh, const Mesh<Dim>& target_mesh);
     111             : 
     112             : /// @{
     113             : /// \brief Performs a perfect hash of the mortars into $2^{d-1}$ slots on the
     114             : /// range $[0, 2^{d-1})$.
     115             : ///
     116             : /// This is particularly useful when hashing into statically-sized maps based
     117             : /// on the number of dimensions.
     118             : template <size_t DimMinusOne>
     119           1 : size_t hash(const std::array<Spectral::SegmentSize, DimMinusOne>& mortar_size);
     120             : 
     121             : template <size_t Dim>
     122           0 : struct MortarSizeHash {
     123             :   template <size_t MaxSize>
     124           0 :   static constexpr bool is_perfect = MaxSize == two_to_the(Dim);
     125             : 
     126           0 :   size_t operator()(
     127             :       const std::array<Spectral::SegmentSize, Dim - 1>& mortar_size);
     128             : };
     129             : /// @}
     130             : }  // namespace Spectral