Enable interpolation for shells using spherical harmonics

docs/References.bib

@@ -1003,6 +1003,19 @@ @ARTICLE{Font1998
   adsurl =       {https://ui.adsabs.harvard.edu/abs/1998ApJ...494..297F},
 }
 
+@article{Fornberg1998,
+  author =       {Fornberg, Bengt},
+  title =        {Classroom Note: Calculation of Weights in Finite
+                  Difference Formulas},
+  journal =      {SIAM Review},
+  volume =       40,
+  number =       3,
+  pages =        {685-691},
+  year =         1998,
+  doi =          {10.1137/S0036144596322507},
+  url =          "https://doi.org/10.1137/S0036144596322507"
+}
+
 @article{Fortunato2019jl,
   title     = "Efficient Operator-Coarsening Multigrid Schemes for Local
                Discontinuous {Galerkin} Methods",

src/NumericalAlgorithms/Interpolation/CMakeLists.txt

@@ -10,8 +10,10 @@ spectre_target_sources(
   PRIVATE
   BarycentricRational.cpp
   BarycentricRationalSpanInterpolator.cpp
+  CardinalInterpolator.cpp
   CubicSpanInterpolator.cpp
   CubicSpline.cpp
+  InterpolationWeights.cpp
   IrregularInterpolant.cpp
   LinearLeastSquares.cpp
   LinearRegression.cpp
@@ -29,8 +31,10 @@ spectre_target_headers(
   HEADERS
   BarycentricRational.hpp
   BarycentricRationalSpanInterpolator.hpp
+  CardinalInterpolator.hpp
   CubicSpanInterpolator.hpp
   CubicSpline.hpp
+  InterpolationWeights.hpp
   IrregularInterpolant.hpp
   LagrangePolynomial.hpp
   LinearLeastSquares.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

src/NumericalAlgorithms/Interpolation/InterpolationWeights.cpp

@@ -0,0 +1,70 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#include "NumericalAlgorithms/Interpolation/InterpolationWeights.hpp"
+
+#include <cmath>
+#include <cstddef>
+
+#include "DataStructures/DataVector.hpp"
+#include "DataStructures/Matrix.hpp"
+
+namespace intrp {
+Matrix fornberg_interpolation_matrix(const DataVector& x_target,
+                                     const DataVector& x_source) {
+  Matrix result{x_target.size(), x_source.size()};
+  const size_t n_source_pts = x_source.size();
+  for (size_t k = 0; k < x_target.size(); ++k) {
+    double c1 = 1.0;
+    double c4 = x_source[0] - x_target[k];
+    result(k, 0) = 1.0;
+    for (size_t i = 1; i < n_source_pts; ++i) {
+      double c2 = 1.0;
+      const double c5 = c4;
+      c4 = x_source[i] - x_target[k];
+      for (size_t j = 0; j < i; ++j) {
+        const double c3 = x_source[i] - x_source[j];
+        c2 *= c3;
+        if (j + 1 == i) {
+          result(k, i) = -c1 * c5 * result(k, i - 1) / c2;
+        }
+        result(k, j) = c4 * result(k, j) / c3;
+      }
+      c1 = c2;
+    }
+  }
+  return result;
+}
+
+Matrix fourier_interpolation_matrix(const DataVector& x_target,
+                                    const size_t n_source_points) {
+  if (n_source_points == 1) {
+    return Matrix{x_target.size(), 1, 1.0};
+  }
+  DataVector x_source{n_source_points,
+                      2.0 * M_PI / static_cast<double>(n_source_points)};
+  for (size_t i = 0; i < n_source_points; ++i) {
+    x_source[i] *= static_cast<double>(i);
+  }
+  Matrix result{x_target.size(), n_source_points,
+                1.0 / static_cast<double>(n_source_points)};
+  const bool n_source_points_is_even = n_source_points % 2 == 0;
+  for (size_t k = 0; k < x_target.size(); ++k) {
+    for (size_t i = 0; i < n_source_points; ++i) {
+      double c0 = 1.0;
+      double c1 = cos(x_target[k] - x_source[i]);
+      const double cdx2 = 2.0 * c1;
+      double sum = c1;
+      for (size_t j = 2; j <= n_source_points / 2; ++j) {
+        const double tmp = c1;
+        c1 = cdx2 * c1 - c0;
+        c0 = tmp;
+        sum += c1;
+      }
+      result(k, i) *=
+          n_source_points_is_even ? 2.0 * sum + 1.0 - c1 : 2.0 * sum + 1.0;
+    }
+  }
+  return result;
+}
+}  // namespace intrp