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   VECCORE_ATT_HOST_DEVICE
0034   static void PrintType() {}
0035
0036   template <typename Stream>
0037   static void PrintType(Stream &st, int transCodeT = translation::kGeneric, int rotCodeT = rotation::kGeneric)
0038   {
0039     st << "SpecializedMultiUnion<" << transCodeT << "," << rotCodeT << ">";
0040   }
0041
0042   template <typename Stream>
0043   static void PrintImplementationType(Stream &st)
0044   {
0045     (void)st;
0046   }
0047
0048   template <typename Stream>
0049   static void PrintUnplacedType(Stream &st)
0050   {
0051     (void)st;
0052   }
0053
0054   template <typename Real_v, typename Bool_v>
0055   VECGEOM_FORCE_INLINE
0056   VECCORE_ATT_HOST_DEVICE
0057   static void Contains(UnplacedStruct_t const &munion, Vector3D<Real_v> const &point, Bool_v &inside)
0058   {
0059     auto containshook = [&](size_t id) {
0060       inside = munion.fVolumes[id]->Contains(point);
0061       return inside;
0062     };
0063
0064     HybridNavigator<> *boxNav = (HybridNavigator<> *)HybridNavigator<>::Instance();
0065     boxNav->BVHContainsLooper(*munion.fNavHelper, point, containshook);
0066   }
0067
0068   template <typename Real_v, typename Inside_v>
0069   VECGEOM_FORCE_INLINE
0070   VECCORE_ATT_HOST_DEVICE
0071   static void Inside(UnplacedStruct_t const &munion, Vector3D<Real_v> const &point, Inside_v &inside)
0072   {
0073     inside          = EInside::kOutside;
0074     auto insidehook = [&](size_t id) {
0075       auto inside_crt = munion.fVolumes[id]->Inside(point);
0076       if (inside_crt == EInside::kInside) {
0077         inside = EInside::kInside;
0078         return true;
0079       }
0080       if (inside_crt == EInside::kSurface) inside = EInside::kSurface;
0081
0082       return false;
0083     };
0084
0085     HybridNavigator<> *boxNav = (HybridNavigator<> *)HybridNavigator<>::Instance();
0086     boxNav->BVHContainsLooper(*munion.fNavHelper, point, insidehook);
0087   }
0088
0089   template <typename Real_v>
0090   VECGEOM_FORCE_INLINE
0091   VECCORE_ATT_HOST_DEVICE
0092   static void InsideComponent(UnplacedStruct_t const &munion, Vector3D<Real_v> const &point, int &component)
0093   {
0094     component       = -1;
0095     auto insidehook = [&](size_t id) {
0096       auto inside_crt = munion.fVolumes[id]->Inside(point);
0097       if (inside_crt != EInside::kOutside) {
0098         component = id;
0099         return true;
0100       }
0101       return false;
0102     };
0103
0104     HybridNavigator<> *boxNav = (HybridNavigator<> *)HybridNavigator<>::Instance();
0105     boxNav->BVHContainsLooper(*munion.fNavHelper, point, insidehook);
0106   }
0107
0108   template <typename Real_v>
0109   VECGEOM_FORCE_INLINE
0110   VECCORE_ATT_HOST_DEVICE
0111   static void InsideCluster(UnplacedStruct_t const &munion, Vector3D<Real_v> const &point, int &component)
0112   {
0113     // loop cluster overlap candidates for current component, then update component
0114     size_t *cluster = munion.fNeighbours[component];
0115     size_t ncluster = munion.fNneighbours[component];
0116     for (size_t i = 0; i < ncluster; ++i) {
0117       if (munion.fVolumes[cluster[i]]->Inside(point) == EInside::kInside) {
0118         component = cluster[i];
0119         return;
0120       }
0121     }
0122     component = -1;
0123   }
0124
0125   template <typename Real_v>
0126   VECGEOM_FORCE_INLINE
0127   VECCORE_ATT_HOST_DEVICE
0128   static void DistanceToIn(UnplacedStruct_t const &munion, Vector3D<Real_v> const &point,
0129                            Vector3D<Real_v> const &direction, Real_v const &stepMax, Real_v &distance)
0130   {
0131     // Check if the bounding box is hit
0132     const Vector3D<Real_v> invdir(Real_v(1.0) / NonZero(direction.x()), Real_v(1.0) / NonZero(direction.y()),
0133                                   Real_v(1.0) / NonZero(direction.z()));
0134     Vector3D<int> sign;
0135     sign[0]  = invdir.x() < 0;
0136     sign[1]  = invdir.y() < 0;
0137     sign[2]  = invdir.z() < 0;
0138     distance = BoxImplementation::IntersectCachedKernel2<Real_v, Real_v>(
0139         &munion.fMinExtent, point, invdir, sign.x(), sign.y(), sign.z(), -kTolerance, InfinityLength<Real_v>());
0140     if (distance >= stepMax) return;
0141     distance = kInfLength;
0142     // Lambda function to be called for each candidate selected by the bounding box navigator
0143     auto userhook = [&](HybridManager2::BoxIdDistancePair_t hitbox) {
0144       // Stop searching if the distance to the current box is bigger than the
0145       // requested limit or than the current distance
0146       if (hitbox.second > vecCore::math::Min(stepMax, distance)) return true;
0147       // Compute distance to the cluster (in both ToIn or ToOut assumptions)
0148       auto distance_crt = munion.fVolumes[hitbox.first]->DistanceToIn(point, direction, stepMax);
0149       if (distance_crt < distance) distance = distance_crt;
0150       return false;
0151     };
0152
0153     HybridNavigator<> *boxNav = (HybridNavigator<> *)HybridNavigator<>::Instance();
0154     // intersect ray with the BVH structure and use hook
0155     boxNav->BVHSortedIntersectionsLooper(*munion.fNavHelper, point, direction, stepMax, userhook);
0156   }
0157
0158   template <typename Real_v>
0159   VECGEOM_FORCE_INLINE
0160   VECCORE_ATT_HOST_DEVICE
0161   static void DistanceToOut(UnplacedStruct_t const &munion, Vector3D<Real_v> const &point,
0162                             Vector3D<Real_v> const &direction, Real_v const &stepMax, Real_v &distance)
0163   {
0164     constexpr Real_v eps = 10 * kTolerance;
0165     distance             = -1.;
0166     // Locate the component containing the point
0167     int comp;
0168     InsideComponent(munion, point, comp);
0169     if (comp < 0) return; // Point not inside
0170     // Compute distance to exit current component
0171     distance              = -eps;
0172     Real_v dstep          = 1.;
0173     Vector3D<Real_v> pnew = point;
0174     Vector3D<Real_v> local, ldir;
0175     while (dstep > kTolerance && comp >= 0) {
0176       size_t component = (size_t)comp;
0177       munion.fVolumes[component]->GetTransformation()->Transform(pnew, local);
0178       munion.fVolumes[component]->GetTransformation()->TransformDirection(direction, ldir);
0179       dstep = munion.fVolumes[component]->DistanceToOut(local, ldir, stepMax);
0180       assert(dstep < kInfLength);
0181       distance += dstep + eps;
0182       // If no neighbours, exit
0183       if (!munion.fNneighbours[component]) return;
0184       // Propagate to exit of current component
0185       pnew += (dstep + eps) * direction;
0186       // Try to relocate inside the cluster of neighbours
0187       MultiUnionImplementation::InsideCluster(munion, pnew, comp);
0188     }
0189   }
0190
0191   template <typename Real_v>
0192   VECGEOM_FORCE_INLINE
0193   VECCORE_ATT_HOST_DEVICE
0194   static void SafetyToInComp(UnplacedStruct_t const &munion, Vector3D<Real_v> const &point, Real_v &safety,
0195                              int &component)
0196   {
0197     safety        = vecgeom::InfinityLength<Real_v>();
0198     component     = -1;
0199     auto userhook = [&](HybridManager2::BoxIdDistancePair_t hitbox) {
0200       // Stop searching if the safety to the current cluster is bigger than the
0201       // current safety
0202       if (hitbox.second > safety * safety) return true;
0203       // Compute distance to the cluster
0204       Real_v safetycrt = munion.fVolumes[hitbox.first]->SafetyToIn(point);
0205       if (safetycrt > 0 && safetycrt < safety) {
0206         safety    = safetycrt;
0207         component = hitbox.first;
0208       }
0209       return false;
0210     };
0211
0212     HybridSafetyEstimator *safEstimator = (HybridSafetyEstimator *)HybridSafetyEstimator::Instance();
0213     // Use the BVH structure and connect hook
0214     safEstimator->BVHSortedSafetyLooper(*munion.fNavHelper, point, userhook, safety);
0215   }
0216
0217   template <typename Real_v>
0218   VECGEOM_FORCE_INLINE
0219   VECCORE_ATT_HOST_DEVICE
0220   static void SafetyToIn(UnplacedStruct_t const &munion, Vector3D<Real_v> const &point, Real_v &safety)
0221   {
0222     int comp;
0223     InsideComponent(munion, point, comp);
0224     if (comp > -1) {
0225       safety = -1.;
0226       return;
0227     }
0228
0229     SafetyToInComp<Real_v>(munion, point, safety, comp);
0230   }
0231
0232   template <typename Real_v>
0233   VECGEOM_FORCE_INLINE
0234   VECCORE_ATT_HOST_DEVICE
0235   static void SafetyToOut(UnplacedStruct_t const &munion, Vector3D<Real_v> const &point, Real_v &safety)
0236   {
0237     // Locate the component containing the point
0238     int comp;
0239     MultiUnionImplementation::InsideComponent(munion, point, comp);
0240     if (comp < 0) {
0241       safety = -1.; // Point not inside
0242       return;
0243     }
0244     // Compute safety to exit current component
0245     Vector3D<Real_v> const local = munion.fVolumes[comp]->GetTransformation()->Transform(point);
0246     safety                       = munion.fVolumes[comp]->SafetyToOut(local);
0247     assert(safety > -kTolerance);
0248     // Loop cluster of neighbours
0249     size_t *cluster = munion.fNeighbours[comp];
0250     size_t ncluster = munion.fNneighbours[comp];
0251     for (size_t i = 0; i < ncluster; ++i) {
0252       Vector3D<Real_v> const local = munion.fVolumes[cluster[i]]->GetTransformation()->Transform(point);
0253       Real_v safetycrt             = munion.fVolumes[cluster[i]]->SafetyToOut(local);
0254       if (safetycrt > 0 && safetycrt < safety) safety = safetycrt;
0255     }
0256   }
0257
0258   template <typename Real_v>
0259   VECGEOM_FORCE_INLINE
0260   VECCORE_ATT_HOST_DEVICE
0261   static Vector3D<Real_v> NormalKernel(UnplacedStruct_t const &munion, Vector3D<Real_v> const &point,
0262                                        typename vecCore::Mask_v<Real_v> &valid)
0263   {
0264     // Locate the component containing the point
0265     int comp;
0266     valid = false;
0267     Vector3D<Real_v> direction;
0268
0269     InsideComponent(munion, point, comp);
0270     // If component not found, locate closest one
0271     Real_v safety;
0272     if (comp < 0) SafetyToInComp(munion, point, safety, comp);
0273     if (comp < 0) return direction;
0274
0275     Vector3D<Real_v> local = munion.fVolumes[comp]->GetTransformation()->Transform(point);
0276     Vector3D<Real_v> ldir;
0277     valid = munion.fVolumes[comp]->Normal(local, ldir);
0278     if (valid) direction = munion.fVolumes[comp]->GetTransformation()->InverseTransformDirection(ldir);
0279     return direction;
0280   }
0281
0282 }; // end MultiUnionImplementation
0283 } // namespace VECGEOM_IMPL_NAMESPACE
0284 } // namespace vecgeom
0285
0286 #endif // VECGEOM_VOLUMES_KERNEL_MULTIUNIONIMPLEMENTATION_H_