volumes/kernel/TetImplementation.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 tetrahedron
0006 /// @file volumes/kernel/TetImplementation.h
0007 /// @author Raman Sehgal, Evgueni Tcherniaev
0008
0009 #ifndef VECGEOM_VOLUMES_KERNEL_TETIMPLEMENTATION_H_
0010 #define VECGEOM_VOLUMES_KERNEL_TETIMPLEMENTATION_H_
0011
0012 #include "VecGeom/base/Vector3D.h"
0013 #include "VecGeom/volumes/TetStruct.h"
0014 #include "VecGeom/volumes/kernel/GenericKernels.h"
0015 #include <VecCore/VecCore>
0016
0017 #include <cstdio>
0018
0019 namespace vecgeom {
0020
0021 VECGEOM_DEVICE_FORWARD_DECLARE(struct TetImplementation;);
0022 VECGEOM_DEVICE_DECLARE_CONV(struct, TetImplementation);
0023
0024 inline namespace VECGEOM_IMPL_NAMESPACE {
0025
0026 class PlacedTet;
0027 template <typename T>
0028 struct TetStruct;
0029 class UnplacedTet;
0030
0031 struct TetImplementation {
0032
0033   using PlacedShape_t    = PlacedTet;
0034   using UnplacedStruct_t = TetStruct<Precision>;
0035   using UnplacedVolume_t = UnplacedTet;
0036
0037   template <typename Real_v, typename Bool_v>
0038   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void Contains(UnplacedStruct_t const &tet,
0039                                                                     Vector3D<Real_v> const &point, Bool_v &inside)
0040   {
0041     Bool_v unused(false), outside(false);
0042     GenericKernelForContainsAndInside<Real_v, Bool_v, false>(tet, point, unused, outside);
0043     inside = !outside;
0044   }
0045
0046   // BIG QUESTION: DO WE WANT TO GIVE ALL 3 TEMPLATE PARAMETERS
0047   // -- OR -- DO WE WANT TO DEDUCE Bool_v, Index_t from Real_v???
0048   template <typename Real_v, typename Inside_t>
0049   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void Inside(UnplacedStruct_t const &tet,
0050                                                                   Vector3D<Real_v> const &point, Inside_t &inside)
0051   {
0052
0053     using Bool_v       = vecCore::Mask_v<Real_v>;
0054     using InsideBool_v = vecCore::Mask_v<Inside_t>;
0055     Bool_v completelyinside, completelyoutside;
0056     GenericKernelForContainsAndInside<Real_v, Bool_v, true>(tet, point, completelyinside, completelyoutside);
0057     inside = EInside::kSurface;
0058     vecCore::MaskedAssign(inside, (InsideBool_v)completelyoutside, Inside_t(EInside::kOutside));
0059     vecCore::MaskedAssign(inside, (InsideBool_v)completelyinside, Inside_t(EInside::kInside));
0060   }
0061
0062   template <typename Real_v, typename Bool_v, bool ForInside>
0063   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void GenericKernelForContainsAndInside(
0064       UnplacedStruct_t const &tet, Vector3D<Real_v> const &localPoint, Bool_v &completelyinside,
0065       Bool_v &completelyoutside)
0066   {
0067     /* Logic to check where the point is inside or not.
0068     **
0069     ** if ForInside is false then it will only check if the point is outside,
0070     ** and is used by Contains function
0071     **
0072     ** if ForInside is true then it will check whether the point is inside or outside,
0073     ** and if neither inside nor outside then it is on the surface.
0074     ** and is used by Inside function
0075     */
0076
0077     Real_v dist[4];
0078     for (int i = 0; i < 4; ++i) {
0079       Vector3D<Real_v> n = tet.fPlane[i].n;
0080       dist[i]            = n.Dot(localPoint) + tet.fPlane[i].d;
0081     }
0082     Real_v safety = vecCore::math::Max(vecCore::math::Max(vecCore::math::Max(dist[0], dist[1]), dist[2]), dist[3]);
0083
0084     completelyoutside = safety > kHalfTolerance;
0085     if (ForInside) completelyinside = safety <= -kHalfTolerance;
0086     return;
0087   }
0088
0089   template <typename Real_v>
0090   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void DistanceToIn(UnplacedStruct_t const &tet,
0091                                                                         Vector3D<Real_v> const &point,
0092                                                                         Vector3D<Real_v> const &direction,
0093                                                                         Real_v const & /*stepMax*/, Real_v &distance)
0094   {
0095     /* Logic to calculate Distance from outside point to the Tet surface */
0096     // using Bool_v       = vecCore::Mask_v<Real_v>;
0097     distance           = -kInfLength;
0098     Real_v distanceOut = kInfLength;
0099     Real_v absSafe     = kInfLength;
0100
0101     Real_v cosa[4];
0102     Real_v safe[4];
0103     Real_v dist[4];
0104     for (int i = 0; i < 4; ++i) {
0105       cosa[i] = NonZero(Vector3D<Real_v>(tet.fPlane[i].n).Dot(direction));
0106       safe[i] = Vector3D<Real_v>(tet.fPlane[i].n).Dot(point) + tet.fPlane[i].d;
0107       dist[i] = -safe[i] / cosa[i];
0108     }
0109
0110     for (int i = 0; i < 4; ++i) {
0111       vecCore__MaskedAssignFunc(distance, (cosa[i] < -kTolerance), vecCore::math::Max(distance, dist[i]));
0112       vecCore__MaskedAssignFunc(distanceOut, (cosa[i] > kTolerance), vecCore::math::Min(distanceOut, dist[i]));
0113       vecCore__MaskedAssignFunc(absSafe, (cosa[i] > kTolerance),
0114                                 vecCore::math::Min(absSafe, vecCore::math::Abs(safe[i])));
0115     }
0116
0117     vecCore::MaskedAssign(distance, distance >= distanceOut || distanceOut <= kTolerance || absSafe <= -kTolerance,
0118                           Real_v(kInfLength));
0119   }
0120
0121   template <typename Real_v>
0122   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void DistanceToOut(UnplacedStruct_t const &tet,
0123                                                                          Vector3D<Real_v> const &point,
0124                                                                          Vector3D<Real_v> const &direction,
0125                                                                          Real_v const & /* stepMax */, Real_v &distance)
0126   {
0127     /* Logic to calculate Distance from inside point to the Tet surface */
0128     distance      = kInfLength;
0129     Real_v safety = -kInfLength;
0130
0131     Real_v cosa[4];
0132     Real_v safe[4];
0133     for (int i = 0; i < 4; ++i) {
0134       cosa[i] = NonZero(Vector3D<Real_v>(tet.fPlane[i].n).Dot(direction));
0135       safe[i] = Vector3D<Real_v>(tet.fPlane[i].n).Dot(point) + tet.fPlane[i].d;
0136       safety  = vecCore::math::Max(safety, safe[i]);
0137     }
0138
0139     for (int i = 0; i < 4; ++i) {
0140       vecCore__MaskedAssignFunc(distance, (cosa[i] > kTolerance), vecCore::math::Min(distance, -safe[i] / cosa[i]));
0141     }
0142     vecCore::MaskedAssign(distance, safety > kHalfTolerance, Real_v(-1.));
0143   }
0144
0145   template <typename Real_v>
0146   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void SafetyToIn(UnplacedStruct_t const &tet,
0147                                                                       Vector3D<Real_v> const &point, Real_v &safety)
0148   {
0149     /* Logic to calculate Safety from outside point to the Tet surface */
0150
0151     Real_v dist[4];
0152     for (int i = 0; i < 4; ++i) {
0153       dist[i] = point.Dot(tet.fPlane[i].n) + tet.fPlane[i].d;
0154     }
0155     safety = vecCore::math::Max(vecCore::math::Max(vecCore::math::Max(dist[0], dist[1]), dist[2]), dist[3]);
0156     vecCore::MaskedAssign(safety, vecCore::math::Abs(safety) <= kHalfTolerance, Real_v(0.));
0157   }
0158
0159   template <typename Real_v>
0160   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void SafetyToOut(UnplacedStruct_t const &tet,
0161                                                                        Vector3D<Real_v> const &point, Real_v &safety)
0162   {
0163     /* Logic to calculate Safety from inside point to the Tet surface */
0164
0165     Real_v dist[4];
0166     for (int i = 0; i < 4; ++i) {
0167       dist[i] = point.Dot(tet.fPlane[i].n) + tet.fPlane[i].d;
0168     }
0169     safety = -vecCore::math::Max(vecCore::math::Max(vecCore::math::Max(dist[0], dist[1]), dist[2]), dist[3]);
0170     vecCore::MaskedAssign(safety, vecCore::math::Abs(safety) <= kHalfTolerance, Real_v(0.));
0171   }
0172
0173   template <typename Real_v>
0174   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static Vector3D<Real_v> NormalKernel(
0175       UnplacedStruct_t const &tet, Vector3D<Real_v> const &point, typename vecCore::Mask_v<Real_v> &valid)
0176   {
0177     Vector3D<Real_v> normal(0.);
0178     valid = true;
0179
0180     Real_v dist[4];
0181     for (int i = 0; i < 4; ++i) {
0182       Vector3D<Real_v> n = tet.fPlane[i].n;
0183       dist[i]            = n.Dot(point) + tet.fPlane[i].d;
0184       vecCore__MaskedAssignFunc(normal, vecCore::math::Abs(dist[i]) <= kHalfTolerance, normal + tet.fPlane[i].n)
0185     }
0186     vecCore::Mask_v<Real_v> done = normal.Mag2() > 1;
0187     vecCore__MaskedAssignFunc(normal, done, normal.Unit());
0188
0189     done = normal.Mag2() > 0.;
0190     if (vecCore::MaskFull(done)) return normal;
0191
0192     // Point is not on the surface - normally, this should never be.
0193     // Return normal of the nearest face.
0194     //
0195     vecCore__MaskedAssignFunc(valid, !done, false);
0196
0197     Real_v safety(-kInfLength);
0198     for (int i = 0; i < 4; ++i) {
0199       vecCore__MaskedAssignFunc(normal, dist[i] > safety && !done, tet.fPlane[i].n);
0200       vecCore__MaskedAssignFunc(safety, dist[i] > safety && !done, dist[i]);
0201     }
0202     return normal;
0203   }
0204 };
0205 } // namespace VECGEOM_IMPL_NAMESPACE
0206 } // namespace vecgeom
0207
0208 #endif // VECGEOM_VOLUMES_KERNEL_TETIMPLEMENTATION_H_