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   VECCORE_ATT_HOST_DEVICE
0038   static void PrintType()
0039   {
0040     //  printf("SpecializedTet<%i, %i>", transCodeT, rotCodeT);
0041   }
0042
0043   template <typename Stream>
0044   static void PrintType(Stream &st, int transCodeT = translation::kGeneric, int rotCodeT = rotation::kGeneric)
0045   {
0046     st << "SpecializedTet<" << transCodeT << "," << rotCodeT << ">";
0047   }
0048
0049   template <typename Stream>
0050   static void PrintImplementationType(Stream &st)
0051   {
0052     (void)st;
0053     // st << "TetImplementation<" << transCodeT << "," << rotCodeT << ">";
0054   }
0055
0056   template <typename Stream>
0057   static void PrintUnplacedType(Stream &st)
0058   {
0059     (void)st;
0060     // TODO: this is wrong
0061     // st << "UnplacedTet";
0062   }
0063
0064   template <typename Real_v, typename Bool_v>
0065   VECGEOM_FORCE_INLINE
0066   VECCORE_ATT_HOST_DEVICE
0067   static void Contains(UnplacedStruct_t const &tet, Vector3D<Real_v> const &point, Bool_v &inside)
0068   {
0069     Bool_v unused, outside;
0070     GenericKernelForContainsAndInside<Real_v, Bool_v, false>(tet, point, unused, outside);
0071     inside = !outside;
0072   }
0073
0074   // BIG QUESTION: DO WE WANT TO GIVE ALL 3 TEMPLATE PARAMETERS
0075   // -- OR -- DO WE WANT TO DEDUCE Bool_v, Index_t from Real_v???
0076   template <typename Real_v, typename Inside_t>
0077   VECGEOM_FORCE_INLINE
0078   VECCORE_ATT_HOST_DEVICE
0079   static void Inside(UnplacedStruct_t const &tet, Vector3D<Real_v> const &point, Inside_t &inside)
0080   {
0081
0082     using Bool_v       = vecCore::Mask_v<Real_v>;
0083     using InsideBool_v = vecCore::Mask_v<Inside_t>;
0084     Bool_v completelyinside, completelyoutside;
0085     GenericKernelForContainsAndInside<Real_v, Bool_v, true>(tet, point, completelyinside, completelyoutside);
0086     inside = EInside::kSurface;
0087     vecCore::MaskedAssign(inside, (InsideBool_v)completelyoutside, Inside_t(EInside::kOutside));
0088     vecCore::MaskedAssign(inside, (InsideBool_v)completelyinside, Inside_t(EInside::kInside));
0089   }
0090
0091   template <typename Real_v, typename Bool_v, bool ForInside>
0092   VECGEOM_FORCE_INLINE
0093   VECCORE_ATT_HOST_DEVICE
0094   static void GenericKernelForContainsAndInside(UnplacedStruct_t const &tet, Vector3D<Real_v> const &localPoint,
0095                                                 Bool_v &completelyinside, Bool_v &completelyoutside)
0096   {
0097     /* Logic to check where the point is inside or not.
0098     **
0099     ** if ForInside is false then it will only check if the point is outside,
0100     ** and is used by Contains function
0101     **
0102     ** if ForInside is true then it will check whether the point is inside or outside,
0103     ** and if neither inside nor outside then it is on the surface.
0104     ** and is used by Inside function
0105     */
0106
0107     Real_v dist[4];
0108     for (int i = 0; i < 4; ++i) {
0109       Vector3D<Real_v> n = tet.fPlane[i].n;
0110       dist[i]            = n.Dot(localPoint) + tet.fPlane[i].d;
0111     }
0112     Real_v safety = vecCore::math::Max(vecCore::math::Max(vecCore::math::Max(dist[0], dist[1]), dist[2]), dist[3]);
0113
0114     completelyoutside = safety > kHalfTolerance;
0115     if (ForInside) completelyinside = safety <= -kHalfTolerance;
0116     return;
0117   }
0118
0119   template <typename Real_v>
0120   VECGEOM_FORCE_INLINE
0121   VECCORE_ATT_HOST_DEVICE
0122   static void DistanceToIn(UnplacedStruct_t const &tet, Vector3D<Real_v> const &point,
0123                            Vector3D<Real_v> const &direction, Real_v const & /*stepMax*/, Real_v &distance)
0124   {
0125     /* Logic to calculate Distance from outside point to the Tet surface */
0126     // using Bool_v       = vecCore::Mask_v<Real_v>;
0127     distance           = -kInfLength;
0128     Real_v distanceOut = kInfLength;
0129     Real_v absSafe     = kInfLength;
0130
0131     Real_v cosa[4];
0132     Real_v safe[4];
0133     Real_v dist[4];
0134     for (int i = 0; i < 4; ++i) {
0135       cosa[i] = NonZero(Vector3D<Real_v>(tet.fPlane[i].n).Dot(direction));
0136       safe[i] = Vector3D<Real_v>(tet.fPlane[i].n).Dot(point) + tet.fPlane[i].d;
0137       dist[i] = -safe[i] / cosa[i];
0138     }
0139
0140     for (int i = 0; i < 4; ++i) {
0141       vecCore__MaskedAssignFunc(distance, (cosa[i] < Real_v(0.)), vecCore::math::Max(distance, dist[i]));
0142       vecCore__MaskedAssignFunc(distanceOut, (cosa[i] > Real_v(0.)), vecCore::math::Min(distanceOut, dist[i]));
0143       vecCore__MaskedAssignFunc(absSafe, (cosa[i] > Real_v(0.)), vecCore::math::Min(absSafe, vecCore::math::Abs(safe[i])));
0144     }
0145
0146     vecCore::MaskedAssign(distance,
0147                           distance >= distanceOut || distanceOut <= kHalfTolerance || absSafe <= -kHalfTolerance,
0148                           Real_v(kInfLength));
0149   }
0150
0151   template <typename Real_v>
0152   VECGEOM_FORCE_INLINE
0153   VECCORE_ATT_HOST_DEVICE
0154   static void DistanceToOut(UnplacedStruct_t const &tet, Vector3D<Real_v> const &point,
0155                             Vector3D<Real_v> const &direction, Real_v const & /* stepMax */, Real_v &distance)
0156   {
0157     /* Logic to calculate Distance from inside point to the Tet surface */
0158     distance      = kInfLength;
0159     Real_v safety = -kInfLength;
0160
0161     Real_v cosa[4];
0162     Real_v safe[4];
0163     for (int i = 0; i < 4; ++i) {
0164       cosa[i] = NonZero(Vector3D<Real_v>(tet.fPlane[i].n).Dot(direction));
0165       safe[i] = Vector3D<Real_v>(tet.fPlane[i].n).Dot(point) + tet.fPlane[i].d;
0166       safety  = vecCore::math::Max(safety, safe[i]);
0167     }
0168
0169     for (int i = 0; i < 4; ++i) {
0170       vecCore__MaskedAssignFunc(distance, (cosa[i] > Real_v(0.)), vecCore::math::Min(distance, -safe[i] / cosa[i]));
0171     }
0172     vecCore::MaskedAssign(distance, safety > kHalfTolerance, Real_v(-1.));
0173   }
0174
0175   template <typename Real_v>
0176   VECGEOM_FORCE_INLINE
0177   VECCORE_ATT_HOST_DEVICE
0178   static void SafetyToIn(UnplacedStruct_t const &tet, Vector3D<Real_v> const &point, Real_v &safety)
0179   {
0180     /* Logic to calculate Safety from outside point to the Tet surface */
0181
0182     Real_v dist[4];
0183     for (int i = 0; i < 4; ++i) {
0184       dist[i] = point.Dot(tet.fPlane[i].n) + tet.fPlane[i].d;
0185     }
0186     safety = vecCore::math::Max(vecCore::math::Max(vecCore::math::Max(dist[0], dist[1]), dist[2]), dist[3]);
0187     vecCore::MaskedAssign(safety, vecCore::math::Abs(safety) <= kHalfTolerance, Real_v(0.));
0188   }
0189
0190   template <typename Real_v>
0191   VECGEOM_FORCE_INLINE
0192   VECCORE_ATT_HOST_DEVICE
0193   static void SafetyToOut(UnplacedStruct_t const &tet, Vector3D<Real_v> const &point, Real_v &safety)
0194   {
0195     /* Logic to calculate Safety from inside point to the Tet surface */
0196
0197     Real_v dist[4];
0198     for (int i = 0; i < 4; ++i) {
0199       dist[i] = point.Dot(tet.fPlane[i].n) + tet.fPlane[i].d;
0200     }
0201     safety = -vecCore::math::Max(vecCore::math::Max(vecCore::math::Max(dist[0], dist[1]), dist[2]), dist[3]);
0202     vecCore::MaskedAssign(safety, vecCore::math::Abs(safety) <= kHalfTolerance, Real_v(0.));
0203   }
0204
0205   template <typename Real_v>
0206   VECGEOM_FORCE_INLINE
0207   VECCORE_ATT_HOST_DEVICE
0208   static Vector3D<Real_v> NormalKernel(UnplacedStruct_t const &tet, Vector3D<Real_v> const &point,
0209                                        typename vecCore::Mask_v<Real_v> &valid)
0210   {
0211     Vector3D<Real_v> normal(0.);
0212     valid = true;
0213
0214     Real_v dist[4];
0215     for (int i = 0; i < 4; ++i) {
0216       Vector3D<Real_v> n = tet.fPlane[i].n;
0217       dist[i]            = n.Dot(point) + tet.fPlane[i].d;
0218       vecCore__MaskedAssignFunc(normal, vecCore::math::Abs(dist[i]) <= kHalfTolerance, normal + tet.fPlane[i].n)
0219     }
0220     vecCore::Mask_v<Real_v> done = normal.Mag2() > 1;
0221     vecCore__MaskedAssignFunc(normal, done, normal.Unit());
0222
0223     done = normal.Mag2() > 0.;
0224     if (vecCore::MaskFull(done)) return normal;
0225
0226     // Point is not on the surface - normally, this should never be.
0227     // Return normal of the nearest face.
0228     //
0229     vecCore__MaskedAssignFunc(valid, !done, false);
0230
0231     Real_v safety(-kInfLength);
0232     for (int i = 0; i < 4; ++i) {
0233       vecCore__MaskedAssignFunc(normal, dist[i] > safety && !done, tet.fPlane[i].n);
0234       vecCore__MaskedAssignFunc(safety, dist[i] > safety && !done, dist[i]);
0235     }
0236     return normal;
0237   }
0238 };
0239 } // namespace VECGEOM_IMPL_NAMESPACE
0240 } // namespace vecgeom
0241
0242 #endif // VECGEOM_VOLUMES_KERNEL_TETIMPLEMENTATION_H_