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