Line data    Source code

       1            0 : // Distributed under the MIT License.
       2              : // See LICENSE.txt for details.
       3              : 
       4              : #pragma once
       5              : 
       6              : #include <algorithm>
       7              : #include <array>
       8              : #include <cmath>
       9              : #include <cstddef>
      10              : 
      11              : #include "DataStructures/ApplyMatrices.hpp"
      12              : #include "DataStructures/DataVector.hpp"
      13              : #include "DataStructures/Matrix.hpp"
      14              : #include "DataStructures/Tensor/Tensor.hpp"
      15              : #include "DataStructures/Variables.hpp"
      16              : #include "NumericalAlgorithms/Interpolation/IrregularInterpolant.hpp"
      17              : #include "NumericalAlgorithms/Spectral/Basis.hpp"
      18              : #include "NumericalAlgorithms/Spectral/InterpolationMatrix.hpp"
      19              : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
      20              : #include "Utilities/ErrorHandling/Assert.hpp"
      21              : #include "Utilities/Gsl.hpp"
      22              : 
      23              : namespace intrp {
      24              : namespace InterpolateOnSphericalShell_detail {
      25              : // Target points whose radial logical coordinates differ by at most this much
      26              : // from the first point of a group are interpolated together. This absorbs the
      27              : // roundoff in the logical coordinates of points at the same radius.
      28              : constexpr double radial_coordinate_tolerance = 1.0e-12;
      29              : }  // namespace InterpolateOnSphericalShell_detail
      30              : 
      31              : /*!
      32              :  * \brief Interpolates `vars` on a spherical-harmonic shell `mesh` to the
      33              :  * `target_points`, first in the radial dimension and then in the angular
      34              :  * dimensions.
      35              :  *
      36              :  * Consecutive target points whose radial logical coordinates agree with that of
      37              :  * the first point of the group to within a tolerance share the
      38              :  * radial interpolation, which leaves an interpolation on the sphere with
      39              :  * `intrp::Irregular<2>`. This is much cheaper than `intrp::Irregular<3>` when
      40              :  * the target points lie on a few spheres of constant radial logical coordinate.
      41              :  */
      42              : template <typename TagsList>
      43            1 : Variables<TagsList> interpolate_on_spherical_shell(
      44              :     const Variables<TagsList>& vars, const Mesh<3>& mesh,
      45              :     const tnsr::I<DataVector, 3, Frame::ElementLogical>& target_points) {
      46              :   ASSERT((mesh.basis(0) == Spectral::Basis::Legendre or
      47              :           mesh.basis(0) == Spectral::Basis::Chebyshev) and
      48              :              mesh.basis(1) == Spectral::Basis::SphericalHarmonic and
      49              :              mesh.basis(2) == Spectral::Basis::SphericalHarmonic,
      50              :          "Expected a spherical-harmonic shell, but the mesh is " << mesh);
      51              :   ASSERT(vars.number_of_grid_points() == mesh.number_of_grid_points(),
      52              :          "The mesh has " << mesh.number_of_grid_points()
      53              :                          << " grid points but the variables have "
      54              :                          << vars.number_of_grid_points() << ".");
      55              :   const size_t number_of_target_points = get<0>(target_points).size();
      56              :   Variables<TagsList> result(number_of_target_points);
      57              :   const Mesh<1> radial_mesh = mesh.slice_through(0);
      58              :   const Mesh<2> angular_mesh = mesh.slice_away(0);
      59              :   size_t begin = 0;
      60              :   while (begin < number_of_target_points) {
      61              :     const double radial_coord = get<0>(target_points)[begin];
      62              :     size_t end = begin + 1;
      63              :     while (
      64              :         end < number_of_target_points and
      65              :         std::abs(get<0>(target_points)[end] - radial_coord) <=
      66              :             InterpolateOnSphericalShell_detail::radial_coordinate_tolerance) {
      67              :       ++end;
      68              :     }
      69              :     const size_t group_size = end - begin;
      70              :     // Interpolate radially to the sphere of the group
      71              :     const Variables<TagsList> vars_on_sphere = apply_matrices(
      72              :         std::array<Matrix, 3>{{Spectral::interpolation_matrix(
      73              :                                    radial_mesh, DataVector{1, radial_coord}),
      74              :                                Matrix{}, Matrix{}}},
      75              :         vars, mesh.extents());
      76              :     // Then on the sphere to the angular coordinates of the group
      77              :     const tnsr::I<DataVector, 2, Frame::ElementLogical> angular_points{};
      78              :     for (size_t d = 0; d < 2; ++d) {
      79              :       make_const_view(make_not_null(&angular_points.get(d)),
      80              :                       target_points.get(d + 1), begin, group_size);
      81              :     }
      82              :     const Variables<TagsList> result_on_sphere =
      83              :         intrp::Irregular<2>(angular_mesh, angular_points)
      84              :             .interpolate(vars_on_sphere);
      85              :     for (size_t component = 0;
      86              :          component < result.number_of_independent_components; ++component) {
      87              :       std::copy_n(result_on_sphere.data() + component * group_size, group_size,
      88              :                   result.data() + component * number_of_target_points + begin);
      89              :     }
      90              :     begin = end;
      91              :   }
      92              :   return result;
      93              : }
      94              : }  // namespace intrp