DROP. Debug test failure

src/Domain/CoordinateMaps/TimeDependent/Shape.cpp

@@ -55,10 +55,8 @@ void Shape::jacobian_helper(
   tnsr::i<DataVector, 3, Frame::Inertial> cartesian_gradient(
       extended_ylm_.physical_size());
 
-  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();
+  std::array<DataVector, 2> collocation_theta_phis =
+      extended_ylm_.theta_phi_points();
 
   const auto& col_thetas = collocation_theta_phis[0];
   const auto& col_phis = collocation_theta_phis[1];
@@ -77,10 +75,10 @@ void Shape::jacobian_helper(
 
   get<2>(cartesian_gradient) = -sin(col_thetas) * get<0>(angular_gradient);
 
-  // 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);
+  // 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);
 
   // interpolate the cartesian gradient to the thetas and phis of the
   // `source_coords`
@@ -108,11 +106,16 @@ 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,
-               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;
+  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];
 
   get<0, 0>(*result) =
       -x_centered * ((x_transition_gradient_over_radius -
@@ -323,7 +326,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);
 
@@ -369,15 +372,15 @@ void Shape::coords_frame_velocity_jacobian(
   center_coordinates(make_not_null(&centered_coords),
                      *source_and_target_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);
+  std::array<DataVector, 2> theta_phis =
+      cartesian_to_spherical(centered_coords);
   const auto interpolation_info =
       extended_ylm_.set_up_interpolation_info(theta_phis);
 
-  const auto [coefs, coef_derivs] =
+  const auto func_and_deriv =
       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.);
 
@@ -413,7 +416,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,7 +158,8 @@ 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
@@ -264,6 +265,12 @@ 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,7 +88,25 @@ template <typename T>
 T SphereTransition::call_impl(const std::array<T, 3>& source_coords) const {
   T mag = magnitude(source_coords);
   check_magnitudes(mag, false);
-  return blaze::clamp(a_ * mag + b_, 0.0, 1.0);
+  // 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_;
 }
 
 template <typename T>

src/Domain/Creators/Sphere.hpp

@@ -391,6 +391,10 @@ 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,14 +161,18 @@ 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;
-    }
+      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};
     // 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,8 +6,16 @@
 #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"
@@ -22,6 +30,7 @@
 #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"
@@ -55,6 +64,60 @@
 #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;
@@ -96,21 +159,24 @@ 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>,
-      ::Events::Tags::ObserverDetInvJacobianCompute<Frame::ElementLogical,
-                                                    Frame::Inertial>,
-      spacetime_quantities_compute, hydro_quantities_compute, error_compute>;
+  using observer_compute_tags =
+      tmpl::list<::Events::Tags::ObserverMeshCompute<volume_dim>,
+                 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>;
+      tmpl::list<domain::Tags::RadiallyCompressedCoordinatesOptions,
+                 Tags::InnerOuterRadius>;
 
   struct factory_creation
       : tt::ConformsTo<Options::protocols::FactoryCreation> {