VecGeom/volumes/Planes.h

0001 /// \file Planes.h
0002 /// \author Johannes de Fine Licht (johannes.definelicht@cern.ch)
0003
0004 #ifndef VECGEOM_VOLUMES_PLANES_H_
0005 #define VECGEOM_VOLUMES_PLANES_H_
0006
0007 #include "VecGeom/base/Global.h"
0008
0009 #include "VecGeom/base/AlignedBase.h"
0010 #include "VecGeom/base/Array.h"
0011 #include "VecGeom/base/SOA3D.h"
0012 #include "VecGeom/volumes/kernel/GenericKernels.h"
0013
0014 #if defined(VECGEOM_VC) && defined(VECGEOM_QUADRILATERALS_VC)
0015 #include <Vc/Vc>
0016 typedef Vc::Vector<vecgeom::Precision> VcPrecision;
0017 typedef Vc::Vector<vecgeom::Precision>::Mask VcBool;
0018 constexpr int kVectorSize = VcPrecision::Size;
0019 #endif
0020
0021 namespace vecgeom {
0022
0023 VECGEOM_DEVICE_FORWARD_DECLARE(class Planes;);
0024 VECGEOM_DEVICE_DECLARE_CONV(class, Planes);
0025
0026 inline namespace VECGEOM_IMPL_NAMESPACE {
0027
0028 class Planes : public AlignedBase {
0029
0030 private:
0031   SOA3D<Precision> fNormals;   ///< Normalized normals of the planes.
0032   Array<Precision> fDistances; ///< Distance from plane to origin (0, 0, 0).
0033   bool fConvex{true};          ///< Convexity of the planes array (drives the inside reduction)
0034
0035 public:
0036   VECCORE_ATT_HOST_DEVICE
0037   Planes(int size, bool convex = true);
0038
0039   VECCORE_ATT_HOST_DEVICE
0040   ~Planes();
0041
0042   VECCORE_ATT_HOST_DEVICE
0043   Planes &operator=(Planes const &rhs);
0044
0045   VECCORE_ATT_HOST_DEVICE
0046   VECGEOM_FORCE_INLINE
0047   Precision const *operator[](int index) const;
0048
0049   VECCORE_ATT_HOST_DEVICE
0050   VECGEOM_FORCE_INLINE
0051   int size() const;
0052
0053   VECCORE_ATT_HOST_DEVICE
0054   VECGEOM_FORCE_INLINE
0055   void reserve(size_t newsize)
0056   {
0057     fNormals.reserve(newsize);
0058     fDistances.Allocate(newsize);
0059   }
0060
0061   VECCORE_ATT_HOST_DEVICE
0062   VECGEOM_FORCE_INLINE
0063   SOA3D<Precision> const &GetNormals() const;
0064
0065   VECCORE_ATT_HOST_DEVICE
0066   VECGEOM_FORCE_INLINE
0067   Vector3D<Precision> GetNormal(int i) const;
0068
0069   VECCORE_ATT_HOST_DEVICE
0070   VECGEOM_FORCE_INLINE
0071   Array<Precision> const &GetDistances() const;
0072
0073   VECCORE_ATT_HOST_DEVICE
0074   VECGEOM_FORCE_INLINE
0075   Precision GetDistance(int i) const;
0076
0077   VECCORE_ATT_HOST_DEVICE
0078   VECGEOM_FORCE_INLINE
0079   bool IsConvex() const { return fConvex; }
0080
0081   VECCORE_ATT_HOST_DEVICE
0082   void Set(int index, Vector3D<Precision> const &normal, Vector3D<Precision> const &origin);
0083
0084   VECCORE_ATT_HOST_DEVICE
0085   void Set(int index, Vector3D<Precision> const &normal, Precision distance);
0086
0087   /// Flip the sign of the normal and distance at the specified index
0088   VECCORE_ATT_HOST_DEVICE
0089   void FlipSign(int index);
0090
0091   template <typename Real_v>
0092   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE vecCore::Mask_v<Real_v> Contains(Vector3D<Real_v> const &point) const;
0093
0094   template <typename Real_v, typename Inside_v>
0095   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE Inside_v Inside(Vector3D<Real_v> const &point) const;
0096
0097   template <typename Real_v, typename Inside_v>
0098   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE Inside_v Inside(Vector3D<Real_v> const &point, int i) const;
0099
0100   template <typename Real_v>
0101   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE Real_v Distance(Vector3D<Real_v> const &point,
0102                                                                Vector3D<Real_v> const &direction) const;
0103
0104   template <bool pointInsideT, typename Real_v>
0105   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE Real_v Distance(Vector3D<Real_v> const &point) const;
0106 };
0107
0108 std::ostream &operator<<(std::ostream &os, Planes const &planes);
0109
0110 VECCORE_ATT_HOST_DEVICE
0111 Precision const *Planes::operator[](int i) const
0112 {
0113   if (i == 0) return fNormals.x();
0114   if (i == 1) return fNormals.y();
0115   if (i == 2) return fNormals.z();
0116   if (i == 3) return &fDistances[0];
0117   return NULL;
0118 }
0119
0120 VECCORE_ATT_HOST_DEVICE
0121 int Planes::size() const
0122 {
0123   return fNormals.size();
0124 }
0125
0126 VECCORE_ATT_HOST_DEVICE
0127 SOA3D<Precision> const &Planes::GetNormals() const
0128 {
0129   return fNormals;
0130 }
0131
0132 VECCORE_ATT_HOST_DEVICE
0133 Vector3D<Precision> Planes::GetNormal(int i) const
0134 {
0135   return fNormals[i];
0136 }
0137
0138 VECCORE_ATT_HOST_DEVICE
0139 Array<Precision> const &Planes::GetDistances() const
0140 {
0141   return fDistances;
0142 }
0143
0144 VECCORE_ATT_HOST_DEVICE
0145 Precision Planes::GetDistance(int i) const
0146 {
0147   return fDistances[i];
0148 }
0149
0150 namespace {
0151
0152 template <typename Real_v, bool = true>
0153 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE void AcceleratedContains(int & i, const int,
0154                                                                       SOA3D<Precision> const &,
0155                                                                       Array<Precision> const &,
0156                                                                       Vector3D<Real_v> const &,
0157                                                                       vecCore::Mask_v<Real_v> &)
0158 {
0159   return;
0160 }
0161
0162 #if defined(VECGEOM_VC) && defined(VECGEOM_QUADRILATERALS_VC)
0163 template <>
0164 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE void AcceleratedContains<Precision, true>(
0165     int &i, const int n, SOA3D<Precision> const &normals, Array<Precision> const &distances,
0166     Vector3D<Precision> const &point, vecCore::Mask_v<double> &result)
0167 {
0168   for (; i < n - kVectorSize; i += kVectorSize) {
0169     VcBool inside = VcPrecision(normals.x() + i) * point[0] + VcPrecision(normals.y() + i) * point[1] +
0170                         VcPrecision(normals.z() + i) * point[2] + VcPrecision(&distances[0] + i) <
0171                     0;
0172     // Early return if not inside all planes (convex case)
0173     result = vecCore::MaskFull(inside);
0174     if (!result) {
0175       i = n;
0176       break;
0177     }
0178   }
0179 }
0180
0181 template <>
0182 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE void AcceleratedContains<Precision, false>(
0183     int &i, const int n, SOA3D<Precision> const &normals, Array<Precision> const &distances,
0184     Vector3D<Precision> const &point, vecCore::Mask_v<double> &result)
0185 {
0186   for (; i < n - kVectorSize; i += kVectorSize) {
0187     VcBool inside = VcPrecision(normals.x() + i) * point[0] + VcPrecision(normals.y() + i) * point[1] +
0188                         VcPrecision(normals.z() + i) * point[2] + VcPrecision(&distances[0] + i) <
0189                     0;
0190     // Early return ifinside any planes (non-convex case)
0191     result = !vecCore::MaskEmpty(inside);
0192     if (result) {
0193       i = n;
0194       break;
0195     }
0196   }
0197 }
0198 #endif
0199
0200 } // End anonymous namespace
0201
0202 template <typename Real_v>
0203 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE vecCore::Mask_v<Real_v> Planes::Contains(
0204     Vector3D<Real_v> const &point) const
0205 {
0206
0207   vecCore::Mask_v<Real_v> result = vecCore::Mask_v<Real_v>(true);
0208
0209   int i       = 0;
0210   const int n = size();
0211   if (fConvex) {
0212     AcceleratedContains<Real_v, true>(i, n, fNormals, fDistances, point, result);
0213     for (; i < n; ++i) {
0214       result &= point.Dot(fNormals[i]) + fDistances[i] <= 0;
0215     }
0216   } else {
0217     AcceleratedContains<Real_v, false>(i, n, fNormals, fDistances, point, result);
0218     for (; i < n; ++i) {
0219       result |= point.Dot(fNormals[i]) + fDistances[i] <= 0;
0220     }
0221   }
0222
0223   return result;
0224 }
0225
0226 namespace {
0227
0228 template <typename Real_v, typename Inside_v, bool = true>
0229 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE void AcceleratedInside(int & /*i*/, const int /*n*/,
0230                                                                     SOA3D<Precision> const & /*normals*/,
0231                                                                     Array<Precision> const & /*distances*/,
0232                                                                     Vector3D<Real_v> const & /*point*/,
0233                                                                     Inside_v & /*result*/)
0234 {
0235   return;
0236 }
0237
0238 #if defined(VECGEOM_VC) and defined(VECGEOM_QUADRILATERALS_VC)
0239 template <>
0240 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE void AcceleratedInside<Precision, Inside_t, true>(
0241     int &i, const int n, SOA3D<Precision> const &normals, Array<Precision> const &distances,
0242     Vector3D<Precision> const &point, Inside_t &result)
0243 {
0244   for (; i < n - kVectorSize; i += kVectorSize) {
0245     VcPrecision distance = VcPrecision(normals.x() + i) * point[0] + VcPrecision(normals.y() + i) * point[1] +
0246                            VcPrecision(normals.z() + i) * point[2] + VcPrecision(&distances[0] + i);
0247     // If point is outside tolerance of any plane, it is safe to return
0248     if (!vecCore::MaskEmpty(distance > kTolerance)) {
0249       result = EInside::kOutside;
0250       i      = n;
0251       break;
0252     }
0253     // If point is inside tolerance of all planes, keep looking
0254     if (vecCore::MaskFull(distance < -kTolerance)) continue;
0255     // Otherwise point must be on a surface, but could still be outside
0256     result = EInside::kSurface;
0257   }
0258 }
0259
0260 template <>
0261 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE void AcceleratedInside<Precision, Inside_t, false>(
0262     int &i, const int n, SOA3D<Precision> const &normals, Array<Precision> const &distances,
0263     Vector3D<Precision> const &point, Inside_t &result)
0264 {
0265   for (; i < n - kVectorSize; i += kVectorSize) {
0266     VcPrecision distance = VcPrecision(normals.x() + i) * point[0] + VcPrecision(normals.y() + i) * point[1] +
0267                            VcPrecision(normals.z() + i) * point[2] + VcPrecision(&distances[0] + i);
0268     // If point is inside tolerance of any plane, it is safe to return
0269     if (!vecCore::MaskEmpty(distance < -kTolerance)) {
0270       result = EInside::kInside;
0271       i      = n;
0272       break;
0273     }
0274     // If point is outside tolerance of all planes, keep looking
0275     if (vecCore::MaskFull(distance > kTolerance)) continue;
0276     // Otherwise point must be on a surface, but could still be outside
0277     result = EInside::kSurface;
0278   }
0279 }
0280 #endif
0281
0282 } // End anonymous namespace
0283
0284 template <typename Real_v, typename Inside_v>
0285 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE Inside_v Planes::Inside(Vector3D<Real_v> const &point) const
0286 {
0287   Inside_v result = fConvex ? Inside_v(EInside::kInside) : Inside_v(EInside::kOutside);
0288
0289   int i       = 0;
0290   const int n = size();
0291   if (fConvex) {
0292     AcceleratedInside<Real_v, Inside_v, true>(i, n, fNormals, fDistances, point, result);
0293     for (; i < n; ++i) {
0294       Real_v distanceResult =
0295           fNormals.x(i) * point[0] + fNormals.y(i) * point[1] + fNormals.z(i) * point[2] + fDistances[i];
0296       vecCore::MaskedAssign(result, distanceResult > Real_v(kTolerance), EInside::kOutside);
0297       vecCore::MaskedAssign(result, result == Inside_v(EInside::kInside) && distanceResult > Real_v(-kTolerance),
0298                             Inside_v(EInside::kSurface));
0299       if (vecCore::MaskFull(result == Inside_v(EInside::kOutside))) break;
0300     }
0301   } else {
0302     AcceleratedInside<Real_v, Inside_v, false>(i, n, fNormals, fDistances, point, result);
0303     for (; i < n; ++i) {
0304       Real_v distanceResult =
0305           fNormals.x(i) * point[0] + fNormals.y(i) * point[1] + fNormals.z(i) * point[2] + fDistances[i];
0306       vecCore::MaskedAssign(result, distanceResult < Real_v(-kTolerance), EInside::kInside);
0307       vecCore::MaskedAssign(result, result == Inside_v(EInside::kOutside) && distanceResult < Real_v(kTolerance),
0308                             Inside_v(EInside::kSurface));
0309       if (vecCore::MaskFull(result == Inside_v(EInside::kInside))) break;
0310     }
0311   }
0312
0313   return result;
0314 }
0315
0316 template <typename Real_v, typename Inside_v>
0317 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE Inside_v Planes::Inside(Vector3D<Real_v> const &point, int i) const
0318 {
0319
0320   Inside_v result = Inside_v(EInside::kInside);
0321   Real_v distanceResult =
0322       fNormals.x(i) * point[0] + fNormals.y(i) * point[1] + fNormals.z(i) * point[2] + fDistances[i];
0323   vecCore::MaskedAssign(result, distanceResult > Real_v(kTolerance), Inside_v(EInside::kOutside));
0324   vecCore::MaskedAssign(result, result == Inside_v(EInside::kInside) && distanceResult > Real_v(-kTolerance),
0325                         Inside_v(EInside::kSurface));
0326
0327   return result;
0328 }
0329
0330 template <typename Real_v>
0331 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE Real_v Planes::Distance(Vector3D<Real_v> const &point,
0332                                                                      Vector3D<Real_v> const &direction) const
0333 {
0334
0335   Real_v bestDistance = InfinityLength<Real_v>();
0336   for (int i = 0, iMax = size(); i < iMax; ++i) {
0337     Vector3D<Precision> normal = fNormals[i];
0338     Real_v distance            = -(point.Dot(normal) + fDistances[i]) / direction.Dot(normal);
0339     vecCore::MaskedAssign(bestDistance, distance >= Real_v(0) && distance < bestDistance, distance);
0340   }
0341
0342   return bestDistance;
0343 }
0344
0345 template <bool pointInsideT, typename Real_v>
0346 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE Real_v Planes::Distance(Vector3D<Real_v> const &point) const
0347 {
0348
0349   Real_v bestDistance = InfinityLength<Real_v>();
0350   for (int i = 0, iMax = size(); i < iMax; ++i) {
0351     Real_v distance = Flip<!pointInsideT>::FlipSign(point.Dot(fNormals[i]) + fDistances[i]);
0352     vecCore::MaskedAssign(bestDistance, distance >= Real_v(0) && distance < bestDistance, distance);
0353   }
0354   return bestDistance;
0355 }
0356
0357 } // namespace VECGEOM_IMPL_NAMESPACE
0358
0359 } // namespace vecgeom
0360
0361 #endif // VECGEOM_VOLUMES_PLANES_H_