volumes/kernel/BooleanIntersectionImplementation.h

0001 /*
0002  * BooleanImplementation.h
0003  */
0004
0005 #ifndef BOOLEANINTERSECTIONIMPLEMENTATION_H_
0006 #define BOOLEANINTERSECTIONIMPLEMENTATION_H_
0007
0008 #include "VecGeom/base/Global.h"
0009 #include "VecGeom/base/Vector3D.h"
0010 #include "VecGeom/volumes/BooleanStruct.h"
0011 #include <VecCore/VecCore>
0012
0013 namespace vecgeom {
0014
0015 inline namespace VECGEOM_IMPL_NAMESPACE {
0016
0017 /**
0018  * partial template specialization for UNION implementation
0019  */
0020 template <>
0021 struct BooleanImplementation<kIntersection> {
0022   using PlacedShape_t    = PlacedBooleanVolume<kIntersection>;
0023   using UnplacedVolume_t = UnplacedBooleanVolume<kIntersection>;
0024   using UnplacedStruct_t = BooleanStruct;
0025
0026   template <typename Real_v, typename Bool_v>
0027   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void Contains(BooleanStruct const &unplaced,
0028                                                                     Vector3D<Real_v> const &point, Bool_v &inside)
0029   {
0030     const auto insideA = unplaced.fLeftVolume->Contains(point);
0031     const auto insideB = unplaced.fRightVolume->Contains(point);
0032     inside             = insideA && insideB;
0033   }
0034
0035   template <typename Real_v, typename Inside_t>
0036   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void Inside(BooleanStruct const &unplaced,
0037                                                                   Vector3D<Real_v> const &point, Inside_t &inside)
0038   {
0039     // now use the Inside functionality of left and right components
0040     // algorithm taken from Geant4 implementation
0041     VPlacedVolume const *const fPtrSolidA = unplaced.fLeftVolume;
0042     VPlacedVolume const *const fPtrSolidB = unplaced.fRightVolume;
0043
0044     const auto positionA = fPtrSolidA->Inside(point);
0045
0046     if (positionA == EInside::kOutside) {
0047       inside = EInside::kOutside;
0048       return;
0049     }
0050
0051     const auto positionB = fPtrSolidB->Inside(point);
0052     if (positionA == EInside::kInside && positionB == EInside::kInside) {
0053       inside = EInside::kInside;
0054       return;
0055     } else {
0056       if ((positionA == EInside::kInside && positionB == EInside::kSurface) ||
0057           (positionB == EInside::kInside && positionA == EInside::kSurface) ||
0058           (positionA == EInside::kSurface && positionB == EInside::kSurface)) {
0059         inside = EInside::kSurface;
0060         return;
0061       } else {
0062         inside = EInside::kOutside;
0063         return;
0064       }
0065     }
0066   }
0067
0068   template <typename Real_v>
0069   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void DistanceToIn(BooleanStruct const &unplaced,
0070                                                                         Vector3D<Real_v> const &point,
0071                                                                         Vector3D<Real_v> const &dir,
0072                                                                         Real_v const &stepMax, Real_v &distance)
0073   {
0074     Vector3D<Real_v> hitpoint = point;
0075
0076     auto inleft  = unplaced.fLeftVolume->Contains(hitpoint);
0077     auto inright = unplaced.fRightVolume->Contains(hitpoint);
0078     Real_v d1    = 0.;
0079     Real_v d2    = 0.;
0080     Real_v snext = 0.0;
0081
0082     // just a pre-check before entering main algorithm
0083     if (inleft && inright) {
0084       d1 = unplaced.fLeftVolume->PlacedDistanceToOut(hitpoint, dir, stepMax);
0085       d2 = unplaced.fRightVolume->PlacedDistanceToOut(hitpoint, dir, stepMax);
0086
0087       // if we are close to a boundary continue
0088       if (d1 < 2 * kTolerance) inleft = false;  // Backend::kFalse;
0089       if (d2 < 2 * kTolerance) inright = false; // Backend::kFalse;
0090
0091       // otherwise exit
0092       if (inleft && inright) {
0093         // TODO: WE are inside both so should return a negative number
0094         distance = 0.0;
0095         return;
0096       }
0097     }
0098
0099     // main loop
0100     while (1) {
0101       d1 = d2 = 0;
0102       if (!inleft) {
0103         d1 = unplaced.fLeftVolume->DistanceToIn(hitpoint, dir);
0104         d1 = Max(d1, kTolerance);
0105         if (d1 > 1E20) {
0106           distance = kInfLength;
0107           return;
0108         }
0109       }
0110       if (!inright) {
0111         d2 = unplaced.fRightVolume->DistanceToIn(hitpoint, dir);
0112         d2 = Max(d2, kTolerance);
0113         if (d2 > 1E20) {
0114           distance = kInfLength;
0115           return;
0116         }
0117       }
0118
0119       if (d1 > d2) {
0120         // propagate to left shape
0121         snext += d1;
0122         inleft = true; // Backend::kTrue;
0123         hitpoint += d1 * dir;
0124
0125         // check if propagated point is inside right shape
0126         // check is done with a little push
0127         inright = unplaced.fRightVolume->Contains(hitpoint + kTolerance * dir);
0128         if (inright) {
0129           distance = snext;
0130           return;
0131         }
0132         // here inleft=true, inright=false
0133       } else {
0134         // propagate to right shape
0135         snext += d2;
0136         inright = true; // Backend::kTrue;
0137         hitpoint += d2 * dir;
0138
0139         // check if propagated point is inside left shape
0140         inleft = unplaced.fLeftVolume->Contains(hitpoint + kTolerance * dir);
0141         if (inleft) {
0142           distance = snext;
0143           return;
0144         }
0145       }
0146       // here inleft=false, inright=true
0147     } // end while loop
0148     distance = snext;
0149     return;
0150   }
0151
0152   template <typename Real_v>
0153   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void DistanceToOut(BooleanStruct const &unplaced,
0154                                                                          Vector3D<Real_v> const &point,
0155                                                                          Vector3D<Real_v> const &direction,
0156                                                                          Real_v const &stepMax, Real_v &distance)
0157   {
0158     distance = Min(unplaced.fLeftVolume->DistanceToOut(point, direction),
0159                    unplaced.fRightVolume->PlacedDistanceToOut(point, direction));
0160   }
0161
0162   template <typename Real_v>
0163   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void SafetyToIn(BooleanStruct const &unplaced,
0164                                                                       Vector3D<Real_v> const &point, Real_v &safety)
0165   {
0166     // This is the Geant4 algorithm
0167     // TODO: ROOT seems to produce better safeties
0168     const auto insideA = unplaced.fLeftVolume->Contains(point);
0169     const auto insideB = unplaced.fRightVolume->Contains(point);
0170
0171     if (!insideA && insideB) {
0172       safety = unplaced.fLeftVolume->SafetyToIn(point);
0173     } else {
0174       if (!insideB && insideA) {
0175         safety = unplaced.fRightVolume->SafetyToIn(point);
0176       } else {
0177         safety = Min(unplaced.fLeftVolume->SafetyToIn(point), unplaced.fRightVolume->SafetyToIn(point));
0178       }
0179     }
0180     return;
0181   }
0182
0183   template <typename Real_v>
0184   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void SafetyToOut(BooleanStruct const &unplaced,
0185                                                                        Vector3D<Real_v> const &point, Real_v &safety)
0186   {
0187     safety = Min(
0188         // TODO: could fail if left volume is placed shape
0189         unplaced.fLeftVolume->SafetyToOut(point),
0190
0191         // TODO: consider introducing PlacedSafetyToOut function
0192         unplaced.fRightVolume->SafetyToOut(unplaced.fRightVolume->GetTransformation()->Transform(point)));
0193     vecCore::MaskedAssign(safety, safety < Real_v(0.), Real_v(0.));
0194   }
0195
0196   template <typename Real_v, typename Bool_v>
0197   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void NormalKernel(BooleanStruct const &unplaced,
0198                                                                         Vector3D<Real_v> const &point,
0199                                                                         Vector3D<Real_v> &normal, Bool_v &valid)
0200   {
0201     Vector3D<Real_v> localNorm;
0202     Vector3D<Real_v> localPoint;
0203     valid = false; // Backend::kFalse;
0204
0205     VPlacedVolume const *const fPtrSolidA = unplaced.fLeftVolume;
0206     VPlacedVolume const *const fPtrSolidB = unplaced.fRightVolume;
0207     Real_v safetyA, safetyB;
0208
0209     if (fPtrSolidA->Contains(point)) {
0210       fPtrSolidA->GetTransformation()->Transform(point, localPoint);
0211       safetyA = fPtrSolidA->SafetyToOut(localPoint);
0212     } else {
0213       safetyA = fPtrSolidA->SafetyToIn(point);
0214     }
0215
0216     if (fPtrSolidB->Contains(point)) {
0217       fPtrSolidB->GetTransformation()->Transform(point, localPoint);
0218       safetyB = fPtrSolidB->SafetyToOut(localPoint);
0219     } else {
0220       safetyB = fPtrSolidB->SafetyToIn(point);
0221     }
0222     const auto onA = safetyA < safetyB;
0223     if (vecCore::MaskFull(onA)) {
0224       fPtrSolidA->GetTransformation()->Transform(point, localPoint);
0225       valid = fPtrSolidA->Normal(localPoint, localNorm);
0226       fPtrSolidA->GetTransformation()->InverseTransformDirection(localNorm, normal);
0227       return;
0228     } else {
0229       //  if (vecCore::MaskEmpty(onA)) {  // to use real mask operation when supporting vectors
0230       fPtrSolidB->GetTransformation()->Transform(point, localPoint);
0231       valid = fPtrSolidB->Normal(localPoint, localNorm);
0232       fPtrSolidB->GetTransformation()->InverseTransformDirection(localNorm, normal);
0233       return;
0234     }
0235     // Some particles are on A, some on B. We never arrive here in the scalar case
0236     // If the interface to Normal will support the vector case, we have to write code here.
0237     return;
0238   }
0239 }; // End struct BooleanImplementation
0240
0241 } // namespace VECGEOM_IMPL_NAMESPACE
0242
0243 } // namespace vecgeom
0244
0245 #endif /* BooleanImplementation_H_ */