Enable Ylm meshes

src/NumericalAlgorithms/Spectral/CMakeLists.txt

@@ -18,6 +18,7 @@ spectre_target_sources(
   Projection.cpp
   Quadrature.cpp
   Spectral.cpp
+  SphericalHarmonics.cpp
   )
 
 spectre_target_headers(

src/NumericalAlgorithms/Spectral/Mesh.cpp

@@ -33,8 +33,6 @@ template <size_t Dim>
 Mesh<Dim>::Mesh(const size_t isotropic_extents, const Spectral::Basis basis,
                 const Spectral::Quadrature quadrature) {
   static_assert(sizeof(Mesh<Dim>) == (Dim == 0 ? 2 : 2 * Dim));
-  ASSERT(basis != Spectral::Basis::SphericalHarmonic,
-         "SphericalHarmonic is not a valid basis for the Mesh");
   ASSERT(isotropic_extents <= 255, "Cannot have more than 255 grid points");
   extents_[0] = isotropic_extents;
   quadrature_and_basis_[0] = combine(basis, quadrature);
@@ -52,8 +50,6 @@ template <size_t Dim>
 Mesh<Dim>::Mesh(const std::array<size_t, Dim>& extents,
                 const Spectral::Basis basis,
                 const Spectral::Quadrature quadrature) {
-  ASSERT(basis != Spectral::Basis::SphericalHarmonic,
-         "SphericalHarmonic is not a valid basis for the Mesh");
   if constexpr (Dim > 0) {
     ASSERT(
         extents[0] <= 255,
@@ -86,10 +82,6 @@ template <size_t Dim>
 Mesh<Dim>::Mesh(const std::array<size_t, Dim>& extents,
                 const std::array<Spectral::Basis, Dim>& bases,
                 const std::array<Spectral::Quadrature, Dim>& quadratures) {
-  for (auto it = bases.begin(); it != bases.end(); it++) {
-    ASSERT(*it != Spectral::Basis::SphericalHarmonic,
-           "SphericalHarmonic is not a valid basis for the Mesh");
-  }
   if constexpr (Dim > 0) {
     ASSERT(
         extents[0] <= 255,

src/NumericalAlgorithms/Spectral/Spectral.cpp

@@ -428,6 +428,8 @@ struct DifferentiationMatrixGenerator {
           diff_matrix(i, j) *= inv_delta;
         }
       }
+    } else if constexpr (BasisType == Basis::SphericalHarmonic) {
+      ERROR("Not implemented");
     } else {
       const DataVector& bary_weights =
           barycentric_weights<BasisType, QuadratureType>(num_points);
@@ -797,9 +799,7 @@ template <typename F>
 decltype(auto) get_spectral_quantity_for_mesh(F&& f, const Mesh<1>& mesh) {
   const auto num_points = mesh.extents(0);
   // Switch on runtime values of basis and quadrature to select
-  // corresponding template specialization. For basis functions spanning
-  // multiple dimensions we can generalize this function to take a
-  // higher-dimensional Mesh.
+  // corresponding template specialization.
   switch (mesh.basis(0)) {
     case Basis::Legendre:
       switch (mesh.quadrature(0)) {
@@ -829,6 +829,20 @@ decltype(auto) get_spectral_quantity_for_mesh(F&& f, const Mesh<1>& mesh) {
         default:
           ERROR("Missing quadrature case for spectral quantity");
       }
+    case Basis::SphericalHarmonic:
+      switch (mesh.quadrature(0)) {
+        case Quadrature::Gauss:  // [0, pi] direction
+          return f(std::integral_constant<Basis, Basis::SphericalHarmonic>{},
+                   std::integral_constant<Quadrature, Quadrature::Gauss>{},
+                   num_points);
+        case Quadrature::Equiangular:  // [0, 2 pi) direction
+          return f(
+              std::integral_constant<Basis, Basis::SphericalHarmonic>{},
+              std::integral_constant<Quadrature, Quadrature::Equiangular>{},
+              num_points);
+        default:
+          ERROR("Missing quadrature case for spectral quantity");
+      }
     case Basis::FiniteDifference:
       switch (mesh.quadrature(0)) {
         case Quadrature::CellCentered:
@@ -932,6 +946,12 @@ GENERATE_INSTANTIATIONS(INSTANTIATE,
 #undef QUAD
 #undef INSTANTIATE
 
+template const DataVector& Spectral::collocation_points<
+    Spectral::Basis::SphericalHarmonic, Spectral::Quadrature::Gauss>(size_t);
+template const DataVector&
+    Spectral::collocation_points<Spectral::Basis::SphericalHarmonic,
+                                 Spectral::Quadrature::Equiangular>(size_t);
+
 template const DataVector&
     Spectral::collocation_points<Spectral::Basis::FiniteDifference,
                                  Spectral::Quadrature::CellCentered>(size_t);

src/NumericalAlgorithms/Spectral/SphericalHarmonics.cpp

@@ -0,0 +1,60 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#include "NumericalAlgorithms/Spectral/Spectral.hpp"
+
+#include <cstddef>
+
+#include "DataStructures/DataVector.hpp"
+#include "DataStructures/Matrix.hpp"
+#include "Utilities/ErrorHandling/Error.hpp"
+
+namespace Spectral {
+
+// Algorithms to compute spherical harmonic basis functions
+// These functions specialize the templates declared in `Spectral.hpp`.
+
+template <>
+std::pair<DataVector, DataVector> compute_collocation_points_and_weights<
+    Basis::SphericalHarmonic, Quadrature::Gauss>(const size_t num_points) {
+  return compute_collocation_points_and_weights<Basis::Legendre,
+                                                Quadrature::Gauss>(num_points);
+}
+
+template <>
+std::pair<DataVector, DataVector> compute_collocation_points_and_weights<
+    Basis::SphericalHarmonic, Quadrature::Equiangular>(
+    const size_t num_points) {
+  return compute_collocation_points_and_weights<Basis::FiniteDifference,
+                                                Quadrature::CellCentered>(
+      num_points);
+}
+
+template <>
+DataVector compute_basis_function_value<Basis::SphericalHarmonic>(
+    const size_t /*k*/, const DataVector& /*x*/) {
+  ERROR("not implemented");
+}
+
+template <>
+DataVector compute_inverse_weight_function_values<Basis::SphericalHarmonic>(
+    const DataVector& /*x*/) {
+  ERROR("not implemented");
+}
+
+template <>
+double compute_basis_function_normalization_square<Basis::SphericalHarmonic>(
+    const size_t /*k*/) {
+  ERROR("not implemented");
+}
+
+template <Basis BasisType>
+Matrix spectral_indefinite_integral_matrix(size_t num_points);
+
+template <>
+Matrix spectral_indefinite_integral_matrix<Basis::SphericalHarmonic>(
+    const size_t /*num_points*/) {
+  ERROR("not implemented");
+}
+
+}  // namespace Spectral