volumes/kernel/ParallelepipedImplementation.h

0001 // This file is part of VecGeom and is distributed under the
0002 // conditions in the file LICENSE.txt in the top directory.
0003 // For the full list of authors see CONTRIBUTORS.txt and `git log`.
0004
0005 /// This file implements the algorithms for Paralleliped
0006 /// @file volumes/kernel/ParallelepipedImplementation.h
0007 /// @author First version by Johannes de Fine Licht
0008 /// @author Revised by Evgueni Tcherniaev
0009
0010 #ifndef VECGEOM_VOLUMES_KERNEL_PARALLELEPIPEDIMPLEMENTATION_H_
0011 #define VECGEOM_VOLUMES_KERNEL_PARALLELEPIPEDIMPLEMENTATION_H_
0012
0013 #include "VecGeom/base/Vector3D.h"
0014 #include "VecGeom/volumes/ParallelepipedStruct.h"
0015 #include "VecGeom/volumes/kernel/GenericKernels.h"
0016 #include <VecCore/VecCore>
0017
0018 #include <cstdio>
0019
0020 namespace vecgeom {
0021
0022 VECGEOM_DEVICE_FORWARD_DECLARE(struct ParallelepipedImplementation;);
0023 VECGEOM_DEVICE_DECLARE_CONV(struct, ParallelepipedImplementation);
0024
0025 inline namespace VECGEOM_IMPL_NAMESPACE {
0026
0027 class PlacedParallelepiped;
0028 template <typename T>
0029 struct ParallelepipedStruct;
0030 class UnplacedParallelepiped;
0031
0032 struct ParallelepipedImplementation {
0033
0034   using PlacedShape_t    = PlacedParallelepiped;
0035   using UnplacedStruct_t = ParallelepipedStruct<Precision>;
0036   using UnplacedVolume_t = UnplacedParallelepiped;
0037
0038   VECCORE_ATT_HOST_DEVICE
0039   static void PrintType()
0040   {
0041     // printf("SpecializedParallelepiped<%i, %i>", transCodeT, rotCodeT);
0042   }
0043
0044   template <typename Stream>
0045   static void PrintType(Stream &s, int transCodeT = translation::kGeneric, int rotCodeT = rotation::kGeneric)
0046   {
0047     s << "SpecializedParallelepiped<" << transCodeT << "," << rotCodeT << ">";
0048   }
0049
0050   template <typename Stream>
0051   static void PrintImplementationType(Stream & /*s*/)
0052   {
0053     // s << "ParallelepipedImplementation<" << transCodeT << "," << rotCodeT << ">";
0054   }
0055
0056   template <typename Stream>
0057   static void PrintUnplacedType(Stream & /*s*/)
0058   {
0059     // s << "UnplacedParallelepiped";
0060   }
0061
0062   template <typename Real_v>
0063   VECGEOM_FORCE_INLINE
0064   VECCORE_ATT_HOST_DEVICE
0065   static void Transform(UnplacedStruct_t const &unplaced, Vector3D<Real_v> &point)
0066   {
0067     point.y() -= unplaced.fTanThetaSinPhi * point.z();
0068     point.x() -= unplaced.fTanThetaCosPhi * point.z() + unplaced.fTanAlpha * point.y();
0069   }
0070
0071   template <typename Real_v>
0072   VECGEOM_FORCE_INLINE
0073   VECCORE_ATT_HOST_DEVICE
0074   static void SafetyVector(UnplacedStruct_t const &unplaced, Vector3D<Real_v> const &localPoint,
0075                            Vector3D<Real_v> &safety)
0076   {
0077     safety = localPoint.Abs() - Vector3D<Real_v>(unplaced.fDimensions);
0078     safety.x() *= unplaced.fCtx;
0079     safety.y() *= unplaced.fCty;
0080   }
0081
0082   template <typename Real_v, typename Bool_v>
0083   VECGEOM_FORCE_INLINE
0084   VECCORE_ATT_HOST_DEVICE
0085   static void Contains(UnplacedStruct_t const &unplaced, Vector3D<Real_v> const &point, Bool_v &inside)
0086   {
0087     Vector3D<Real_v> localPoint(point);
0088     Vector3D<Real_v> safetyVector;
0089     Transform<Real_v>(unplaced, localPoint);
0090     SafetyVector<Real_v>(unplaced, localPoint, safetyVector);
0091
0092     inside = safetyVector.Max() < Real_v(0.0);
0093   }
0094
0095   template <typename Real_v, typename Inside_v>
0096   VECGEOM_FORCE_INLINE
0097   VECCORE_ATT_HOST_DEVICE
0098   static void Inside(UnplacedStruct_t const &unplaced, Vector3D<Real_v> const &point, Inside_v &inside)
0099   {
0100     Vector3D<Real_v> localPoint(point);
0101     Vector3D<Real_v> safetyVector;
0102     Transform<Real_v>(unplaced, localPoint);
0103     SafetyVector<Real_v>(unplaced, localPoint, safetyVector);
0104
0105     Real_v safety = safetyVector.Max();
0106     inside        = vecCore::Blend(safety < Real_v(0.0), Inside_v(kInside), Inside_v(kOutside));
0107     vecCore__MaskedAssignFunc(inside, vecCore::math::Abs(safety) < Real_v(kHalfTolerance), Inside_v(kSurface));
0108   }
0109
0110   template <typename Real_v>
0111   VECGEOM_FORCE_INLINE
0112   VECCORE_ATT_HOST_DEVICE
0113   static void SafetyToIn(UnplacedStruct_t const &unplaced, Vector3D<Real_v> const &point, Real_v &safety)
0114   {
0115     Vector3D<Real_v> localPoint(point);
0116     Vector3D<Real_v> safetyVector;
0117     Transform<Real_v>(unplaced, localPoint);
0118     SafetyVector<Real_v>(unplaced, localPoint, safetyVector);
0119
0120     safety = safetyVector.Max();
0121     vecCore::MaskedAssign(safety, vecCore::math::Abs(safety) < Real_v(kHalfTolerance), Real_v(0.0));
0122   }
0123
0124   template <typename Real_v>
0125   VECGEOM_FORCE_INLINE
0126   VECCORE_ATT_HOST_DEVICE
0127   static void SafetyToOut(UnplacedStruct_t const &unplaced, Vector3D<Real_v> const &point, Real_v &safety)
0128   {
0129     Vector3D<Real_v> localPoint(point);
0130     Vector3D<Real_v> safetyVector;
0131     Transform<Real_v>(unplaced, localPoint);
0132     SafetyVector<Real_v>(unplaced, localPoint, safetyVector);
0133
0134     safety = -safetyVector.Max();
0135     vecCore::MaskedAssign(safety, vecCore::math::Abs(safety) < Real_v(kHalfTolerance), Real_v(0.0));
0136   }
0137
0138   template <typename Real_v>
0139   VECGEOM_FORCE_INLINE
0140   VECCORE_ATT_HOST_DEVICE
0141   static void DistanceToIn(UnplacedStruct_t const &unplaced, Vector3D<Real_v> const &point,
0142                            Vector3D<Real_v> const &direction, Real_v const & /* stepMax */, Real_v &distance)
0143   {
0144     using Bool_v = vecCore::Mask_v<Real_v>;
0145
0146     // Transform point and direction to local (oblique) system of coordinates,
0147     // compute safety vector
0148     Vector3D<Real_v> p(point);
0149     Vector3D<Real_v> v(direction);
0150     Vector3D<Real_v> safetyVector;
0151
0152     Transform<Real_v>(unplaced, p);
0153     Transform<Real_v>(unplaced, v);
0154     SafetyVector<Real_v>(unplaced, p, safetyVector);
0155
0156     // Check if point is leaving shape
0157     Bool_v leaving(false);
0158     leaving |= (safetyVector.x() >= -kHalfTolerance && p.x() * v.x() >= Real_v(0.));
0159     leaving |= (safetyVector.y() >= -kHalfTolerance && p.y() * v.y() >= Real_v(0.));
0160     leaving |= (safetyVector.z() >= -kHalfTolerance && p.z() * v.z() >= Real_v(0.));
0161
0162     // Compute distances
0163     const Vector3D<Real_v> invDir(Real_v(1.) / NonZero(v.x()), Real_v(1.) / NonZero(v.y()),
0164                                   Real_v(1.) / NonZero(v.z()));
0165     const Vector3D<Real_v> signDir(Sign(invDir.x()), Sign(invDir.y()), Sign(invDir.z()));
0166     const Vector3D<Real_v> temp = signDir * unplaced.fDimensions;
0167     const Real_v distIn         = ((-temp - p) * invDir).Max();
0168     const Real_v distOut        = ((temp - p) * invDir).Min();
0169
0170     // Set distance to in
0171     distance = Real_v(kInfLength);
0172     vecCore__MaskedAssignFunc(distance, !leaving && distOut > distIn + kHalfTolerance, distIn);
0173   }
0174
0175   template <typename Real_v>
0176   VECGEOM_FORCE_INLINE
0177   VECCORE_ATT_HOST_DEVICE
0178   static void DistanceToOut(UnplacedStruct_t const &unplaced, Vector3D<Real_v> const &point,
0179                             Vector3D<Real_v> const &direction, Real_v const & /* stepMax */, Real_v &distance)
0180   {
0181     // Transform point and direction to local (oblique) system of coordinates,
0182     // compute safety vector
0183     Vector3D<Real_v> p(point);
0184     Vector3D<Real_v> v(direction);
0185     Vector3D<Real_v> safetyVector;
0186
0187     Transform<Real_v>(unplaced, p);
0188     Transform<Real_v>(unplaced, v);
0189     SafetyVector<Real_v>(unplaced, p, safetyVector);
0190
0191     // Compute distance to out
0192     const Vector3D<Real_v> dir(NonZero(v.x()), NonZero(v.y()), NonZero(v.z()));
0193     const Vector3D<Real_v> signDir(Sign(dir.x()), Sign(dir.y()), Sign(dir.z()));
0194     distance = ((signDir * unplaced.fDimensions - p) / dir).Min();
0195
0196     // Set distance to out
0197     vecCore__MaskedAssignFunc(distance, safetyVector.Max() > kHalfTolerance, Real_v(-1.0));
0198   }
0199
0200   template <typename Real_v>
0201   VECGEOM_FORCE_INLINE
0202   VECCORE_ATT_HOST_DEVICE
0203   static Vector3D<Real_v> NormalKernel(UnplacedStruct_t const &unplaced, Vector3D<Real_v> const &point,
0204                                        typename vecCore::Mask_v<Real_v> &valid)
0205   {
0206     // Compute normal at the point on the surface.
0207     // In case the point is not on the surface, set valid = false.
0208     // Must return a valid vector (even if the point is not on the surface).
0209     // On edge or corner, provide an average normal of all facets within the tolerance.
0210     Vector3D<Real_v> normal(0.);
0211     valid = true;
0212
0213     // Transform point to local (oblique) system of coordinates and compute safety vector
0214     Vector3D<Real_v> p(point);
0215     Vector3D<Real_v> safetyVector;
0216     Transform<Real_v>(unplaced, p);
0217     SafetyVector<Real_v>(unplaced, p, safetyVector);
0218
0219     // Set normal
0220     const Vector3D<Real_v> signs(Sign(p.x()), Sign(p.y()), Sign(p.z()));
0221     vecCore__MaskedAssignFunc(normal, Abs(safetyVector.z()) <= kHalfTolerance, Vector3D<Real_v>(0., 0., signs.z()));
0222     vecCore__MaskedAssignFunc(normal, Abs(safetyVector.y()) <= kHalfTolerance,
0223                               normal + signs.y() * unplaced.fNormals[1]);
0224     vecCore__MaskedAssignFunc(normal, Abs(safetyVector.x()) <= kHalfTolerance,
0225                               normal + signs.x() * unplaced.fNormals[0]);
0226
0227     Real_v mag2 = normal.Mag2();
0228     vecCore__MaskedAssignFunc(normal, mag2 > 1., normal.Unit());
0229     if (vecCore::MaskFull(mag2 > Real_v(0.))) return normal;
0230
0231     // Point is not on the surface - normally, this should never be.
0232     // Return normal of the nearest face.
0233     vecCore__MaskedAssignFunc(valid, mag2 == Real_v(0.), false);
0234     Real_v safety = safetyVector.Max();
0235     normal        = signs.x() * unplaced.fNormals[0];
0236     vecCore__MaskedAssignFunc(normal, safetyVector.y() == safety, signs.y() * unplaced.fNormals[1]);
0237     vecCore__MaskedAssignFunc(normal, safetyVector.z() == safety, signs.z() * unplaced.fNormals[2]);
0238     return normal;
0239   }
0240 };
0241 } // namespace VECGEOM_IMPL_NAMESPACE
0242 } // namespace vecgeom
0243
0244 #endif // VECGEOM_VOLUMES_KERNEL_PARALLELEPIPEDIMPLEMENTATION_H_