Add function to write strahlkorper coords to a text file

src/NumericalAlgorithms/SphericalHarmonics/IO/CMakeLists.txt

@@ -10,6 +10,7 @@ spectre_target_sources(
   PRIVATE
   FillYlmLegendAndData.cpp
   ReadSurfaceYlm.cpp
+  StrahlkorperCoordsToTextFile.cpp
   )
 
 spectre_target_headers(
@@ -18,6 +19,7 @@ spectre_target_headers(
   HEADERS
   FillYlmLegendAndData.hpp
   ReadSurfaceYlm.hpp
+  StrahlkorperCoordsToTextFile.hpp
   )
 
 target_link_libraries(

src/NumericalAlgorithms/SphericalHarmonics/IO/StrahlkorperCoordsToTextFile.cpp

@@ -0,0 +1,94 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#include "NumericalAlgorithms/SphericalHarmonics/IO/StrahlkorperCoordsToTextFile.hpp"
+
+#include <fstream>
+#include <iomanip>
+#include <limits>
+#include <ostream>
+#include <string>
+
+#include "DataStructures/DataVector.hpp"
+#include "DataStructures/Tensor/Tensor.hpp"
+#include "DataStructures/Transpose.hpp"
+#include "NumericalAlgorithms/SphericalHarmonics/AngularOrdering.hpp"
+#include "NumericalAlgorithms/SphericalHarmonics/Strahlkorper.hpp"
+#include "NumericalAlgorithms/SphericalHarmonics/StrahlkorperFunctions.hpp"
+#include "Utilities/ErrorHandling/Error.hpp"
+#include "Utilities/FileSystem.hpp"
+#include "Utilities/GenerateInstantiations.hpp"
+
+namespace Frame {
+struct Inertial;
+struct Distorted;
+struct Grid;
+}  // namespace Frame
+
+namespace ylm {
+template <typename Frame>
+void write_strahlkorper_coords_to_text_file(
+    const Strahlkorper<Frame>& strahlkorper,
+    const std::string& output_file_name, const AngularOrdering ordering,
+    const bool overwrite) {
+  if (not overwrite and file_system::check_if_file_exists(output_file_name)) {
+    ERROR_NO_TRACE("The output file " << output_file_name
+                                      << " already exists.");
+  }
+
+  tnsr::I<DataVector, 3, Frame> cartesian_coords =
+      ylm::cartesian_coords(strahlkorper);
+
+  // Cce expects coordinates in a different order than a typical Strahlkorper
+  if (ordering == AngularOrdering::Cce) {
+    const auto physical_extents =
+        strahlkorper.ylm_spherepack().physical_extents();
+    auto transposed_coords =
+        tnsr::I<DataVector, 3, Frame>(get<0>(cartesian_coords).size());
+    for (size_t i = 0; i < 3; ++i) {
+      transpose(make_not_null(&transposed_coords.get(i)),
+                cartesian_coords.get(i), physical_extents[0],
+                physical_extents[1]);
+    }
+
+    cartesian_coords = std::move(transposed_coords);
+  }
+
+  std::ofstream output_file(output_file_name);
+  output_file << std::fixed
+              << std::setprecision(std::numeric_limits<double>::digits10 + 4)
+              << std::scientific;
+
+  const size_t num_points = get<0>(cartesian_coords).size();
+  for (size_t i = 0; i < num_points; i++) {
+    output_file << get<0>(cartesian_coords)[i] << " "
+                << get<1>(cartesian_coords)[i] << " "
+                << get<2>(cartesian_coords)[i] << std::endl;
+  }
+}
+
+void write_strahlkorper_coords_to_text_file(const double radius,
+                                            const size_t l_max,
+                                            const std::array<double, 3>& center,
+                                            const std::string& output_file_name,
+                                            const AngularOrdering ordering,
+                                            const bool overwrite) {
+  const Strahlkorper<Frame::Inertial> strahlkorper{l_max, radius, center};
+  write_strahlkorper_coords_to_text_file(strahlkorper, output_file_name,
+                                         ordering, overwrite);
+}
+
+#define FRAME(data) BOOST_PP_TUPLE_ELEM(0, data)
+
+#define INSTANTIATE(_, data)                                               \
+  template void write_strahlkorper_coords_to_text_file(                    \
+      const Strahlkorper<FRAME(data)>& strahlkorper,                       \
+      const std::string& output_file_name, const AngularOrdering ordering, \
+      const bool overwrite);
+
+GENERATE_INSTANTIATIONS(INSTANTIATE,
+                        (Frame::Grid, Frame::Distorted, Frame::Inertial))
+
+#undef INSTANTIATE
+#undef FRAME
+}  // namespace ylm

src/NumericalAlgorithms/SphericalHarmonics/IO/StrahlkorperCoordsToTextFile.hpp

@@ -0,0 +1,39 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#pragma once
+
+#include <array>
+#include <cstddef>
+#include <string>
+
+#include "NumericalAlgorithms/SphericalHarmonics/AngularOrdering.hpp"
+#include "NumericalAlgorithms/SphericalHarmonics/Strahlkorper.hpp"
+
+namespace ylm {
+/// @{
+/*!
+ * \brief Writes the collocation points of a `ylm::Strahlkorper` to an output
+ * text file.
+ *
+ * \details The ordering of the points can be either the typical
+ * `ylm::Spherepack` ordering or CCE ordering that works with libsharp. Also, an
+ * error will occur if the output file already exists, but the output file can
+ * be overwritten with the \p overwrite option.
+ *
+ * The second overload will construct a spherical `ylm::Strahlkorper` with the
+ * given \p radius, \p l_max, and \p center.
+ */
+template <typename Frame>
+void write_strahlkorper_coords_to_text_file(
+    const Strahlkorper<Frame>& strahlkorper,
+    const std::string& output_file_name, AngularOrdering ordering,
+    bool overwrite = false);
+
+void write_strahlkorper_coords_to_text_file(double radius, size_t l_max,
+                                            const std::array<double, 3>& center,
+                                            const std::string& output_file_name,
+                                            AngularOrdering ordering,
+                                            bool overwrite = false);
+/// @}
+}  // namespace ylm

tests/Unit/NumericalAlgorithms/SphericalHarmonics/IO/CMakeLists.txt

@@ -6,6 +6,7 @@ set(LIBRARY "Test_SphericalHarmonicsIO")
 set(LIBRARY_SOURCES
   Test_FillYlmLegendAndData.cpp
   Test_ReadSurfaceYlm.cpp
+  Test_StrahlkorperCoordsToTextFile.cpp
   )
 
 add_test_library(${LIBRARY} "${LIBRARY_SOURCES}")

tests/Unit/NumericalAlgorithms/SphericalHarmonics/IO/Test_StrahlkorperCoordsToTextFile.cpp

@@ -0,0 +1,124 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#include "Framework/TestingFramework.hpp"
+
+#include <array>
+#include <fstream>
+#include <istream>
+#include <sstream>
+#include <string>
+
+#include "DataStructures/DataVector.hpp"
+#include "DataStructures/Tensor/Tensor.hpp"
+#include "DataStructures/Transpose.hpp"
+#include "NumericalAlgorithms/SphericalHarmonics/AngularOrdering.hpp"
+#include "NumericalAlgorithms/SphericalHarmonics/IO/StrahlkorperCoordsToTextFile.hpp"
+#include "NumericalAlgorithms/SphericalHarmonics/Strahlkorper.hpp"
+#include "NumericalAlgorithms/SphericalHarmonics/StrahlkorperFunctions.hpp"
+#include "Utilities/FileSystem.hpp"
+
+namespace {
+std::array<DataVector, 3> read_text_file(const std::string& filename) {
+  std::array<std::vector<double>, 3> vector_result{};
+  std::ifstream file(filename);
+  if (not file.is_open()) {
+    ERROR("Unable to open text file " << filename);
+  }
+  std::string line{};
+  double value = 0.0;
+
+  while (std::getline(file, line)) {
+    std::stringstream ss(line);
+    for (size_t i = 0; i < 3; i++) {
+      ss >> value;
+      gsl::at(vector_result, i).push_back(value);
+    }
+  }
+
+  const size_t num_points = vector_result[0].size();
+
+  std::array<DataVector, 3> result{
+      DataVector{num_points}, DataVector{num_points}, DataVector{num_points}};
+
+  for (size_t i = 0; i < 3; i++) {
+    for (size_t j = 0; j < num_points; j++) {
+      gsl::at(result, i)[j] = gsl::at(vector_result, i)[j];
+    }
+  }
+
+  return result;
+}
+
+void test(const ylm::AngularOrdering ordering) {
+  const std::string filename{"StrahlkorperCoords.txt"};
+  const double radius = 1.5;
+  const size_t l_max = 16;
+  const std::array<double, 3> center{-0.1, -0.2, -0.3};
+  const ylm::Strahlkorper<Frame::Inertial> strahlkorper{l_max, radius, center};
+
+  tnsr::I<DataVector, 3, Frame::Inertial> expected_points =
+      ylm::cartesian_coords(strahlkorper);
+  if (ordering == ylm::AngularOrdering::Cce) {
+    const auto physical_extents =
+        strahlkorper.ylm_spherepack().physical_extents();
+    auto transpose_expected_points =
+        tnsr::I<DataVector, 3, Frame::Inertial>(get<0>(expected_points).size());
+    for (size_t i = 0; i < 3; ++i) {
+      transpose(make_not_null(&transpose_expected_points.get(i)),
+                expected_points.get(i), physical_extents[0],
+                physical_extents[1]);
+    }
+
+    expected_points = std::move(transpose_expected_points);
+  }
+
+  if (file_system::check_if_file_exists(filename)) {
+    file_system::rm(filename, true);
+  }
+
+  {
+    ylm::write_strahlkorper_coords_to_text_file(strahlkorper, filename,
+                                                ordering);
+
+    std::array<DataVector, 3> points_from_file = read_text_file(filename);
+
+    for (size_t i = 0; i < 3; i++) {
+      CHECK(expected_points.get(i) == gsl::at(points_from_file, i));
+    }
+
+    CHECK_THROWS_WITH((ylm::write_strahlkorper_coords_to_text_file(
+                          strahlkorper, filename, ordering)),
+                      Catch::Matchers::ContainsSubstring(
+                          "The output file " + filename + " already exists"));
+
+    ylm::write_strahlkorper_coords_to_text_file(strahlkorper, filename,
+                                                ordering, true);
+
+    for (size_t i = 0; i < 3; i++) {
+      CHECK(expected_points.get(i) == gsl::at(points_from_file, i));
+    }
+  }
+
+  {
+    ylm::write_strahlkorper_coords_to_text_file(radius, l_max, center, filename,
+                                                ordering, true);
+
+    const std::array<DataVector, 3> points_from_file = read_text_file(filename);
+
+    for (size_t i = 0; i < 3; i++) {
+      CHECK(expected_points.get(i) == gsl::at(points_from_file, i));
+    }
+  }
+
+  if (file_system::check_if_file_exists(filename)) {
+    file_system::rm(filename, true);
+  }
+}
+}  // namespace
+
+SPECTRE_TEST_CASE("Unit.SphericalHarmonics.StrahlkorperCoordsToTextFile",
+                  "[NumericalAlgorithms][Unit]") {
+  test(ylm::AngularOrdering::Strahlkorper);
+  test(ylm::AngularOrdering::Cce);
+}