volumes/kernel/BooleanUnionImplementation.h

0001 /*
0002  * BooleanImplementation.h
0003  */
0004
0005 #ifndef BOOLEANUNIONIMPLEMENTATION_H_
0006 #define BOOLEANUNIONIMPLEMENTATION_H_
0007
0008 #include "VecGeom/base/Global.h"
0009 #include "VecGeom/base/Vector3D.h"
0010 #include "VecGeom/volumes/BooleanStruct.h"
0011
0012 namespace vecgeom {
0013
0014 inline namespace VECGEOM_IMPL_NAMESPACE {
0015
0016 /**
0017  * partial template specialization for UNION implementation
0018  */
0019 template <>
0020 struct BooleanImplementation<kUnion> {
0021   using PlacedShape_t    = PlacedBooleanVolume<kUnion>;
0022   using UnplacedVolume_t = UnplacedBooleanVolume<kUnion>;
0023   using UnplacedStruct_t = BooleanStruct;
0024
0025   template <typename Real_v, typename Bool_v>
0026   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void Contains(BooleanStruct const &unplaced,
0027                                                                     Vector3D<Real_v> const &point, Bool_v &inside)
0028   {
0029     inside = unplaced.fLeftVolume->Contains(point);
0030     if (vecCore::MaskFull(inside)) return;
0031     inside |= unplaced.fRightVolume->Contains(point);
0032   }
0033
0034   template <typename Real_v, typename Inside_t>
0035   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void Inside(BooleanStruct const &unplaced,
0036                                                                   Vector3D<Real_v> const &point, Inside_t &inside)
0037   {
0038     // now use the Inside functionality of left and right components
0039     // algorithm taken from Geant4 implementation
0040     VPlacedVolume const *const fPtrSolidA = unplaced.fLeftVolume;
0041     VPlacedVolume const *const fPtrSolidB = unplaced.fRightVolume;
0042
0043     const auto positionA = fPtrSolidA->Inside(point);
0044     if (positionA == EInside::kInside) {
0045       inside = EInside::kInside;
0046       return;
0047     }
0048
0049     const auto positionB = fPtrSolidB->Inside(point);
0050     if (positionB == EInside::kInside) {
0051       inside = EInside::kInside;
0052       return;
0053     }
0054
0055     if ((positionA == EInside::kSurface) && (positionB == EInside::kSurface)) {
0056       Vector3D<Precision> normalA, normalB, localPoint, localNorm;
0057       fPtrSolidA->GetTransformation()->Transform(point, localPoint);
0058       fPtrSolidA->Normal(localPoint, localNorm);
0059       fPtrSolidA->GetTransformation()->InverseTransformDirection(localNorm, normalA);
0060
0061       fPtrSolidB->GetTransformation()->Transform(point, localPoint);
0062       fPtrSolidB->Normal(localPoint, localNorm);
0063       fPtrSolidB->GetTransformation()->InverseTransformDirection(localNorm, normalB);
0064
0065       if (normalA.Dot(normalB) < 0)
0066         inside = EInside::kInside; // touching solids -)(-
0067       else
0068         inside = EInside::kSurface; // overlapping solids =))
0069       return;
0070     } else {
0071       if ((positionB == EInside::kSurface) || (positionA == EInside::kSurface)) {
0072         inside = EInside::kSurface;
0073         return;
0074       } else {
0075         inside = EInside::kOutside;
0076         return;
0077       }
0078     }
0079   }
0080
0081   template <typename Real_v>
0082   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void DistanceToIn(BooleanStruct const &unplaced,
0083                                                                         Vector3D<Real_v> const &point,
0084                                                                         Vector3D<Real_v> const &direction,
0085                                                                         Real_v const &stepMax, Real_v &distance)
0086   {
0087     const auto d1 = unplaced.fLeftVolume->DistanceToIn(point, direction, stepMax);
0088     const auto d2 = unplaced.fRightVolume->DistanceToIn(point, direction, stepMax);
0089     distance      = Min(d1, d2);
0090   }
0091
0092   template <typename Real_v>
0093   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void DistanceToOut(BooleanStruct const &unplaced,
0094                                                                          Vector3D<Real_v> const &point,
0095                                                                          Vector3D<Real_v> const &dir,
0096                                                                          Real_v const &stepMax, Real_v &distance)
0097   {
0098     VPlacedVolume const *const ptrSolidA = unplaced.fLeftVolume;
0099     VPlacedVolume const *const ptrSolidB = unplaced.fRightVolume;
0100
0101     Real_v dist = 0.;
0102     Real_v pushdist(kPushTolerance);
0103     // size_t push          = 0;
0104     const auto positionA = ptrSolidA->Inside(point);
0105     Vector3D<Real_v> nextp(point);
0106     bool connectingstep(false);
0107
0108     // reusable kernel as lambda
0109     auto kernel = [&](VPlacedVolume const *A, VPlacedVolume const *B) {
0110       do {
0111         connectingstep    = false;
0112         const auto disTmp = A->PlacedDistanceToOut(nextp, dir);
0113         dist += (disTmp >= 0. && disTmp < kInfLength) ? disTmp : 0;
0114         // give a push
0115         dist += pushdist;
0116         // push++;
0117         nextp = point + dist * dir;
0118         // B could be overlapping with A -- and/or connecting A to another part of A
0119         // if (B->Contains(nextp)) {
0120         if (B->Inside(nextp) != vecgeom::kOutside) {
0121           const auto disTmp = B->PlacedDistanceToOut(nextp, dir);
0122           dist += (disTmp >= 0. && disTmp < kInfLength) ? disTmp : 0;
0123           dist += pushdist;
0124           // push++;
0125           nextp          = point + dist * dir;
0126           connectingstep = true;
0127         }
0128       } while (connectingstep && (A->Inside(nextp) != kOutside));
0129     };
0130
0131     if (positionA != kOutside) { // initially in A
0132       kernel(ptrSolidA, ptrSolidB);
0133     }
0134     // if( positionB != kOutside )
0135     else {
0136       kernel(ptrSolidB, ptrSolidA);
0137     }
0138     // At the end we need to subtract just one push distance, since intermediate distances
0139     // from pushed points are smaller than the real distance with the push value
0140     distance = dist - pushdist;
0141     if (distance < kTolerance && positionA == kOutside && ptrSolidB->Inside(point) == kOutside) distance = -kTolerance;
0142     return;
0143   }
0144
0145   template <typename Real_v>
0146   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void SafetyToIn(BooleanStruct const &unplaced,
0147                                                                       Vector3D<Real_v> const &point, Real_v &safety)
0148   {
0149     VPlacedVolume const *const fPtrSolidA = unplaced.fLeftVolume;
0150     VPlacedVolume const *const fPtrSolidB = unplaced.fRightVolume;
0151     const auto distA                      = fPtrSolidA->SafetyToIn(point);
0152     const auto distB                      = fPtrSolidB->SafetyToIn(point);
0153     safety                                = Min(distA, distB);
0154     // If safety is negative it should not be made 0 (convention)
0155     // vecCore::MaskedAssign(safety, safety < 0.0, 0.0);
0156   }
0157
0158   template <typename Real_v>
0159   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void SafetyToOut(BooleanStruct const &unplaced,
0160                                                                        Vector3D<Real_v> const &point, Real_v &safety)
0161   {
0162
0163     safety                                = -kTolerance; // invalid side
0164     VPlacedVolume const *const fPtrSolidA = unplaced.fLeftVolume;
0165     VPlacedVolume const *const fPtrSolidB = unplaced.fRightVolume;
0166
0167     const auto insideA = fPtrSolidA->Inside(point);
0168     const auto insideB = fPtrSolidB->Inside(point);
0169
0170     // Is point already outside?
0171     if (insideA == kOutside && insideB == kOutside) return;
0172
0173     if (insideA != kOutside && insideB != kOutside) /* in both */
0174     {
0175       safety = Max(fPtrSolidA->SafetyToOut(point),
0176                    fPtrSolidB->SafetyToOut(fPtrSolidB->GetTransformation()->Transform(point)));
0177     } else {
0178       if (insideA == kSurface || insideB == kSurface) return;
0179       /* only contained in B */
0180       if (insideA == kOutside) {
0181         safety = fPtrSolidB->SafetyToOut(fPtrSolidB->GetTransformation()->Transform(point));
0182       } else {
0183         safety = fPtrSolidA->SafetyToOut(point);
0184       }
0185     }
0186   }
0187
0188   template <typename Real_v, typename Bool_v>
0189   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void NormalKernel(BooleanStruct const &unplaced,
0190                                                                         Vector3D<Real_v> const &point,
0191                                                                         Vector3D<Real_v> &normal, Bool_v &valid)
0192   {
0193     Vector3D<Real_v> localNorm;
0194     Vector3D<Real_v> localPoint;
0195     valid = false; // Backend::kFalse;
0196
0197     VPlacedVolume const *const fPtrSolidA = unplaced.fLeftVolume;
0198     VPlacedVolume const *const fPtrSolidB = unplaced.fRightVolume;
0199
0200     // If point is inside A, then it must be on a surface of A (points on the
0201     // intersection between A and B cannot be on surface, or if they are they
0202     // are on a common surface and the normal can be computer for A or B)
0203     if (fPtrSolidA->Contains(point)) {
0204       fPtrSolidA->GetTransformation()->Transform(point, localPoint);
0205       valid = fPtrSolidA->Normal(localPoint, localNorm);
0206       fPtrSolidA->GetTransformation()->InverseTransformDirection(localNorm, normal);
0207       return;
0208     }
0209     // Same for points inside B
0210     if (fPtrSolidB->Contains(point)) {
0211       fPtrSolidB->GetTransformation()->Transform(point, localPoint);
0212       valid = fPtrSolidB->Normal(localPoint, localNorm);
0213       fPtrSolidB->GetTransformation()->InverseTransformDirection(localNorm, normal);
0214       return;
0215     }
0216     // Points outside both A and B can be on any surface. We use the safety.
0217     const auto safetyA = fPtrSolidA->SafetyToIn(point);
0218     const auto safetyB = fPtrSolidB->SafetyToIn(point);
0219     auto onA           = safetyA < safetyB;
0220     if (vecCore::MaskFull(onA)) {
0221       fPtrSolidA->GetTransformation()->Transform(point, localPoint);
0222       valid = fPtrSolidA->Normal(localPoint, localNorm);
0223       fPtrSolidA->GetTransformation()->InverseTransformDirection(localNorm, normal);
0224       return;
0225     } else {
0226       //  if (vecCore::MaskEmpty(onA)) {  // to use real mask operation when supporting vectors
0227       fPtrSolidB->GetTransformation()->Transform(point, localPoint);
0228       valid = fPtrSolidB->Normal(localPoint, localNorm);
0229       fPtrSolidB->GetTransformation()->InverseTransformDirection(localNorm, normal);
0230       return;
0231     }
0232     // Some particles are on A, some on B. We never arrive here in the scalar case
0233     // If the interface to Normal will support the vector case, we have to write code here.
0234     return;
0235   }
0236
0237 }; // End struct BooleanImplementation
0238
0239 } // namespace VECGEOM_IMPL_NAMESPACE
0240
0241 } // namespace vecgeom
0242
0243 #endif /* BooleanImplementation_H_ */