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 "VecGeom/volumes/kernel/BoxImplementation.h"
0017 #include <VecCore/VecCore>
0018
0019 #include <cstdio>
0020
0021 namespace vecgeom {
0022
0023 VECGEOM_DEVICE_FORWARD_DECLARE(struct ParallelepipedImplementation;);
0024 VECGEOM_DEVICE_DECLARE_CONV(struct, ParallelepipedImplementation);
0025
0026 inline namespace VECGEOM_IMPL_NAMESPACE {
0027
0028 class PlacedParallelepiped;
0029 template <typename T>
0030 struct ParallelepipedStruct;
0031 class UnplacedParallelepiped;
0032
0033 struct ParallelepipedImplementation {
0034
0035   using PlacedShape_t    = PlacedParallelepiped;
0036   using UnplacedStruct_t = ParallelepipedStruct<Precision>;
0037   using UnplacedVolume_t = UnplacedParallelepiped;
0038
0039   template <typename Real_v>
0040   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void Transform(UnplacedStruct_t const &unplaced,
0041                                                                      Vector3D<Real_v> &point)
0042   {
0043     point.y() -= unplaced.fTanThetaSinPhi * point.z();
0044     point.x() -= unplaced.fTanThetaCosPhi * point.z() + unplaced.fTanAlpha * point.y();
0045   }
0046
0047   template <typename Real_v>
0048   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void SafetyVector(UnplacedStruct_t const &unplaced,
0049                                                                         Vector3D<Real_v> const &localPoint,
0050                                                                         Vector3D<Real_v> &safety)
0051   {
0052     safety = localPoint.Abs() - Vector3D<Real_v>(unplaced.fDimensions);
0053     safety.x() *= unplaced.fCtx;
0054     safety.y() *= unplaced.fCty;
0055   }
0056
0057   template <typename Real_v, typename Bool_v>
0058   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void Contains(UnplacedStruct_t const &unplaced,
0059                                                                     Vector3D<Real_v> const &point, Bool_v &inside)
0060   {
0061     Vector3D<Real_v> localPoint(point);
0062     Vector3D<Real_v> safetyVector;
0063     Transform<Real_v>(unplaced, localPoint);
0064     SafetyVector<Real_v>(unplaced, localPoint, safetyVector);
0065
0066     inside = safetyVector.Max() < Real_v(0.0);
0067   }
0068
0069   template <typename Real_v, typename Inside_v>
0070   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void Inside(UnplacedStruct_t const &unplaced,
0071                                                                   Vector3D<Real_v> const &point, Inside_v &inside)
0072   {
0073     Vector3D<Real_v> localPoint(point);
0074     Vector3D<Real_v> safetyVector;
0075     Transform<Real_v>(unplaced, localPoint);
0076     SafetyVector<Real_v>(unplaced, localPoint, safetyVector);
0077
0078     Real_v safety = safetyVector.Max();
0079     inside        = vecCore::Blend(safety < Real_v(0.0), Inside_v(kInside), Inside_v(kOutside));
0080     vecCore__MaskedAssignFunc(inside, vecCore::math::Abs(safety) < Real_v(kHalfTolerance), Inside_v(kSurface));
0081   }
0082
0083   template <typename Real_v>
0084   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void SafetyToIn(UnplacedStruct_t const &unplaced,
0085                                                                       Vector3D<Real_v> const &point, Real_v &safety)
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     safety = safetyVector.Max();
0093     vecCore::MaskedAssign(safety, vecCore::math::Abs(safety) < Real_v(kHalfTolerance), Real_v(0.0));
0094   }
0095
0096   template <typename Real_v>
0097   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void SafetyToOut(UnplacedStruct_t const &unplaced,
0098                                                                        Vector3D<Real_v> const &point, Real_v &safety)
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     safety = -safetyVector.Max();
0106     vecCore::MaskedAssign(safety, vecCore::math::Abs(safety) < Real_v(kHalfTolerance), Real_v(0.0));
0107   }
0108
0109   template <typename Real_v>
0110   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void DistanceToIn(UnplacedStruct_t const &unplaced,
0111                                                                         Vector3D<Real_v> const &point,
0112                                                                         Vector3D<Real_v> const &direction,
0113                                                                         Real_v const &stepMax, Real_v &distance)
0114   {
0115     // Transform point and direction to local (oblique) system of coordinates,
0116     // compute safety vector
0117     Vector3D<Real_v> p(point);
0118     Vector3D<Real_v> v(direction);
0119
0120     Transform<Real_v>(unplaced, p);
0121     Transform<Real_v>(unplaced, v);
0122
0123     BoxImplementation::DistanceToIn(BoxStruct<Precision>(unplaced.fDimensions), p, v, stepMax, distance);
0124   }
0125
0126   template <typename Real_v>
0127   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void DistanceToOut(UnplacedStruct_t const &unplaced,
0128                                                                          Vector3D<Real_v> const &point,
0129                                                                          Vector3D<Real_v> const &direction,
0130                                                                          Real_v const &stepMax, Real_v &distance)
0131   {
0132     // Transform point and direction to local (oblique) system of coordinates,
0133     // compute safety vector
0134     Vector3D<Real_v> p(point);
0135     Vector3D<Real_v> v(direction);
0136
0137     Transform<Real_v>(unplaced, p);
0138     Transform<Real_v>(unplaced, v);
0139
0140     BoxImplementation::DistanceToOut(BoxStruct<Precision>(unplaced.fDimensions), p, v, stepMax, distance);
0141   }
0142
0143   template <typename Real_v>
0144   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static Vector3D<Real_v> NormalKernel(
0145       UnplacedStruct_t const &unplaced, Vector3D<Real_v> const &point, typename vecCore::Mask_v<Real_v> &valid)
0146   {
0147     // Compute normal at the point on the surface.
0148     // In case the point is not on the surface, set valid = false.
0149     // Must return a valid vector (even if the point is not on the surface).
0150     // On edge or corner, provide an average normal of all facets within the tolerance.
0151     Vector3D<Real_v> normal(0.);
0152     valid = true;
0153
0154     // Transform point to local (oblique) system of coordinates and compute safety vector
0155     Vector3D<Real_v> p(point);
0156     Vector3D<Real_v> safetyVector;
0157     Transform<Real_v>(unplaced, p);
0158     SafetyVector<Real_v>(unplaced, p, safetyVector);
0159
0160     // Set normal
0161     const Vector3D<Real_v> signs(Sign(p.x()), Sign(p.y()), Sign(p.z()));
0162     vecCore__MaskedAssignFunc(normal, Abs(safetyVector.z()) <= kHalfTolerance, Vector3D<Real_v>(0., 0., signs.z()));
0163     vecCore__MaskedAssignFunc(normal, Abs(safetyVector.y()) <= kHalfTolerance,
0164                               normal + signs.y() * unplaced.fNormals[1]);
0165     vecCore__MaskedAssignFunc(normal, Abs(safetyVector.x()) <= kHalfTolerance,
0166                               normal + signs.x() * unplaced.fNormals[0]);
0167
0168     Real_v mag2 = normal.Mag2();
0169     vecCore__MaskedAssignFunc(normal, mag2 > 1., normal.Unit());
0170     if (vecCore::MaskFull(mag2 > Real_v(0.))) return normal;
0171
0172     // Point is not on the surface - normally, this should never be.
0173     // Return normal of the nearest face.
0174     vecCore__MaskedAssignFunc(valid, mag2 == Real_v(0.), false);
0175     Real_v safety = safetyVector.Max();
0176     normal        = signs.x() * unplaced.fNormals[0];
0177     vecCore__MaskedAssignFunc(normal, safetyVector.y() == safety, signs.y() * unplaced.fNormals[1]);
0178     vecCore__MaskedAssignFunc(normal, safetyVector.z() == safety, signs.z() * unplaced.fNormals[2]);
0179     return normal;
0180   }
0181 };
0182 } // namespace VECGEOM_IMPL_NAMESPACE
0183 } // namespace vecgeom
0184
0185 #endif // VECGEOM_VOLUMES_KERNEL_PARALLELEPIPEDIMPLEMENTATION_H_