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   VECCORE_ATT_HOST_DEVICE
0046   static void PrintType()
0047   { /*printf("SpecializedBooleanVolume<%i, %i, %i>", kSubtraction, transCodeT, rotCodeT);*/
0048   }
0049
0050   template <typename Stream>
0051   static void PrintType(Stream &s)
0052   {
0053     /*s << "SpecializedBooleanVolume<kSubtraction"
0054       << "," << transCodeT << "," << rotCodeT << ">";
0055     */
0056   }
0057
0058   template <typename Stream>
0059   static void PrintType(Stream &st, int transCodeT = translation::kGeneric, int rotCodeT = rotation::kGeneric)
0060   {
0061     st << "SpecializedBooleanVolume<kSubtraction" << transCodeT << "," << rotCodeT << ">";
0062   }
0063
0064   template <typename Stream>
0065   static void PrintImplementationType(Stream &s)
0066   {
0067     /*
0068     s << "BooleanImplementation<kSubtraction"
0069       << "," << transCodeT << "," << rotCodeT << ">";
0070     */
0071   }
0072
0073   template <typename Stream>
0074   static void PrintUnplacedType(Stream &s)
0075   {
0076     s << "BooleanStruct";
0077   }
0078
0079   template <typename Real_v, typename Bool_v>
0080   VECGEOM_FORCE_INLINE
0081   VECCORE_ATT_HOST_DEVICE
0082   static void Contains(BooleanStruct const &unplaced, Vector3D<Real_v> const &point, Bool_v &inside)
0083   {
0084     Vector3D<Real_v> tmp;
0085     inside = unplaced.fLeftVolume->Contains(point);
0086     if (vecCore::MaskEmpty(inside)) return;
0087
0088     auto rightInside = unplaced.fRightVolume->Contains(point);
0089     inside &= !rightInside;
0090   }
0091
0092   template <typename Real_v, typename Inside_t>
0093   VECGEOM_FORCE_INLINE
0094   VECCORE_ATT_HOST_DEVICE
0095   static void Inside(BooleanStruct const &unplaced, Vector3D<Real_v> const &p, Inside_t &inside)
0096   {
0097
0098     // now use the Inside functionality of left and right components
0099     // algorithm taken from Geant4 implementation
0100     VPlacedVolume const *const fPtrSolidA = unplaced.fLeftVolume;
0101     VPlacedVolume const *const fPtrSolidB = unplaced.fRightVolume;
0102
0103     const auto positionA = fPtrSolidA->Inside(p);
0104     if (positionA == EInside::kOutside) {
0105       inside = EInside::kOutside;
0106       return;
0107     }
0108
0109     const auto positionB = fPtrSolidB->Inside(p);
0110
0111     if (positionA == EInside::kInside && positionB == EInside::kOutside) {
0112       inside = EInside::kInside;
0113       return;
0114     } else {
0115       if ((positionA == EInside::kInside && positionB == EInside::kSurface) ||
0116           (positionB == EInside::kOutside && positionA == EInside::kSurface)
0117           /*
0118            ||( positionA == EInside::kSurface && positionB == EInside::kSurface &&
0119              (   fPtrSolidA->Normal(p) -
0120                fPtrSolidB->Normal(p) ).mag2() >
0121              1000.0*G4GeometryTolerance::GetInstance()->GetRadialTolerance() ) )
0122           */) {
0123         inside = EInside::kSurface;
0124         return;
0125       } else {
0126         inside = EInside::kOutside;
0127         return;
0128       }
0129     }
0130     // going to be a bit more complicated due to Surface states
0131   }
0132
0133   template <typename Real_v>
0134   VECGEOM_FORCE_INLINE
0135   VECCORE_ATT_HOST_DEVICE
0136   static void DistanceToIn(BooleanStruct const &unplaced, Vector3D<Real_v> const &p, Vector3D<Real_v> const &dir,
0137                            Real_v const &stepMax, Real_v &distance)
0138   {
0139
0140     // TOBEDONE: ASK Andrei about the while loop
0141     // Compute distance from a given point outside to the shape.
0142     Real_v d1, d2, snxt = 0.;
0143     Vector3D<Real_v> hitpoint = p;
0144     // check if inside '-'
0145     auto insideRight = unplaced.fRightVolume->Contains(p);
0146     //  // epsilon is used to push across boundaries
0147     Precision epsil(kPushTolerance);
0148     //
0149     //  // we should never subtract a volume such that B - A > 0
0150     //
0151     while (1) {
0152       if (insideRight) {
0153         //    // propagate to outside of '- / RightShape'
0154         d1 = unplaced.fRightVolume->PlacedDistanceToOut(hitpoint, dir, stepMax);
0155         snxt += (d1 >= 0. && d1 < kInfLength) ? (d1 + epsil) : 0.;
0156         hitpoint += (d1 >= 0. && d1 < kInfLength) ? (d1 + epsil) * dir : 0. * dir;
0157
0158         // now master outside 'B'; check if inside 'A'
0159         //    Bool_t insideLeft =
0160         if (unplaced.fLeftVolume->Contains(hitpoint)) {
0161           auto check = unplaced.fLeftVolume->PlacedDistanceToOut(hitpoint, dir);
0162           if (check > epsil) {
0163             distance = snxt;
0164             //  std::cerr << "hitting  " << distance << "\n";
0165             return;
0166           }
0167         }
0168       }
0169
0170       // if outside of both we do a max operation
0171       // master outside '-' and outside '+' ;  find distances to both
0172       //        fLeftMat->MasterToLocal(&master[0], &local[0]);
0173       d2 = unplaced.fLeftVolume->DistanceToIn(hitpoint, dir, stepMax);
0174       d2 = Max(d2, Precision(0.));
0175       if (d2 == kInfLength) {
0176         distance = kInfLength;
0177         // std::cerr << "missing A " << d2 << "\n";
0178         return;
0179       }
0180
0181       d1 = unplaced.fRightVolume->DistanceToIn(hitpoint, dir, stepMax);
0182       if (d2 < d1 - kTolerance) {
0183         snxt += d2 + epsil;
0184         // std::cerr << "returning  " << snxt << "\n";
0185         distance = snxt;
0186         return;
0187       }
0188
0189       //        // propagate to '-'
0190       snxt += (d1 >= 0. && d1 < kInfLength) ? d1 + epsil : 0.;
0191       hitpoint += (d1 >= 0. && d1 < kInfLength) ? (d1 + epsil) * dir : epsil * dir;
0192       insideRight = true;
0193     } // end while
0194   }
0195
0196   template <typename Real_v>
0197   VECGEOM_FORCE_INLINE
0198   VECCORE_ATT_HOST_DEVICE
0199   static void DistanceToOut(BooleanStruct const &unplaced, Vector3D<Real_v> const &point,
0200                             Vector3D<Real_v> const &direction, Real_v const &stepMax, Real_v &distance)
0201   {
0202     const auto distancel  = unplaced.fLeftVolume->PlacedDistanceToOut(point, direction, stepMax);
0203     const Real_v dinright = unplaced.fRightVolume->DistanceToIn(point, direction, stepMax);
0204     distance              = Min(distancel, dinright);
0205     return;
0206   }
0207
0208   template <typename Real_v>
0209   VECGEOM_FORCE_INLINE
0210   VECCORE_ATT_HOST_DEVICE
0211   static void SafetyToIn(BooleanStruct const &unplaced, Vector3D<Real_v> const &point, Real_v &safety)
0212   {
0213     VPlacedVolume const *const fPtrSolidA = unplaced.fLeftVolume;
0214     VPlacedVolume const *const fPtrSolidB = unplaced.fRightVolume;
0215
0216     // very approximate
0217     if ((fPtrSolidA->Contains(point)) && // case 1
0218         (fPtrSolidB->Contains(point))) {
0219       safety = fPtrSolidB->SafetyToOut(fPtrSolidB->GetTransformation()->Transform(point));
0220     } else {
0221       // po
0222       safety = fPtrSolidA->SafetyToIn(point);
0223     }
0224   }
0225
0226   template <typename Real_v>
0227   VECGEOM_FORCE_INLINE
0228   VECCORE_ATT_HOST_DEVICE
0229   static void SafetyToOut(BooleanStruct const &unplaced, Vector3D<Real_v> const &point, Real_v &safety)
0230   {
0231     const auto safetyleft  = unplaced.fLeftVolume->SafetyToOut(point);
0232     const auto safetyright = unplaced.fRightVolume->SafetyToIn(point);
0233     safety                 = Min(safetyleft, safetyright);
0234   }
0235
0236   template <typename Real_v, typename Bool_v>
0237   VECGEOM_FORCE_INLINE
0238   VECCORE_ATT_HOST_DEVICE
0239   static void NormalKernel(BooleanStruct const &unplaced, Vector3D<Real_v> const &point, Vector3D<Real_v> &normal,
0240                            Bool_v &valid)
0241   {
0242     Vector3D<Real_v> localNorm;
0243     Vector3D<Real_v> localPoint;
0244     valid = false; // Backend::kFalse;
0245
0246     VPlacedVolume const *const fPtrSolidA = unplaced.fLeftVolume;
0247     VPlacedVolume const *const fPtrSolidB = unplaced.fRightVolume;
0248
0249     // If point is inside B, then it must be on a surface of B
0250     if (fPtrSolidB->Contains(point)) {
0251       fPtrSolidB->GetTransformation()->Transform(point, localPoint);
0252       valid = fPtrSolidB->Normal(localPoint, localNorm);
0253       // The normal to the subtracted solid has to be inverted and transformed back
0254       localNorm *= -1.;
0255       fPtrSolidB->GetTransformation()->InverseTransformDirection(localNorm, normal);
0256       return;
0257     }
0258
0259     // If point is outside A, then it must be on a surface of A
0260     if (!fPtrSolidA->Contains(point)) {
0261       fPtrSolidA->GetTransformation()->Transform(point, localPoint);
0262       valid = fPtrSolidA->Normal(localPoint, localNorm);
0263       fPtrSolidA->GetTransformation()->InverseTransformDirection(localNorm, normal);
0264       return;
0265     }
0266
0267     // Point is inside A and outside B, check safety
0268     fPtrSolidA->GetTransformation()->Transform(point, localPoint);
0269     Real_v safetyA = fPtrSolidA->SafetyToOut(localPoint);
0270     Real_v safetyB = fPtrSolidB->SafetyToIn(point);
0271     Bool_v onA     = safetyA < safetyB;
0272     if (vecCore::MaskFull(onA)) {
0273       valid = fPtrSolidA->Normal(localPoint, localNorm);
0274       fPtrSolidA->GetTransformation()->InverseTransformDirection(localNorm, normal);
0275       return;
0276     } else {
0277       //  if (vecCore::MaskEmpty(onA)) {  // to use real mask operation when supporting vectors
0278       fPtrSolidB->GetTransformation()->Transform(point, localPoint);
0279       valid = fPtrSolidB->Normal(localPoint, localNorm);
0280       // The normal to the subtracted solid has to be inverted and transformed back
0281       localNorm *= -1.;
0282       fPtrSolidB->GetTransformation()->InverseTransformDirection(localNorm, normal);
0283       return;
0284     }
0285     // Some particles are on A, some on B. We never arrive here in the scalar case
0286     // If the interface to Normal will support the vector case, we have to write code here.
0287     return;
0288   }
0289
0290 }; // End struct BooleanImplementation
0291
0292 } // namespace VECGEOM_IMPL_NAMESPACE
0293
0294 } // namespace vecgeom
0295
0296 // include stuff for boolean union
0297 #include "BooleanUnionImplementation.h"
0298
0299 // include stuff for boolean intersection
0300 #include "BooleanIntersectionImplementation.h"
0301
0302 #endif /* BOOLEANIMPLEMENTATION_H_ */