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