VecGeom/volumes/PolygonalShell.h

0001 #ifndef VECGEOM_POLYGONAL_SHELL_H
0002 #define VECGEOM_POLYGONAL_SHELL_H
0003
0004 #include "VecGeom/base/Global.h"
0005 #include "VecGeom/volumes/PlanarPolygon.h"
0006
0007 namespace vecgeom {
0008
0009 VECGEOM_DEVICE_FORWARD_DECLARE(class PolygonalShell;);
0010 VECGEOM_DEVICE_DECLARE_CONV(class, PolygonalShell);
0011
0012 inline namespace VECGEOM_IMPL_NAMESPACE {
0013
0014 // a set of z-axis aligned rectangles
0015 // looking from the z - direction the rectangles form a convex or concave polygon
0016 class PolygonalShell : AlignedBase {
0017
0018 private:
0019   // the polygon (with friend access)
0020   PlanarPolygon fPolygon;
0021   Precision fLowerZ; // lower z plane
0022   Precision fUpperZ; // upper z plane
0023
0024   friend class SimpleExtruPolygon;
0025   friend struct SExtruImplementation;
0026   friend class UnplacedSExtruVolume;
0027
0028 public:
0029   VECCORE_ATT_HOST_DEVICE
0030   PolygonalShell() : fPolygon() {}
0031
0032   VECCORE_ATT_HOST_DEVICE
0033   PolygonalShell(int nvertices, Precision *x, Precision *y, Precision lowerz, Precision upperz)
0034       : fPolygon(nvertices, x, y), fLowerZ(lowerz), fUpperZ(upperz)
0035   {
0036   }
0037
0038   VECCORE_ATT_HOST_DEVICE
0039   void Init(int nvertices, Precision *x, Precision *y, Precision lowerz, Precision upperz)
0040   {
0041     fPolygon.Init(nvertices, x, y);
0042     fLowerZ = lowerz;
0043     fUpperZ = upperz;
0044   }
0045
0046   // the area of the shell ( does not include the area of the planar polygon )
0047   VECCORE_ATT_HOST_DEVICE
0048   Precision SurfaceArea() const
0049   {
0050     const auto kS = fPolygon.fVertices.size();
0051     Precision area(0.);
0052     for (size_t i = 0; i < kS; ++i) {
0053       // vertex lengh x (fUpperZ - fLowerZ)
0054       area += fPolygon.fLengthSqr[i];
0055     }
0056     return std::sqrt(area) * (fUpperZ - fLowerZ);
0057   }
0058
0059   VECCORE_ATT_HOST_DEVICE
0060   PlanarPolygon const &GetPolygon() const { return fPolygon; }
0061
0062   VECCORE_ATT_HOST_DEVICE
0063   Precision GetLowerZ() const { return fLowerZ; }
0064
0065   VECCORE_ATT_HOST_DEVICE
0066   Precision GetUpperZ() const { return fUpperZ; }
0067
0068   template <typename Real_v>
0069   VECCORE_ATT_HOST_DEVICE void Extent(Vector3D<Real_v> &aMin, Vector3D<Real_v> &aMax) const
0070   {
0071     aMin[0] = Real_v(fPolygon.GetMinX());
0072     aMin[1] = Real_v(fPolygon.GetMinY());
0073     aMin[2] = Real_v(fLowerZ);
0074
0075     aMax[0] = Real_v(fPolygon.GetMaxX());
0076     aMax[1] = Real_v(fPolygon.GetMaxY());
0077     aMax[2] = Real_v(fUpperZ);
0078   }
0079
0080   template <typename Real_v>
0081   VECCORE_ATT_HOST_DEVICE Real_v DistanceToIn(Vector3D<Real_v> const &point, Vector3D<Real_v> const &dir) const
0082   {
0083     return fPolygon.IsConvex() ? DistanceToInConvex(point, dir) : DistanceToInConcave(point, dir);
0084   }
0085
0086   template <typename Real_v>
0087   VECCORE_ATT_HOST_DEVICE Real_v DistanceToInConvex(Vector3D<Real_v> const &point, Vector3D<Real_v> const &dir) const
0088   {
0089     using Bool_v = vecCore::Mask_v<Real_v>;
0090     Bool_v done(false);
0091     Real_v result(kInfLength);
0092     const auto S = fPolygon.fVertices.size();
0093
0094     for (size_t i = 0; i < S; ++i) { // side/rectangle index
0095       // approaching from right side?
0096       // under the assumption that surface normals points "inwards"
0097       const Real_v proj        = fPolygon.fA[i] * dir.x() + fPolygon.fB[i] * dir.y();
0098       const Bool_v sidecorrect = proj >= -kTolerance;
0099       if (vecCore::MaskEmpty(sidecorrect)) {
0100         continue;
0101       }
0102
0103       // the distance to the plane (specialized for fNormalsZ == 0)
0104       const Real_v pdist = fPolygon.fA[i] * point.x() + fPolygon.fB[i] * point.y() + fPolygon.fD[i];
0105
0106       const Bool_v moving_away = pdist > kTolerance;
0107       if (vecCore::MaskFull(moving_away)) {
0108         continue;
0109       }
0110
0111       const Real_v dist = -pdist / NonZero(proj);
0112
0113       // propagate to plane (first just z)
0114       const Real_v zInters(point.z() + dist * dir.z());
0115       const Bool_v zRangeOk = (zInters <= fUpperZ) && (zInters >= fLowerZ);
0116       if (!vecCore::MaskEmpty(zRangeOk)) {
0117         // check intersection with rest of rectangle
0118         const Real_v xInters(point.x() + dist * dir.x());
0119         const Real_v yInters(point.y() + dist * dir.y());
0120
0121         // we could already check if intersection within the known extent
0122         const Bool_v intersects = fPolygon.OnSegment<Real_v, Precision, Bool_v>(i, xInters, yInters);
0123
0124         vecCore::MaskedAssign(result, !done && intersects, dist);
0125         done |= intersects;
0126       }
0127       if (vecCore::MaskFull(done)) {
0128         return result;
0129       }
0130     }
0131     return result;
0132   }
0133
0134   template <typename Real_v>
0135   VECCORE_ATT_HOST_DEVICE Real_v DistanceToInConcave(Vector3D<Real_v> const &point, Vector3D<Real_v> const &dir) const
0136   {
0137     using Bool_v = vecCore::Mask_v<Real_v>;
0138     Real_v result(kInfLength);
0139     const auto S = fPolygon.fVertices.size();
0140
0141     for (size_t i = 0; i < S; ++i) { // side/rectangle index
0142       // approaching from right side?
0143       // under the assumption that surface normals points "inwards"
0144       const Real_v proj        = fPolygon.fA[i] * dir.x() + fPolygon.fB[i] * dir.y();
0145       const Bool_v sidecorrect = proj >= -kTolerance;
0146       if (vecCore::MaskEmpty(sidecorrect)) {
0147         continue;
0148       }
0149
0150       // the distance to the plane (specialized for fNormalsZ == 0)
0151       const Real_v pdist = fPolygon.fA[i] * point.x() + fPolygon.fB[i] * point.y() + fPolygon.fD[i];
0152
0153       const Bool_v moving_away = pdist > kTolerance;
0154       if (vecCore::MaskFull(moving_away)) {
0155         continue;
0156       }
0157
0158       const Real_v dist = -pdist / NonZero(proj);
0159
0160       // propagate to plane (first just z)
0161       const Real_v zInters(point.z() + dist * dir.z());
0162       const Bool_v zRangeOk = (zInters <= fUpperZ) && (zInters >= fLowerZ);
0163       if (!vecCore::MaskEmpty(zRangeOk)) {
0164         // check intersection with rest of rectangle
0165         const Real_v xInters(point.x() + dist * dir.x());
0166         const Real_v yInters(point.y() + dist * dir.y());
0167
0168         // we could already check if intersection within the known extent
0169         const Bool_v intersects = fPolygon.OnSegment<Real_v, Precision, Bool_v>(i, xInters, yInters);
0170
0171         vecCore__MaskedAssignFunc(result, intersects, Min(dist, result));
0172       }
0173       // if (vecCore::MaskFull(done)) {
0174       //        return result;
0175       //      }
0176     }
0177     return result;
0178   }
0179
0180   // -- DistanceToOut --
0181
0182   template <typename Real_v>
0183   VECCORE_ATT_HOST_DEVICE Real_v DistanceToOut(Vector3D<Real_v> const &point, Vector3D<Real_v> const &dir) const
0184   {
0185     return fPolygon.IsConvex() ? DistanceToOutConvex(point, dir) : DistanceToOutConcave(point, dir);
0186   }
0187
0188   // convex distance to out; checks for hits and aborts loop if hit found
0189   // NOTE: this kernel is the same as DistanceToIn apart from the comparisons for early return
0190   // these could become a template parameter
0191   template <typename Real_v>
0192   VECCORE_ATT_HOST_DEVICE Real_v DistanceToOutConvex(Vector3D<Real_v> const &point, Vector3D<Real_v> const &dir) const
0193   {
0194     using Bool_v = vecCore::Mask_v<Real_v>;
0195     Bool_v done(false);
0196     Real_v result(kInfLength);
0197     const auto S = fPolygon.fVertices.size();
0198
0199     for (size_t i = 0; i < S; ++i) { // side/rectangle index
0200       // approaching from right side?
0201       // under the assumption that surface normals points "inwards"
0202       const Real_v proj        = fPolygon.fA[i] * dir.x() + fPolygon.fB[i] * dir.y();
0203       const Bool_v sidecorrect = proj <= kTolerance;
0204       if (vecCore::MaskEmpty(sidecorrect)) {
0205         continue;
0206       }
0207
0208       // the distance to the plane (specialized for fNormalsZ == 0)
0209       const Real_v pdist = fPolygon.fA[i] * point.x() + fPolygon.fB[i] * point.y() + fPolygon.fD[i];
0210
0211       const Bool_v moving_away = pdist < -kTolerance;
0212       if (vecCore::MaskFull(moving_away)) {
0213         continue;
0214       }
0215
0216       const Real_v dist = -pdist / NonZero(proj);
0217
0218       // propagate to plane (first just z)
0219       const Real_v zInters(point.z() + dist * dir.z());
0220       const Bool_v zRangeOk = (zInters <= fUpperZ) && (zInters >= fLowerZ) && sidecorrect && !moving_away;
0221       if (!vecCore::MaskEmpty(zRangeOk)) {
0222         // check intersection with rest of rectangle
0223         const Real_v xInters(point.x() + dist * dir.x());
0224         const Real_v yInters(point.y() + dist * dir.y());
0225
0226         // we could already check if intersection within the known extent
0227         const Bool_v intersects = fPolygon.OnSegment<Real_v, Precision, Bool_v>(i, xInters, yInters) && zRangeOk &&
0228                                   (dist >= -Real_v(kTolerance));
0229
0230         vecCore::MaskedAssign(result, !done && intersects, dist);
0231         done |= intersects;
0232       }
0233       if (vecCore::MaskFull(done)) {
0234         return result;
0235       }
0236     }
0237     return result;
0238   }
0239
0240   // DistanceToOut for the concave case
0241   // we should ideally combine this with the other kernel
0242   template <typename Real_v>
0243   VECCORE_ATT_HOST_DEVICE Real_v DistanceToOutConcave(Vector3D<Real_v> const &point, Vector3D<Real_v> const &dir) const
0244   {
0245     using Bool_v = vecCore::Mask_v<Real_v>;
0246     Real_v result(kInfLength);
0247     const auto S = fPolygon.fVertices.size();
0248
0249     for (size_t i = 0; i < S; ++i) { // side/rectangle index
0250       // approaching from right side?
0251       // under the assumption that surface normals points "inwards"
0252       const Real_v proj        = fPolygon.fA[i] * dir.x() + fPolygon.fB[i] * dir.y();
0253       const Bool_v sidecorrect = proj <= kTolerance;
0254       if (vecCore::MaskEmpty(sidecorrect)) {
0255         continue;
0256       }
0257
0258       // the distance to the plane (specialized for fNormalsZ == 0)
0259       const Real_v pdist = fPolygon.fA[i] * point.x() + fPolygon.fB[i] * point.y() + fPolygon.fD[i];
0260
0261       const Bool_v moving_away = pdist < -kTolerance;
0262       if (vecCore::MaskFull(moving_away)) {
0263         continue;
0264       }
0265
0266       const Real_v dist = -pdist / NonZero(proj);
0267
0268       // propagate to plane (first just z)
0269       const Real_v zInters(point.z() + dist * dir.z());
0270       const Bool_v zRangeOk = (zInters <= fUpperZ) && (zInters >= fLowerZ) && sidecorrect && !moving_away;
0271       if (!vecCore::MaskEmpty(zRangeOk)) {
0272         // check intersection with rest of rectangle
0273         const Real_v xInters(point.x() + dist * dir.x());
0274         const Real_v yInters(point.y() + dist * dir.y());
0275
0276         // we could already check if intersection within the known extent
0277         const Bool_v intersects = fPolygon.OnSegment<Real_v, Precision, Bool_v>(i, xInters, yInters) && zRangeOk &&
0278                                   (dist >= -Real_v(kTolerance));
0279
0280         vecCore__MaskedAssignFunc(result, intersects, Min(dist, result));
0281         // done |= intersects;
0282       }
0283       // if (vecCore::MaskFull(done)) {
0284       //  return result;
0285       //}
0286     }
0287     return result;
0288   }
0289
0290 }; // end class
0291
0292 #define SPECIALIZATION
0293 #ifdef SPECIALIZATION
0294 // template specialization for Distance functions
0295 template <>
0296 VECCORE_ATT_HOST_DEVICE inline Precision PolygonalShell::DistanceToOutConvex(Vector3D<Precision> const &point,
0297                                                                              Vector3D<Precision> const &dir) const
0298 {
0299   Precision dz         = 0.5 * (fUpperZ - fLowerZ);
0300   Precision pz         = point.z() - 0.5 * (fLowerZ + fUpperZ);
0301   const Precision safz = vecCore::math::Abs(pz) - dz;
0302   if (safz > kTolerance) return -kTolerance;
0303
0304   Precision vz   = dir.z();
0305   Precision tmax = (vecCore::math::CopySign(dz, vz) - pz) / NonZero(vz);
0306   const auto S   = fPolygon.fVertices.size();
0307   for (size_t i = 0; i < S; ++i) { // side/rectangle index
0308
0309     const Precision proj = -(fPolygon.fA[i] * dir.x() + fPolygon.fB[i] * dir.y());
0310     // normals pointing inwards
0311     const Precision pdist = -(fPolygon.fA[i] * point.x() + fPolygon.fB[i] * point.y() + fPolygon.fD[i]);
0312     if (pdist > kTolerance) return -kTolerance;
0313     if (proj > 0) {
0314       const Precision dist = -pdist / NonZero(proj);
0315       if (tmax > dist) tmax = dist;
0316     }
0317   }
0318   return tmax;
0319 }
0320
0321 // template specialization for Distance functions
0322 template <>
0323 VECCORE_ATT_HOST_DEVICE inline Precision PolygonalShell::DistanceToInConvex(Vector3D<Precision> const &point,
0324                                                                             Vector3D<Precision> const &dir) const
0325 {
0326   Precision dz = 0.5 * (fUpperZ - fLowerZ);
0327   Precision pz = point.z() - 0.5 * (fLowerZ + fUpperZ);
0328   if ((vecCore::math::Abs(pz) - dz) > -kTolerance && pz * dir.z() >= 0) return kInfLength;
0329   const Precision invz = -1. / NonZero(dir.z());
0330   const Precision ddz  = (invz < 0) ? dz : -dz;
0331   Precision tmin       = (pz + ddz) * invz;
0332   Precision tmax       = (pz - ddz) * invz;
0333   const auto S         = fPolygon.fVertices.size();
0334   for (size_t i = 0; i < S; ++i) { // side/rectangle index
0335
0336     const Precision proj = -(fPolygon.fA[i] * dir.x() + fPolygon.fB[i] * dir.y());
0337     // normals pointing inwards
0338     const bool moving_away = proj > -kTolerance;
0339     // the distance to the plane (specialized for fNormalsZ == 0)
0340     const Precision pdist   = -(fPolygon.fA[i] * point.x() + fPolygon.fB[i] * point.y() + fPolygon.fD[i]);
0341     const bool side_correct = pdist > -kTolerance;
0342     if (side_correct) {
0343       if (moving_away) return kInfLength;
0344       const Precision dist = -pdist / NonZero(proj);
0345       if (dist > tmin) tmin = dist;
0346     } else if (moving_away) {
0347       const Precision dist = -pdist / NonZero(proj);
0348       if (dist < tmax) tmax = dist;
0349     }
0350   }
0351   if (tmax < tmin + kTolerance) return kInfLength;
0352   return tmin;
0353 }
0354
0355 #endif
0356
0357 } // namespace VECGEOM_IMPL_NAMESPACE
0358 } // namespace vecgeom
0359
0360 #endif