VecGeom/volumes/Tile.h

0001 /// \file QuadrilateralFacet.h
0002 /// \author Mihaela Gheata (mihaela.gheata@cern.ch)
0003
0004 #ifndef VECGEOM_VOLUMES_TILE_H_
0005 #define VECGEOM_VOLUMES_TILE_H_
0006
0007 #include "VecGeom/base/Vector3D.h"
0008 #include "kernel/GenericKernels.h"
0009
0010 namespace vecgeom {
0011
0012 enum TileType { kTriangle = 3, kQuadrilateral = 4 };
0013
0014 VECGEOM_DEVICE_DECLARE_CONV_TEMPLATE_1v_1t(struct, Tile, size_t, typename);
0015
0016 inline namespace VECGEOM_IMPL_NAMESPACE {
0017
0018 template <size_t, typename>
0019 struct Tile;
0020
0021 template <typename T>
0022 using TriangleFacet = Tile<3, T>;
0023
0024 template <typename T>
0025 using QuadrilateralFacet = Tile<4, T>;
0026
0027 //______________________________________________________________________________
0028 // Basic facet tile structure having NVERT vertices making a convex polygon.
0029 // The vertices making the tile have to be given in anti-clockwise
0030 // order looking from the outsider of the solid where it belongs.
0031 //______________________________________________________________________________
0032 template <size_t NVERT, typename T = double>
0033 struct Tile {
0034   size_t fNvert = 0;               ///< the tile is fully defined after adding the last vertex
0035   Vector3D<T> fVertices[NVERT];    ///< vertices of the tile
0036   Vector3D<T> fSideVectors[NVERT]; ///< side vectors perpendicular to edges
0037   Vector3D<T> fNormal;             ///< normal vector pointing outside
0038   Vector3D<T> fCenter;             ///< Center of the tile
0039   size_t fIndices[NVERT] = {0};    ///< indices for 3 distinct vertices
0040   T fSurfaceArea         = 0;      ///< surface area
0041   bool fConvex           = false;  ///< convexity of the facet with respect to the solid
0042   T fDistance            = 0;      ///< distance between the origin and the triangle plane
0043
0044   VECCORE_ATT_HOST_DEVICE
0045   Tile() {}
0046
0047   VECCORE_ATT_HOST_DEVICE
0048   VECGEOM_FORCE_INLINE
0049   bool SetVertices(Vector3D<T> const &vtx0, Vector3D<T> const &vtx1, Vector3D<T> const &vtx2, size_t ind0 = 0,
0050                    size_t ind1 = 0, size_t ind2 = 0)
0051   {
0052     assert(NVERT == 3);
0053     AddVertex(vtx0, ind0);
0054     AddVertex(vtx1, ind1);
0055     return AddVertex(vtx2, ind2);
0056   }
0057
0058   VECCORE_ATT_HOST_DEVICE
0059   VECGEOM_FORCE_INLINE
0060   bool SetVertices(Vector3D<T> const &vtx0, Vector3D<T> const &vtx1, Vector3D<T> const &vtx2, Vector3D<T> const &vtx3,
0061                    size_t ind0 = 0, size_t ind1 = 0, size_t ind2 = 0, size_t ind3 = 0)
0062   {
0063     assert(NVERT == 4);
0064     AddVertex(vtx0, ind0);
0065     AddVertex(vtx1, ind1);
0066     AddVertex(vtx2, ind2);
0067     return AddVertex(vtx3, ind3);
0068   }
0069
0070   VECCORE_ATT_HOST_DEVICE
0071   VECGEOM_FORCE_INLINE
0072   bool AddVertex(Vector3D<T> const &vtx, size_t ind = 0)
0073   {
0074     fVertices[fNvert] = vtx;
0075     fIndices[fNvert]  = ind;
0076     fNvert++;
0077     if (fNvert < NVERT) return true;
0078     // Check validity
0079     // Get number of different vertices
0080     size_t nvert = NVERT;
0081     for (size_t i = 0; i < NVERT; ++i) {
0082       const Vector3D<T> vi = fVertices[(i + 1) % NVERT] - fVertices[i];
0083       if (vi.Mag2() < kTolerance) {
0084         nvert--;
0085       }
0086     }
0087
0088     if (nvert < 3) {
0089       std::cout << "Tile degenerated: Length of sides of facet are too small." << std::endl;
0090       return false;
0091     }
0092
0093     // Compute normal using non-zero segments
0094
0095     bool degenerated = true;
0096     for (size_t i = 0; i < NVERT - 1; ++i) {
0097       Vector3D<T> e1 = fVertices[i + 1] - fVertices[i];
0098       if (e1.Mag2() < kTolerance) continue;
0099       for (size_t j = i + 1; j < NVERT; ++j) {
0100         Vector3D<T> e2 = fVertices[(j + 1) % NVERT] - fVertices[j];
0101         if (e2.Mag2() < kTolerance) continue;
0102         fNormal = e1.Cross(e2);
0103         // e1 and e2 may be colinear
0104         if (fNormal.Mag2() < kTolerance) continue;
0105         fNormal.Normalize();
0106         degenerated = false;
0107         break;
0108       }
0109       if (!degenerated) break;
0110     }
0111
0112     if (degenerated) {
0113       std::cout << "Tile degenerated 2: Length of sides of facet are too small." << std::endl;
0114       return false;
0115     }
0116
0117     // Compute side vectors
0118     for (size_t i = 0; i < NVERT; ++i) {
0119       Vector3D<T> e1 = fVertices[(i + 1) % NVERT] - fVertices[i];
0120       if (e1.Mag2() < kTolerance) continue;
0121       fSideVectors[i] = fNormal.Cross(e1).Normalized();
0122       fDistance       = -fNormal.Dot(fVertices[i]);
0123       for (size_t j = i + 1; j < i + NVERT; ++j) {
0124         Vector3D<T> e2 = fVertices[(j + 1) % NVERT] - fVertices[j % NVERT];
0125         if (e2.Mag2() < kTolerance)
0126           fSideVectors[j % NVERT] = fSideVectors[(j - 1) % NVERT];
0127         else
0128           fSideVectors[j % NVERT] = fNormal.Cross(e2).Normalized();
0129       }
0130       break;
0131     }
0132
0133     // Compute surface area
0134     fSurfaceArea = 0.;
0135     for (size_t i = 1; i < NVERT - 1; ++i) {
0136       Vector3D<T> e1 = fVertices[i] - fVertices[0];
0137       Vector3D<T> e2 = fVertices[i + 1] - fVertices[0];
0138       fSurfaceArea += 0.5 * (e1.Cross(e2)).Mag();
0139     }
0140     assert(fSurfaceArea > kTolerance * kTolerance);
0141
0142     // Center of the tile
0143     for (size_t i = 0; i < NVERT; ++i)
0144       fCenter += fVertices[i];
0145     fCenter /= NVERT;
0146     return true;
0147   }
0148
0149   VECCORE_ATT_HOST_DEVICE
0150   VECGEOM_FORCE_INLINE
0151   Vector3D<T> const &GetNormal() const { return fNormal; }
0152
0153   VECCORE_ATT_HOST_DEVICE
0154   VECGEOM_FORCE_INLINE
0155   size_t IsNeighbor(Tile<NVERT, T> const &other)
0156   {
0157     // Check if a segment is common
0158     size_t ncommon = 0;
0159     for (size_t ind1 = 0; ind1 < NVERT; ++ind1) {
0160       for (size_t ind2 = 0; ind2 < NVERT; ++ind2) {
0161         if (fIndices[ind1] == other.fIndices[ind2]) ncommon++;
0162       }
0163     }
0164     return ncommon;
0165   }
0166
0167   VECGEOM_FORCE_INLINE
0168   VECCORE_ATT_HOST_DEVICE
0169   bool Contains(Vector3D<T> const &point) const
0170   {
0171     // Check id point within the triangle plane is inside the triangle.
0172     bool inside = true;
0173     for (size_t i = 0; i < NVERT; ++i) {
0174       T saf = (point - fVertices[i]).Dot(fSideVectors[i]);
0175       inside &= saf > -kTolerance;
0176     }
0177     return inside;
0178   }
0179
0180   VECGEOM_FORCE_INLINE
0181   VECCORE_ATT_HOST_DEVICE
0182   T DistPlane(Vector3D<T> const &point) const
0183   {
0184     // Returns distance from point to plane. This is positive if the point is on
0185     // the outside halfspace, negative otherwise.
0186     return (point.Dot(fNormal) + fDistance);
0187   }
0188
0189   VECCORE_ATT_HOST_DEVICE
0190   T DistanceToIn(Vector3D<T> const &point, Vector3D<T> const &direction) const
0191   {
0192     T ndd      = NonZero(direction.Dot(fNormal));
0193     T saf      = DistPlane(point);
0194     bool valid = ndd < 0. && saf > -kTolerance;
0195     if (!valid) return InfinityLength<T>();
0196     T distance = -saf / ndd;
0197     // Propagate the point with the distance to the plane.
0198     Vector3D<T> point_prop = point + distance * direction;
0199     // Check if propagated points hit the triangle
0200     if (!Contains(point_prop)) return InfinityLength<T>();
0201     return distance;
0202   }
0203
0204   VECCORE_ATT_HOST_DEVICE
0205   Precision DistanceToOut(Vector3D<T> const &point, Vector3D<T> const &direction) const
0206   {
0207     T ndd      = NonZero(direction.Dot(fNormal));
0208     T saf      = DistPlane(point);
0209     bool valid = ndd > 0. && saf < kTolerance;
0210     if (!valid) return InfinityLength<T>();
0211     T distance = -saf / ndd;
0212     // Propagate the point with the distance to the plane.
0213     Vector3D<T> point_prop = point + distance * direction;
0214     // Check if propagated points hit the triangle
0215     if (!Contains(point_prop)) return InfinityLength<T>();
0216     return distance;
0217   }
0218
0219   template <bool ToIn>
0220   VECCORE_ATT_HOST_DEVICE
0221   T SafetySq(Vector3D<T> const &point) const
0222   {
0223     T safety = DistPlane(point);
0224     // Find the projection of the point on each plane
0225     Vector3D<T> intersection = point - safety * fNormal;
0226     bool withinBound         = Contains(intersection);
0227     if (ToIn)
0228       withinBound &= safety > -kTolerance;
0229     else
0230       withinBound &= safety < kTolerance;
0231     safety *= safety;
0232     if (withinBound) return safety;
0233
0234     Vector3D<T> safety_outbound = InfinityLength<T>();
0235     for (size_t ivert = 0; ivert < NVERT; ++ivert) {
0236       safety_outbound[ivert] =
0237           DistanceToLineSegmentSquared<kScalar>(fVertices[ivert], fVertices[(ivert + 1) % NVERT], point);
0238     }
0239     return (safety_outbound.Min());
0240   }
0241 };
0242
0243 #ifndef VECCORE_CUDA
0244 std::ostream &operator<<(std::ostream &os, TriangleFacet<double> const &facet);
0245 std::ostream &operator<<(std::ostream &os, QuadrilateralFacet<double> const &facet);
0246 #endif
0247 } // namespace VECGEOM_IMPL_NAMESPACE
0248 } // end namespace vecgeom
0249
0250 #endif // VECGEOM_VOLUMES_TILE_H_