Add HalfFourier basis functions

Represending functions on a half-circle could be useful to Cartoon evolutions;
in paticular for outer boundaries in GH where an angular basis could be more
suitable to constraint-preserving BC

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

Author: michaeldmurphy1

src/NumericalAlgorithms/Spectral/BasisFunctions/CMakeLists.txt

@@ -6,6 +6,7 @@ spectre_target_sources(
   PRIVATE
   Chebyshev.cpp
   Fourier.cpp
+  HalfFourier.cpp
   Jacobi.cpp
   Legendre.cpp
   Zernike.cpp
@@ -17,6 +18,7 @@ spectre_target_headers(
   HEADERS
   Chebyshev.hpp
   Fourier.hpp
+  HalfFourier.hpp
   Jacobi.hpp
   Zernike.hpp
   )

src/NumericalAlgorithms/Spectral/BasisFunctions/HalfFourier.cpp

@@ -0,0 +1,202 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#include "NumericalAlgorithms/Spectral/BasisFunctions/HalfFourier.hpp"
+
+#include <cmath>
+#include <cstddef>
+#include <numbers>
+
+#include "DataStructures/DataVector.hpp"
+#include "DataStructures/Matrix.hpp"
+#include "NumericalAlgorithms/Spectral/Basis.hpp"
+#include "NumericalAlgorithms/Spectral/BasisFunctionNormalizationSquare.hpp"
+#include "NumericalAlgorithms/Spectral/BasisFunctionValue.hpp"
+#include "NumericalAlgorithms/Spectral/CollocationPointsAndWeights.hpp"
+#include "NumericalAlgorithms/Spectral/Parity.hpp"
+#include "NumericalAlgorithms/Spectral/Quadrature.hpp"
+#include "Utilities/ContainerHelpers.hpp"
+#include "Utilities/ErrorHandling/Assert.hpp"
+#include "Utilities/GenerateInstantiations.hpp"
+
+namespace Spectral {
+
+DataVector HalfFourier::collocation_points(const size_t num_points) {
+  const double pi_over_n = std::numbers::pi / static_cast<double>(num_points);
+  DataVector result(num_points);
+  for (size_t j = 0; j < num_points; ++j) {
+    result[j] = pi_over_n * (static_cast<double>(j) + 0.5);
+  }
+  return result;
+}
+
+DataVector HalfFourier::quadrature_weights(const size_t num_points) {
+  return DataVector{num_points,
+                    std::numbers::pi / static_cast<double>(num_points)};
+}
+
+Matrix HalfFourier::even_differentiation_matrix(const size_t num_points) {
+  const double pi_over_n = std::numbers::pi / static_cast<double>(num_points);
+  const double two_over_n = 2.0 / static_cast<double>(num_points);
+  Matrix result(num_points, num_points, 0.0);
+  for (size_t i = 0; i < num_points; ++i) {
+    const double phi_i = pi_over_n * (static_cast<double>(i) + 0.5);
+    for (size_t j = 0; j < num_points; ++j) {
+      const double phi_j = pi_over_n * (static_cast<double>(j) + 0.5);
+      double val = 0.0;
+      for (size_t n = 1; n < num_points; ++n) {
+        val += static_cast<double>(n) * sin(static_cast<double>(n) * phi_i) *
+               cos(static_cast<double>(n) * phi_j);
+      }
+      result(i, j) = -two_over_n * val;
+    }
+  }
+  return result;
+}
+
+Matrix HalfFourier::odd_differentiation_matrix(const size_t num_points) {
+  const double pi_over_n = std::numbers::pi / static_cast<double>(num_points);
+  const double two_over_n = 2.0 / static_cast<double>(num_points);
+  Matrix result(num_points, num_points, 0.0);
+  for (size_t i = 0; i < num_points; ++i) {
+    const double phi_i = pi_over_n * (static_cast<double>(i) + 0.5);
+    for (size_t j = 0; j < num_points; ++j) {
+      const double phi_j = pi_over_n * (static_cast<double>(j) + 0.5);
+      double val = 0.0;
+      // Note: the n=num_points term vanishes because cos(num_points * phi_i) =
+      // cos(pi*(i+0.5)) = 0, so only n=1,...,num_points-1 contribute.
+      for (size_t n = 1; n < num_points; ++n) {
+        val += static_cast<double>(n) * cos(static_cast<double>(n) * phi_i) *
+               sin(static_cast<double>(n) * phi_j);
+      }
+      result(i, j) = two_over_n * val;
+    }
+  }
+  return result;
+}
+
+template <typename T>
+Matrix HalfFourier::even_interpolation_matrix(const size_t num_points,
+                                              const T& target_points) {
+  const double pi_over_n = std::numbers::pi / static_cast<double>(num_points);
+  const double inv_n = 1.0 / static_cast<double>(num_points);
+  const size_t num_target_points = get_size(target_points);
+  Matrix result(num_target_points, num_points);
+  for (size_t i = 0; i < num_target_points; ++i) {
+    const double x = get_element(target_points, i);
+    for (size_t j = 0; j < num_points; ++j) {
+      const double phi_j = pi_over_n * (static_cast<double>(j) + 0.5);
+      double val = 1.0;
+      for (size_t n = 1; n < num_points; ++n) {
+        val += 2.0 * cos(static_cast<double>(n) * x) *
+               cos(static_cast<double>(n) * phi_j);
+      }
+      result(i, j) = inv_n * val;
+    }
+  }
+  return result;
+}
+
+template <typename T>
+Matrix HalfFourier::odd_interpolation_matrix(const size_t num_points,
+                                             const T& target_points) {
+  const double pi_over_n = std::numbers::pi / static_cast<double>(num_points);
+  const double two_over_n = 2.0 / static_cast<double>(num_points);
+  const double inv_n = 1.0 / static_cast<double>(num_points);
+  const size_t num_target_points = get_size(target_points);
+  Matrix result(num_target_points, num_points);
+  for (size_t i = 0; i < num_target_points; ++i) {
+    const double x = get_element(target_points, i);
+    for (size_t j = 0; j < num_points; ++j) {
+      const double phi_j = pi_over_n * (static_cast<double>(j) + 0.5);
+      double val = 0.0;
+      // Modes n=1,...,N-1 have discrete norm N/2 → coefficient 2/N.
+      // Mode n=N (Nyquist) has discrete norm N → coefficient 1/N.
+      for (size_t n = 1; n < num_points; ++n) {
+        val += sin(static_cast<double>(n) * x) *
+               sin(static_cast<double>(n) * phi_j);
+      }
+      result(i, j) =
+          two_over_n * val + inv_n * sin(static_cast<double>(num_points) * x) *
+                                 sin(static_cast<double>(num_points) * phi_j);
+    }
+  }
+  return result;
+}
+
+template <typename T>
+Matrix HalfFourier::interpolation_matrix(const size_t num_points,
+                                         const T& target_points,
+                                         const Parity parity) {
+  ASSERT(parity != Parity::Uninitialized,
+         "Parity must be set to either Even or Odd");
+  if (parity == Parity::Even) {
+    return even_interpolation_matrix(num_points, target_points);
+  } else {
+    return odd_interpolation_matrix(num_points, target_points);
+  }
+}
+
+// Specializations of function templates defined in the Spectral directory
+
+template <>
+std::pair<DataVector, DataVector> compute_collocation_points_and_weights<
+    Basis::HalfFourier, Quadrature::Equiangular>(const size_t num_points) {
+  return std::make_pair(HalfFourier::collocation_points(num_points),
+                        HalfFourier::quadrature_weights(num_points));
+}
+
+template <Basis BasisType>
+Matrix spectral_indefinite_integral_matrix(size_t num_points);
+
+#if defined(__GNUC__) && !defined(__clang__)
+#pragma GCC diagnostic push
+#pragma GCC diagnostic ignored "-Wsuggest-attribute=noreturn"
+#endif
+
+template <>
+Matrix spectral_indefinite_integral_matrix<Basis::HalfFourier>(
+    size_t /*num_points*/) {
+  ERROR("Indefinite integral matrix is not implemented for HalfFourier basis");
+}
+
+template <>
+DataVector compute_basis_function_value<Basis::HalfFourier>(
+    const size_t /*k*/, const DataVector& /*x*/) {
+  ERROR("HalfFourier basis function value requires a parity argument");
+}
+
+template <>
+double compute_basis_function_value<Basis::HalfFourier>(const size_t /*k*/,
+                                                        const double& /*x*/) {
+  ERROR("HalfFourier basis function value requires a parity argument");
+}
+
+template <>
+double compute_basis_function_normalization_square<Basis::HalfFourier>(
+    const size_t /*k*/) {
+  ERROR(
+      "HalfFourier normalization square requires a parity argument; "
+      "norms are pi (even k=0), pi/2 (even k>=1 or odd k=1,...,N).");
+}
+
+#if defined(__GNUC__) && !defined(__clang__)
+#pragma GCC diagnostic pop
+#endif
+
+#define GET_TYPE(data) BOOST_PP_TUPLE_ELEM(0, data)
+
+#define INSTANTIATE_INTERPOLATION(r, data)                \
+  template Matrix HalfFourier::even_interpolation_matrix( \
+      size_t num_points, const GET_TYPE(data)&);          \
+  template Matrix HalfFourier::odd_interpolation_matrix(  \
+      size_t num_points, const GET_TYPE(data)&);          \
+  template Matrix HalfFourier::interpolation_matrix(      \
+      size_t, const GET_TYPE(data)&, Parity);
+
+GENERATE_INSTANTIATIONS(INSTANTIATE_INTERPOLATION,
+                        (double, DataVector, std::vector<double>))
+
+#undef INSTANTIATE_INTERP
+#undef GET_TYPE
+}  // namespace Spectral

src/NumericalAlgorithms/Spectral/BasisFunctions/HalfFourier.hpp

@@ -0,0 +1,144 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#pragma once
+
+#include <cstddef>
+
+/// \cond
+class DataVector;
+class Matrix;
+namespace Spectral {
+enum class Parity : uint8_t;
+}  // namespace Spectral
+/// \endcond
+
+namespace Spectral {
+
+/*!
+ * \ingroup SpectralGroup
+ *
+ * \brief A collection of helper functions for the half-Fourier spectral basis
+ *
+ * \details The half-Fourier basis represents functions on the interval
+ * \f$\phi \in [0, \pi)\f$ using \f$N\f$ equispaced interior collocation
+ * points \f$\phi_j = (j + \tfrac{1}{2})\pi/N\f$ for \f$j = 0, \ldots, N-1\f$.
+ *
+ * Functions of even parity under \f$\phi \to -\phi\f$ (i.e. those satisfying
+ * \f$f(-\phi)=f(\phi)\f$) are expanded in cosines:
+ * \f[
+ * f(\phi) = \sum_{n=0}^{N-1} a_n \cos(n\phi)
+ * \f]
+ *
+ * Functions of odd parity under \f$\phi \to -\phi\f$ (i.e. those satisfying
+ * \f$f(-\phi)=-f(\phi)\f$) are expanded in sines:
+ * \f[
+ * f(\phi) = \sum_{n=1}^{N} b_n \sin(n\phi)
+ * \f]
+ *
+ * The derivative \f$\partial / \partial \phi\f$ maps even-parity functions to
+ * odd-parity functions and vice versa.
+ *
+ * This basis is intended for use in the Cartoon method for axisymmetric
+ * problems on a cylinder, where the azimuthal direction covers only half a
+ * circle due to the reflection symmetry, and the parity boundary conditions
+ * are internal to the spectral representation.
+ */
+class HalfFourier {
+ public:
+  /*!
+   * \brief Collocation points \f$\{\phi_j\}\f$
+   *
+   * \details The collocation points on the interval \f$(0, \pi)\f$ are given
+   * by
+   * \f[
+   * \phi_j = \frac{(j + \tfrac{1}{2})\pi}{N}
+   * \f]
+   */
+  static DataVector collocation_points(size_t num_points);
+
+  /*!
+   * \brief Quadrature weights \f$\{w_j\}\f$
+   *
+   * \details The quadrature weights are uniform:
+   * \f[
+   * w_j = \frac{\pi}{N}
+   * \f]
+   */
+  static DataVector quadrature_weights(size_t num_points);
+
+  /*!
+   * \brief Differentiation matrix \f$D^{\text{even}}_{ij}\f$ for even-parity
+   * functions
+   *
+   * \details Maps an even-parity function (expanded in cosines) to its
+   * derivative, which is an odd-parity function (expanded in sines).
+   * Explicitly:
+   * \f[
+   * D^{\text{even}}_{ij} = \frac{2}{N} \sum_{n=1}^{N-1}
+   * (-n) \sin(n\phi_i) \cos(n\phi_j)
+   * \f]
+   */
+  static Matrix even_differentiation_matrix(size_t num_points);
+
+  /*!
+   * \brief Differentiation matrix \f$D^{\text{odd}}_{ij}\f$ for odd-parity
+   * functions
+   *
+   * \details Maps an odd-parity function (expanded in sines) to its
+   * derivative, which is an even-parity function (expanded in cosines).
+   * Explicitly:
+   * \f[
+   * D^{\text{odd}}_{ij} = \frac{2}{N} \sum_{n=1}^{N-1}
+   * n \cos(n\phi_i) \sin(n\phi_j)
+   * \f]
+   *
+   * Note that \f$D^{\text{even}} = -(D^{\text{odd}})^T\f$.
+   */
+  static Matrix odd_differentiation_matrix(size_t num_points);
+
+  /*!
+   * \brief Interpolation matrix for even-parity functions to
+   * \p target_points.
+   *
+   * \details Using the discrete cosine transform (DCT-II) representation, the
+   * interpolation weights at a target point \f$x\f$ are:
+   * \f[
+   * I^{\text{even}}_j(x) = \frac{1}{N}\left[1 + 2\sum_{n=1}^{N-1}
+   * \cos(nx)\cos(n\phi_j)\right]
+   * \f]
+   */
+  template <typename T>
+  static Matrix even_interpolation_matrix(size_t num_points,
+                                          const T& target_points);
+
+  /*!
+   * \brief Interpolation matrix for odd-parity functions to
+   * \p target_points.
+   *
+   * \details Using the discrete sine transform (DST-II) representation, the
+   * interpolation weights at a target point \f$x\f$ are:
+   * \f[
+   * I^{\text{odd}}_j(x) = \frac{2}{N}\sum_{n=1}^{N-1}
+   * \sin(nx)\sin(n\phi_j) + \frac{1}{N}\sin(Nx)\sin(N\phi_j)
+   * \f]
+   * where the Nyquist mode \f$n=N\f$ carries half the weight of the other
+   * modes.
+   */
+  template <typename T>
+  static Matrix odd_interpolation_matrix(size_t num_points,
+                                         const T& target_points);
+
+  /*!
+   * \brief %Matrix used to interpolate to the \p target_points for a given
+   * \p parity.
+   *
+   * \details Dispatches to `even_interpolation_matrix` when \p parity is
+   * `Parity::Even` and to `odd_interpolation_matrix` when it is `Parity::Odd`.
+   * This mirrors the interface of `Zernike<1>::interpolation_matrix`.
+   */
+  template <typename T>
+  static Matrix interpolation_matrix(size_t num_points, const T& target_points,
+                                     Parity parity);
+};
+}  // namespace Spectral

src/NumericalAlgorithms/Spectral/CollocationPoints.cpp

@@ -47,6 +47,8 @@ template const DataVector& collocation_points<Basis::FiniteDifference,
                                               Quadrature::FaceCentered>(size_t);
 template const DataVector&
     collocation_points<Basis::Fourier, Quadrature::Equiangular>(size_t);
+template const DataVector&
+    collocation_points<Basis::HalfFourier, Quadrature::Equiangular>(size_t);
 template const DataVector&
     collocation_points<Basis::ZernikeB1, Quadrature::GaussRadauUpper>(size_t);
 template const DataVector&

src/NumericalAlgorithms/Spectral/CollocationPointsAndWeights.cpp

@@ -37,6 +37,8 @@ template struct CollocationPointsAndWeightsGenerator<Basis::FiniteDifference,
                                                      Quadrature::FaceCentered>;
 template struct CollocationPointsAndWeightsGenerator<Basis::Fourier,
                                                      Quadrature::Equiangular>;
+template struct CollocationPointsAndWeightsGenerator<Basis::HalfFourier,
+                                                     Quadrature::Equiangular>;
 template struct CollocationPointsAndWeightsGenerator<
     Basis::ZernikeB1, Quadrature::GaussRadauUpper>;
 template struct CollocationPointsAndWeightsGenerator<