Support surface-based VecGeom 2.x navigator

src/geocel/g4vg/SolidConverter.cc

@@ -129,6 +129,17 @@ make_unplaced_boolean(VPlacedVolume const* left, VPlacedVolume const* right)
     return GeoManager::MakeInstance<UnplacedBooleanVolume<Op>>(Op, left, right);
 }
 
+//---------------------------------------------------------------------------//
+VUnplacedVolume* make_sphere(double radius)
+{
+    CELER_EXPECT(radius > 0);
+#ifdef VECGEOM_USE_SURF
+    return GeoManager::MakeInstance<UnplacedSphere>(
+        0, radius, 0, vecgeom::kPi * 2, 0, vecgeom::kPi);
+#endif
+    return GeoManager::MakeInstance<UnplacedOrb>(radius);
+}
+
 //---------------------------------------------------------------------------//
 }  // namespace
 
@@ -156,8 +167,7 @@ auto SolidConverter::operator()(arg_type solid_base) -> result_type
 auto SolidConverter::to_sphere(arg_type solid_base) const -> result_type
 {
     double vol = this->calc_capacity(solid_base);
-    double radius = std::cbrt(vol / (4.0 / 3.0 * constants::pi));
-    return GeoManager::MakeInstance<UnplacedOrb>(radius);
+    return make_sphere(std::cbrt(vol / (4.0 / 3.0 * constants::pi)));
 }
 
 //---------------------------------------------------------------------------//
@@ -177,17 +187,25 @@ auto SolidConverter::convert_impl(arg_type solid_base) -> result_type
         VGSC_TYPE_FUNC(Box              , box),
         VGSC_TYPE_FUNC(Cons             , cons),
         VGSC_TYPE_FUNC(CutTubs          , cuttubs),
+#ifndef VECGEOM_USE_SURF
         VGSC_TYPE_FUNC(Ellipsoid        , ellipsoid),
         VGSC_TYPE_FUNC(EllipticalCone   , ellipticalcone),
         VGSC_TYPE_FUNC(EllipticalTube   , ellipticaltube),
+#endif
         VGSC_TYPE_FUNC(ExtrudedSolid    , extrudedsolid),
+#ifndef VECGEOM_USE_SURF
         VGSC_TYPE_FUNC(GenericPolycone  , genericpolycone),
+#endif
         VGSC_TYPE_FUNC(GenericTrap      , generictrap),
+#ifndef VECGEOM_USE_SURF
         VGSC_TYPE_FUNC(Hype             , hype),
+#endif
         VGSC_TYPE_FUNC(IntersectionSolid, intersectionsolid),
         VGSC_TYPE_FUNC(Orb              , orb),
         VGSC_TYPE_FUNC(Para             , para),
+#ifndef VECGEOM_USE_SURF
         VGSC_TYPE_FUNC(Paraboloid       , paraboloid),
+#endif
         VGSC_TYPE_FUNC(Polycone         , polycone),
         VGSC_TYPE_FUNC(Polyhedra        , polyhedra),
         VGSC_TYPE_FUNC(ReflectedSolid   , reflectedsolid),
@@ -416,7 +434,12 @@ auto SolidConverter::intersectionsolid(arg_type solid_base) -> result_type
 auto SolidConverter::orb(arg_type solid_base) -> result_type
 {
     auto const& solid = dynamic_cast<G4Orb const&>(solid_base);
-    return GeoManager::MakeInstance<UnplacedOrb>(scale_(solid.GetRadius()));
+    auto radius = scale_(solid.GetRadius());
+#ifdef VECGEOM_USE_SURF
+    CELER_LOG(info) << "Replacing orb with equivalent sphere for surface "
+                       "geometry";
+#endif
+    return make_sphere(radius);
 }
 
 //---------------------------------------------------------------------------//

src/geocel/vg/VecgeomParams.cc

@@ -11,6 +11,7 @@
 #include <vector>
 #include <VecGeom/base/Config.h>
 #include <VecGeom/base/Cuda.h>
+#include <VecGeom/base/Version.h>
 #include <VecGeom/management/ABBoxManager.h>
 #include <VecGeom/management/BVHManager.h>
 #include <VecGeom/management/GeoManager.h>
@@ -410,7 +411,11 @@ void VecgeomParams::build_volume_tracking()
 
     {
         ScopedTimeAndRedirect time_and_output_("vecgeom::ABBoxManager");
+#if VECGEOM_VERSION < 0x020000
         vecgeom::ABBoxManager::Instance().InitABBoxesForCompleteGeometry();
+#else
+        vecgeom::ABBoxManager<real_type>::Instance().InitABBoxesForCompleteGeometry();
+#endif
     }
 
     // Init the bounding volume hierarchy structure
