1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163
|
// This file is part of VecGeom and is distributed under the
// conditions in the file LICENSE.txt in the top directory.
// For the full list of authors see CONTRIBUTORS.txt and `git log`.
/// Declaration of a struct with data members for the UnplacedTet class
/// @file volumes/TetStruct.h
/// @author Raman Sehgal, Evgueni Tcherniaev
#ifndef VECGEOM_VOLUMES_TETSTRUCT_H_
#define VECGEOM_VOLUMES_TETSTRUCT_H_
#include "VecGeom/base/Global.h"
#include "VecGeom/base/Vector3D.h"
namespace vecgeom {
inline namespace VECGEOM_IMPL_NAMESPACE {
/// Struct encapsulating data members of the unplaced tetrahedron
template <typename T = double>
struct TetStruct {
Vector3D<T> fVertex[4]; ///< Array of the tetrahedron vertices
struct {
Vector3D<T> n; ///< Normal of the plane
T d; ///< Distance from origin to the plane
} fPlane[4]; ///< The tetrahedron face planes
Precision fCubicVolume; ///< Volume of the tetrahedron
Precision fSurfaceArea; ///< Surface area of the tetrahedron
/// Empty constructor
VECCORE_ATT_HOST_DEVICE
TetStruct() {}
/// Constructor from four points
/// @param p0 Point given as array
/// @param p1 Point given as array
/// @param p2 Point given as array
/// @param p3 Point given as array
VECCORE_ATT_HOST_DEVICE
TetStruct(const T p0[], const T p1[], const T p2[], const T p3[])
{
Vector3D<T> vertices[4];
vertices[0].Set(p0[0], p0[1], p0[2]);
vertices[1].Set(p1[0], p1[1], p1[2]);
vertices[2].Set(p2[0], p2[1], p2[2]);
vertices[3].Set(p3[0], p3[1], p3[2]);
CalculateCached(vertices[0], vertices[1], vertices[2], vertices[3]);
}
/// Constructor from four points
/// @param p0 Point given as 3D vector
/// @param p1 Point given as 3D vector
/// @param p2 Point given as 3D vector
/// @param p3 Point given as 3D vector
VECCORE_ATT_HOST_DEVICE
TetStruct(const Vector3D<T> p0, const Vector3D<T> p1, const Vector3D<T> p2, const Vector3D<T> p3)
: fCubicVolume(0.), fSurfaceArea(0.)
{
CalculateCached(p0, p1, p2, p3);
}
/// Set the tetrahedron data members
/// @param p0 Point given as array
/// @param p1 Point given as array
/// @param p2 Point given as array
/// @param p3 Point given as array
VECCORE_ATT_HOST_DEVICE
void CalculateCached(const Vector3D<T> p0, const Vector3D<T> p1, const Vector3D<T> p2, const Vector3D<T> p3)
{
// Fill all the cached values
fVertex[0] = p0;
fVertex[1] = p1;
fVertex[2] = p2;
fVertex[3] = p3;
// if (CheckDegeneracy()) std::cout << "DeGenerate Tetrahedron not allowed" << std::endl;
CheckDegeneracy();
Vector3D<Precision> n0 = (fVertex[1] - fVertex[0]).Cross(fVertex[2] - fVertex[0]).Unit();
Vector3D<Precision> n1 = (fVertex[2] - fVertex[1]).Cross(fVertex[3] - fVertex[1]).Unit();
Vector3D<Precision> n2 = (fVertex[3] - fVertex[2]).Cross(fVertex[0] - fVertex[2]).Unit();
Vector3D<Precision> n3 = (fVertex[0] - fVertex[3]).Cross(fVertex[1] - fVertex[3]).Unit();
if (n0.Dot(fVertex[3] - fVertex[0]) > 0) n0 = -n0;
if (n1.Dot(fVertex[0] - fVertex[1]) > 0) n1 = -n1;
if (n2.Dot(fVertex[1] - fVertex[2]) > 0) n2 = -n2;
if (n3.Dot(fVertex[2] - fVertex[3]) > 0) n3 = -n3;
fPlane[0].n = n0;
fPlane[0].d = -n0.Dot(fVertex[0]);
// fPlane[0].d = n0.Dot(fVertex[0]);
fPlane[1].n = n1;
fPlane[1].d = -n1.Dot(fVertex[1]);
// fPlane[1].d = n1.Dot(fVertex[1]);
fPlane[2].n = n2;
fPlane[2].d = -n2.Dot(fVertex[2]);
// fPlane[2].d = n2.Dot(fVertex[2]);
fPlane[3].n = n3;
fPlane[3].d = -n3.Dot(fVertex[3]);
// fPlane[3].d = n3.Dot(fVertex[3]);
for (int i = 0; i < 4; i++) {
// std::cout << "Plane[" << i << "] = " << fPlane[i].n << " " << fPlane[i].d << std::endl;
}
}
/// Set volume of the tetrahedron
VECCORE_ATT_HOST_DEVICE
void CalcCapacity()
{
fCubicVolume =
vecCore::math::Abs((fVertex[1] - fVertex[0]).Dot((fVertex[2] - fVertex[0]).Cross(fVertex[3] - fVertex[0]))) /
6.;
}
/// Set surface area of the tetrahedron
VECCORE_ATT_HOST_DEVICE
void CalcSurfaceArea()
{
fSurfaceArea = ((fVertex[1] - fVertex[0]).Cross(fVertex[2] - fVertex[0]).Mag() +
(fVertex[2] - fVertex[1]).Cross(fVertex[3] - fVertex[1]).Mag() +
(fVertex[3] - fVertex[2]).Cross(fVertex[0] - fVertex[2]).Mag() +
(fVertex[0] - fVertex[3]).Cross(fVertex[1] - fVertex[3]).Mag()) *
0.5;
}
/// Set the tetrahedron
/// @param p0 Point given as array
/// @param p1 Point given as array
/// @param p2 Point given as array
/// @param p3 Point given as array
VECCORE_ATT_HOST_DEVICE
void SetParameters(const Vector3D<T> p0, const Vector3D<T> p1, const Vector3D<T> p2, const Vector3D<T> p3)
{
CalculateCached(p0, p1, p2, p3);
}
/// Check correctness of the tetrahedron data
/// @return false if tetrahedron is degenerate
VECCORE_ATT_HOST_DEVICE
bool CheckDegeneracy()
{
CalcCapacity();
CalcSurfaceArea();
if (fCubicVolume < kTolerance * kTolerance * kTolerance) return true;
Precision s0 = (fVertex[1] - fVertex[0]).Cross(fVertex[2] - fVertex[0]).Mag() * 0.5;
Precision s1 = (fVertex[2] - fVertex[1]).Cross(fVertex[3] - fVertex[1]).Mag() * 0.5;
Precision s2 = (fVertex[3] - fVertex[2]).Cross(fVertex[0] - fVertex[2]).Mag() * 0.5;
Precision s3 = (fVertex[0] - fVertex[3]).Cross(fVertex[1] - fVertex[3]).Mag() * 0.5;
return (vecCore::math::Max(vecCore::math::Max(vecCore::math::Max(s0, s1), s2), s3) < 2 * kTolerance);
}
};
} // namespace VECGEOM_IMPL_NAMESPACE
} // namespace vecgeom
#endif
|