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             : 
      10             : #include "DataStructures/DataVector.hpp"
      11             : #include "DataStructures/Matrix.hpp"
      12             : #include "DataStructures/Variables.hpp"
      13             : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
      14             : #include "NumericalAlgorithms/Spectral/SegmentSize.hpp"
      15             : #include "Utilities/ConstantExpressions.hpp"
      16             : #include "Utilities/Gsl.hpp"
      17             : 
      18             : namespace Spectral {
      19             : /// Determine whether data needs to be projected between a child mesh and its
      20             : /// parent mesh. If no projection is necessary the data may be used as-is.
      21             : /// Projection is necessary if the child is either p-refined or h-refined
      22             : /// relative to its parent, or both. This operation is symmetric, i.e. it is
      23             : /// irrelevant in which order the child and the parent mesh are passed in.
      24             : template <size_t Dim>
      25           1 : bool needs_projection(const Mesh<Dim>& mesh1, const Mesh<Dim>& mesh2,
      26             :                       const std::array<SegmentSize, Dim>& child_sizes);
      27             : 
      28             : /*!
      29             :  * \brief The projection matrix from a child mesh to its parent.
      30             :  *
      31             :  * The projection matrices returned by this function (and by
      32             :  * projection_matrix_parent_to_child()) define orthogonal projection operators
      33             :  * between the spaces of functions on a parent mesh and its children. These
      34             :  * projections are usually the correct way to transfer data between meshes in
      35             :  * a mesh-refinement hierarchy, as well as between an element face and its
      36             :  * adjacent mortars.
      37             :  *
      38             :  * These functions assume that the `child_mesh` is at least as fine as the
      39             :  * `parent_mesh`, i.e. functions on the `parent_mesh` can be represented exactly
      40             :  * on the `child_mesh`. In practice this means that functions can be projected
      41             :  * to a mortar (the `child_mesh`) from both adjacent element faces (the
      42             :  * `parent_mesh`) without losing accuracy. Similarly, functions in a
      43             :  * mesh-refinement hierarchy don't lose accuracy when an element is split
      44             :  * (h-refined). For this reason, the `projection_matrix_child_to_parent` is
      45             :  * sometimes referred to as a "restriction operator" and the
      46             :  * `projection_matrix_parent_to_child` as a "prolongation operator".
      47             :  *
      48             :  * \par Massive quantities
      49             :  * If the quantity that should be projected is not a function over the
      50             :  * computational grid but a "massive" residual, i.e. a quantity
      51             :  * \f$\int_{\Omega_k} f(x) \psi_p(x) \mathrm{d}V\f$ where \f$\psi_p\f$ are the
      52             :  * basis functions on the mesh, then pass `true` for the parameter
      53             :  * `operand_is_massive` (default is `false`). The restriction operator for this
      54             :  * case is just the transpose of the prolongation operator, i.e. just an
      55             :  * interpolation matrix transpose. Note that the "massive" residual already
      56             :  * takes the difference in element size between parent and children into account
      57             :  * by including a Jacobian in the volume element of the integral.
      58             :  *
      59             :  * \par Implementation details
      60             :  * The half-interval projections are based on an equation derived by
      61             :  * Saul.  This shows that the projection from the spectral basis for
      62             :  * the entire interval to the spectral basis for the upper half
      63             :  * interval is
      64             :  * \f{equation*}
      65             :  * T_{jk} = \frac{2 j + 1}{2} 2^j \sum_{n=0}^{j-k} \binom{j}{k+n}
      66             :  * \binom{(j + k + n - 1)/2}{j} \frac{(k + n)!^2}{(2 k + n + 1)! n!}
      67             :  * \f}
      68             :  *
      69             :  * \note This and the other matrix-returning projection functions in this file
      70             :  * (projection_matrix_parent_to_child(), projection_matrices()) only support
      71             :  * tensor-product bases (Legendre, Fourier, ZernikeB2). They cannot be used for
      72             :  * spherical-shell meshes that use the `SphericalHarmonic` basis, because the
      73             :  * angular projection couples the \f$\theta\f$ and \f$\phi\f$ directions and so
      74             :  * cannot be expressed as the per-dimension matrices these functions return. Use
      75             :  * Spectral::project() instead whenever spherical-shell meshes must be
      76             :  * supported; it dispatches to these matrices for tensor-product meshes and to a
      77             :  * Spherepack-based angular projection for spherical-shell meshes.
      78             :  */
      79           1 : const Matrix& projection_matrix_child_to_parent(
      80             :     const Mesh<1>& child_mesh, const Mesh<1>& parent_mesh, SegmentSize size,
      81             :     bool operand_is_massive = false);
      82             : 
      83             : /// The projection matrix from a child mesh to its parent, in `Dim` dimensions.
      84             : template <size_t Dim>
      85             : std::array<std::reference_wrapper<const Matrix>, Dim>
      86           1 : projection_matrix_child_to_parent(
      87             :     const Mesh<Dim>& child_mesh, const Mesh<Dim>& parent_mesh,
      88             :     const std::array<SegmentSize, Dim>& child_sizes,
      89             :     bool operand_is_massive = false);
      90             : 
      91             : /// The projection matrix from a parent mesh to one of its children.
      92             : ///
      93             : /// \see projection_matrix_child_to_parent()
      94           1 : const Matrix& projection_matrix_parent_to_child(const Mesh<1>& parent_mesh,
      95             :                                                 const Mesh<1>& child_mesh,
      96             :                                                 SegmentSize size);
      97             : 
      98             : /// The projection matrix from a parent mesh to one of its children, in `Dim`
      99             : /// dimensions
     100             : template <size_t Dim>
     101             : std::array<std::reference_wrapper<const Matrix>, Dim>
     102           1 : projection_matrix_parent_to_child(
     103             :     const Mesh<Dim>& parent_mesh, const Mesh<Dim>& child_mesh,
     104             :     const std::array<SegmentSize, Dim>& child_sizes);
     105             : 
     106             : /// The projection matrices from a source mesh to a target mesh
     107             : /// covering given portions of an element
     108             : template <size_t Dim>
     109           1 : std::array<std::reference_wrapper<const Matrix>, Dim> projection_matrices(
     110             :     const Mesh<Dim>& source_mesh, const Mesh<Dim>& target_mesh,
     111             :     const std::array<SegmentSize, Dim>& source_sizes,
     112             :     const std::array<SegmentSize, Dim>& target_sizes,
     113             :     bool operand_is_massive = false);
     114             : 
     115             : /// Change the angular resolution (`l_max`, with `m_max == l_max`) of volume
     116             : /// data on a spherical shell from `source_data` to `result_data`, for
     117             : /// `num_components` components each laid out with the radial dimension varying
     118             : /// fastest. Uses Spherepack prolong/restrict. \see Spectral::project() for the
     119             : /// higher-level interface to call this.
     120             : ///
     121             : /// For non-massive operands this is the L2 (Galerkin) projection,
     122             : /// $P_\mathrm{L2}$. For massive operands (`operand_is_massive == true`) the
     123             : /// restriction is the transpose of the prolongation (interpolation) operator,
     124             : /// $I^T$, which is the L2 projection conjugated by the diagonal matrices $W$ of
     125             : /// angular quadrature weights:
     126             : /// $I^T = W_\mathrm{target} P_\mathrm{L2} W_\mathrm{source}^{-1}$.
     127             : /// This means we divide by the source weights before projecting and multiply by
     128             : /// the target weights afterwards. See projection_matrix_child_to_parent() for
     129             : /// details on massive operands.
     130           1 : void project_spherical_harmonics(gsl::not_null<double*> result_data,
     131             :                                  const double* source_data,
     132             :                                  size_t num_components,
     133             :                                  size_t num_radial_points, size_t l_max_source,
     134             :                                  size_t l_max_target, bool operand_is_massive);
     135             : 
     136             : /// @{
     137             : /*!
     138             :  * \brief Project volume data from `source_mesh` to `target_mesh`, writing the
     139             :  * (resized) result into `*result`.
     140             :  *
     141             :  * This is the unified entry point for projection regardless of basis. It
     142             :  * handles both tensor-product meshes and spherical-shell meshes (Legendre
     143             :  * radial dimension plus `SphericalHarmonic` angular dimensions with
     144             :  * `m_max == l_max`):
     145             :  * - For tensor-product meshes it delegates to the per-dimension projection
     146             :  *   matrices above.
     147             :  * - For spherical-shell meshes it projects the radial dimension with the 1D
     148             :  *   projection matrix and changes the angular `l_max` with a Spherepack
     149             :  *   prolong/restrict. The angular dimensions cannot be h-refined, so their
     150             :  *   segment sizes must be `Full`.
     151             :  *
     152             :  * `source_sizes` and `target_sizes` carry h-refinement information (which
     153             :  * portion of an element the source/target cover); pass `Full` in every
     154             :  * dimension for pure p-refinement.
     155             :  *
     156             :  * \warning `*result` must not point to `source` (the result is resized to the
     157             :  * target number of grid points).
     158             :  */
     159             : template <typename VectorType, size_t Dim>
     160           1 : void project(gsl::not_null<VectorType*> result, const VectorType& source,
     161             :              const Mesh<Dim>& source_mesh, const Mesh<Dim>& target_mesh,
     162             :              const std::array<SegmentSize, Dim>& source_sizes,
     163             :              const std::array<SegmentSize, Dim>& target_sizes,
     164             :              bool operand_is_massive = false);
     165             : 
     166             : template <size_t Dim, typename TagList>
     167           1 : void project(const gsl::not_null<Variables<TagList>*> result,
     168             :              const Variables<TagList>& source, const Mesh<Dim>& source_mesh,
     169             :              const Mesh<Dim>& target_mesh,
     170             :              const std::array<SegmentSize, Dim>& source_sizes,
     171             :              const std::array<SegmentSize, Dim>& target_sizes,
     172             :              const bool operand_is_massive = false) {
     173             :   // Type-erase to the vector implementation with multiple components
     174             :   using VectorType = typename Variables<TagList>::vector_type;
     175             :   using ValueType = typename Variables<TagList>::value_type;
     176             :   result->initialize(target_mesh.number_of_grid_points());
     177             :   VectorType result_view(result->data(), result->size());
     178             :   // NOLINTNEXTLINE(cppcoreguidelines-pro-type-const-cast)
     179             :   const VectorType source_view(const_cast<ValueType*>(source.data()),
     180             :                                source.size());
     181             :   project(make_not_null(&result_view), source_view, source_mesh, target_mesh,
     182             :           source_sizes, target_sizes, operand_is_massive);
     183             : }
     184             : 
     185             : template <size_t Dim, typename T>
     186           1 : T project(const T& source, const Mesh<Dim>& source_mesh,
     187             :           const Mesh<Dim>& target_mesh,
     188             :           const std::array<SegmentSize, Dim>& source_sizes,
     189             :           const std::array<SegmentSize, Dim>& target_sizes,
     190             :           const bool operand_is_massive = false) {
     191             :   T result{};
     192             :   project(make_not_null(&result), source, source_mesh, target_mesh,
     193             :           source_sizes, target_sizes, operand_is_massive);
     194             :   return result;
     195             : }
     196             : /// @}
     197             : 
     198             : /// @{
     199             : /// \brief Performs a perfect hash of the mortars into $2^{d-1}$ slots on the
     200             : /// range $[0, 2^{d-1})$.
     201             : ///
     202             : /// This is particularly useful when hashing into statically-sized maps based
     203             : /// on the number of dimensions.
     204             : template <size_t DimMinusOne>
     205           1 : size_t hash(const std::array<Spectral::SegmentSize, DimMinusOne>& mortar_size);
     206             : 
     207             : template <size_t Dim>
     208           0 : struct MortarSizeHash {
     209             :   template <size_t MaxSize>
     210           0 :   static constexpr bool is_perfect = MaxSize == two_to_the(Dim);
     211             : 
     212           0 :   size_t operator()(
     213             :       const std::array<Spectral::SegmentSize, Dim - 1>& mortar_size);
     214             : };
     215             : /// @}
     216             : }  // namespace Spectral