diff --git a/src/NumericalAlgorithms/Interpolation/IrregularInterpolant.cpp b/src/NumericalAlgorithms/Interpolation/IrregularInterpolant.cpp index f7282b53bf26..e0dd0d9f3cbb 100644 --- a/src/NumericalAlgorithms/Interpolation/IrregularInterpolant.cpp +++ b/src/NumericalAlgorithms/Interpolation/IrregularInterpolant.cpp @@ -10,6 +10,7 @@ #include "DataStructures/DataVector.hpp" #include "DataStructures/Tensor/Tensor.hpp" +#include "NumericalAlgorithms/Interpolation/CardinalInterpolator.hpp" #include "NumericalAlgorithms/Spectral/LogicalCoordinates.hpp" #include "NumericalAlgorithms/Spectral/Mesh.hpp" #include "NumericalAlgorithms/Spectral/Spectral.hpp" @@ -105,9 +106,26 @@ Matrix interpolation_matrix( } } return result; + } else if (mesh.basis()[0] == Spectral::Basis::SphericalHarmonic) { + ASSERT(mesh.basis()[1] == Spectral::Basis::SphericalHarmonic, + "Expected both dimensions to have spherical harmonic basis. Mesh = " + << mesh); + const intrp::Cardinal cardinal_interpolator(mesh, points); + const auto [n_th, n_ph] = mesh.extents().indices(); + const auto& [m_th, m_ph] = cardinal_interpolator.interpolation_matrices(); + for (size_t i_ph = 0, s = 0; i_ph < n_ph; ++i_ph) { + for (size_t i_th = 0; i_th < n_th; ++i_th) { + for (size_t k = 0; k < number_of_target_points; ++k) { + result(k, s) = m_th(k, i_th) * m_ph(2 * k, i_ph) + + m_th(k, i_th + n_th) * m_ph(2 * k + 1, i_ph); + } + ++s; + } + } + return result; } - // Not FD, so use spectral interpolation + // Not FD or special basis, so use 1D spectral interpolation matrices const std::array matrices{ {Spectral::interpolation_matrix(mesh.slice_through(0), get<0>(points)), Spectral::interpolation_matrix(mesh.slice_through(1), get<1>(points))}}; @@ -175,9 +193,31 @@ Matrix interpolation_matrix( } } return result; + } else if (mesh.basis()[1] == Spectral::Basis::SphericalHarmonic) { + ASSERT(mesh.basis()[2] == Spectral::Basis::SphericalHarmonic, + "Expected last two dimensions to each have spherical harmonic " + "basis. Mesh = " + << mesh); + const intrp::Cardinal cardinal_interpolator(mesh, points); + const auto [n_r, n_th, n_ph] = mesh.extents().indices(); + const auto& [m_r, m_th, m_ph] = + cardinal_interpolator.interpolation_matrices(); + for (size_t i_ph = 0, s = 0; i_ph < n_ph; ++i_ph) { + for (size_t i_th = 0; i_th < n_th; ++i_th) { + for (size_t i_r = 0; i_r < n_r; ++i_r) { + for (size_t k = 0; k < number_of_target_points; ++k) { + result(k, s) = + m_r(k, i_r) * (m_th(k, i_th) * m_ph(2 * k, i_ph) + + m_th(k, i_th + n_th) * m_ph(2 * k + 1, i_ph)); + } + ++s; + } + } + } + return result; } - // Not FD, so use spectral interpolation + // Not FD or special basis, so use 1D spectral interpolation matrices const std::array matrices{ {Spectral::interpolation_matrix(mesh.slice_through(0), get<0>(points)), Spectral::interpolation_matrix(mesh.slice_through(1), get<1>(points)), diff --git a/tests/Unit/NumericalAlgorithms/Interpolation/Test_IrregularInterpolant.cpp b/tests/Unit/NumericalAlgorithms/Interpolation/Test_IrregularInterpolant.cpp index ee850bd7b1f4..143588e6a102 100644 --- a/tests/Unit/NumericalAlgorithms/Interpolation/Test_IrregularInterpolant.cpp +++ b/tests/Unit/NumericalAlgorithms/Interpolation/Test_IrregularInterpolant.cpp @@ -30,6 +30,7 @@ #include "Domain/Structure/ElementId.hpp" #include "Framework/TestHelpers.hpp" #include "Helpers/DataStructures/MakeWithRandomValues.hpp" +#include "Helpers/NumericalAlgorithms/SphericalHarmonics/YlmTestFunctions.hpp" #include "NumericalAlgorithms/Interpolation/IrregularInterpolant.hpp" #include "NumericalAlgorithms/Spectral/Basis.hpp" #include "NumericalAlgorithms/Spectral/LogicalCoordinates.hpp" @@ -47,6 +48,26 @@ #include "Utilities/TaggedTuple.hpp" namespace { +// Polynomial of given degree with leading coefficinet \f$a_0\f$, and with +// \f$a_{n+1} = a_n / falloff +class Polynomial { + public: + Polynomial(const size_t degree, const double a_0, const double falloff) + : coefficients_(degree + 1) { + double n = falloff * a_0; + std::generate(coefficients_.begin(), coefficients_.end(), [&falloff, &n]() { + n /= falloff; + return n; + }); + } + Polynomial() = default; + DataVector operator()(const DataVector& x) const { + return evaluate_polynomial(coefficients_, x); + } + + private: + std::vector coefficients_; +}; using Affine = domain::CoordinateMaps::Affine; using Affine2D = domain::CoordinateMaps::ProductOf2Maps; @@ -411,6 +432,81 @@ void test_tov() { } } +void test_2d_spherical(const gsl::not_null generator) { + std::uniform_real_distribution<> xi_distribution(-1.0, 1.0); + std::uniform_real_distribution<> phi_distribution(0.0, 2.0 * M_PI); + for (size_t n_target_points = 1; n_target_points < 13; + n_target_points += 11) { + tnsr::I xi_target{n_target_points}; + get<0>(xi_target) = acos(make_with_random_values( + generator, make_not_null(&xi_distribution), xi_target)); + get<1>(xi_target) = make_with_random_values( + generator, make_not_null(&phi_distribution), xi_target); + for (size_t n_z = 0; n_z < 4; ++n_z) { + for (size_t n_y = 0; n_y < 4; ++n_y) { + for (size_t n_x = 0; n_x < 4; ++n_x) { + const Mesh<2> source_mesh{ + std::array{n_x + n_y + n_z + 2, 2 * (n_x + n_y) + 3}, + std::array{Spectral::Basis::SphericalHarmonic, + Spectral::Basis::SphericalHarmonic}, + std::array{Spectral::Quadrature::Gauss, + Spectral::Quadrature::Equiangular}}; + const YlmTestFunctions::ProductOfPolynomials f(n_x, n_y, n_z); + const auto xi_source = logical_coordinates(source_mesh); + const DataVector f_source = f(xi_source); + const DataVector f_expected = f(xi_target); + const intrp::Irregular<2> interpolator(source_mesh, xi_target); + const DataVector f_interpolated = interpolator.interpolate(f_source); + CHECK_ITERABLE_APPROX(f_interpolated, f_expected); + } + } + } + } +} + +void test_3d_spherical(const gsl::not_null generator) { + std::uniform_real_distribution<> xi_distribution(-1.0, 1.0); + std::uniform_real_distribution<> phi_distribution(0.0, 2.0 * M_PI); + for (size_t n_target_points = 1; n_target_points < 13; + n_target_points += 11) { + tnsr::I xi_target{n_target_points}; + get<0>(xi_target) = make_with_random_values( + generator, make_not_null(&xi_distribution), xi_target); + get<1>(xi_target) = acos(make_with_random_values( + generator, make_not_null(&xi_distribution), xi_target)); + get<2>(xi_target) = make_with_random_values( + generator, make_not_null(&phi_distribution), xi_target); + for (size_t n_r = 2; n_r < 4; ++n_r) { + for (size_t n_z = 0; n_z < 4; ++n_z) { + for (size_t n_y = 0; n_y < 4; ++n_y) { + for (size_t n_x = 0; n_x < 4; ++n_x) { + const Mesh<3> source_mesh{ + std::array{n_r, n_x + n_y + n_z + 2, 2 * (n_x + n_y) + 3}, + std::array{Spectral::Basis::Legendre, + Spectral::Basis::SphericalHarmonic, + Spectral::Basis::SphericalHarmonic}, + std::array{Spectral::Quadrature::GaussLobatto, + Spectral::Quadrature::Gauss, + Spectral::Quadrature::Equiangular}}; + const Polynomial f_r{n_r - 1, 1.5, 2.0}; + const YlmTestFunctions::ProductOfPolynomials f_a(n_x, n_y, n_z); + const auto xi_source = logical_coordinates(source_mesh); + const DataVector f_source = + f_r(get<0>(xi_source)) * + f_a(get<1>(xi_source), get<2>(xi_source)); + const DataVector f_expected = + f_r(get<0>(xi_target)) * + f_a(get<1>(xi_target), get<2>(xi_target)); + const intrp::Irregular<3> interpolator(source_mesh, xi_target); + const DataVector f_interpolated = + interpolator.interpolate(f_source); + CHECK_ITERABLE_APPROX(f_interpolated, f_expected); + } + } + } + } + } +} } // namespace SPECTRE_TEST_CASE("Unit.Numerical.Interpolation.IrregularInterpolant", @@ -429,4 +525,7 @@ SPECTRE_TEST_CASE("Unit.Numerical.Interpolation.IrregularInterpolant", test_polynomial_interpolant<3, 1>({{11, 9, 9}}); test_polynomial_interpolant<3, 1>({{11, 9, 13}}); test_tov(); + MAKE_GENERATOR(generator); + test_2d_spherical(make_not_null(&generator)); + test_3d_spherical(make_not_null(&generator)); }