volumes/kernel/BooleanImplementation.h

0001 /*
0002  * BooleanImplementation.h
0003  */
0004
0005 #ifndef BOOLEANIMPLEMENTATION_H_
0006 #define BOOLEANIMPLEMENTATION_H_
0007
0008 #include "VecGeom/base/Vector3D.h"
0009 #include "VecGeom/volumes/BooleanStruct.h"
0010 #include <VecCore/VecCore>
0011
0012 namespace vecgeom {
0013
0014 VECGEOM_DEVICE_DECLARE_CONV_TEMPLATE_1v(struct, BooleanImplementation, BooleanOperation, Arg1);
0015
0016 inline namespace VECGEOM_IMPL_NAMESPACE {
0017
0018 template <BooleanOperation Op>
0019 class PlacedBooleanVolume;
0020 template <BooleanOperation Op>
0021 class UnplacedBooleanVolume;
0022
0023 template <BooleanOperation boolOp>
0024 struct BooleanImplementation {
0025   using PlacedShape_t    = PlacedBooleanVolume<boolOp>;
0026   using UnplacedVolume_t = UnplacedBooleanVolume<boolOp>;
0027   using UnplacedStruct_t = BooleanStruct;
0028
0029   // empty since functionality will be implemented in
0030   // partially template specialized structs
0031 };
0032
0033 /**
0034  * an ordinary (non-templated) implementation of a Boolean solid
0035  * using the virtual function interface of its constituents
0036  *
0037  * TEMPLATE SPECIALIZATION FOR SUBTRACTION
0038  */
0039 template <>
0040 struct BooleanImplementation<kSubtraction> {
0041   using PlacedShape_t    = PlacedBooleanVolume<kSubtraction>;
0042   using UnplacedVolume_t = UnplacedBooleanVolume<kSubtraction>;
0043   using UnplacedStruct_t = BooleanStruct;
0044
0045   template <typename Real_v, typename Bool_v>
0046   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void Contains(BooleanStruct const &unplaced,
0047                                                                     Vector3D<Real_v> const &point, Bool_v &inside)
0048   {
0049     Vector3D<Real_v> tmp;
0050     inside = unplaced.fLeftVolume->Contains(point);
0051     if (vecCore::MaskEmpty(inside)) return;
0052
0053     auto rightInside = unplaced.fRightVolume->Contains(point);
0054     inside &= !rightInside;
0055   }
0056
0057   template <typename Real_v, typename Inside_t>
0058   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void Inside(BooleanStruct const &unplaced,
0059                                                                   Vector3D<Real_v> const &p, Inside_t &inside)
0060   {
0061
0062     // now use the Inside functionality of left and right components
0063     // algorithm taken from Geant4 implementation
0064     VPlacedVolume const *const fPtrSolidA = unplaced.fLeftVolume;
0065     VPlacedVolume const *const fPtrSolidB = unplaced.fRightVolume;
0066
0067     const auto positionA = fPtrSolidA->Inside(p);
0068     if (positionA == EInside::kOutside) {
0069       inside = EInside::kOutside;
0070       return;
0071     }
0072
0073     const auto positionB = fPtrSolidB->Inside(p);
0074
0075     if (positionA == EInside::kInside && positionB == EInside::kOutside) {
0076       inside = EInside::kInside;
0077       return;
0078     } else {
0079       if ((positionA == EInside::kInside && positionB == EInside::kSurface) ||
0080           (positionB == EInside::kOutside && positionA == EInside::kSurface)
0081           /*
0082            ||( positionA == EInside::kSurface && positionB == EInside::kSurface &&
0083              (   fPtrSolidA->Normal(p) -
0084                fPtrSolidB->Normal(p) ).mag2() >
0085              1000.0*G4GeometryTolerance::GetInstance()->GetRadialTolerance() ) )
0086           */) {
0087         inside = EInside::kSurface;
0088         return;
0089       } else {
0090         inside = EInside::kOutside;
0091         return;
0092       }
0093     }
0094     // going to be a bit more complicated due to Surface states
0095   }
0096
0097   template <typename Real_v>
0098   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void DistanceToIn(BooleanStruct const &unplaced,
0099                                                                         Vector3D<Real_v> const &p,
0100                                                                         Vector3D<Real_v> const &dir,
0101                                                                         Real_v const &stepMax, Real_v &distance)
0102   {
0103     Real_v dist_right, dist_left, advance = 0.;
0104     Real_v limit = stepMax;
0105     Vector3D<Real_v> hitpoint(p);
0106     // check if inside '-'
0107     auto insideRight = unplaced.fRightVolume->Inside(p) != kOutside;
0108     // epsil is used to push across boundaries, push records the actual push
0109     Precision epsil(0.), push(0.);
0110     while (1) {
0111       if (insideRight) {
0112         //    // propagate to outside of '- / RightShape'
0113         dist_right = unplaced.fRightVolume->PlacedDistanceToOut(hitpoint, dir, limit);
0114         if (dist_right >= 0.) {
0115           advance += dist_right + push;
0116           limit = stepMax - advance;
0117           epsil = kRelTolerance(hitpoint + dist_right * dir);
0118           // Push point across the boundary and record push
0119           hitpoint += (dist_right + epsil) * dir;
0120           push = epsil;
0121         } else {
0122           push = 0.;
0123         }
0124
0125         // now master outside 'B'; check if inside 'A'
0126         if (unplaced.fLeftVolume->Inside(hitpoint) != kOutside) {
0127           auto check = unplaced.fLeftVolume->PlacedDistanceToOut(hitpoint, dir);
0128           if (check > epsil) {
0129             distance = advance;
0130             return;
0131           }
0132         }
0133       }
0134
0135       // master outside '-' and outside '+' ;  find distances to both
0136       dist_left = unplaced.fLeftVolume->DistanceToIn(hitpoint, dir, limit);
0137       dist_left = vecCore::math::Max(dist_left, 0.);
0138       if (dist_left >= limit) {
0139         distance = kInfLength;
0140         return;
0141       }
0142
0143       dist_right = unplaced.fRightVolume->DistanceToIn(hitpoint, dir, limit);
0144       if (dist_left < dist_right - kTolerance) {
0145         advance += dist_left + push;
0146         distance = advance;
0147         return;
0148       }
0149
0150       //        // propagate to '-'
0151       if (dist_right >= 0. && dist_right < kInfLength) {
0152         advance += dist_right + push;
0153         limit = stepMax - advance;
0154         epsil = kRelTolerance(hitpoint + dist_right * dir);
0155         hitpoint += (dist_right + epsil) * dir;
0156         push = epsil;
0157       } else {
0158         push = 0.;
0159       }
0160       insideRight = true;
0161     } // end while
0162   }
0163
0164   template <typename Real_v>
0165   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void DistanceToOut(BooleanStruct const &unplaced,
0166                                                                          Vector3D<Real_v> const &point,
0167                                                                          Vector3D<Real_v> const &direction,
0168                                                                          Real_v const &stepMax, Real_v &distance)
0169   {
0170     const auto distancel  = unplaced.fLeftVolume->PlacedDistanceToOut(point, direction, stepMax);
0171     const Real_v dinright = unplaced.fRightVolume->DistanceToIn(point, direction, stepMax);
0172     distance              = Min(distancel, dinright);
0173     return;
0174   }
0175
0176   template <typename Real_v>
0177   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void SafetyToIn(BooleanStruct const &unplaced,
0178                                                                       Vector3D<Real_v> const &point, Real_v &safety)
0179   {
0180     VPlacedVolume const *const fPtrSolidA = unplaced.fLeftVolume;
0181     VPlacedVolume const *const fPtrSolidB = unplaced.fRightVolume;
0182
0183     // very approximate
0184     if ((fPtrSolidA->Contains(point)) && // case 1
0185         (fPtrSolidB->Contains(point))) {
0186       safety = fPtrSolidB->SafetyToOut(fPtrSolidB->GetTransformation()->Transform(point));
0187     } else {
0188       // po
0189       safety = fPtrSolidA->SafetyToIn(point);
0190     }
0191   }
0192
0193   template <typename Real_v>
0194   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void SafetyToOut(BooleanStruct const &unplaced,
0195                                                                        Vector3D<Real_v> const &point, Real_v &safety)
0196   {
0197     const auto safetyleft  = unplaced.fLeftVolume->SafetyToOut(point);
0198     const auto safetyright = unplaced.fRightVolume->SafetyToIn(point);
0199     safety                 = Min(safetyleft, safetyright);
0200   }
0201
0202   template <typename Real_v, typename Bool_v>
0203   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void NormalKernel(BooleanStruct const &unplaced,
0204                                                                         Vector3D<Real_v> const &point,
0205                                                                         Vector3D<Real_v> &normal, Bool_v &valid)
0206   {
0207     Vector3D<Real_v> localNorm;
0208     Vector3D<Real_v> localPoint;
0209     valid = false; // Backend::kFalse;
0210
0211     VPlacedVolume const *const fPtrSolidA = unplaced.fLeftVolume;
0212     VPlacedVolume const *const fPtrSolidB = unplaced.fRightVolume;
0213
0214     // If point is inside B, then it must be on a surface of B
0215     if (fPtrSolidB->Contains(point)) {
0216       fPtrSolidB->GetTransformation()->Transform(point, localPoint);
0217       valid = fPtrSolidB->Normal(localPoint, localNorm);
0218       // The normal to the subtracted solid has to be inverted and transformed back
0219       localNorm *= -1.;
0220       fPtrSolidB->GetTransformation()->InverseTransformDirection(localNorm, normal);
0221       return;
0222     }
0223
0224     // If point is outside A, then it must be on a surface of A
0225     if (!fPtrSolidA->Contains(point)) {
0226       fPtrSolidA->GetTransformation()->Transform(point, localPoint);
0227       valid = fPtrSolidA->Normal(localPoint, localNorm);
0228       fPtrSolidA->GetTransformation()->InverseTransformDirection(localNorm, normal);
0229       return;
0230     }
0231
0232     // Point is inside A and outside B, check safety
0233     fPtrSolidA->GetTransformation()->Transform(point, localPoint);
0234     Real_v safetyA = fPtrSolidA->SafetyToOut(localPoint);
0235     Real_v safetyB = fPtrSolidB->SafetyToIn(point);
0236     Bool_v onA     = safetyA < safetyB;
0237     if (vecCore::MaskFull(onA)) {
0238       valid = fPtrSolidA->Normal(localPoint, localNorm);
0239       fPtrSolidA->GetTransformation()->InverseTransformDirection(localNorm, normal);
0240       return;
0241     } else {
0242       //  if (vecCore::MaskEmpty(onA)) {  // to use real mask operation when supporting vectors
0243       fPtrSolidB->GetTransformation()->Transform(point, localPoint);
0244       valid = fPtrSolidB->Normal(localPoint, localNorm);
0245       // The normal to the subtracted solid has to be inverted and transformed back
0246       localNorm *= -1.;
0247       fPtrSolidB->GetTransformation()->InverseTransformDirection(localNorm, normal);
0248       return;
0249     }
0250     // Some particles are on A, some on B. We never arrive here in the scalar case
0251     // If the interface to Normal will support the vector case, we have to write code here.
0252     return;
0253   }
0254
0255 }; // End struct BooleanImplementation
0256
0257 } // namespace VECGEOM_IMPL_NAMESPACE
0258
0259 } // namespace vecgeom
0260
0261 // include stuff for boolean union
0262 #include "BooleanUnionImplementation.h"
0263
0264 // include stuff for boolean intersection
0265 #include "BooleanIntersectionImplementation.h"
0266
0267 #endif /* BOOLEANIMPLEMENTATION_H_ */