volumes/kernel/MultiUnionImplementation.h

0001 //===-- kernel/MultiUnionImplementation.h ---------------------------*- C++ -*-===//
0002 //===--------------------------------------------------------------------------===//
0003 /// @file MultiUnionImplementation.h
0004 /// @author Mihaela Gheata (mihaela.gheata@cern.ch)
0005
0006 #ifndef VECGEOM_VOLUMES_KERNEL_MULTIUNIONIMPLEMENTATION_H_
0007 #define VECGEOM_VOLUMES_KERNEL_MULTIUNIONIMPLEMENTATION_H_
0008
0009 #include "VecGeom/base/Vector3D.h"
0010 #include "VecGeom/volumes/MultiUnionStruct.h"
0011 #include "VecGeom/volumes/kernel/GenericKernels.h"
0012 #include <VecCore/VecCore>
0013
0014 #include <cstdio>
0015
0016 namespace vecgeom {
0017
0018 VECGEOM_DEVICE_FORWARD_DECLARE(struct MultiUnionImplementation;);
0019 VECGEOM_DEVICE_DECLARE_CONV(struct, MultiUnionImplementation);
0020
0021 inline namespace VECGEOM_IMPL_NAMESPACE {
0022
0023 class PlacedMultiUnion;
0024 struct MultiUnionStruct;
0025 class UnplacedMultiUnion;
0026
0027 struct MultiUnionImplementation {
0028
0029   using PlacedShape_t    = PlacedMultiUnion;
0030   using UnplacedStruct_t = MultiUnionStruct;
0031   using UnplacedVolume_t = UnplacedMultiUnion;
0032
0033   template <typename Real_v, typename Bool_v>
0034   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void Contains(UnplacedStruct_t const &munion,
0035                                                                     Vector3D<Real_v> const &point, Bool_v &inside)
0036   {
0037     auto containshook = [&](size_t id) {
0038       inside = munion.fVolumes[id]->Contains(point);
0039       return inside;
0040     };
0041
0042     HybridNavigator<> *boxNav = (HybridNavigator<> *)HybridNavigator<>::Instance();
0043     boxNav->BVHContainsLooper(*munion.fNavHelper, point, containshook);
0044   }
0045
0046   template <typename Real_v, typename Inside_v>
0047   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void Inside(UnplacedStruct_t const &munion,
0048                                                                   Vector3D<Real_v> const &point, Inside_v &inside)
0049   {
0050     inside          = EInside::kOutside;
0051     auto insidehook = [&](size_t id) {
0052       auto inside_crt = munion.fVolumes[id]->Inside(point);
0053       if (inside_crt == EInside::kInside) {
0054         inside = EInside::kInside;
0055         return true;
0056       }
0057       if (inside_crt == EInside::kSurface) inside = EInside::kSurface;
0058
0059       return false;
0060     };
0061
0062     HybridNavigator<> *boxNav = (HybridNavigator<> *)HybridNavigator<>::Instance();
0063     boxNav->BVHContainsLooper(*munion.fNavHelper, point, insidehook);
0064   }
0065
0066   template <typename Real_v>
0067   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void InsideComponent(UnplacedStruct_t const &munion,
0068                                                                            Vector3D<Real_v> const &point,
0069                                                                            int &component)
0070   {
0071     component       = -1;
0072     auto insidehook = [&](size_t id) {
0073       auto inside_crt = munion.fVolumes[id]->Inside(point);
0074       if (inside_crt != EInside::kOutside) {
0075         component = id;
0076         return true;
0077       }
0078       return false;
0079     };
0080
0081     HybridNavigator<> *boxNav = (HybridNavigator<> *)HybridNavigator<>::Instance();
0082     boxNav->BVHContainsLooper(*munion.fNavHelper, point, insidehook);
0083   }
0084
0085   template <typename Real_v>
0086   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void InsideCluster(UnplacedStruct_t const &munion,
0087                                                                          Vector3D<Real_v> const &point, int &component)
0088   {
0089     // loop cluster overlap candidates for current component, then update component
0090     size_t *cluster = munion.fNeighbours[component];
0091     size_t ncluster = munion.fNneighbours[component];
0092     for (size_t i = 0; i < ncluster; ++i) {
0093       if (munion.fVolumes[cluster[i]]->Inside(point) == EInside::kInside) {
0094         component = cluster[i];
0095         return;
0096       }
0097     }
0098     component = -1;
0099   }
0100
0101   template <typename Real_v>
0102   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void DistanceToIn(UnplacedStruct_t const &munion,
0103                                                                         Vector3D<Real_v> const &point,
0104                                                                         Vector3D<Real_v> const &direction,
0105                                                                         Real_v const &stepMax, Real_v &distance)
0106   {
0107     // Check if the bounding box is hit
0108     const Vector3D<Real_v> invdir(Real_v(1.0) / NonZero(direction.x()), Real_v(1.0) / NonZero(direction.y()),
0109                                   Real_v(1.0) / NonZero(direction.z()));
0110     Vector3D<int> sign;
0111     sign[0]  = invdir.x() < 0;
0112     sign[1]  = invdir.y() < 0;
0113     sign[2]  = invdir.z() < 0;
0114     distance = BoxImplementation::IntersectCachedKernel2<Real_v, Real_v>(
0115         &munion.fMinExtent, point, invdir, sign.x(), sign.y(), sign.z(), -kTolerance, InfinityLength<Real_v>());
0116     if (distance >= stepMax) return;
0117     distance = kInfLength;
0118     // Lambda function to be called for each candidate selected by the bounding box navigator
0119     auto userhook = [&](HybridManager2::BoxIdDistancePair_t hitbox) {
0120       // Stop searching if the distance to the current box is bigger than the
0121       // requested limit or than the current distance
0122       if (hitbox.second > vecCore::math::Min(stepMax, distance)) return true;
0123       // Compute distance to the cluster (in both ToIn or ToOut assumptions)
0124       auto distance_crt = munion.fVolumes[hitbox.first]->DistanceToIn(point, direction, stepMax);
0125       if (distance_crt < distance) distance = distance_crt;
0126       return false;
0127     };
0128
0129     HybridNavigator<> *boxNav = (HybridNavigator<> *)HybridNavigator<>::Instance();
0130     // intersect ray with the BVH structure and use hook
0131     boxNav->BVHSortedIntersectionsLooper(*munion.fNavHelper, point, direction, stepMax, userhook);
0132   }
0133
0134   template <typename Real_v>
0135   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void DistanceToOut(UnplacedStruct_t const &munion,
0136                                                                          Vector3D<Real_v> const &point,
0137                                                                          Vector3D<Real_v> const &direction,
0138                                                                          Real_v const &stepMax, Real_v &distance)
0139   {
0140     constexpr Real_v eps = 10 * kTolerance;
0141     distance             = -1.;
0142     // Locate the component containing the point
0143     int comp;
0144     InsideComponent(munion, point, comp);
0145     if (comp < 0) return; // Point not inside
0146     // Compute distance to exit current component
0147     distance              = -eps;
0148     Real_v dstep          = 1.;
0149     Vector3D<Real_v> pnew = point;
0150     Vector3D<Real_v> local, ldir;
0151     while (dstep > kTolerance && comp >= 0) {
0152       size_t component = (size_t)comp;
0153       munion.fVolumes[component]->GetTransformation()->Transform(pnew, local);
0154       munion.fVolumes[component]->GetTransformation()->TransformDirection(direction, ldir);
0155       dstep = munion.fVolumes[component]->DistanceToOut(local, ldir, stepMax);
0156       VECGEOM_ASSERT(dstep < kInfLength);
0157       distance += dstep + eps;
0158       // If no neighbours, exit
0159       if (!munion.fNneighbours[component]) return;
0160       // Propagate to exit of current component
0161       pnew += (dstep + eps) * direction;
0162       // Try to relocate inside the cluster of neighbours
0163       MultiUnionImplementation::InsideCluster(munion, pnew, comp);
0164     }
0165   }
0166
0167   template <typename Real_v>
0168   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void SafetyToInComp(UnplacedStruct_t const &munion,
0169                                                                           Vector3D<Real_v> const &point, Real_v &safety,
0170                                                                           int &component)
0171   {
0172     safety        = vecgeom::InfinityLength<Real_v>();
0173     component     = -1;
0174     auto userhook = [&](HybridManager2::BoxIdDistancePair_t hitbox) {
0175       // Stop searching if the safety to the current cluster is bigger than the
0176       // current safety
0177       if (hitbox.second > safety * safety) return true;
0178       // Compute distance to the cluster
0179       Real_v safetycrt = munion.fVolumes[hitbox.first]->SafetyToIn(point);
0180       if (safetycrt > 0 && safetycrt < safety) {
0181         safety    = safetycrt;
0182         component = hitbox.first;
0183       }
0184       return false;
0185     };
0186
0187     HybridSafetyEstimator *safEstimator = (HybridSafetyEstimator *)HybridSafetyEstimator::Instance();
0188     // Use the BVH structure and connect hook
0189     safEstimator->BVHSortedSafetyLooper(*munion.fNavHelper, point, userhook, safety);
0190   }
0191
0192   template <typename Real_v>
0193   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void SafetyToIn(UnplacedStruct_t const &munion,
0194                                                                       Vector3D<Real_v> const &point, Real_v &safety)
0195   {
0196     int comp;
0197     InsideComponent(munion, point, comp);
0198     if (comp > -1) {
0199       safety = -1.;
0200       return;
0201     }
0202
0203     SafetyToInComp<Real_v>(munion, point, safety, comp);
0204   }
0205
0206   template <typename Real_v>
0207   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void SafetyToOut(UnplacedStruct_t const &munion,
0208                                                                        Vector3D<Real_v> const &point, Real_v &safety)
0209   {
0210     // Locate the component containing the point
0211     int comp;
0212     MultiUnionImplementation::InsideComponent(munion, point, comp);
0213     if (comp < 0) {
0214       safety = -1.; // Point not inside
0215       return;
0216     }
0217     // Compute safety to exit current component
0218     Vector3D<Real_v> const local = munion.fVolumes[comp]->GetTransformation()->Transform(point);
0219     safety                       = munion.fVolumes[comp]->SafetyToOut(local);
0220     VECGEOM_ASSERT(safety > -kTolerance);
0221     // Loop cluster of neighbours
0222     size_t *cluster = munion.fNeighbours[comp];
0223     size_t ncluster = munion.fNneighbours[comp];
0224     for (size_t i = 0; i < ncluster; ++i) {
0225       Vector3D<Real_v> const local = munion.fVolumes[cluster[i]]->GetTransformation()->Transform(point);
0226       Real_v safetycrt             = munion.fVolumes[cluster[i]]->SafetyToOut(local);
0227       if (safetycrt > 0 && safetycrt < safety) safety = safetycrt;
0228     }
0229   }
0230
0231   template <typename Real_v>
0232   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static Vector3D<Real_v> NormalKernel(
0233       UnplacedStruct_t const &munion, Vector3D<Real_v> const &point, typename vecCore::Mask_v<Real_v> &valid)
0234   {
0235     // Locate the component containing the point
0236     int comp;
0237     valid = false;
0238     Vector3D<Real_v> direction;
0239
0240     InsideComponent(munion, point, comp);
0241     // If component not found, locate closest one
0242     Real_v safety;
0243     if (comp < 0) SafetyToInComp(munion, point, safety, comp);
0244     if (comp < 0) return direction;
0245
0246     Vector3D<Real_v> local = munion.fVolumes[comp]->GetTransformation()->Transform(point);
0247     Vector3D<Real_v> ldir;
0248     valid = munion.fVolumes[comp]->Normal(local, ldir);
0249     if (valid) direction = munion.fVolumes[comp]->GetTransformation()->InverseTransformDirection(ldir);
0250     return direction;
0251   }
0252
0253 }; // end MultiUnionImplementation
0254 } // namespace VECGEOM_IMPL_NAMESPACE
0255 } // namespace vecgeom
0256
0257 #endif // VECGEOM_VOLUMES_KERNEL_MULTIUNIONIMPLEMENTATION_H_