@@ -570,9 +575,13 @@ void VecgeomParams::build_metadata()
         VPlacedVolume const* pv = GeoManager::Instance().GetWorld();
 
         // Calculate bounding box
+#if VECGEOM_VERSION < 0x020000
+        auto bbox_mgr = ABBoxManager::Instance();
+#else
+        auto bbox_mgr = ABBoxManager<real_type>::Instance();
+#endif
         Vector3D<real_type> lower, upper;
-        ABBoxManager::Instance().ComputeABBox(pv, &lower, &upper);
-
+        bbox_mgr.ComputeABBox(pv, &lower, &upper);
         return BBox{detail::to_array(lower), detail::to_array(upper)};
     }();
 }

src/geocel/vg/VecgeomTrackView.hh

@@ -25,12 +25,10 @@
 
 #include "detail/VecgeomCompatibility.hh"
 
-#if VECGEOM_VERSION < 0x020000
-#    include "detail/BVHNavigator.hh"
-#elif defined(VECGEOM_USE_SURF)
-#    include "detail/SurfNavigator.hh"
+#if defined(VECGEOM_USE_SURF)
+#    include "geocel/vg/detail/SurfNavigator.hh"
 #else
-#    include <VecGeom/navigation/BVHNavigator.h>
+#    include "geocel/vg/detail/BVHNavigator.hh"
 #endif
 
 namespace celeritas
@@ -57,12 +55,10 @@ class VecgeomTrackView
     using Initializer_t = GeoTrackInitializer;
     using ParamsRef = NativeCRef<VecgeomParamsData>;
     using StateRef = NativeRef<VecgeomStateData>;
-#if VECGEOM_VERSION < 0x020000
-    using Navigator = celeritas::detail::BVHNavigator;
-#elif defined(VECGEOM_USE_SURF)
+#if defined(VECGEOM_USE_SURF)
     using Navigator = celeritas::detail::SurfNavigator;
 #else
-    using Navigator = vecgeom::BVHNavigator;
+    using Navigator = celeritas::detail::BVHNavigator;
 #endif
     //!@}
 
@@ -215,13 +211,17 @@ VecgeomTrackView::operator=(Initializer_t const& init)
 
     // Set up current state and locate daughter volume.
     vgstate_.Clear();
-    vecgeom::VPlacedVolume const* worldvol = params_.world_volume;
     bool const contains_point = true;
-
+#ifdef VECGEOM_USE_SURF
+    auto world_id = vecgeom::NavigationState::WorldId();
     // LocatePointIn sets `vgstate_`
-    Navigator::LocatePointIn(
-        worldvol, detail::to_vector(pos_), vgstate_, contains_point);
-
+    Navigator::LocatePointIn(world_id, detail::to_vector(pos_), vgstate_, contains_point);
+#else
+    vecgeom::VPlacedVolume const* worldvol = params_.world_volume;
+    // LocatePointIn sets `vgstate_`
+    Navigator::LocatePointIn(worldvol,
+        detail::to_vector(pos_), vgstate_, contains_point);
+#endif
     return *this;
 }
 
@@ -421,10 +421,16 @@ CELER_FUNCTION void VecgeomTrackView::cross_boundary()
     // Relocate to next tracking volume (maybe across multiple boundaries)
     if (vgnext_.Top() != nullptr)
     {
+    #ifdef VECGEOM_USE_SURF
+        long hitsurf_id = -1;
+        Navigator::RelocateToNextVolume(detail::to_vector(this->pos_),
+            detail::to_vector(this->dir_), hitsurf_id, vgnext_);
+    #else
         // Some navigators require an lvalue temp_pos
         auto temp_pos = detail::to_vector(this->pos_);
         Navigator::RelocateToNextVolume(
             temp_pos, detail::to_vector(this->dir_), vgnext_);
+    #endif
     }
 
     vgstate_ = vgnext_;

src/geocel/vg/detail/SurfNavigator.hh

@@ -14,8 +14,8 @@
 #include <VecGeom/base/Config.h>
 #include <VecGeom/base/Global.h>
 #include <VecGeom/base/Vector3D.h>
-#include <VecGeom/navigation/NavStateIndex.h>
-#include <VecGeom/surfaces/Navigator.h>
+#include <VecGeom/navigation/NavigationState.h>
+#include <VecGeom/surfaces/BVHSurfNavigator.h>
 
 #ifdef VECGEOM_ENABLE_CUDA
 #    include <VecGeom/backend/cuda/Interface.h>
@@ -32,29 +32,39 @@ class SurfNavigator
   public:
     using Precision = vecgeom::Precision;
     using Vector3D = vecgeom::Vector3D<vecgeom::Precision>;
-    using VPlacedVolumePtr_t = vecgeom::VPlacedVolume const*;
     using SurfData = vgbrep::SurfData<Precision>;
