volumes/kernel/SExtruImplementation.h

0001 #ifndef VECGEOM_VOLUMES_KERNEL_SEXTRUIMPLEMENTATION_H_
0002 #define VECGEOM_VOLUMES_KERNEL_SEXTRUIMPLEMENTATION_H_
0003
0004 #include "VecGeom/base/Vector3D.h"
0005 #include "VecGeom/volumes/PolygonalShell.h"
0006 #include "VecGeom/volumes/kernel/BoxImplementation.h"
0007
0008 namespace vecgeom {
0009
0010 VECGEOM_DEVICE_FORWARD_DECLARE(struct SExtruImplementation;);
0011 VECGEOM_DEVICE_DECLARE_CONV(struct, SExtruImplementation);
0012
0013 inline namespace VECGEOM_IMPL_NAMESPACE {
0014
0015 class PlacedSExtru;
0016 class PolygonalShell;
0017 class UnplacedSExtruVolume;
0018
0019 struct SExtruImplementation {
0020
0021   using PlacedShape_t    = PlacedSExtru;
0022   using UnplacedStruct_t = PolygonalShell;
0023   using UnplacedVolume_t = UnplacedSExtruVolume;
0024
0025   template <typename Real_v, typename Bool_v>
0026   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void Contains(UnplacedStruct_t const &unplaced,
0027                                                                     Vector3D<Real_v> const &p, Bool_v &inside)
0028   {
0029     inside    = Bool_v(false);
0030     auto done = p.z() > Real_v(unplaced.fUpperZ);
0031     done |= p.z() < Real_v(unplaced.fLowerZ);
0032     if (vecCore::MaskFull(done)) return;
0033     if (unplaced.fPolygon.IsConvex())
0034       inside = !done && unplaced.fPolygon.ContainsConvex(p);
0035     else
0036       inside = !done && unplaced.fPolygon.Contains(p);
0037   }
0038
0039   template <typename Real_v, typename Inside_t>
0040   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void Inside(UnplacedStruct_t const &unplaced,
0041                                                                   Vector3D<Real_v> const &point, Inside_t &inside)
0042   {
0043     // this is a quick / non-optimized ans scalar only implementation:
0044     if (point.z() > unplaced.fUpperZ + kTolerance) {
0045       inside = vecgeom::kOutside;
0046       return;
0047     }
0048     if (point.z() < unplaced.fLowerZ - kTolerance) {
0049       inside = vecgeom::kOutside;
0050       return;
0051     }
0052
0053     // check conditions for surface first
0054     using Bool_v = vecCore::Mask_v<Real_v>;
0055     Bool_v onZ   = Abs(point.z() - unplaced.fUpperZ) < kTolerance;
0056     onZ |= Abs(point.z() - unplaced.fLowerZ) < kTolerance;
0057
0058     if (unplaced.fPolygon.IsConvex()) {
0059       inside = unplaced.fPolygon.InsideConvex(point);
0060       if (onZ && inside != vecgeom::kOutside) inside = vecgeom::kSurface;
0061       return;
0062     }
0063
0064     if (onZ) {
0065       if (unplaced.fPolygon.Contains(point)) {
0066         inside = vecgeom::kSurface;
0067         return;
0068       }
0069     }
0070
0071     // not on z-surface --> check other surface with safety for moment
0072     if (unplaced.fLowerZ <= point.z() && point.z() <= unplaced.fUpperZ) {
0073       int unused;
0074       auto s = unplaced.fPolygon.SafetySqr(point, unused);
0075       if (s < kTolerance * kTolerance) {
0076         inside = vecgeom::kSurface;
0077         return;
0078       }
0079     }
0080
0081     Bool_v c;
0082     Contains(unplaced, point, c);
0083
0084     if (c)
0085       inside = vecgeom::kInside;
0086     else
0087       inside = vecgeom::kOutside;
0088     return;
0089   }
0090
0091   template <typename Real_v>
0092   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void DistanceToIn(UnplacedStruct_t const &polyshell,
0093                                                                         Vector3D<Real_v> const &p,
0094                                                                         Vector3D<Real_v> const &dir,
0095                                                                         Real_v const & /*stepMax*/, Real_v &distance)
0096   {
0097     if (polyshell.fPolygon.IsConvex()) {
0098       distance = polyshell.DistanceToInConvex(p, dir);
0099       return;
0100     }
0101     distance = Real_v(kInfLength);
0102     // consider adding bounding box check
0103
0104     // check collision with +z or -z
0105     const auto s = vecCore::Blend(dir.z() > Real_v(0.), p.z() - polyshell.fLowerZ, polyshell.fUpperZ - p.z());
0106
0107     const auto canhit = s < Real_v(kTolerance);
0108     if (!vecCore::MaskEmpty(canhit)) {
0109       const auto dist = -s / Abs(dir.z());
0110       // propagate
0111       const auto xInters = p.x() + dist * dir.x();
0112       const auto yInters = p.y() + dist * dir.y();
0113
0114       const auto hits = polyshell.fPolygon.Contains(Vector3D<Real_v>(xInters, yInters, Real_v(0.)));
0115
0116       vecCore::MaskedAssign(distance, hits, dist);
0117       if (vecCore::MaskFull(hits)) {
0118         return;
0119       }
0120     }
0121
0122     // check collision with polyshell
0123     vecCore__MaskedAssignFunc(distance, distance == Real_v(kInfLength), polyshell.DistanceToIn(p, dir));
0124     return;
0125   }
0126
0127   template <typename Real_v>
0128   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void DistanceToOut(UnplacedStruct_t const &polyshell,
0129                                                                          Vector3D<Real_v> const &p,
0130                                                                          Vector3D<Real_v> const &dir,
0131                                                                          Real_v const & /* stepMax */, Real_v &distance)
0132   {
0133     if (polyshell.fPolygon.IsConvex()) {
0134       distance = polyshell.DistanceToOutConvex(p, dir);
0135       return;
0136     }
0137     distance = Real_v(-1.);
0138     // or do a hit check; if not then it has to be the z planes
0139     const auto dshell   = polyshell.DistanceToOut(p, dir);
0140     const auto hitshell = dshell < Real_v(kInfLength);
0141     if (vecCore::MaskFull(hitshell)) {
0142       distance = dshell;
0143       return;
0144     }
0145     const auto correctZ = vecCore::Blend(dir.z() > Real_v(0.), Real_v(polyshell.fUpperZ), Real_v(polyshell.fLowerZ));
0146     distance            = (correctZ - p.z()) / dir.z();
0147     return;
0148   }
0149
0150   template <typename Real_v>
0151   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void SafetyToIn(UnplacedStruct_t const &polyshell,
0152                                                                       Vector3D<Real_v> const &point, Real_v &safety)
0153   {
0154     if (polyshell.fPolygon.IsConvex()) {
0155       Real_v safeZ = vecCore::math::Max(polyshell.fLowerZ - point.z(), point.z() - polyshell.fUpperZ);
0156       safety       = vecCore::math::Max(safeZ, polyshell.fPolygon.SafetyConvex(point, false));
0157       return;
0158     }
0159
0160     Vector3D<Precision> aMin, aMax;
0161     polyshell.Extent(aMin, aMax);
0162
0163     using Bool_v = vecCore::Mask_v<Real_v>;
0164     Bool_v isInExtent;
0165     ABBoxImplementation::ABBoxContainsKernelGeneric(aMin, aMax, point, isInExtent);
0166
0167     // no one is in --> return precise safety to box
0168     if (vecCore::MaskEmpty(isInExtent)) {
0169       const auto ssqr = ABBoxImplementation::ABBoxSafetySqr(aMin, aMax, point);
0170       if (ssqr <= 0.) {
0171         safety = 0.;
0172         return;
0173       }
0174       safety = std::sqrt(ssqr);
0175       return;
0176     }
0177
0178     const auto zSafety1 = polyshell.fLowerZ - point.z();
0179     const auto zSafety2 = polyshell.fUpperZ - point.z();
0180     if (Abs(zSafety1) < kTolerance || Abs(zSafety2) < kTolerance) {
0181       // on the z - entering surface:
0182       // need more careful treatment
0183       bool c;
0184       Contains(polyshell, point, c);
0185       if (c) {
0186         safety = 0.;
0187         return;
0188       }
0189     }
0190     int unused;
0191     safety = std::sqrt(polyshell.fPolygon.SafetySqr(point, unused));
0192   }
0193
0194   template <typename Real_v>
0195   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void SafetyToOut(UnplacedStruct_t const &polyshell,
0196                                                                        Vector3D<Real_v> const &point, Real_v &safety)
0197   {
0198     int unused;
0199     if (polyshell.fPolygon.IsConvex()) {
0200       Real_v safeZ = vecCore::math::Min(point.z() - polyshell.fLowerZ, polyshell.fUpperZ - point.z());
0201       safety       = vecCore::math::Min(safeZ, polyshell.fPolygon.SafetyConvex(point, true));
0202       return;
0203     }
0204     safety = std::sqrt(polyshell.fPolygon.SafetySqr(point, unused));
0205     safety = Min(safety, polyshell.fUpperZ - point.z());
0206     safety = Min(safety, point.z() - polyshell.fLowerZ);
0207   }
0208
0209   template <typename Real_v>
0210   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static Vector3D<Real_v> NormalKernel(
0211       UnplacedStruct_t const &unplaced, Vector3D<Real_v> const &point, typename vecCore::Mask_v<Real_v> &valid)
0212   {
0213     // very rough implementation
0214     // not doing any sort of vector addition for normals on corners etc.
0215     valid = false;
0216     Vector3D<Real_v> normal(0., 0., 0.);
0217
0218     // check conditions for surface first
0219     using Bool_v    = vecCore::Mask_v<Real_v>;
0220     Bool_v onUpperZ = Abs(point.z() - unplaced.fUpperZ) < kTolerance;
0221     Bool_v onLowerZ = Abs(point.z() - unplaced.fLowerZ) < kTolerance;
0222
0223     if (onUpperZ || onLowerZ) {
0224       if (unplaced.fPolygon.Contains(point)) {
0225         valid = true;
0226         if (onUpperZ)
0227           normal = Vector3D<Real_v>(0., 0., 1);
0228         else {
0229           normal = Vector3D<Real_v>(0., 0., -1.);
0230         }
0231         return normal;
0232       }
0233     }
0234
0235     // not on z-surface --> check other surface with safety for moment
0236     if (unplaced.fLowerZ <= point.z() && point.z() <= unplaced.fUpperZ) {
0237       int surfaceindex;
0238       auto s = unplaced.fPolygon.SafetySqr(point, surfaceindex);
0239       normal = Vector3D<Real_v>(-unplaced.fPolygon.fA[surfaceindex], -unplaced.fPolygon.fB[surfaceindex], 0.);
0240       if (s < kTolerance * kTolerance) {
0241         valid = true;
0242       }
0243     }
0244     return normal;
0245   }
0246
0247 }; // End struct SExtruImplementation
0248 } // namespace VECGEOM_IMPL_NAMESPACE
0249 } // namespace vecgeom
0250
0251 #endif