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     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       assert(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, ind2, ind3;
0155       FacetInd(int i1, int i2, int i3)
0156       {
0157         ind1 = i1;
0158         ind2 = i2;
0159         ind3 = i3;
0160       }
0161     };
0162
0163     // Store sections
0164     for (size_t isect = 0; isect < nsections; ++isect)
0165       fSections.push_back(sections[isect]);
0166
0167     // Create the polygon
0168     Precision *vx = new Precision[nvertices];
0169     Precision *vy = new Precision[nvertices];
0170     for (size_t i = 0; i < nvertices; ++i) {
0171       vx[i] = vertices[i].x;
0172       vy[i] = vertices[i].y;
0173     }
0174     fPolygon.Init(nvertices, vx, vy);
0175 #ifndef VECGEOM_ENABLE_CUDA
0176     fUseTslSections &= fPolygon.IsConvex();
0177     if (fUseTslSections) {
0178       // Create tessellated sections
0179       fTslSections.reserve(nsections);
0180       for (size_t i = 0; i < nsections - 1; ++i) {
0181         fTslSections[i] =
0182             new TessellatedSection<Precision>(nvertices, sections[i].fOrigin.z(), sections[i + 1].fOrigin.z());
0183       }
0184     }
0185 #endif
0186     // TRIANGULATE POLYGON
0187
0188     VectorBase<FacetInd> facets(nvertices);
0189     // Fill a vector of vertex indices
0190     vector_t<size_t> vtx;
0191     for (size_t i = 0; i < nvertices; ++i)
0192       vtx.push_back(i);
0193
0194     size_t i1 = 0;
0195     size_t i2 = 1;
0196     size_t i3 = 2;
0197
0198     while (vtx.size() > 2) {
0199       // Find convex parts of the polygon (ears)
0200       size_t counter = 0;
0201       while (!IsConvexSide(vtx[i1], vtx[i2], vtx[i3])) {
0202         i1++;
0203         i2++;
0204         i3 = (i3 + 1) % vtx.size();
0205         counter++;
0206         assert(counter < nvertices && "Triangulation failed");
0207         (void)counter; // silence unused variable warnings in release builds
0208       }
0209       bool good = true;
0210       // Check if any of the remaining vertices are in the ear
0211       for (auto i : vtx) {
0212         if (i == vtx[i1] || i == vtx[i2] || i == vtx[i3]) continue;
0213         if (IsPointInside(vtx[i1], vtx[i2], vtx[i3], i)) {
0214           good = false;
0215           i1++;
0216           i2++;
0217           i3 = (i3 + 1) % vtx.size();
0218           break;
0219         }
0220       }
0221
0222       if (good) {
0223         // Make triangle
0224         facets.push_back(FacetInd(vtx[i1], vtx[i2], vtx[i3]));
0225         // Remove the middle vertex of the ear and restart
0226         vtx.erase(vtx.begin() + i2);
0227         i1 = 0;
0228         i2 = 1;
0229         i3 = 2;
0230       }
0231     }
0232     // We have all index facets, create now the real facets
0233     // Bottom (normals pointing down)
0234     for (size_t i = 0; i < facets.size(); ++i) {
0235       i1 = facets[i].ind1;
0236       i2 = facets[i].ind2;
0237       i3 = facets[i].ind3;
0238       fTslHelper.AddTriangularFacet(VertexToSection(i1, 0), VertexToSection(i2, 0), VertexToSection(i3, 0));
0239     }
0240     // Sections
0241     for (size_t isect = 0; isect < nsections - 1; ++isect) {
0242       for (size_t i = 0; i < (size_t)nvertices; ++i) {
0243         size_t j = (i + 1) % nvertices;
0244         // Quadrilateral isect:(j, i)  isect+1: (i, j)
0245         fTslHelper.AddQuadrilateralFacet(VertexToSection(j, isect), VertexToSection(i, isect),
0246                                          VertexToSection(i, isect + 1), VertexToSection(j, isect + 1));
0247 #ifndef VECGEOM_ENABLE_CUDA
0248         if (fUseTslSections)
0249           fTslSections[isect]->AddQuadrilateralFacet(VertexToSection(j, isect), VertexToSection(i, isect),
0250                                                      VertexToSection(i, isect + 1), VertexToSection(j, isect + 1));
0251 #endif
0252       }
0253     }
0254     // Top (normals pointing up)
0255     for (size_t i = 0; i < facets.size(); ++i) {
0256       i1 = facets[i].ind1;
0257       i2 = facets[i].ind2;
0258       i3 = facets[i].ind3;
0259       fTslHelper.AddTriangularFacet(VertexToSection(i1, nsections - 1), VertexToSection(i3, nsections - 1),
0260                                     VertexToSection(i2, nsections - 1));
0261     }
0262     // Now close the tessellated structure
0263     fTslHelper.Close();
0264   }
0265
0266   /** @brief Get the number of sections */
0267   VECCORE_ATT_HOST_DEVICE
0268   VECGEOM_FORCE_INLINE
0269   size_t GetNSections() const { return fSections.size(); }
0270
0271   /** @brief Get the number of planes */
0272   VECCORE_ATT_HOST_DEVICE
0273   VECGEOM_FORCE_INLINE
0274   size_t GetNSegments() const { return (fSections.size() - 1); }
0275
0276   /** @brief Get section i */
0277   VECCORE_ATT_HOST_DEVICE
0278   VECGEOM_FORCE_INLINE
0279   XtruSection GetSection(int i) const { return fSections[i]; }
0280
0281   /** @brief Get the number of vertices */
0282   VECCORE_ATT_HOST_DEVICE
0283   VECGEOM_FORCE_INLINE
0284   size_t GetNVertices() const { return fPolygon.GetNVertices(); }
0285
0286   /** @brief Get the polygone vertex i */
0287   VECCORE_ATT_HOST_DEVICE
0288   VECGEOM_FORCE_INLINE
0289   void GetVertex(int i, Precision &x, Precision &y) const
0290   {
0291     x = fPolygon.GetVertices().x()[i];
0292     y = fPolygon.GetVertices().y()[i];
0293   }
0294
0295   /** Return true if i is on the line through i1, i2 */
0296   VECCORE_ATT_HOST_DEVICE
0297   VECGEOM_FORCE_INLINE
0298   bool IsSameLine(size_t i, size_t i1, size_t i2) const
0299   {
0300     const Precision *x = fPolygon.GetVertices().x();
0301     const Precision *y = fPolygon.GetVertices().y();
0302     if (x[i1] == x[i2]) return std::fabs(x[i] - x[i1]) < kTolerance * 0.5;
0303
0304     Precision slope = (y[i2] - y[i1]) / (x[i2] - x[i1]);
0305     Precision predy = y[i1] + slope * (x[i] - x[i1]);
0306     Precision dy    = y[i] - predy;
0307
0308     // Check perpendicular distance vs tolerance 'directly'
0309     const Precision tol = 0.5 * kTolerance;
0310     bool squareComp     = (dy * dy < (1 + slope * slope) * tol * tol);
0311     return squareComp;
0312   }
0313
0314   /** @brief Return true if point i is on the line through i1, i2 and lies between i1 and i2 */
0315   VECCORE_ATT_HOST_DEVICE
0316   VECGEOM_FORCE_INLINE
0317   bool IsSameLineSegment(size_t i, size_t i1, size_t i2) const
0318   {
0319     const Precision *x = fPolygon.GetVertices().x();
0320     const Precision *y = fPolygon.GetVertices().y();
0321     if (x[i] < std::min(x[i1], x[i2]) - kTolerance * 0.5 || x[i] > std::max(x[i1], x[i2]) + kTolerance * 0.5 ||
0322         y[i] < std::min(y[i1], y[i2]) - kTolerance * 0.5 || y[i] > std::max(y[i1], y[i2]) + kTolerance * 0.5)
0323       return false;
0324
0325     return IsSameLine(i, i1, i2);
0326   }
0327
0328   /** @brief Return true if i and j are on the same side of the line through i1, i2 */
0329   VECCORE_ATT_HOST_DEVICE
0330   VECGEOM_FORCE_INLINE
0331   bool IsSameSide(size_t i, size_t j, size_t i1, size_t i2) const
0332   {
0333     const Precision *x = fPolygon.GetVertices().x();
0334     const Precision *y = fPolygon.GetVertices().y();
0335
0336     return ((x[i] - x[i1]) * (y[i2] - y[i1]) - (x[i2] - x[i1]) * (y[i] - y[i1])) *
0337                ((x[j] - x[i1]) * (y[i2] - y[i1]) - (x[i2] - x[i1]) * (y[j] - y[i1])) >
0338            0;
0339   }
0340
0341   /** Return true if i is inside of triangle (i1, i2, i3) or on its edges, else returns false */
0342   VECCORE_ATT_HOST_DEVICE
0343   VECGEOM_FORCE_INLINE
0344   bool IsPointInside(size_t i1, size_t i2, size_t i3, size_t i) const
0345   {
0346     const Precision *x = fPolygon.GetVertices().x();
0347     const Precision *y = fPolygon.GetVertices().y();
0348
0349     // Check extent first
0350     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]) ||
0351         (y[i] < y[i1] && y[i] < y[i2] && y[i] < y[i3]) || (y[i] > y[i1] && y[i] > y[i2] && y[i] > y[i3]))
0352       return false;
0353
0354     bool inside = IsSameSide(i, i1, i2, i3) && IsSameSide(i, i2, i1, i3) && IsSameSide(i, i3, i1, i2);
0355
0356     bool onEdge = IsSameLineSegment(i, i1, i2) || IsSameLineSegment(i, i2, i3) || IsSameLineSegment(i, i3, i1);
0357
0358     return inside || onEdge;
0359   }
0360
0361   /** @brief Check if the polygone segments (i0, i1) and (i1, i2) make a convex side */
0362   VECCORE_ATT_HOST_DEVICE
0363   VECGEOM_FORCE_INLINE
0364   bool IsConvexSide(size_t i0, size_t i1, size_t i2)
0365   {
0366     const Precision *x = fPolygon.GetVertices().x();
0367     const Precision *y = fPolygon.GetVertices().y();
0368     Precision cross    = (x[i1] - x[i0]) * (y[i2] - y[i1]) - (x[i2] - x[i1]) * (y[i1] - y[i0]);
0369     return cross < 0.;
0370   }
0371
0372   /** @brief Returns convexity of polygon */
0373   VECCORE_ATT_HOST_DEVICE
0374   VECGEOM_FORCE_INLINE
0375   bool IsConvexPolygon() const { return fPolygon.IsConvex(); }
0376
0377   /** @brief Returns the coordinates for a given vertex index at a given section */
0378   VECCORE_ATT_HOST_DEVICE
0379   VECGEOM_FORCE_INLINE
0380   Vector3D<Precision> VertexToSection(size_t ivert, size_t isect) const
0381   {
0382     const Precision *x = fPolygon.GetVertices().x();
0383     const Precision *y = fPolygon.GetVertices().y();
0384     Vector3D<Precision> vert(fSections[isect].fOrigin[0] + fSections[isect].fScale * x[ivert],
0385                              fSections[isect].fOrigin[1] + fSections[isect].fScale * y[ivert],
0386                              fSections[isect].fOrigin[2]);
0387     return vert;
0388   }
0389 };
0390
0391 } // namespace VECGEOM_IMPL_NAMESPACE
0392 } // namespace vecgeom
0393
0394 #endif