Revert "DROP. Debug test failure"

This reverts commit d8dadb5.

src/Domain/CoordinateMaps/TimeDependent/Shape.cpp

@@ -55,8 +55,10 @@ void Shape::jacobian_helper(
   tnsr::i<DataVector, 3, Frame::Inertial> cartesian_gradient(
       extended_ylm_.physical_size());
 
-  std::array<DataVector, 2> collocation_theta_phis =
-      extended_ylm_.theta_phi_points();
+  std::array<DataVector, 2> collocation_theta_phis{};
+  collocation_theta_phis[0].set_data_ref(&get<2>(cartesian_gradient));
+  collocation_theta_phis[1].set_data_ref(&get<1>(cartesian_gradient));
+  collocation_theta_phis = extended_ylm_.theta_phi_points();
 
   const auto& col_thetas = collocation_theta_phis[0];
   const auto& col_phis = collocation_theta_phis[1];
@@ -75,10 +77,10 @@ void Shape::jacobian_helper(
 
   get<2>(cartesian_gradient) = -sin(col_thetas) * get<0>(angular_gradient);
 
-  // don't re-use allocation
-  auto target_gradient_x = get<2, 0>(*result);
-  auto target_gradient_y = get<2, 1>(*result);
-  auto target_gradient_z = get<2, 2>(*result);
+  // re-use allocation
+  auto& target_gradient_x = get<2, 0>(*result);
+  auto& target_gradient_y = get<2, 1>(*result);
+  auto& target_gradient_z = get<2, 2>(*result);
 
   // interpolate the cartesian gradient to the thetas and phis of the
   // `source_coords`
@@ -106,16 +108,11 @@ void Shape::jacobian_helper(
   target_gradient_y_times_spatial_part *= transition_func_over_square_radius;
   target_gradient_z_times_spatial_part *= transition_func_over_square_radius;
 
-  const auto& x_transition_gradient_over_radius =
-      transition_func_gradient_over_radius[0];
-  const auto& y_transition_gradient_over_radius =
-      transition_func_gradient_over_radius[1];
-  const auto& z_transition_gradient_over_radius =
-      transition_func_gradient_over_radius[2];
-  ;
-  const auto& x_centered = centered_coords[0];
-  const auto& y_centered = centered_coords[1];
-  const auto& z_centered = centered_coords[2];
+  const auto& [x_transition_gradient_over_radius,
+               y_transition_gradient_over_radius,
+               z_transition_gradient_over_radius] =
+      transition_func_gradient_over_radius;
+  const auto& [x_centered, y_centered, z_centered] = centered_coords;
 
   get<0, 0>(*result) =
       -x_centered * ((x_transition_gradient_over_radius -
@@ -326,7 +323,7 @@ tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame> Shape::jacobian(
   check_size(make_not_null(&extended_coefs), functions_of_time, time, false);
 
   // Re-use allocation
-  auto radial_distortion = get<0>(theta_phis);
+  auto& radial_distortion = get<0>(theta_phis);
   extended_ylm_.interpolate_from_coefs(make_not_null(&radial_distortion),
                                        extended_coefs, interpolation_info);
 
@@ -372,15 +369,15 @@ void Shape::coords_frame_velocity_jacobian(
   center_coordinates(make_not_null(&centered_coords),
                      *source_and_target_coords);
 
-  std::array<DataVector, 2> theta_phis =
-      cartesian_to_spherical(centered_coords);
+  std::array<DataVector, 2> theta_phis{};
+  theta_phis[0].set_data_ref(&get<0, 0>(*jac));
+  theta_phis[1].set_data_ref(&get<0, 1>(*jac));
+  cartesian_to_spherical(make_not_null(&theta_phis), centered_coords);
   const auto interpolation_info =
       extended_ylm_.set_up_interpolation_info(theta_phis);
 
-  const auto func_and_deriv =
+  const auto [coefs, coef_derivs] =
       functions_of_time.at(shape_f_of_t_name_)->func_and_deriv(time);
-  const auto& coefs = func_and_deriv[0];
-  const auto& coef_derivs = func_and_deriv[1];
   DataVector extended_coefs_derivs(extended_ylm_.spectral_size(), 0.);
   DataVector extended_coefs(extended_ylm_.spectral_size(), 0.);
 
@@ -416,7 +413,7 @@ void Shape::coords_frame_velocity_jacobian(
       center_ +
       centered_coords * (1. - radial_distortion * transition_func_over_radius);
 
-  auto radii_velocities = get<0, 1>(*jac);
+  auto& radii_velocities = get<0, 1>(*jac);
   extended_ylm_.interpolate_from_coefs(make_not_null(&radii_velocities),
                                        extended_coefs_derivs,
                                        interpolation_info);

src/Domain/CoordinateMaps/TimeDependent/Shape.hpp

@@ -158,8 +158,7 @@ namespace domain::CoordinateMaps::TimeDependent {
  * \f}
  *
  * For more details, see
- * \link
- * domain::CoordinateMaps::ShapeMapTransitionFunctions::ShapeMapTransitionFunction::original_radius_over_radius
+ * \link domain::CoordinateMaps::ShapeMapTransitionFunctions::ShapeMapTransitionFunction::original_radius_over_radius
  * ShapeMapTransitionFunction::original_radius_over_radius \endlink.
  *
  * ### Frame velocity
@@ -265,12 +264,6 @@ class Shape {
     return f_of_t_names_;
   }
 
-  const std::unique_ptr<
-      ShapeMapTransitionFunctions::ShapeMapTransitionFunction>&
-  transition_function() {
-    return transition_func_;
-  }
-
  private:
   std::string shape_f_of_t_name_;
   std::optional<std::string> size_f_of_t_name_;

src/Domain/CoordinateMaps/TimeDependent/ShapeMapTransitionFunctions/SphereTransition.cpp

@@ -88,25 +88,7 @@ template <typename T>
 T SphereTransition::call_impl(const std::array<T, 3>& source_coords) const {
   T mag = magnitude(source_coords);
   check_magnitudes(mag, false);
-  // T tmp = a_ * mag + b_;
-  // for (size_t i=0; i<get_size(mag); i++) {
-  //   if (get_element(mag, i) /*+ eps_*/ < r_min_) {
-  //     get_element(tmp, i) = a_ > 0.0 ? 0.0 : 1.0;
-  //   } else if (get_element(mag, i) /*- eps_*/ > r_max_) {
-  //     get_element(tmp, i) = a_ > 0.0 ? 1.0 : 0.0;
-  //   }
-  // }
-  // for (size_t i=0; i<get_size(mag); i++) {
-  //   if (get_element(tmp, i) /*+ eps_*/ > 1.0) {
-  //     get_element(tmp, i) = 1.0;
-  //   } else if (get_element(mag, i) /*- eps_*/ < 0.0) {
-  //     get_element(tmp, i) = 0.0;
-  //   }
-  // }
-  // return tmp;
-  return blaze::clamp(T(a_ * mag + b_), 0.0, 1.0);
-  // return blaze::clamp(a_ * mag + b_, 0.0, 1.0);
-  // return a_ * mag + b_;
+  return blaze::clamp(a_ * mag + b_, 0.0, 1.0);
 }
 
 template <typename T>

src/Domain/Creators/Sphere.hpp

@@ -391,10 +391,6 @@ class Sphere : public DomainCreator<3> {
     return grid_anchors_;
   }
 
-  std::array<double, 2> inner_outer_radius() const {
-    return {inner_radius_, outer_radius_};
-  }
-
   std::vector<DirectionMap<
       3, std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>>>
   external_boundary_conditions() const override;

src/Elliptic/DiscontinuousGalerkin/Initialization.hpp

@@ -161,18 +161,14 @@ struct ProjectGeometry : tt::ConformsTo<::amr::protocols::Projector> {
       const gsl::not_null<InverseJacobian<
           DataVector, Dim, Frame::ElementLogical, Frame::Inertial>*>
           det_times_inv_jacobian,
-      const Mesh<Dim>& mesh, const Element<Dim>& element,
-      const Domain<Dim>& domain,
+      const Mesh<Dim>& mesh, const Element<Dim>& /*element*/,
+      const Domain<Dim>& /*domain*/,
       const domain::FunctionsOfTimeMap& functions_of_time,
-      const std::pair<Mesh<Dim>, Element<Dim>>& /*old_mesh_and_element*/) {
-    // if (mesh == old_mesh_and_element.first) {
-    //   // Neighbors may have changed, but this element hasn't. Nothing to do.
-    //   return;
-    // }
-    // QUESTION: Just try it out
-    const auto& element_id = element.id();
-    const auto& block = domain.blocks()[element_id.block_id()];
-    *element_map = ElementMap<Dim, Frame::Inertial>{element_id, block};
+      const std::pair<Mesh<Dim>, Element<Dim>>& old_mesh_and_element) {
+    if (mesh == old_mesh_and_element.first) {
+      // Neighbors may have changed, but this element hasn't. Nothing to do.
+      return;
+    }
     // The element map doesn't change with p-refinement
     detail::initialize_coords_and_jacobians(
         logical_coords, inertial_coords, inv_jacobian, det_inv_jacobian,

src/Elliptic/Executables/Xcts/SolveXcts.hpp

@@ -6,16 +6,8 @@
 #include <cstddef>
 
 #include "DataStructures/DataBox/PrefixHelpers.hpp"
-#include "DataStructures/DataBox/Tag.hpp"
-#include "DataStructures/DataVector.hpp"
 #include "DataStructures/Tensor/EagerMath/Magnitude.hpp"
-#include "DataStructures/Tensor/IndexType.hpp"
-#include "Domain/CoordinateMaps/TimeDependent/ShapeMapTransitionFunctions/SphereTransition.hpp"
-#include "Domain/Creators/DomainCreator.hpp"
 #include "Domain/Creators/Factory3D.hpp"
-#include "Domain/Creators/Sphere.hpp"
-#include "Domain/Creators/Tags/Domain.hpp"
-#include "Domain/ElementMap.hpp"
 #include "Domain/RadiallyCompressedCoordinates.hpp"
 #include "Domain/Tags.hpp"
 #include "Elliptic/Actions/RunEventsAndTriggers.hpp"
@@ -30,7 +22,6 @@
 #include "IO/Observer/Actions/RegisterEvents.hpp"
 #include "IO/Observer/Helpers.hpp"
 #include "IO/Observer/ObserverComponent.hpp"
-#include "NumericalAlgorithms/Spectral/LogicalCoordinates.hpp"
 #include "Options/Protocols/FactoryCreation.hpp"
 #include "Options/String.hpp"
 #include "Parallel/GlobalCache.hpp"
@@ -64,60 +55,6 @@
 #include "Utilities/ProtocolHelpers.hpp"
 #include "Utilities/TMPL.hpp"
 
-namespace Tags {
-struct InnerOuterRadius : db::SimpleTag {
-  using type = std::array<double, 2>;
-  static constexpr bool pass_metavariables = false;
-  using option_tags = tmpl::list<domain::OptionTags::DomainCreator<3>>;
-  static type create_from_options(
-      const std::unique_ptr<DomainCreator<3>>& creator) {
-    const auto* sphere =
-        dynamic_cast<const domain::creators::Sphere*>(creator.get());
-    if (sphere == nullptr) {
-      ERROR_NO_TRACE("Why have you forsaken me");
-    }
-
-    return sphere->inner_outer_radius();
-  }
-};
-
-struct TransitionFunction : db::SimpleTag {
-  using type = Scalar<DataVector>;
-};
-
-struct TransitionFunctionCompute : TransitionFunction, db::ComputeTag {
-  using base = TransitionFunction;
-  using return_type = typename base::type;
-  using argument_tags =
-      tmpl::list<Tags::InnerOuterRadius, domain::Tags::Domain<3>,
-                 domain::Tags::Element<3>,
-                 domain::Tags::Coordinates<3, Frame::ElementLogical>>;
-
-  static void function(const gsl::not_null<type*> transition_function,
-                       const std::array<double, 2> inner_outer_radius,
-                       const Domain<3>& domain, const Element<3>& element,
-                       const tnsr::I<DataVector, 3, Frame::ElementLogical>&
-                           element_logical_coords) {
-    const domain::CoordinateMaps::ShapeMapTransitionFunctions::SphereTransition
-        sphere_transition{inner_outer_radius[0], inner_outer_radius[1]};
-
-    const ElementId<3>& element_id = element.id();
-    const Block<3>& block = domain.blocks()[element_id.block_id()];
-
-    ElementMap<3, Frame::Grid> element_map(element_id, block);
-
-    auto grid_coords = element_map(element_logical_coords);
-
-    std::array<DataVector, 3> array_grid_coords{};
-    for (size_t i = 0; i < array_grid_coords.size(); i++) {
-      gsl::at(array_grid_coords, i) = grid_coords.get(i);
-    }
-
-    get(*transition_function) = sphere_transition(array_grid_coords);
-  }
-};
-}  // namespace Tags
-
 /// \cond
 struct Metavariables {
   static constexpr size_t volume_dim = 3;
@@ -159,24 +96,21 @@ struct Metavariables {
       typename hydro_quantities_compute::tags_list, error_tags,
       typename solver::observe_fields,
       tmpl::list<domain::Tags::Coordinates<volume_dim, Frame::Inertial>,
-                 ::Events::Tags::ObserverDetInvJacobianCompute<
-                     Frame::ElementLogical, Frame::Inertial>,
-                 Tags::TransitionFunctionCompute,
                  domain::Tags::RadiallyCompressedCoordinatesCompute<
                      volume_dim, Frame::Inertial>,
                  ::Tags::NonEuclideanMagnitude<
                      Xcts::Tags::ShiftExcess<DataVector, 3, Frame::Inertial>,
                      gr::Tags::SpatialMetric<DataVector, 3>>,
                  hydro::Tags::LowerSpatialFourVelocityCompute>>;
-  using observer_compute_tags =
-      tmpl::list<::Events::Tags::ObserverMeshCompute<volume_dim>,
-                 spacetime_quantities_compute, hydro_quantities_compute,
-                 error_compute>;
+  using observer_compute_tags = tmpl::list<
+      ::Events::Tags::ObserverMeshCompute<volume_dim>,
+      ::Events::Tags::ObserverDetInvJacobianCompute<Frame::ElementLogical,
+                                                    Frame::Inertial>,
+      spacetime_quantities_compute, hydro_quantities_compute, error_compute>;
 
   // Collect all items to store in the cache.
   using const_global_cache_tags =
-      tmpl::list<domain::Tags::RadiallyCompressedCoordinatesOptions,
-                 Tags::InnerOuterRadius>;
+      tmpl::list<domain::Tags::RadiallyCompressedCoordinatesOptions>;
 
   struct factory_creation
       : tt::ConformsTo<Options::protocols::FactoryCreation> {