VecGeom/volumes/TetStruct.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 /// Declaration of a struct with data members for the UnplacedTet class
0006 /// @file volumes/TetStruct.h
0007 /// @author Raman Sehgal, Evgueni Tcherniaev
0008
0009 #ifndef VECGEOM_VOLUMES_TETSTRUCT_H_
0010 #define VECGEOM_VOLUMES_TETSTRUCT_H_
0011 #include "VecGeom/base/Global.h"
0012 #include "VecGeom/base/Vector3D.h"
0013
0014 namespace vecgeom {
0015
0016 inline namespace VECGEOM_IMPL_NAMESPACE {
0017
0018 /// Struct encapsulating data members of the unplaced tetrahedron
0019 template <typename T = double>
0020 struct TetStruct {
0021
0022   Vector3D<T> fVertex[4]; ///< Array of the tetrahedron vertices
0023   struct {
0024     Vector3D<T> n; ///< Normal of the plane
0025     T d;           ///< Distance from  origin to the plane
0026   } fPlane[4];     ///< The tetrahedron face planes
0027
0028   Precision fCubicVolume; ///< Volume of the tetrahedron
0029   Precision fSurfaceArea; ///< Surface area of the tetrahedron
0030
0031   /// Empty constructor
0032   VECCORE_ATT_HOST_DEVICE
0033   TetStruct() {}
0034
0035   /// Constructor from four points
0036   /// @param p0 Point given as array
0037   /// @param p1 Point given as array
0038   /// @param p2 Point given as array
0039   /// @param p3 Point given as array
0040   VECCORE_ATT_HOST_DEVICE
0041   TetStruct(const T p0[], const T p1[], const T p2[], const T p3[])
0042   {
0043     Vector3D<T> vertices[4];
0044     vertices[0].Set(p0[0], p0[1], p0[2]);
0045     vertices[1].Set(p1[0], p1[1], p1[2]);
0046     vertices[2].Set(p2[0], p2[1], p2[2]);
0047     vertices[3].Set(p3[0], p3[1], p3[2]);
0048
0049     CalculateCached(vertices[0], vertices[1], vertices[2], vertices[3]);
0050   }
0051
0052   /// Constructor from four points
0053   /// @param p0 Point given as 3D vector
0054   /// @param p1 Point given as 3D vector
0055   /// @param p2 Point given as 3D vector
0056   /// @param p3 Point given as 3D vector
0057   VECCORE_ATT_HOST_DEVICE
0058   TetStruct(const Vector3D<T> p0, const Vector3D<T> p1, const Vector3D<T> p2, const Vector3D<T> p3)
0059       : fCubicVolume(0.), fSurfaceArea(0.)
0060   {
0061
0062     CalculateCached(p0, p1, p2, p3);
0063   }
0064
0065   /// Set the tetrahedron data members
0066   /// @param p0 Point given as array
0067   /// @param p1 Point given as array
0068   /// @param p2 Point given as array
0069   /// @param p3 Point given as array
0070   VECCORE_ATT_HOST_DEVICE
0071   void CalculateCached(const Vector3D<T> p0, const Vector3D<T> p1, const Vector3D<T> p2, const Vector3D<T> p3)
0072   {
0073     // Fill all the cached values
0074     fVertex[0] = p0;
0075     fVertex[1] = p1;
0076     fVertex[2] = p2;
0077     fVertex[3] = p3;
0078
0079     // if (CheckDegeneracy()) std::cout << "DeGenerate Tetrahedron not allowed" << std::endl;
0080     CheckDegeneracy();
0081
0082     Vector3D<Precision> n0 = (fVertex[1] - fVertex[0]).Cross(fVertex[2] - fVertex[0]).Unit();
0083     Vector3D<Precision> n1 = (fVertex[2] - fVertex[1]).Cross(fVertex[3] - fVertex[1]).Unit();
0084     Vector3D<Precision> n2 = (fVertex[3] - fVertex[2]).Cross(fVertex[0] - fVertex[2]).Unit();
0085     Vector3D<Precision> n3 = (fVertex[0] - fVertex[3]).Cross(fVertex[1] - fVertex[3]).Unit();
0086
0087     if (n0.Dot(fVertex[3] - fVertex[0]) > 0) n0 = -n0;
0088     if (n1.Dot(fVertex[0] - fVertex[1]) > 0) n1 = -n1;
0089     if (n2.Dot(fVertex[1] - fVertex[2]) > 0) n2 = -n2;
0090     if (n3.Dot(fVertex[2] - fVertex[3]) > 0) n3 = -n3;
0091
0092     fPlane[0].n = n0;
0093     fPlane[0].d = -n0.Dot(fVertex[0]);
0094     // fPlane[0].d = n0.Dot(fVertex[0]);
0095
0096     fPlane[1].n = n1;
0097     fPlane[1].d = -n1.Dot(fVertex[1]);
0098     // fPlane[1].d = n1.Dot(fVertex[1]);
0099
0100     fPlane[2].n = n2;
0101     fPlane[2].d = -n2.Dot(fVertex[2]);
0102     // fPlane[2].d = n2.Dot(fVertex[2]);
0103
0104     fPlane[3].n = n3;
0105     fPlane[3].d = -n3.Dot(fVertex[3]);
0106     // fPlane[3].d = n3.Dot(fVertex[3]);
0107
0108     for (int i = 0; i < 4; i++) {
0109       // std::cout << "Plane[" << i << "] = " << fPlane[i].n << "  " << fPlane[i].d << std::endl;
0110     }
0111   }
0112
0113   /// Set volume of the tetrahedron
0114   VECCORE_ATT_HOST_DEVICE
0115   void CalcCapacity()
0116   {
0117     fCubicVolume =
0118         vecCore::math::Abs((fVertex[1] - fVertex[0]).Dot((fVertex[2] - fVertex[0]).Cross(fVertex[3] - fVertex[0]))) /
0119         6.;
0120   }
0121
0122   /// Set surface area of the tetrahedron
0123   VECCORE_ATT_HOST_DEVICE
0124   void CalcSurfaceArea()
0125   {
0126     fSurfaceArea = ((fVertex[1] - fVertex[0]).Cross(fVertex[2] - fVertex[0]).Mag() +
0127                     (fVertex[2] - fVertex[1]).Cross(fVertex[3] - fVertex[1]).Mag() +
0128                     (fVertex[3] - fVertex[2]).Cross(fVertex[0] - fVertex[2]).Mag() +
0129                     (fVertex[0] - fVertex[3]).Cross(fVertex[1] - fVertex[3]).Mag()) *
0130                    0.5;
0131   }
0132
0133   /// Set the tetrahedron
0134   /// @param p0 Point given as array
0135   /// @param p1 Point given as array
0136   /// @param p2 Point given as array
0137   /// @param p3 Point given as array
0138   VECCORE_ATT_HOST_DEVICE
0139   void SetParameters(const Vector3D<T> p0, const Vector3D<T> p1, const Vector3D<T> p2, const Vector3D<T> p3)
0140   {
0141     CalculateCached(p0, p1, p2, p3);
0142   }
0143
0144   /// Check correctness of the tetrahedron data
0145   /// @return false if tetrahedron is degenerate
0146   VECCORE_ATT_HOST_DEVICE
0147   bool CheckDegeneracy()
0148   {
0149     CalcCapacity();
0150     CalcSurfaceArea();
0151     if (fCubicVolume < kTolerance * kTolerance * kTolerance) return true;
0152     Precision s0 = (fVertex[1] - fVertex[0]).Cross(fVertex[2] - fVertex[0]).Mag() * 0.5;
0153     Precision s1 = (fVertex[2] - fVertex[1]).Cross(fVertex[3] - fVertex[1]).Mag() * 0.5;
0154     Precision s2 = (fVertex[3] - fVertex[2]).Cross(fVertex[0] - fVertex[2]).Mag() * 0.5;
0155     Precision s3 = (fVertex[0] - fVertex[3]).Cross(fVertex[1] - fVertex[3]).Mag() * 0.5;
0156     return (vecCore::math::Max(vecCore::math::Max(vecCore::math::Max(s0, s1), s2), s3) < 2 * kTolerance);
0157   }
0158 };
0159
0160 } // namespace VECGEOM_IMPL_NAMESPACE
0161 } // namespace vecgeom
0162
0163 #endif