BOOLEAN scale_by_resolution "Scale error estimate by geometric mean of grid spacing"

KEYWORD region_shape "Shape of refined region" STEERABLE=always
{
  "sphere" :: ""
  "cube" :: ""
} "sphere"
BOOLEAN scale_by_resolution "Scale error estimate by geometric mean of grid spacing" STEERABLE=always

CCTK_REAL epsilon "Epsilon to avoid infinity when dividing by zero"

CCTK_REAL epsilon "Epsilon to avoid infinity when dividing by zero" STEERABLE=always

template <typename T> constexpr T pow2(const T x) { return x * x; }

namespace {
template <typename T> T pow2(T x) { return x * x; }

template <typename T> constexpr T fmax3(const T x, const T y, const T z) {
  return fmax(fmax(x, y), z);
}
enum class shape_t { sphere, cube };
template <typename T>
constexpr T radius_sphere(const T x, const T y, const T z) {
  return sqrt(pow2(x) + pow2(y) + pow2(z));
}
template <typename T> constexpr T radius_cube(const T x, const T y, const T z) {
  return fmax3(fabs(x), fabs(y), fabs(z));
}

} // namespace

  typedef CCTK_REAL radius_t(CCTK_REAL x, CCTK_REAL y, CCTK_REAL z);
  radius_t *const radius = [](const char *const sh) {
    if (CCTK_EQUALS(sh, "sphere"))
      return (radius_t *)radius_sphere;
    else if (CCTK_EQUALS(sh, "cube"))
      return (radius_t *)radius_cube;
    else
      assert(0);
  }(region_shape);

          const auto r = sqrt(pow2(p.x + i * dx / 2) + pow2(p.y + j * dy / 2) +
                              pow2(p.z + k * dz / 2));

          const auto r =
              radius(p.x + i * dx / 2, p.y + j * dy / 2, p.z + k * dz / 2);