+    using Real_b = typename SurfData::Real_b;
+    using VPlacedVolumePtr_t = vecgeom::VPlacedVolume const*;
 
     static constexpr Precision kBoundaryPush = 10 * vecgeom::kTolerance;
 
-    // Locates the point in the geometry volume tree
-    CELER_FUNCTION static VPlacedVolumePtr_t
-    LocatePointIn(VPlacedVolumePtr_t vol,
+    /// @brief Locates the point in the geometry volume tree
+    /// @param pvol_id Placed volume id to be checked first
+    /// @param point Point to be checked, in the local frame of pvol
+    /// @param path Path to a parent of pvol that must contain the point
+    /// @param top Check first if pvol contains the point
+    /// @param exclude Placed volume id to exclude from the search
+    /// @return Index of the placed volume that contains the point
+    CELER_FUNCTION static int
+    LocatePointIn(int pvol_id,
                   Vector3D const& point,
-                  vecgeom::NavStateIndex& path,
-                  bool top)
+                  vecgeom::NavigationState& path,
+                  bool top,
+		  int *exclude = nullptr)
     {
-        return vgbrep::protonav::LocatePointIn(vol, point, path, top);
+      return vgbrep::protonav::BVHSurfNavigator<Precision>::LocatePointIn(pvol_id, point, path, top, exclude);
     }
 
-    // Computes the isotropic safety from the globalpoint.
+    /// @brief Computes the isotropic safety from the globalpoint.
+    /// @param globalpoint Point in global coordinates
+    /// @param state Path where to compute safety
+    /// @return Isotropic safe distance
     CELER_FUNCTION static Precision
     ComputeSafety(Vector3D const& globalpoint,
-                  vecgeom::NavStateIndex const& state)
+                  vecgeom::NavigationState const& state)
     {
-        int closest_surf = 0;
-        return vgbrep::protonav::ComputeSafety(
-            globalpoint, state, closest_surf);
+        auto safety = vgbrep::protonav::BVHSurfNavigator<Precision>::ComputeSafety(globalpoint, state);
+	return safety;
     }
 
     // Computes a step from the globalpoint (which must be in the current
@@ -68,9 +78,8 @@ class SurfNavigator
     ComputeStepAndNextVolume(Vector3D const& globalpoint,
                              Vector3D const& globaldir,
                              Precision step_limit,
-                             vecgeom::NavStateIndex const& in_state,
-                             vecgeom::NavStateIndex& out_state,
-                             [[maybe_unused]] Precision push = 0)
+                             vecgeom::NavigationState const& in_state,
+                             vecgeom::NavigationState& out_state)
     {
         if (step_limit <= 0)
         {
@@ -79,9 +88,9 @@ class SurfNavigator
             return step_limit;
         }
 
-        int exit_surf = 0;
-        auto step = vgbrep::protonav::ComputeStepAndHit(
-            globalpoint, globaldir, in_state, out_state, exit_surf, step_limit);
+	long hitsurf_id = 0;
+        auto step = vgbrep::protonav::BVHSurfNavigator<Precision>::ComputeStepAndNextSurface(
+            globalpoint, globaldir, in_state, out_state, hitsurf_id, step_limit);
         return step;
     }
 
@@ -93,22 +102,29 @@ class SurfNavigator
     ComputeStepAndPropagatedState(Vector3D const& globalpoint,
                                   Vector3D const& globaldir,
                                   Precision step_limit,
-                                  vecgeom::NavStateIndex const& in_state,
-                                  vecgeom::NavStateIndex& out_state,
-                                  Precision push = 0)
+                                  vecgeom::NavigationState const& in,
+                                  vecgeom::NavigationState& out)
     {
-        return ComputeStepAndNextVolume(
-            globalpoint, globaldir, step_limit, in_state, out_state, push);
+      return ComputeStepAndNextVolume(globalpoint, globaldir, step_limit, in, out);
     }
 
     // Relocate a state that was returned from ComputeStepAndNextVolume: the
     // surface model does this computation within ComputeStepAndNextVolume, so
     // the relocation does nothing
     CELER_FUNCTION static void
-    RelocateToNextVolume(Vector3D const& /*globalpoint*/,
-                         Vector3D const& /*globaldir*/,
-                         vecgeom::NavStateIndex& /*state*/)
+    RelocateToNextVolume(Vector3D const& globalpoint,
+                         Vector3D const& globaldir,
+			 long hitsurf_index,
+                         vecgeom::NavigationState& out_state)
     {
+        vgbrep::CrossedSurface crossed_surf;
+        vgbrep::protonav::BVHSurfNavigator<Precision>::RelocateToNextVolume(
+	          globalpoint,
+	          globaldir,
+	          Precision(0),
+	          hitsurf_index,
+	          out_state,
+	          crossed_surf);
     }
 };