volumes/kernel/TessellatedImplementation.h

0001 //===-- kernel/TessellatedImplementation.h ----------------------------------*- C++ -*-===//
0002 //===--------------------------------------------------------------------------===//
0003 /// @file TessellatedImplementation.h
0004 /// @author mihaela.gheata@cern.ch
0005
0006 #ifndef VECGEOM_VOLUMES_KERNEL_TESSELLATEDIMPLEMENTATION_H_
0007 #define VECGEOM_VOLUMES_KERNEL_TESSELLATEDIMPLEMENTATION_H_
0008
0009 #include "VecGeom/base/Config.h"
0010 #include "VecGeom/base/Vector3D.h"
0011 #include "VecGeom/volumes/TessellatedStruct.h"
0012 #include "VecGeom/volumes/kernel/GenericKernels.h"
0013 #include <VecCore/VecCore>
0014
0015 #include <cstdio>
0016
0017 namespace vecgeom {
0018
0019 VECGEOM_DEVICE_FORWARD_DECLARE(struct TessellatedImplementation;);
0020 VECGEOM_DEVICE_DECLARE_CONV(struct, TessellatedImplementation);
0021
0022 inline namespace VECGEOM_IMPL_NAMESPACE {
0023
0024 class PlacedTessellated;
0025 template <size_t NVERT, typename T>
0026 class TessellatedStruct;
0027 class UnplacedTessellated;
0028
0029 struct TessellatedImplementation {
0030
0031   using PlacedShape_t    = PlacedTessellated;
0032   using UnplacedStruct_t = TessellatedStruct<3, Precision>;
0033   using UnplacedVolume_t = UnplacedTessellated;
0034
0035   template <typename Real_v, typename Bool_v>
0036   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void Contains(UnplacedStruct_t const &tessellated,
0037                                                                     Vector3D<Real_v> const &point, Bool_v &inside)
0038   {
0039     inside = Bool_v(false);
0040     int isurfOut, isurfIn;
0041     Real_v distOut, distIn;
0042     DistanceToSolid<Real_v, false>(tessellated, point, tessellated.fTestDir, InfinityLength<Real_v>(), distOut,
0043                                    isurfOut, distIn, isurfIn);
0044     if (isurfOut >= 0) inside = Bool_v(true);
0045     /*
0046         DistanceToSolid<Real_v, true>(tessellated, point, tessellated.fTestDir, stepMax, distIn, isurf);
0047         // If distance to out is finite and less than distance to in, the point is inside
0048         if (distOut < distIn) inside = Bool_v(true);
0049     */
0050   }
0051
0052   template <typename Real_v, typename Inside_v>
0053   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void Inside(UnplacedStruct_t const &tessellated,
0054                                                                   Vector3D<Real_v> const &point, Inside_v &inside)
0055   {
0056     inside = Inside_v(kOutside);
0057     int isurfOut, isurfIn;
0058     Real_v distOut, distIn;
0059     DistanceToSolid<Real_v, false>(tessellated, point, tessellated.fTestDir, InfinityLength<Real_v>(), distOut,
0060                                    isurfOut, distIn, isurfIn);
0061     // If no surface is hit then the point is outside
0062     if (isurfOut < 0) return;
0063     if (distOut < 0 || distOut * tessellated.fTestDir.Dot(tessellated.fFacets[isurfOut]->fNormal) < kTolerance) {
0064       inside = Inside_v(kSurface);
0065       return;
0066     }
0067
0068     // DistanceToSolid<Real_v, true>(tessellated, point, tessellated.fTestDir, stepMax, distIn, isurf);
0069     // If distance to out is finite and less than distance to in, the point is inside
0070     if (isurfIn < 0 || distOut < distIn) {
0071       inside = Inside_v(kInside);
0072       return;
0073     }
0074     if (distIn < 0 || distIn * tessellated.fTestDir.Dot(tessellated.fFacets[isurfIn]->fNormal) > -kTolerance)
0075       inside = Inside_v(kSurface);
0076   }
0077
0078   template <typename Real_v>
0079   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void DistanceToIn(UnplacedStruct_t const &tessellated,
0080                                                                         Vector3D<Real_v> const &point,
0081                                                                         Vector3D<Real_v> const &direction,
0082                                                                         Real_v const &stepMax, Real_v &distance)
0083   {
0084     int isurf, isurfOut;
0085     Real_v distOut;
0086     DistanceToSolid<Real_v, true>(tessellated, point, direction, stepMax, distance, isurf, distOut, isurfOut);
0087   }
0088
0089   template <typename Real_v>
0090   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void DistanceToOut(UnplacedStruct_t const &tessellated,
0091                                                                          Vector3D<Real_v> const &point,
0092                                                                          Vector3D<Real_v> const &direction,
0093                                                                          Real_v const &stepMax, Real_v &distance)
0094   {
0095     int isurf, isurfIn;
0096     Real_v distIn;
0097     DistanceToSolid<Real_v, false>(tessellated, point, direction, stepMax, distance, isurf, distIn, isurfIn);
0098   }
0099
0100   template <typename Real_v>
0101   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void SafetyToIn(UnplacedStruct_t const &tessellated,
0102                                                                       Vector3D<Real_v> const &point, Real_v &safety)
0103   {
0104     using Bool_v = vecCore::Mask_v<Real_v>;
0105     Bool_v inside;
0106     TessellatedImplementation::Contains<Real_v, Bool_v>(tessellated, point, inside);
0107     if (inside) {
0108       safety = -1.;
0109       return;
0110     }
0111     int isurf;
0112     Real_v safetysq = SafetySq<Real_v, true>(tessellated, point, isurf);
0113     safety          = vecCore::math::Sqrt(safetysq);
0114   }
0115
0116   template <typename Real_v>
0117   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void SafetyToOut(UnplacedStruct_t const &tessellated,
0118                                                                        Vector3D<Real_v> const &point, Real_v &safety)
0119   {
0120     using Bool_v = vecCore::Mask_v<Real_v>;
0121     Bool_v inside;
0122     TessellatedImplementation::Contains<Real_v, Bool_v>(tessellated, point, inside);
0123     if (!inside) {
0124       safety = -1.;
0125       return;
0126     }
0127     int isurf;
0128     Real_v safetysq = SafetySq<Real_v, false>(tessellated, point, isurf);
0129     safety          = vecCore::math::Sqrt(safetysq);
0130   }
0131
0132   template <typename Real_v>
0133   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static Vector3D<Real_v> NormalKernel(
0134       UnplacedStruct_t const &tessellated, Vector3D<Real_v> const &point, typename vecCore::Mask_v<Real_v> &valid)
0135   {
0136     // Computes the normal on a surface and returns it as a unit vector
0137     valid = true;
0138     int isurf;
0139     // We may need to check the value of safety to declare the validity of the normal
0140     SafetySq<Real_v, false>(tessellated, point, isurf);
0141     return tessellated.fFacets[isurf]->fNormal;
0142   }
0143
0144   template <typename Real_v, bool ToIn>
0145   VECCORE_ATT_HOST_DEVICE static void DistanceToSolid(UnplacedStruct_t const &tessellated,
0146                                                       Vector3D<Real_v> const &point, Vector3D<Real_v> const &direction,
0147                                                       Real_v const &stepMax, Real_v &distance, int &isurf,
0148                                                       Real_v &distother, int &isurfother)
0149   {
0150 // Common method providing DistanceToIn/Out functionality
0151 // Real_v here is scalar, we need to pass vector point/direction
0152 #ifndef VECGEOM_ENABLE_CUDA
0153     using Float_v = vecgeom::VectorBackend::Real_v;
0154 #else
0155     using Float_v             = vecgeom::ScalarBackend::Real_v;
0156 #endif
0157     isurf      = -1;
0158     isurfother = -1;
0159     if (ToIn) {
0160       // Check if the bounding box is hit
0161       const Vector3D<Real_v> invdir(Real_v(1.0) / NonZero(direction.x()), Real_v(1.0) / NonZero(direction.y()),
0162                                     Real_v(1.0) / NonZero(direction.z()));
0163       Vector3D<int> sign;
0164       sign[0]  = invdir.x() < 0;
0165       sign[1]  = invdir.y() < 0;
0166       sign[2]  = invdir.z() < 0;
0167       distance = BoxImplementation::IntersectCachedKernel2<Real_v, Real_v>(
0168           &tessellated.fMinExtent, point, invdir, sign.x(), sign.y(), sign.z(), -kTolerance, InfinityLength<Real_v>());
0169       if (distance >= stepMax) return;
0170     }
0171
0172     // Define the user hook calling DistanceToIn for the cluster with the same
0173     // index as the bounding box
0174     Vector3D<Float_v> pointv(point);
0175     Vector3D<Float_v> dirv(direction);
0176     distance             = InfinityLength<Real_v>();
0177     distother            = InfinityLength<Real_v>();
0178     Real_v distanceToIn  = InfinityLength<Real_v>();
0179     Real_v distanceToOut = InfinityLength<Real_v>();
0180     int isurfToIn        = -1;
0181     int isurfToOut       = -1;
0182     auto userhook        = [&](HybridManager2::BoxIdDistancePair_t hitbox) {
0183       // Stop searching if the distance to the current box is bigger than the
0184       // requested limit or than the current distance
0185       if (hitbox.second > vecCore::math::Min(stepMax, distance)) return true;
0186       // Compute distance to the cluster (in both ToIn or ToOut assumptions)
0187       Real_v clusterToIn, clusterToOut;
0188       int icrtToIn, icrtToOut;
0189       tessellated.fClusters[hitbox.first]->DistanceToCluster(pointv, dirv, clusterToIn, clusterToOut, icrtToIn,
0190                                                                     icrtToOut);
0191
0192       // Update distanceToIn/Out
0193       if (icrtToIn >= 0 && clusterToIn < distanceToIn) {
0194         distanceToIn = clusterToIn;
0195         isurfToIn    = icrtToIn;
0196         if (ToIn) {
0197           isurf    = isurfToIn;
0198           distance = distanceToIn;
0199         } else {
0200           isurfother = isurfToIn;
0201           distother  = distanceToIn;
0202         }
0203       }
0204
0205       if (icrtToOut >= 0 && clusterToOut < distanceToOut) {
0206         distanceToOut = clusterToOut;
0207         isurfToOut    = icrtToOut;
0208         if (!ToIn) {
0209           isurf    = isurfToOut;
0210           distance = distanceToOut;
0211         } else {
0212           isurfother = isurfToOut;
0213           distother  = distanceToOut;
0214         }
0215       }
0216       return false;
0217     };
0218
0219 #ifdef USEEMBREE
0220     EmbreeNavigator<> *boxNav = (EmbreeNavigator<> *)EmbreeNavigator<>::Instance();
0221     // intersect ray with the BVH structure and use hook
0222     boxNav->BVHSortedIntersectionsLooper(*tessellated.fNavHelper2, point, direction, 1E20, userhook);
0223 #else
0224     HybridNavigator<> *boxNav = (HybridNavigator<> *)HybridNavigator<>::Instance();
0225     boxNav->BVHSortedIntersectionsLooper(*tessellated.fNavHelper2, point, direction, stepMax, userhook);
0226 #endif
0227
0228     // Treat special cases
0229     if (ToIn) {
0230       if (isurfToIn < 0) {
0231         if (isurfToOut >= 0 && distanceToOut * direction.Dot(tessellated.fFacets[isurfToOut]->fNormal) > kTolerance)
0232           distance = -1.; // point inside or on boundary
0233         // else not hitting, distance already inf
0234       } else {
0235         if (isurfToOut >= 0 && distanceToOut > kTolerance && distanceToOut < distanceToIn)
0236           distance = -1.; // point inside exiting first then re-entering
0237         // else valid entry point, distance already set
0238       }
0239     } else {
0240       if (isurfToOut < 0)
0241         distance = -1.; // point outside
0242       else {
0243         if (isurfToIn >= 0 && distanceToIn < distanceToOut &&
0244             distanceToIn * direction.Dot(tessellated.fFacets[isurfToIn]->fNormal) < -kTolerance) {
0245           distance = -1.; // point outside (first entering then exiting)
0246           isurf    = -1;
0247         }
0248       }
0249     }
0250   }
0251
0252   template <typename Real_v, bool ToIn>
0253   VECCORE_ATT_HOST_DEVICE static Real_v SafetySq(UnplacedStruct_t const &tessellated, Vector3D<Real_v> const &point,
0254                                                  int &isurf)
0255   {
0256 #ifndef VECGEOM_ENABLE_CUDA
0257     using Float_v = vecgeom::VectorBackend::Real_v;
0258 #else
0259     using Float_v = vecgeom::ScalarBackend::Real_v;
0260 #endif
0261     Real_v safetysq = InfinityLength<Real_v>();
0262     isurf           = -1;
0263     Vector3D<Float_v> pointv(point);
0264
0265     auto userhook = [&](HybridManager2::BoxIdDistancePair_t hitbox) {
0266       // Stop searching if the safety to the current cluster is bigger than the
0267       // current safety
0268       if (hitbox.second > safetysq) return true;
0269       // Compute distance to the cluster
0270       int isurfcrt;
0271       Real_v safetycrt = tessellated.fClusters[hitbox.first]->template SafetySq<ToIn>(pointv, isurfcrt);
0272       if (safetycrt < safetysq) {
0273         safetysq = safetycrt;
0274         isurf    = isurfcrt;
0275       }
0276       return false;
0277     };
0278
0279     HybridSafetyEstimator *safEstimator = (HybridSafetyEstimator *)HybridSafetyEstimator::Instance();
0280     // Use the BVH structure and connect hook
0281     safEstimator->BVHSortedSafetyLooper(*tessellated.fNavHelper, point, userhook, safetysq);
0282     return safetysq;
0283   }
0284
0285 }; // end TessellatedImplementation
0286 } // namespace VECGEOM_IMPL_NAMESPACE
0287 } // namespace vecgeom
0288
0289 #endif // VECGEOM_VOLUMES_KERNEL_TESSELLATEDIMPLEMENTATION_H_