VecGeom/volumes/GenTrapStruct.h

0001 /*
0002  * GenTrapStruct.h
0003  *
0004  *  Created on: 17.07.2016
0005  *      Author: mgheata
0006  */
0007
0008 #ifndef VECGEOM_VOLUMES_GENTRAPSTRUCT_H_
0009 #define VECGEOM_VOLUMES_GENTRAPSTRUCT_H_
0010 #include "VecGeom/base/Global.h"
0011 #include "VecGeom/volumes/SecondOrderSurfaceShell.h"
0012
0013 #ifndef VECCORE_CUDA
0014 #include "VecGeom/volumes/TessellatedSection.h"
0015 #endif
0016
0017 namespace vecgeom {
0018
0019 inline namespace VECGEOM_IMPL_NAMESPACE {
0020
0021 // A plain struct without member functions to encapsulate just the parameters
0022 // of a generic trapezoid
0023 template <typename T = double>
0024 struct GenTrapStruct {
0025   using Vertex_t = Vector3D<T>;
0026
0027   Vertex_t fBBdimensions; /** Bounding box dimensions */
0028   Vertex_t fBBorigin;     /** Bounding box origin */
0029   Vertex_t fVertices[8];  /** The eight points that define the Arb8 */
0030
0031   // we also store this in SOA form
0032   T fVerticesX[8]; /** Backed-up X positions of vertices */
0033   T fVerticesY[8]; /** Backed-up Y positions of vertices */
0034   T fTwist[4];     /** Twist angles */
0035
0036   T fDz;            /** The half-height of the GenTrap */
0037   T fInverseDz;     /** Pre-computed 1/fDz */
0038   T fHalfInverseDz; /** Pre-computed 0.5/fDz */
0039   bool fIsTwisted;  /** Twisted flag */
0040
0041   // we store the connecting vectors in SOA Form
0042   // these vectors are used to calculate the polygon at a certain z-height
0043   // moreover: they can be precomputed !!
0044   // Compute intersection between Z plane containing point and the shape
0045   //
0046   T fConnectingComponentsX[4]; /** X components of connecting bottom-top vectors vi */
0047   T fConnectingComponentsY[4]; /** Y components of connecting bottom-top vectors vi */
0048
0049   T fDeltaX[8]; /** X components of connecting horizontal vectors hij */
0050   T fDeltaY[8]; /** Y components of connecting horizontal vectors hij */
0051
0052   bool fDegenerated[4]; /** Flags for each top-bottom edge marking that this is degenerated */
0053
0054   SecondOrderSurfaceShell<4> fSurfaceShell; /** Utility class for twisted surface algorithms */
0055
0056 #ifndef VECCORE_CUDA
0057   TessellatedSection<T> *fTslHelper = nullptr; /** SIMD helper using tessellated clusters */
0058 #endif
0059
0060   VECCORE_ATT_HOST_DEVICE
0061   GenTrapStruct()
0062       : fBBdimensions(), fBBorigin(), fVertices(), fVerticesX(), fVerticesY(), fDz(0.), fInverseDz(0.),
0063         fHalfInverseDz(0.), fIsTwisted(false), fConnectingComponentsX(), fConnectingComponentsY(), fDeltaX(), fDeltaY(),
0064         fSurfaceShell()
0065   {
0066     // Dummy constructor
0067   }
0068
0069   VECCORE_ATT_HOST_DEVICE
0070   GenTrapStruct(const Precision verticesx[], const Precision verticesy[], Precision halfzheight)
0071       : fBBdimensions(), fBBorigin(), fVertices(), fVerticesX(), fVerticesY(), fDz(0.), fInverseDz(0.),
0072         fHalfInverseDz(0.), fIsTwisted(false), fConnectingComponentsX(), fConnectingComponentsY(), fDeltaX(), fDeltaY(),
0073         fSurfaceShell()
0074   {
0075     // Constructor
0076     Initialize(verticesx, verticesy, halfzheight);
0077   }
0078
0079   VECCORE_ATT_HOST_DEVICE
0080   bool Initialize(const Precision verticesx[], const Precision verticesy[], Precision halfzheight)
0081   {
0082     // Initialization based on vertices and half length
0083     fDz            = halfzheight;
0084     fInverseDz     = 1. / halfzheight;
0085     fHalfInverseDz = 0.5 / halfzheight;
0086     fIsTwisted     = false;
0087     fSurfaceShell.Initialize(verticesx, verticesy, halfzheight);
0088
0089     // Set vertices in Vector3D form
0090     for (int i = 0; i < 4; ++i) {
0091       fVertices[i].operator[](0) = verticesx[i];
0092       fVertices[i].operator[](1) = verticesy[i];
0093       fVertices[i].operator[](2) = -halfzheight;
0094     }
0095     for (int i = 4; i < 8; ++i) {
0096       fVertices[i].operator[](0) = verticesx[i];
0097       fVertices[i].operator[](1) = verticesy[i];
0098       fVertices[i].operator[](2) = halfzheight;
0099     }
0100
0101     for (int i = 0; i < 4; ++i) {
0102       int j           = (i + 1) % 4;
0103       fDegenerated[i] = (Abs(verticesx[i] - verticesx[j]) < kTolerance) &&
0104                         (Abs(verticesy[i] - verticesy[j]) < kTolerance) &&
0105                         (Abs(verticesx[i + 4] - verticesx[j + 4]) < kTolerance) &&
0106                         (Abs(verticesy[i + 4] - verticesy[j + 4]) < kTolerance);
0107     }
0108
0109     // Make sure vertices are defined clockwise
0110     Precision sum1 = 0.;
0111     Precision sum2 = 0.;
0112     for (int i = 0; i < 4; ++i) {
0113       int j = (i + 1) % 4;
0114       sum1 += fVertices[i].x() * fVertices[j].y() - fVertices[j].x() * fVertices[i].y();
0115       sum2 += fVertices[i + 4].x() * fVertices[j + 4].y() - fVertices[j + 4].x() * fVertices[i + 4].y();
0116     }
0117
0118     // Me should generate an exception here
0119     if (sum1 * sum2 < -kTolerance) {
0120       printf("ERROR: Unplaced generic trap defined with opposite clockwise\n");
0121       Print();
0122       return false;
0123     }
0124
0125     // Revert sequence of vertices to have them clockwise
0126     if (sum1 > kTolerance) {
0127       printf("INFO: Reverting to clockwise vertices of GenTrap shape:\n");
0128       Print();
0129       Vertex_t vtemp;
0130       vtemp        = fVertices[1];
0131       fVertices[1] = fVertices[3];
0132       fVertices[3] = vtemp;
0133       vtemp        = fVertices[5];
0134       fVertices[5] = fVertices[7];
0135       fVertices[7] = vtemp;
0136     }
0137
0138     // Initialize the vertices components and connecting components
0139     for (int i = 0; i < 4; ++i) {
0140       fConnectingComponentsX[i] = (fVertices[i] - fVertices[i + 4]).x();
0141       fConnectingComponentsY[i] = (fVertices[i] - fVertices[i + 4]).y();
0142       fVerticesX[i]             = fVertices[i].x();
0143       fVerticesX[i + 4]         = fVertices[i + 4].x();
0144       fVerticesY[i]             = fVertices[i].y();
0145       fVerticesY[i + 4]         = fVertices[i + 4].y();
0146     }
0147
0148     // Initialize components of horizontal connecting vectors
0149     for (int i = 0; i < 4; ++i) {
0150       int j          = (i + 1) % 4;
0151       fDeltaX[i]     = fVerticesX[j] - fVerticesX[i];
0152       fDeltaX[i + 4] = fVerticesX[j + 4] - fVerticesX[i + 4];
0153       fDeltaY[i]     = fVerticesY[j] - fVerticesY[i];
0154       fDeltaY[i + 4] = fVerticesY[j + 4] - fVerticesY[i + 4];
0155     }
0156
0157     // Check that opposite segments are not crossing -> fatal exception
0158     if (SegmentsCrossing(fVertices[0], fVertices[1], fVertices[3], fVertices[2]) ||
0159         SegmentsCrossing(fVertices[1], fVertices[2], fVertices[0], fVertices[3]) ||
0160         SegmentsCrossing(fVertices[4], fVertices[5], fVertices[7], fVertices[6]) ||
0161         SegmentsCrossing(fVertices[5], fVertices[6], fVertices[4], fVertices[7])) {
0162       printf("ERROR: Unplaced generic trap defined with crossing opposite segments\n");
0163       Print();
0164       return false;
0165     }
0166
0167     // Check that top and bottom quadrilaterals are convex
0168     if (!ComputeIsConvexQuadrilaterals()) {
0169       printf("ERROR: Unplaced generic trap defined with top/bottom quadrilaterals not convex\n");
0170       Print();
0171       return false;
0172     }
0173
0174     fIsTwisted = ComputeIsTwisted();
0175
0176 #ifndef VECCORE_CUDA
0177     // Create the tessellated helper if  the faces are planar
0178     if (IsPlanar()) {
0179       fTslHelper = new TessellatedSection<T>(4, -fDz, fDz);
0180       fTslHelper->AddQuadrilateralFacet(fVertices[0], fVertices[4], fVertices[5], fVertices[1]);
0181       fTslHelper->AddQuadrilateralFacet(fVertices[1], fVertices[5], fVertices[6], fVertices[2]);
0182       fTslHelper->AddQuadrilateralFacet(fVertices[2], fVertices[6], fVertices[7], fVertices[3]);
0183       fTslHelper->AddQuadrilateralFacet(fVertices[3], fVertices[7], fVertices[4], fVertices[0]);
0184     }
0185 #endif
0186     ComputeBoundingBox();
0187     return true;
0188   }
0189
0190   VECCORE_ATT_HOST_DEVICE
0191   VECGEOM_FORCE_INLINE
0192   bool IsPlanar() const { return (!fIsTwisted); }
0193
0194   VECCORE_ATT_HOST_DEVICE
0195   void ComputeBoundingBox()
0196   {
0197     // Computes bounding box parameters
0198     Vertex_t aMin, aMax;
0199     Extent(aMin, aMax);
0200     fBBorigin     = 0.5 * (aMin + aMax);
0201     fBBdimensions = 0.5 * (aMax - aMin);
0202   }
0203
0204   VECCORE_ATT_HOST_DEVICE
0205   void Extent(Vertex_t &aMin, Vertex_t &aMax) const
0206   {
0207     // Returns the full 3D cartesian extent of the solid.
0208     aMin = aMax = fVertices[0];
0209     aMin[2]     = -fDz;
0210     aMax[2]     = fDz;
0211     for (int i = 0; i < 4; ++i) {
0212       // lower -fDz vertices
0213       if (aMin[0] > fVertices[i].x()) aMin[0] = fVertices[i].x();
0214       if (aMax[0] < fVertices[i].x()) aMax[0] = fVertices[i].x();
0215       if (aMin[1] > fVertices[i].y()) aMin[1] = fVertices[i].y();
0216       if (aMax[1] < fVertices[i].y()) aMax[1] = fVertices[i].y();
0217       // upper fDz vertices
0218       if (aMin[0] > fVertices[i + 4].x()) aMin[0] = fVertices[i + 4].x();
0219       if (aMax[0] < fVertices[i + 4].x()) aMax[0] = fVertices[i + 4].x();
0220       if (aMin[1] > fVertices[i + 4].y()) aMin[1] = fVertices[i + 4].y();
0221       if (aMax[1] < fVertices[i + 4].y()) aMax[1] = fVertices[i + 4].y();
0222     }
0223   }
0224
0225   VECCORE_ATT_HOST_DEVICE
0226   bool ComputeIsConvexQuadrilaterals()
0227   {
0228     // Computes if this gentrap top and bottom quadrilaterals are convex. The vertices
0229     // have to be pre-ordered clockwise in the XY plane.
0230
0231     // The cross product of all vector pairs corresponding to ordered consecutive
0232     // segments has to be positive.
0233     for (int i = 0; i < 4; ++i) {
0234       int j = (i + 1) % 4;
0235       // Bottom face
0236       Precision crossij = fDeltaX[i] * fDeltaY[j] - fDeltaY[i] * fDeltaX[j];
0237       if (crossij > kTolerance) return false;
0238       // Top face
0239       crossij = fDeltaX[i + 4] * fDeltaY[j + 4] - fDeltaY[i + 4] * fDeltaX[j + 4];
0240       if (crossij > kTolerance) return false;
0241     }
0242     return true;
0243   }
0244
0245   VECCORE_ATT_HOST_DEVICE
0246   bool ComputeIsTwisted()
0247   {
0248     // Check if the trapezoid is twisted. A lateral face is twisted if the top and
0249     // bottom segments are not parallel (cross product not null)
0250
0251     bool twisted = false;
0252     Precision dx1, dy1, dx2, dy2;
0253     const int nv = 4; // half the number of verices
0254
0255     for (int i = 0; i < 4; ++i) {
0256       dx1 = fVertices[(i + 1) % nv].x() - fVertices[i].x();
0257       dy1 = fVertices[(i + 1) % nv].y() - fVertices[i].y();
0258       if ((dx1 == 0) && (dy1 == 0)) {
0259         continue;
0260       }
0261
0262       dx2 = fVertices[nv + (i + 1) % nv].x() - fVertices[nv + i].x();
0263       dy2 = fVertices[nv + (i + 1) % nv].y() - fVertices[nv + i].y();
0264
0265       if ((dx2 == 0 && dy2 == 0)) {
0266         continue;
0267       }
0268       fTwist[i] = std::fabs(dy1 * dx2 - dx1 * dy2);
0269       if (fTwist[i] < kTolerance) {
0270         fTwist[i] = 0.;
0271         continue;
0272       }
0273       twisted = true;
0274     }
0275     return twisted;
0276   }
0277
0278   VECCORE_ATT_HOST_DEVICE
0279   void Print() const
0280   {
0281     printf("--------------------------------------------------------\n");
0282     printf("    =================================================== \n");
0283     printf(" Solid type: UnplacedGenTrap \n");
0284     printf("   half length Z: %f mm \n", fDz);
0285     printf("   list of vertices:\n");
0286
0287     for (int i = 0; i < 8; ++i) {
0288       printf("#%d", i);
0289       printf("   vx = %f mm", fVertices[i].x());
0290       printf("   vy = %f mm\n", fVertices[i].y());
0291     }
0292     printf("   planar: %s\n", IsPlanar() ? "true" : "false");
0293   }
0294
0295   VECCORE_ATT_HOST_DEVICE
0296   bool SegmentsCrossing(Vertex_t p, Vertex_t p1, Vertex_t q, Vertex_t q1) const
0297   {
0298     // Check if 2 segments defined by (p,p1) and (q,q1) are crossing.
0299     // See: http://stackoverflow.com/questions/563198/how-do-you-detect-where-two-line-segments-intersect
0300     using Vector     = Vertex_t;
0301     Vector r         = p1 - p; // p1 = p+r
0302     Vector s         = q1 - q; // q1 = q+s
0303     Vector r_cross_s = Vector::Cross(r, s);
0304     if (r_cross_s.Mag2() < kTolerance) // parallel, colinear or degenerated - ignore crossing
0305       return false;
0306     // The segments are crossing if the vector equations:
0307     //   t * (r x s) = (q-p) x s
0308     //   u * (r x s) = (q-p) x r
0309     // give t and u values in the range (0,1).
0310     // We allow crossing within kTolerance, so the interval becomes (kTolerance, 1-kTolerance)
0311     Precision t = Vector::Cross(q - p, s).z() / r_cross_s.z();
0312     if (t < kTolerance || t > 1. - kTolerance) return false;
0313     Precision u = Vector::Cross(q - p, r).z() / r_cross_s.z();
0314     if (u < kTolerance || u > 1. - kTolerance) return false;
0315     return true;
0316   }
0317 };
0318 } // namespace VECGEOM_IMPL_NAMESPACE
0319 } // namespace vecgeom
0320
0321 #endif