Add cardinal interpolator

This interpolator interpolates dimension by dimension and can handle
spherical harmonic interpolation.

src/NumericalAlgorithms/Interpolation/CMakeLists.txt

@@ -10,6 +10,7 @@ spectre_target_sources(
   PRIVATE
   BarycentricRational.cpp
   BarycentricRationalSpanInterpolator.cpp
+  CardinalInterpolator.cpp
   CubicSpanInterpolator.cpp
   CubicSpline.cpp
   InterpolationWeights.cpp
@@ -30,6 +31,7 @@ spectre_target_headers(
   HEADERS
   BarycentricRational.hpp
   BarycentricRationalSpanInterpolator.hpp
+  CardinalInterpolator.hpp
   CubicSpanInterpolator.hpp
   CubicSpline.hpp
   InterpolationWeights.hpp

src/NumericalAlgorithms/Interpolation/CardinalInterpolator.cpp

@@ -0,0 +1,211 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#include "CardinalInterpolator.hpp"
+
+#include <array>
+#include <cstddef>
+
+#include "DataStructures/DataVector.hpp"
+#include "DataStructures/Matrix.hpp"
+#include "DataStructures/Tensor/Tensor.hpp"
+#include "NumericalAlgorithms/Interpolation/InterpolationWeights.hpp"
+#include "NumericalAlgorithms/Spectral/LogicalCoordinates.hpp"
+#include "NumericalAlgorithms/Spectral/Mesh.hpp"
+#include "NumericalAlgorithms/Spectral/Spectral.hpp"
+#include "Utilities/Blas.hpp"
+#include "Utilities/ErrorHandling/Assert.hpp"
+#include "Utilities/Gsl.hpp"
+
+namespace intrp {
+template <size_t Dim>
+Cardinal<Dim>::Cardinal(
+    const Mesh<Dim>& source_mesh,
+    const tnsr::I<DataVector, Dim, Frame::ElementLogical>& target_points)
+    : n_target_points_(get<0>(target_points).size()),
+      source_mesh_(source_mesh) {
+  for (size_t d = 0; d < Dim; ++d) {
+    switch (source_mesh_.basis(d)) {
+      case Spectral::Basis::Chebyshev:
+        [[fallthrough]];
+      case Spectral::Basis::Legendre: {
+        const DataVector xi =
+            get<0>(logical_coordinates(source_mesh_.slice_through(d)));
+        gsl::at(interpolation_matrices_, d) =
+            fornberg_interpolation_matrix(target_points.get(d), xi);
+        break;
+      }
+      case Spectral::Basis::SphericalHarmonic: {
+        using_spherical_harmonics_ = true;
+        switch (source_mesh_.quadrature(d)) {
+          case Spectral::Quadrature::Gauss: {
+            const DataVector theta =
+                get<0>(logical_coordinates(source_mesh_.slice_through(d)));
+            const DataVector cos_theta_source = cos(theta);
+            const DataVector cos_theta_target = cos(target_points.get(d));
+            const Matrix theta_matrix = fornberg_interpolation_matrix(
+                cos_theta_target, cos_theta_source);
+            const size_t n_th = source_mesh_.extents(d);
+            gsl::at(interpolation_matrices_, d)
+                .resize(n_target_points_, 2 * n_th);
+            const DataVector csc_theta_source = 1.0 / sin(theta);
+            const DataVector sin_theta_target = sin(target_points.get(d));
+            for (size_t k = 0; k < n_target_points_; ++k) {
+              for (size_t i_th = 0; i_th < n_th; ++i_th) {
+                const double factor =
+                    0.5 * sin_theta_target[k] * csc_theta_source[i_th];
+                gsl::at(interpolation_matrices_, d)(k, i_th) =
+                    theta_matrix(k, i_th) * (0.5 + factor);
+                gsl::at(interpolation_matrices_, d)(k, n_th + i_th) =
+                    theta_matrix(k, i_th) * (0.5 - factor);
+              }
+            }
+            break;
+          }
+          case Spectral::Quadrature::Equiangular: {
+            const DataVector phi =
+                get<0>(logical_coordinates(source_mesh_.slice_through(d)));
+            const DataVector& phi_target = target_points.get(d);
+            DataVector extended_phi_target{2 * n_target_points_};
+            for (size_t k = 0; k < n_target_points_; ++k) {
+              extended_phi_target[2 * k] = phi_target[k];
+              extended_phi_target[2 * k + 1] = phi_target[k] + M_PI;
+            }
+            gsl::at(interpolation_matrices_, d) = fourier_interpolation_matrix(
+                extended_phi_target, source_mesh_.extents(d));
+            break;
+          }
+          default:
+            ERROR(
+                "Quadrature must be Gauss or Equiangular for Basis "
+                "SphericalHarmonic, not "
+                << source_mesh_.quadrature(d));
+        }
+        break;
+      }
+      default:
+        ERROR("Basis " << source_mesh_.basis(d)
+                       << " is not supported by the Cardinal interpolator.");
+    }
+  }
+}
+
+template <>
+DataVector Cardinal<1>::interpolate(const DataVector& f_source) const {
+  const size_t n_source_points = source_mesh_.number_of_grid_points();
+  ASSERT(f_source.size() == n_source_points,
+         "Size of source data ("
+             << f_source.size() << ") does not match size of source mesh ("
+             << n_source_points
+             << ") that this interpolator was constructed with.");
+  DataVector result{n_target_points_};
+  dgemv_('N', n_target_points_, n_source_points, 1.0,
+         interpolation_matrices_[0].data(),
+         interpolation_matrices_[0].spacing(), f_source.data(), 1, 0.0,
+         result.data(), 1);
+  return result;
+}
+
+template <>
+DataVector Cardinal<2>::interpolate(const DataVector& f_source) const {
+  ASSERT(f_source.size() == source_mesh_.number_of_grid_points(),
+         "Size of source data ("
+             << f_source.size() << ") does not match size of source mesh ("
+             << source_mesh_.number_of_grid_points()
+             << ") that this interpolator was constructed with.");
+  DataVector result{n_target_points_};
+  if (using_spherical_harmonics_) {
+    const auto [n_th, n_ph] = source_mesh_.extents().indices();
+    DataVector intermediate_result{2 * n_th};
+    for (size_t k = 0; k < n_target_points_; ++k) {
+      dgemv_('N', n_th, n_ph, 1.0, f_source.data(), n_th,
+             interpolation_matrices_[1].data() + 2 * k, 2 * n_target_points_,
+             0.0, intermediate_result.data(), 1);
+      dgemv_('N', n_th, n_ph, 1.0, f_source.data(), n_th,
+             interpolation_matrices_[1].data() + (2 * k + 1),
+             2 * n_target_points_, 0.0, intermediate_result.data() + n_th, 1);
+      result[k] = ddot_(2 * n_th, interpolation_matrices_[0].data() + k,
+                        n_target_points_, intermediate_result.data(), 1);
+    }
+    return result;
+  }
+  const auto [n_xi, n_eta] = source_mesh_.extents().indices();
+  Matrix intermediate_result{n_target_points_, n_eta};
+  dgemm_('N', 'N', n_target_points_, n_eta, n_xi, 1.0,
+         interpolation_matrices_[0].data(),
+         interpolation_matrices_[0].spacing(), f_source.data(), n_xi, 0.0,
+         intermediate_result.data(), n_target_points_);
+  for (size_t k = 0; k < n_target_points_; ++k) {
+    result[k] =
+        ddot_(n_eta, interpolation_matrices_[1].data() + k, n_target_points_,
+              intermediate_result.data() + k, n_target_points_);
+  }
+  return result;
+}
+
+template <>
+DataVector Cardinal<3>::interpolate(const DataVector& f_source) const {
+  ASSERT(f_source.size() == source_mesh_.number_of_grid_points(),
+         "Size of source data ("
+             << f_source.size() << ") does not match size of source mesh ("
+             << source_mesh_.number_of_grid_points()
+             << ") that this interpolator was constructed with.");
+  const auto [n0, n1, n2] = source_mesh_.extents().indices();
+  Matrix intermediate_matrix{n_target_points_, n1 * n2};
+  dgemm_('N', 'N', n_target_points_, n1 * n2, n0, 1.0,
+         interpolation_matrices_[0].data(),
+         interpolation_matrices_[0].spacing(), f_source.data(), n0, 0.0,
+         intermediate_matrix.data(), n_target_points_);
+  auto transpose = intermediate_matrix.transpose();
+  DataVector result{n_target_points_};
+  if (using_spherical_harmonics_) {
+    DataVector intermediate_vector{2 * n1};
+    for (size_t k = 0; k < n_target_points_; ++k) {
+      dgemv_('N', n1, n2, 1.0, transpose.data() + k * n1 * n2, n1,
+             interpolation_matrices_[2].data() + 2 * k, 2 * n_target_points_,
+             0.0, intermediate_vector.data(), 1);
+      dgemv_('N', n1, n2, 1.0, transpose.data() + k * n1 * n2, n1,
+             interpolation_matrices_[2].data() + (2 * k + 1),
+             2 * n_target_points_, 0.0, intermediate_vector.data() + n1, 1);
+      result[k] = ddot_(2 * n1, interpolation_matrices_[1].data() + k,
+                        n_target_points_, intermediate_vector.data(), 1);
+    }
+    return result;
+  }
+  DataVector intermediate_vector{n2};
+  for (size_t k = 0; k < n_target_points_; ++k) {
+    dgemv_('T', n1, n2, 1.0, transpose.data() + k * n1 * n2, n1,
+           interpolation_matrices_[1].data() + k, n_target_points_, 0.0,
+           intermediate_vector.data(), 1);
+    result[k] = ddot_(n2, interpolation_matrices_[2].data() + k,
+                      n_target_points_, intermediate_vector.data(), 1);
+  }
+  return result;
+}
+
+template <size_t Dim>
+const std::array<Matrix, Dim>& Cardinal<Dim>::interpolation_matrices() const {
+  return interpolation_matrices_;
+}
+
+template <size_t Dim>
+bool operator==(const Cardinal<Dim>& lhs, const Cardinal<Dim>& rhs) {
+  return lhs.interpolation_matrices() == rhs.interpolation_matrices();
+}
+
+template <size_t Dim>
+bool operator!=(const Cardinal<Dim>& lhs, const Cardinal<Dim>& rhs) {
+  return not(lhs == rhs);
+}
+
+template class Cardinal<1>;
+template class Cardinal<2>;
+template class Cardinal<3>;
+template bool operator==(const Cardinal<1>& lhs, const Cardinal<1>& rhs);
+template bool operator==(const Cardinal<2>& lhs, const Cardinal<2>& rhs);
+template bool operator==(const Cardinal<3>& lhs, const Cardinal<3>& rhs);
+template bool operator!=(const Cardinal<1>& lhs, const Cardinal<1>& rhs);
+template bool operator!=(const Cardinal<2>& lhs, const Cardinal<2>& rhs);
+template bool operator!=(const Cardinal<3>& lhs, const Cardinal<3>& rhs);
+
+}  // namespace intrp

src/NumericalAlgorithms/Interpolation/CardinalInterpolator.hpp

@@ -0,0 +1,68 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#pragma once
+
+#include <array>
+#include <cstddef>
+
+#include "DataStructures/Matrix.hpp"
+#include "DataStructures/Tensor/TypeAliases.hpp"
+
+/// \cond
+class DataVector;
+template <size_t Dim>
+class Mesh;
+namespace PUP {
+class er;
+}  // namespace PUP
+/// \endcond
+
+namespace intrp {
+/*!
+ * \brief Interpolates by doing partial summation in each dimension using
+ * one-dimensional interpolation
+ *
+ * \details The one-dimensional matrices used to do the interpolation depend
+ * upon the Spectral::Basis used in each dimension:
+ * - For a Chebyshev or Legendre basis, the matrices are given by
+ * intrp::fornberg_interpolation_matrix at the quadrature points of the
+ * source_mesh.  (These are equivalent to those returned by
+ * Spectral::interpolation_matrix.)
+ * - For a Fourier basis, the matrix is given by
+ * intrp::fourier_interpolation_matrix at the quadrature points of the
+ * source_mesh
+ *
+ */
+template <size_t Dim>
+class Cardinal {
+ public:
+  Cardinal(
+      const Mesh<Dim>& source_mesh,
+      const tnsr::I<DataVector, Dim, Frame::ElementLogical>& target_points);
+
+  Cardinal();
+
+  /// Interpolates the function `f` provided on the `source_mesh` to the
+  /// `target_points` with which the interpolator was constructed.
+  DataVector interpolate(const DataVector& f) const;
+
+  /// The one-dimensional interpolation matrices used to do the interpolation
+  const std::array<Matrix, Dim>& interpolation_matrices() const;
+
+  // NOLINTNEXTLINE(google-runtime-references)
+  void pup(PUP::er& p);
+
+ private:
+  size_t n_target_points_;
+  Mesh<Dim> source_mesh_;
+  std::array<Matrix, Dim> interpolation_matrices_;
+  bool using_spherical_harmonics_{false};
+};
+
+template <size_t Dim>
+bool operator==(const Cardinal<Dim>& lhs, const Cardinal<Dim>& rhs);
+
+template <size_t Dim>
+bool operator!=(const Cardinal<Dim>& lhs, const Cardinal<Dim>& rhs);
+}  // namespace intrp

tests/Unit/NumericalAlgorithms/Interpolation/CMakeLists.txt

@@ -7,6 +7,7 @@ set(LIBRARY "Test_NumericalInterpolation")
 
 set(LIBRARY_SOURCES
   Test_BarycentricRational.cpp
+  Test_CardinalInterpolator.cpp
   Test_CubicSpline.cpp
   Test_InterpolationWeights.cpp
   Test_IrregularInterpolant.cpp
@@ -40,6 +41,7 @@ target_link_libraries(
   ParallelInterpolation
   RelativisticEulerSolutions
   Spectral
+  SphericalHarmonicsHelpers
   SpinWeightedSphericalHarmonics
   Utilities
   )