StrahlkorperCoords src can be Frame::Distorted.

Formerly only Frame::Grid was allowed.

src/ApparentHorizons/StrahlkorperCoordsInDifferentFrame.cpp

@@ -8,13 +8,15 @@
 #include "DataStructures/Tensor/Tensor.hpp"
 #include "DataStructures/Variables.hpp"
 #include "Domain/Block.hpp"
+#include "Domain/BlockLogicalCoordinates.hpp"
 #include "Domain/CoordinateMaps/CoordinateMap.hpp"
 #include "Domain/Domain.hpp"
 #include "Domain/FunctionsOfTime/FunctionOfTime.hpp"
 #include "NumericalAlgorithms/SphericalHarmonics/Strahlkorper.hpp"
 #include "NumericalAlgorithms/SphericalHarmonics/Tags.hpp"
 #include "NumericalAlgorithms/SphericalHarmonics/YlmSpherepack.hpp"
 #include "Utilities/ContainerHelpers.hpp"
+#include "Utilities/ErrorHandling/Assert.hpp"
 #include "Utilities/GenerateInstantiations.hpp"
 #include "Utilities/Gsl.hpp"
 
@@ -27,6 +29,8 @@ void strahlkorper_coords_in_different_frame(
         std::string, std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>&
         functions_of_time,
     const double time) {
+  static_assert(std::is_same_v<DestFrame, ::Frame::Inertial>,
+                "Destination frame must currently be Inertial frame");
   destructive_resize_components(
       dest_cartesian_coords, src_strahlkorper.ylm_spherepack().physical_size());
 
@@ -53,41 +57,34 @@ void strahlkorper_coords_in_different_frame(
 
   // We now wish to map src_cartesian_coords to the destination frame.
   // Each Block will have a different map, so the mapping must be done
-  // Block by Block.
-  for (const auto& block : domain.blocks()) {
-    // Once there are more possible source frames than the grid frame
-    // (i.e. the distorted frame), then this static_assert will change,
-    // and there will be an `if constexpr` below to treat the different
-    // possible source frames.
-    static_assert(std::is_same_v<SrcFrame, ::Frame::Grid>,
-                  "Source frame must currently be Grid frame");
-    static_assert(std::is_same_v<DestFrame, ::Frame::Inertial>,
-                  "Destination frame must currently be Inertial frame");
-    const auto& grid_to_inertial_map = block.moving_mesh_grid_to_inertial_map();
-    const auto& logical_to_grid_map = block.moving_mesh_logical_to_grid_map();
-    // Fill only those dest_cartesian_coords that are in this Block.
-    // Determine which coords are in this Block by checking if the
-    // inverse grid-to-logical map yields logical coords between -1 and 1.
-    for (size_t s = 0; s < get<0>(src_cartesian_coords).size(); ++s) {
-      const tnsr::I<double, 3, SrcFrame> x_src{
-          {get<0>(src_cartesian_coords)[s], get<1>(src_cartesian_coords)[s],
-           get<2>(src_cartesian_coords)[s]}};
-      const auto x_logical = logical_to_grid_map.inverse(x_src);
-      // x_logical might be an empty std::optional.
-      if (x_logical.has_value() and get<0>(x_logical.value()) <= 1.0 and
-          get<0>(x_logical.value()) >= -1.0 and
-          get<1>(x_logical.value()) <= 1.0 and
-          get<1>(x_logical.value()) >= -1.0 and
-          get<2>(x_logical.value()) <= 1.0 and
-          get<2>(x_logical.value()) >= -1.0) {
-        const auto x_dest =
-            grid_to_inertial_map(x_src, time, functions_of_time);
-        get<0>(*dest_cartesian_coords)[s] = get<0>(x_dest);
-        get<1>(*dest_cartesian_coords)[s] = get<1>(x_dest);
-        get<2>(*dest_cartesian_coords)[s] = get<2>(x_dest);
-        // Note that if a point is on the boundary of two or more
-        // Blocks, it might get filled twice, but that's ok.
-      }
+  // Block by Block.  Since each point in general has a different
+  // block (given by block_logical_coordinates), this means we must
+  // map point by point.
+  const auto block_logical_coords = block_logical_coordinates(
+      domain, src_cartesian_coords, time, functions_of_time);
+  for (size_t s = 0; s < block_logical_coords.size(); ++s) {
+    const tnsr::I<double, 3, SrcFrame> x_src{{get<0>(src_cartesian_coords)[s],
+                                              get<1>(src_cartesian_coords)[s],
+                                              get<2>(src_cartesian_coords)[s]}};
+    ASSERT(block_logical_coords[s].has_value(),
+           "Found a point (source coords " << x_src
+                                           << ") that is not in any Block.");
+    const auto& block =
+        domain.blocks()[block_logical_coords[s].value().id.get_index()];
+    if constexpr (std::is_same_v<SrcFrame, ::Frame::Grid>) {
+      const auto x_dest = block.moving_mesh_grid_to_inertial_map()(
+          x_src, time, functions_of_time);
+      get<0>(*dest_cartesian_coords)[s] = get<0>(x_dest);
+      get<1>(*dest_cartesian_coords)[s] = get<1>(x_dest);
+      get<2>(*dest_cartesian_coords)[s] = get<2>(x_dest);
+    } else {
+      static_assert(std::is_same_v<SrcFrame, ::Frame::Distorted>,
+                    "Src frame must be Distorted if it is not Grid");
+      const auto x_dest = block.moving_mesh_distorted_to_inertial_map()(
+          x_src, time, functions_of_time);
+      get<0>(*dest_cartesian_coords)[s] = get<0>(x_dest);
+      get<1>(*dest_cartesian_coords)[s] = get<1>(x_dest);
+      get<2>(*dest_cartesian_coords)[s] = get<2>(x_dest);
     }
   }
 }
@@ -107,7 +104,8 @@ void strahlkorper_coords_in_different_frame(
           functions_of_time,                                         \
       const double time);
 
-GENERATE_INSTANTIATIONS(INSTANTIATE, (::Frame::Grid), (::Frame::Inertial))
+GENERATE_INSTANTIATIONS(INSTANTIATE, (::Frame::Grid, ::Frame::Distorted),
+                        (::Frame::Inertial))
 
 #undef INSTANTIATE
 #undef DESTFRAME

tests/Unit/ApparentHorizons/Test_StrahlkorperCoordsInDifferentFrame.cpp

@@ -23,17 +23,18 @@
 
 namespace {
 
+template<typename SrcFrame>
 void test_strahlkorper_coords_in_different_frame() {
   const size_t grid_points_each_dimension = 5;
 
   // Set up a Strahlkorper corresponding to a Schwarzschild hole of
-  // mass 1, in the grid frame.
+  // mass 1, in the source frame.
   // Center the Strahlkorper at (0.03,0.02,0.01) so that we test a
   // nonzero center.
-  const std::array<double, 3> strahlkorper_grid_center = {0.03, 0.02, 0.01};
+  const std::array<double, 3> strahlkorper_src_center = {0.03, 0.02, 0.01};
   const size_t l_max = 8;
-  const Strahlkorper<Frame::Grid> strahlkorper_grid(l_max, 2.0,
-                                                    strahlkorper_grid_center);
+  const Strahlkorper<SrcFrame> strahlkorper_src(l_max, 2.0,
+                                                strahlkorper_src_center);
 
   // Create a Domain.
   // We choose a spherical shell domain extending from radius 1.9M to
@@ -43,15 +44,32 @@ void test_strahlkorper_coords_in_different_frame() {
   std::vector<double> radial_partitioning{};
   std::vector<domain::CoordinateMaps::Distribution> radial_distribution{
       domain::CoordinateMaps::Distribution::Linear};
+
+  // In computing the time_dependence, make sure that src-to-inertial
+  // velocity is (0.01,0.02,0.03) to agree with the analytic checks
+  // below.  If src is distorted frame, then grid-to-distorted
+  // velocity doesn't matter for the value of the check (but does matter
+  // in terms of which points are in which blocks).
+  std::unique_ptr<domain::creators::time_dependence::TimeDependence<3>>
+      time_dependence;
+  if constexpr (std::is_same_v<SrcFrame, ::Frame::Grid>) {
+    time_dependence = std::make_unique<
+        domain::creators::time_dependence::UniformTranslation<3>>(
+        0.0, std::array<double, 3>({{0.01, 0.02, 0.03}}));
+  } else {
+    static_assert(std::is_same_v<SrcFrame, ::Frame::Distorted>,
+                  "Src frame must be Distorted if it is not Grid");
+    time_dependence = std::make_unique<
+        domain::creators::time_dependence::UniformTranslation<3>>(
+        0.0, std::array<double, 3>({{-0.02, -0.01, -0.01}}),
+        std::array<double, 3>({{0.01, 0.02, 0.03}}));
+  }
   domain::creators::Shell domain_creator(
       1.9, 2.9, 1,
       std::array<size_t, 2>{grid_points_each_dimension,
                             grid_points_each_dimension},
       false, {{1.0, 2}}, radial_partitioning, radial_distribution,
-      ShellWedges::All,
-      std::make_unique<
-          domain::creators::time_dependence::UniformTranslation<3>>(
-          0.0, std::array<double, 3>({{0.01, 0.02, 0.03}})));
+      ShellWedges::All, std::move(time_dependence));
   Domain<3> domain = domain_creator.create_domain();
   const auto functions_of_time = domain_creator.functions_of_time();
 
@@ -60,26 +78,27 @@ void test_strahlkorper_coords_in_different_frame() {
   tnsr::I<DataVector, 3, Frame::Inertial> inertial_coords{};
 
   strahlkorper_coords_in_different_frame(make_not_null(&inertial_coords),
-                                         strahlkorper_grid, domain,
+                                         strahlkorper_src, domain,
                                          functions_of_time, time);
 
-  // Now compare with expected result, which is the grid-frame coords of
+  // Now compare with expected result, which is the src-frame coords of
   // the Strahlkorper translated by (0.005,0.01,0.015).
-  const auto grid_coords = StrahlkorperFunctions::cartesian_coords(
-      strahlkorper_grid, StrahlkorperFunctions::radius(strahlkorper_grid),
+  const auto src_coords = StrahlkorperFunctions::cartesian_coords(
+      strahlkorper_src, StrahlkorperFunctions::radius(strahlkorper_src),
       StrahlkorperFunctions::rhat(
-          StrahlkorperFunctions::theta_phi(strahlkorper_grid)));
-  CHECK_ITERABLE_APPROX(get<0>(grid_coords) + 0.005, get<0>(inertial_coords));
-  CHECK_ITERABLE_APPROX(get<1>(grid_coords) + 0.01, get<1>(inertial_coords));
-  CHECK_ITERABLE_APPROX(get<2>(grid_coords) + 0.015, get<2>(inertial_coords));
+          StrahlkorperFunctions::theta_phi(strahlkorper_src)));
+  CHECK_ITERABLE_APPROX(get<0>(src_coords) + 0.005, get<0>(inertial_coords));
+  CHECK_ITERABLE_APPROX(get<1>(src_coords) + 0.01, get<1>(inertial_coords));
+  CHECK_ITERABLE_APPROX(get<2>(src_coords) + 0.015, get<2>(inertial_coords));
 }
 
 SPECTRE_TEST_CASE("Unit.ApparentHorizons.StrahlkorperCoordsInDifferentFrame",
                   "[Unit]") {
   domain::creators::register_derived_with_charm();
   domain::creators::time_dependence::register_derived_with_charm();
   domain::FunctionsOfTime::register_derived_with_charm();
-  test_strahlkorper_coords_in_different_frame();
+  test_strahlkorper_coords_in_different_frame<Frame::Grid>();
+  test_strahlkorper_coords_in_different_frame<Frame::Distorted>();
 }
 
 }  // namespace