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