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