VecGeom/volumes/ExtrudedStruct.h

0001 /// @file ExtrudedStruct.h
0002 /// @author Mihaela Gheata (mihaela.gheata@cern.ch)
0003
0004 #ifndef VECGEOM_EXTRUDED_STRUCT_H
0005 #define VECGEOM_EXTRUDED_STRUCT_H
0006
0007 #include "VecGeom/base/Config.h"
0008
0009 #include "VecGeom/volumes/PolygonalShell.h"
0010 #include "VecGeom/volumes/TessellatedStruct.h"
0011
0012 #ifndef VECGEOM_ENABLE_CUDA
0013 #include "VecGeom/volumes/TessellatedSection.h"
0014 #endif
0015
0016 namespace vecgeom {
0017
0018 VECGEOM_DEVICE_FORWARD_DECLARE(class ExtrudedStruct;);
0019 VECGEOM_DEVICE_DECLARE_CONV(class, ExtrudedStruct);
0020 VECGEOM_DEVICE_DECLARE_CONV(struct, XtruVertex2);
0021 VECGEOM_DEVICE_DECLARE_CONV(struct, XtruSection);
0022
0023 inline namespace VECGEOM_IMPL_NAMESPACE {
0024
0025 // Structure wrapping either a polygonal shell helper in case of two
0026 // extruded sections or a tessellated structure in case of more
0027
0028 struct XtruVertex2 {
0029   Precision x;
0030   Precision y;
0031 };
0032
0033 struct XtruSection {
0034   Vector3D<Precision> fOrigin; // Origin of the section
0035   Precision fScale;
0036 };
0037
0038 class ExtrudedStruct {
0039
0040   // template <typename U>
0041   // using vector_t = vecgeom::Vector<U>;
0042   template <typename U>
0043   using vector_t = vecgeom::Vector<U>;
0044
0045 public:
0046   bool fIsSxtru                  = false;     ///< Flag for sxtru representation
0047   bool fInitialized              = false;     ///< Flag for initialization
0048   Precision *fZPlanes            = nullptr;   ///< Z position of planes
0049   mutable Precision fCubicVolume = 0.;        ///< Cubic volume
0050   mutable Precision fSurfaceArea = 0.;        ///< Surface area
0051   PolygonalShell fSxtruHelper;                ///< Sxtru helper
0052   TessellatedStruct<3, Precision> fTslHelper; ///< Tessellated helper
0053 #ifndef VECGEOM_ENABLE_CUDA
0054   bool fUseTslSections = false;                           ///< Use tessellated section helper
0055   vector_t<TessellatedSection<Precision> *> fTslSections; ///< Tessellated sections
0056 #endif
0057   vector_t<XtruVertex2> fVertices; ///< Polygone vertices
0058   vector_t<XtruSection> fSections; ///< Vector of sections
0059   PlanarPolygon fPolygon;          ///< Planar polygon
0060
0061 public:
0062   /** @brief Dummy constructor */
0063   VECCORE_ATT_HOST_DEVICE
0064   ExtrudedStruct() {}
0065
0066   /** @brief Constructor providing polygone vertices and sections */
0067   VECCORE_ATT_HOST_DEVICE
0068   ExtrudedStruct(int nvertices, XtruVertex2 const *vertices, int nsections, XtruSection const *sections)
0069   {
0070     Initialize(nvertices, vertices, nsections, sections);
0071   }
0072
0073   // Constructor used during Specialization for nsections == 2
0074   VECCORE_ATT_HOST_DEVICE
0075   ExtrudedStruct(size_t nvertices, const Precision *x, const Precision *y, Precision zmin, Precision zmax)
0076   {
0077     XtruVertex2 *vertices = new XtruVertex2[nvertices];
0078     XtruSection *sections = new XtruSection[2];
0079     for (size_t i = 0; i < nvertices; ++i) {
0080       vertices[i].x = x[i];
0081       vertices[i].y = y[i];
0082     }
0083
0084     sections[0].fScale = 1.;
0085     sections[0].fOrigin.Set(0., 0., zmin);
0086
0087     sections[1].fScale = 1.;
0088     sections[1].fOrigin.Set(0., 0., zmax);
0089
0090     Initialize(nvertices, vertices, 2, sections);
0091     delete[] vertices;
0092     delete[] sections;
0093   }
0094
0095   VECGEOM_FORCE_INLINE
0096   VECCORE_ATT_HOST_DEVICE
0097   int FindZSegment(Precision const &pointZ) const
0098   {
0099     int index              = -1;
0100     Precision const *begin = fZPlanes;
0101     Precision const *end   = fZPlanes + fSections.size() + 1;
0102     while (begin < end - 1 && pointZ - kTolerance > *begin) {
0103       ++index;
0104       ++begin;
0105     }
0106     if (pointZ + kTolerance > *begin) return (index + 1);
0107     return index;
0108   }
0109
0110   /** @brief Initialize based on vertices and sections */
0111   void Initialize(int nvertices, XtruVertex2 const *vertices, int nsections, XtruSection const *sections)
0112   {
0113     if (fInitialized) return;
0114     VECGEOM_ASSERT(nsections > 1 && nvertices > 2);
0115     fZPlanes         = new Precision[nsections];
0116     fZPlanes[0]      = sections[0].fOrigin.z();
0117     bool degenerated = false;
0118     for (size_t i = 1; i < (size_t)nsections; ++i) {
0119       fZPlanes[i] = sections[i].fOrigin.z();
0120       // Make sure sections are defined in increasing order
0121       VECGEOM_VALIDATE(fZPlanes[i] >= fZPlanes[i - 1], << "Extruded sections not defined in increasing Z order");
0122       if (fZPlanes[i] - fZPlanes[i - 1] < kTolerance) degenerated = true;
0123     }
0124 #ifndef VECGEOM_ENABLE_CUDA
0125     if (!degenerated) fUseTslSections = true;
0126 #endif
0127     (void)degenerated; // silence the compiler
0128     // Check if this is an SXtru
0129     if (nsections == 2 && (sections[0].fOrigin - sections[1].fOrigin).Perp2() < kTolerance &&
0130         vecCore::math::Abs(sections[0].fScale - sections[1].fScale) < kTolerance)
0131       fIsSxtru = true;
0132     if (fIsSxtru) {
0133       // Put vertices in arrays
0134       Precision *x = new Precision[nvertices];
0135       Precision *y = new Precision[nvertices];
0136       for (size_t i = 0; i < (size_t)nvertices; ++i) {
0137         x[i] = sections[0].fOrigin.x() + sections[0].fScale * vertices[i].x;
0138         y[i] = sections[0].fOrigin.y() + sections[0].fScale * vertices[i].y;
0139       }
0140       fSxtruHelper.Init(nvertices, x, y, sections[0].fOrigin[2], sections[1].fOrigin[2]);
0141       delete[] x;
0142       delete[] y;
0143     }
0144     // Create the tessellated structure in all cases
0145     CreateTessellated(nvertices, vertices, nsections, sections);
0146     fInitialized = true;
0147   }
0148
0149   /** @brief Construct facets based on vertices and sections */
0150   VECCORE_ATT_HOST_DEVICE
0151   void CreateTessellated(size_t nvertices, XtruVertex2 const *vertices, size_t nsections, XtruSection const *sections)
0152   {
0153     struct FacetInd {
0154       size_t ind1{0}, ind2{0}, ind3{0};
0155
0156       FacetInd() = default;
0157       FacetInd(int i1, int i2, int i3)
0158       {
0159         ind1 = i1;
0160         ind2 = i2;
0161         ind3 = i3;
0162       }
0163     };
0164
0165     // Store sections
0166     for (size_t isect = 0; isect < nsections; ++isect)
0167       fSections.push_back(sections[isect]);
0168
0169     // Create the polygon
0170     Precision *vx = new Precision[nvertices];
0171     Precision *vy = new Precision[nvertices];
0172     for (size_t i = 0; i < nvertices; ++i) {
0173       vx[i] = vertices[i].x;
0174       vy[i] = vertices[i].y;
0175     }
0176     fPolygon.Init(nvertices, vx, vy);
0177 #ifndef VECGEOM_ENABLE_CUDA
0178     fUseTslSections &= fPolygon.IsConvex();
0179     if (fUseTslSections) {
0180       // Create tessellated sections
0181       fTslSections.reserve(nsections);
0182       for (size_t i = 0; i < nsections - 1; ++i) {
0183         fTslSections[i] =
0184             new TessellatedSection<Precision>(nvertices, sections[i].fOrigin.z(), sections[i + 1].fOrigin.z());
0185       }
0186     }
0187 #endif
0188     // TRIANGULATE POLYGON
0189
0190     VectorBase<FacetInd> facets(nvertices);
0191     // Fill a vector of vertex indices
0192     vector_t<size_t> vtx;
0193     for (size_t i = 0; i < nvertices; ++i)
0194       vtx.push_back(i);
0195
0196     size_t i1 = 0;
0197     size_t i2 = 1;
0198     size_t i3 = 2;
0199
0200     while (vtx.size() > 2) {
0201       // Find convex parts of the polygon (ears)
0202       size_t counter = 0;
0203       while (!IsConvexSide(vtx[i1], vtx[i2], vtx[i3])) {
0204         i1++;
0205         i2++;
0206         i3 = (i3 + 1) % vtx.size();
0207         counter++;
0208         VECGEOM_VALIDATE(counter < nvertices, << "Triangulation failed");
0209         (void)counter; // silence unused variable warnings in release builds
0210       }
0211       bool good = true;
0212       // Check if any of the remaining vertices are in the ear
0213       for (auto i : vtx) {
0214         if (i == vtx[i1] || i == vtx[i2] || i == vtx[i3]) continue;
0215         if (IsPointInside(vtx[i1], vtx[i2], vtx[i3], i)) {
0216           good = false;
0217           i1++;
0218           i2++;
0219           i3 = (i3 + 1) % vtx.size();
0220           break;
0221         }
0222       }
0223
0224       if (good) {
0225         // Make triangle
0226         facets.push_back(FacetInd(vtx[i1], vtx[i2], vtx[i3]));
0227         // Remove the middle vertex of the ear and restart
0228         vtx.erase(vtx.begin() + i2);
0229         i1 = 0;
0230         i2 = 1;
0231         i3 = 2;
0232       }
0233     }
0234     // We have all index facets, create now the real facets
0235     // Bottom (normals pointing down)
0236     for (size_t i = 0; i < facets.size(); ++i) {
0237       i1 = facets[i].ind1;
0238       i2 = facets[i].ind2;
0239       i3 = facets[i].ind3;
0240       fTslHelper.AddTriangularFacet(VertexToSection(i1, 0), VertexToSection(i2, 0), VertexToSection(i3, 0));
0241     }
0242     // Sections
0243     for (size_t isect = 0; isect < nsections - 1; ++isect) {
0244       for (size_t i = 0; i < (size_t)nvertices; ++i) {
0245         size_t j = (i + 1) % nvertices;
0246         // Quadrilateral isect:(j, i)  isect+1: (i, j)
0247         fTslHelper.AddQuadrilateralFacet(VertexToSection(j, isect), VertexToSection(i, isect),
0248                                          VertexToSection(i, isect + 1), VertexToSection(j, isect + 1));
0249 #ifndef VECGEOM_ENABLE_CUDA
0250         if (fUseTslSections)
0251           fTslSections[isect]->AddQuadrilateralFacet(VertexToSection(j, isect), VertexToSection(i, isect),
0252                                                      VertexToSection(i, isect + 1), VertexToSection(j, isect + 1));
0253 #endif
0254       }
0255     }
0256     // Top (normals pointing up)
0257     for (size_t i = 0; i < facets.size(); ++i) {
0258       i1 = facets[i].ind1;
0259       i2 = facets[i].ind2;
0260       i3 = facets[i].ind3;
0261       fTslHelper.AddTriangularFacet(VertexToSection(i1, nsections - 1), VertexToSection(i3, nsections - 1),
0262                                     VertexToSection(i2, nsections - 1));
0263     }
0264     // Now close the tessellated structure
0265     fTslHelper.Close();
0266   }
0267
0268   /** @brief Get the number of sections */
0269   VECCORE_ATT_HOST_DEVICE
0270   VECGEOM_FORCE_INLINE
0271   size_t GetNSections() const { return fSections.size(); }
0272
0273   /** @brief Get the number of planes */
0274   VECCORE_ATT_HOST_DEVICE
0275   VECGEOM_FORCE_INLINE
0276   size_t GetNSegments() const { return (fSections.size() - 1); }
0277
0278   /** @brief Get section i */
0279   VECCORE_ATT_HOST_DEVICE
0280   VECGEOM_FORCE_INLINE
0281   XtruSection GetSection(int i) const { return fSections[i]; }
0282
0283   /** @brief Get the number of vertices */
0284   VECCORE_ATT_HOST_DEVICE
0285   VECGEOM_FORCE_INLINE
0286   size_t GetNVertices() const { return fPolygon.GetNVertices(); }
0287
0288   /** @brief Get the polygone vertex i */
0289   VECCORE_ATT_HOST_DEVICE
0290   VECGEOM_FORCE_INLINE
0291   void GetVertex(int i, Precision &x, Precision &y) const
0292   {
0293     x = fPolygon.GetVertices().x()[i];
0294     y = fPolygon.GetVertices().y()[i];
0295   }
0296
0297   /** Return true if i is on the line through i1, i2 */
0298   VECCORE_ATT_HOST_DEVICE
0299   VECGEOM_FORCE_INLINE
0300   bool IsSameLine(size_t i, size_t i1, size_t i2) const
0301   {
0302     const Precision *x = fPolygon.GetVertices().x();
0303     const Precision *y = fPolygon.GetVertices().y();
0304     if (x[i1] == x[i2]) return std::fabs(x[i] - x[i1]) < kTolerance * 0.5;
0305
0306     Precision slope = (y[i2] - y[i1]) / (x[i2] - x[i1]);
0307     Precision predy = y[i1] + slope * (x[i] - x[i1]);
0308     Precision dy    = y[i] - predy;
0309
0310     // Check perpendicular distance vs tolerance 'directly'
0311     const Precision tol = 0.5 * kTolerance;
0312     bool squareComp     = (dy * dy < (1 + slope * slope) * tol * tol);
0313     return squareComp;
0314   }
0315
0316   /** @brief Return true if point i is on the line through i1, i2 and lies between i1 and i2 */
0317   VECCORE_ATT_HOST_DEVICE
0318   VECGEOM_FORCE_INLINE
0319   bool IsSameLineSegment(size_t i, size_t i1, size_t i2) const
0320   {
0321     const Precision *x = fPolygon.GetVertices().x();
0322     const Precision *y = fPolygon.GetVertices().y();
0323     if (x[i] < std::min(x[i1], x[i2]) - kTolerance * 0.5 || x[i] > std::max(x[i1], x[i2]) + kTolerance * 0.5 ||
0324         y[i] < std::min(y[i1], y[i2]) - kTolerance * 0.5 || y[i] > std::max(y[i1], y[i2]) + kTolerance * 0.5)
0325       return false;
0326
0327     return IsSameLine(i, i1, i2);
0328   }
0329
0330   /** @brief Return true if i and j are on the same side of the line through i1, i2 */
0331   VECCORE_ATT_HOST_DEVICE
0332   VECGEOM_FORCE_INLINE
0333   bool IsSameSide(size_t i, size_t j, size_t i1, size_t i2) const
0334   {
0335     const Precision *x = fPolygon.GetVertices().x();
0336     const Precision *y = fPolygon.GetVertices().y();
0337
0338     return ((x[i] - x[i1]) * (y[i2] - y[i1]) - (x[i2] - x[i1]) * (y[i] - y[i1])) *
0339                ((x[j] - x[i1]) * (y[i2] - y[i1]) - (x[i2] - x[i1]) * (y[j] - y[i1])) >
0340            0;
0341   }
0342
0343   /** Return true if i is inside of triangle (i1, i2, i3) or on its edges, else returns false */
0344   VECCORE_ATT_HOST_DEVICE
0345   VECGEOM_FORCE_INLINE
0346   bool IsPointInside(size_t i1, size_t i2, size_t i3, size_t i) const
0347   {
0348     const Precision *x = fPolygon.GetVertices().x();
0349     const Precision *y = fPolygon.GetVertices().y();
0350
0351     // Check extent first
0352     if ((x[i] < x[i1] && x[i] < x[i2] && x[i] < x[i3]) || (x[i] > x[i1] && x[i] > x[i2] && x[i] > x[i3]) ||
0353         (y[i] < y[i1] && y[i] < y[i2] && y[i] < y[i3]) || (y[i] > y[i1] && y[i] > y[i2] && y[i] > y[i3]))
0354       return false;
0355
0356     bool inside = IsSameSide(i, i1, i2, i3) && IsSameSide(i, i2, i1, i3) && IsSameSide(i, i3, i1, i2);
0357
0358     bool onEdge = IsSameLineSegment(i, i1, i2) || IsSameLineSegment(i, i2, i3) || IsSameLineSegment(i, i3, i1);
0359
0360     return inside || onEdge;
0361   }
0362
0363   /** @brief Check if the polygone segments (i0, i1) and (i1, i2) make a convex side */
0364   VECCORE_ATT_HOST_DEVICE
0365   VECGEOM_FORCE_INLINE
0366   bool IsConvexSide(size_t i0, size_t i1, size_t i2)
0367   {
0368     const Precision *x = fPolygon.GetVertices().x();
0369     const Precision *y = fPolygon.GetVertices().y();
0370     Precision cross    = (x[i1] - x[i0]) * (y[i2] - y[i1]) - (x[i2] - x[i1]) * (y[i1] - y[i0]);
0371     return cross < 0.;
0372   }
0373
0374   /** @brief Returns convexity of polygon */
0375   VECCORE_ATT_HOST_DEVICE
0376   VECGEOM_FORCE_INLINE
0377   bool IsConvexPolygon() const { return fPolygon.IsConvex(); }
0378
0379   /** @brief Returns the coordinates for a given vertex index at a given section */
0380   VECCORE_ATT_HOST_DEVICE
0381   VECGEOM_FORCE_INLINE
0382   Vector3D<Precision> VertexToSection(size_t ivert, size_t isect) const
0383   {
0384     const Precision *x = fPolygon.GetVertices().x();
0385     const Precision *y = fPolygon.GetVertices().y();
0386     Vector3D<Precision> vert(fSections[isect].fOrigin[0] + fSections[isect].fScale * x[ivert],
0387                              fSections[isect].fOrigin[1] + fSections[isect].fScale * y[ivert],
0388                              fSections[isect].fOrigin[2]);
0389     return vert;
0390   }
0391 };
0392
0393 } // namespace VECGEOM_IMPL_NAMESPACE
0394 } // namespace vecgeom
0395
0396 #endif