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