surfaces/conv/PolyhedronConverter.h

0001 #ifndef VECGEOM_SURFACE_POLYHEDRONCONVERTER_H_
0002 #define VECGEOM_SURFACE_POLYHEDRONCONVERTER_H_
0003
0004 #include <VecGeom/surfaces/conv/ConvHelper.h>
0005
0006 #include <VecGeom/volumes/Polyhedron.h>
0007
0008 namespace vgbrep {
0009 namespace conv {
0010
0011 /// @brief Converter for Polyhedron
0012 /// @tparam Real_t Precision type
0013 /// @param upoly Polyhedron solid to be converted
0014 /// @param logical_id Id of the logical volume
0015 /// @return Conversion success
0016 template <typename Real_t>
0017 bool CreatePolyhedronSurfaces(vecgeom::UnplacedPolyhedron const &upoly, int logical_id, bool intersection = false)
0018 {
0019   using Vector3D = vecgeom::Vector3D<vecgeom::Precision>;
0020
0021   int isurf;
0022   std::vector<Vector3D> vertices;
0023   int isurfZlast = -1;
0024   LogicExpressionCPU logic; // OR logic: section0 | section1 | ... | sectionN
0025   auto const &poly    = upoly.GetStruct();
0026   size_t nplanes      = poly.fZPlanes.size();
0027   size_t nseg         = nplanes - 1;
0028   size_t sideCount    = poly.fSideCount;
0029   auto phiStart       = poly.fPhiStart;
0030   auto phiDelta       = poly.fPhiDelta;
0031   bool smallerPi      = phiDelta < (vecgeom::kPi - vecgeom::kTolerance);
0032   auto const &rMin    = poly.fRMin;
0033   auto const &rMax    = poly.fRMax;
0034   auto const &zPlanes = poly.fZPlanes;
0035   vecgeom::Precision csphi, ssphi, cephi, sephi;
0036
0037   if (sideCount == 1) VECGEOM_VALIDATE(smallerPi, << "Polyhedron with one segment cannot have angle larger than pi!");
0038
0039   auto sidePhi            = phiDelta / sideCount;
0040   auto cosHalfDeltaPhi    = vecCore::math::Cos(0.5 * sidePhi);
0041   vecgeom::Precision conv = 1. / cosHalfDeltaPhi;
0042
0043   // lambda to get the phi angle of the current side
0044   auto getPhi = [&](size_t side) {
0045     if (!poly.fHasPhiCutout && side == sideCount) {
0046       side = 0;
0047     }
0048     return vecgeom::NormalizeAngle<vecgeom::kScalar>(phiStart + side * sidePhi);
0049   };
0050
0051   auto getSinCos = [&](vecgeom::Precision sphi, vecgeom::Precision ephi) {
0052     csphi = vecCore::math::Cos(sphi);
0053     ssphi = vecCore::math::Sin(sphi);
0054     cephi = vecCore::math::Cos(ephi);
0055     sephi = vecCore::math::Sin(ephi);
0056   };
0057
0058   // lambda to convert a full Z segment
0059   auto convertZsegment = [&](size_t iseg) {
0060     auto dz             = zPlanes[iseg + 1] - zPlanes[iseg];
0061     auto realSeg        = dz > 0;
0062     auto drI            = rMin[iseg + 1] - rMin[iseg];
0063     auto drO            = rMax[iseg + 1] - rMax[iseg];
0064     bool hasInnerRadius = rMin[iseg] > 0 || rMin[iseg + 1] > 0;
0065     bool hasOuter       = realSeg || drO != 0;
0066     bool hasInner       = hasInnerRadius && (realSeg || drI != 0);
0067     bool hasPhi         = poly.fHasPhiCutout && realSeg;
0068     if (!hasOuter && !hasInner && !hasPhi) return;
0069
0070     auto z1 = zPlanes[iseg];
0071     auto z2 = zPlanes[iseg + 1];
0072
0073     if (realSeg) {
0074       if (iseg > 0) logic.push_back(lor); // `OR` between sections
0075       logic.push_back(lplus);             // '(' begin section logic
0076     }
0077
0078     // Add section Z planes as limiters
0079     // the logic is as follows: the first/last segments have a real surface on the bottom/top
0080     // only for the middle segments that are real, virtual Z planes as limiters are needed.
0081     isurf = isurfZlast;
0082     if (realSeg) {
0083       if (iseg > 0) {
0084         logic.push_back(lnot);
0085         logic.push_back(isurf);
0086       }
0087       if (iseg != (nseg - 1)) {
0088         isurf = builder::CreateLocalSurface<Real_t>(
0089             builder::CreateUnplacedSurface<Real_t>(SurfaceType::kPlanar), Frame{FrameType::kNoFrame},
0090             vecgeom::Transformation3DMP<Precision>(0., 0., zPlanes[iseg + 1], 0., 0., 0.));
0091         builder::AddSurfaceToShell<Real_t>(logical_id, isurf);
0092         if (intersection) builder::GetSurface<Real_t>(isurf).fSkipConvexity = true;
0093         if (iseg > 0) logic.push_back(land);
0094         logic.push_back(isurf);
0095         isurfZlast = isurf;
0096       }
0097     }
0098
0099     // outer surfaces
0100     if (hasOuter) {
0101       for (size_t iside = 0; iside < sideCount; iside++) {
0102         getSinCos(getPhi(iside), getPhi(iside + 1));
0103         // corners of the outer surfaces
0104         vertices = {{rMax[iseg] * conv * csphi, rMax[iseg] * conv * ssphi, z1},
0105                     {rMax[iseg] * conv * cephi, rMax[iseg] * conv * sephi, z1},
0106                     {rMax[iseg + 1] * conv * cephi, rMax[iseg + 1] * conv * sephi, z2},
0107                     {rMax[iseg + 1] * conv * csphi, rMax[iseg + 1] * conv * ssphi, z2}};
0108         isurf    = builder::CreateLocalSurfaceFromVertices<Real_t>(vertices, logical_id);
0109         if (intersection) builder::GetSurface<Real_t>(isurf).fSkipConvexity = true;
0110         if (realSeg) {
0111           if (nseg > 1 || iside > 0) logic.push_back(land);
0112           logic.push_back(isurf);
0113           // Only a concave outer surface is embedding, since this is practically a boolean union.
0114           auto is_concave = IsConvexConcave<vecgeom::Precision>(iseg, zPlanes, rMax, /*convex_check=*/0, nseg + 1);
0115           if (!(is_concave)) builder::GetSurface<Real_t>(isurf).fEmbedding = false;
0116         } else {
0117           builder::GetSurface<Real_t>(isurf).fEmbedding = false;
0118         }
0119       }
0120     }
0121
0122     // inner surfaces
0123     if (hasInner) {
0124       for (size_t iside = 0; iside < sideCount; iside++) {
0125         getSinCos(getPhi(iside), getPhi(iside + 1));
0126         vertices = {{rMin[iseg + 1] * conv * csphi, rMin[iseg + 1] * conv * ssphi, z2},
0127                     {rMin[iseg + 1] * conv * cephi, rMin[iseg + 1] * conv * sephi, z2},
0128                     {rMin[iseg] * conv * cephi, rMin[iseg] * conv * sephi, z1},
0129                     {rMin[iseg] * conv * csphi, rMin[iseg] * conv * ssphi, z1}};
0130         isurf    = builder::CreateLocalSurfaceFromVertices<Real_t>(vertices, logical_id);
0131         if (intersection) builder::GetSurface<Real_t>(isurf).fSkipConvexity = true;
0132         if (realSeg) {
0133           if (iside == 0) {
0134             logic.push_back(land);
0135             logic.push_back(lplus);
0136           } else {
0137             logic.push_back(lor);
0138           }
0139           logic.push_back(isurf);
0140           if (iside == sideCount - 1) {
0141             logic.push_back(lminus);
0142           }
0143           // Only a convex inner surface is embedding, since this is practically a boolean union.
0144           auto is_convex = IsConvexConcave<vecgeom::Precision>(iseg, zPlanes, rMin, /*convex_check=*/1, nseg + 1);
0145           if (!(is_convex)) builder::GetSurface<Real_t>(isurf).fEmbedding = false;
0146         } else {
0147           builder::GetSurface<Real_t>(isurf).fEmbedding = false;
0148         }
0149       }
0150     }
0151
0152     if (hasPhi) {
0153       getSinCos(phiStart, phiStart + phiDelta);
0154       vertices = {{rMax[iseg + 1] * conv * csphi, rMax[iseg + 1] * conv * ssphi, z2},
0155                   {rMin[iseg + 1] * conv * csphi, rMin[iseg + 1] * conv * ssphi, z2},
0156                   {rMin[iseg] * conv * csphi, rMin[iseg] * conv * ssphi, z1},
0157                   {rMax[iseg] * conv * csphi, rMax[iseg] * conv * ssphi, z1}};
0158       isurf    = builder::CreateLocalSurfaceFromVertices<Real_t>(vertices, logical_id);
0159       if (intersection) builder::GetSurface<Real_t>(isurf).fSkipConvexity = true;
0160       builder::GetSurface<Real_t>(isurf).fEmbedding = false;
0161       logic.push_back(land);
0162       logic.push_back(lplus); // '('
0163       logic.push_back(isurf);
0164
0165       vertices = {{rMax[iseg] * conv * cephi, rMax[iseg] * conv * sephi, z1},
0166                   {rMin[iseg] * conv * cephi, rMin[iseg] * conv * sephi, z1},
0167                   {rMin[iseg + 1] * conv * cephi, rMin[iseg + 1] * conv * sephi, z2},
0168                   {rMax[iseg + 1] * conv * cephi, rMax[iseg + 1] * conv * sephi, z2}};
0169       isurf    = builder::CreateLocalSurfaceFromVertices<Real_t>(vertices, logical_id);
0170       if (intersection) builder::GetSurface<Real_t>(isurf).fSkipConvexity = true;
0171       builder::GetSurface<Real_t>(isurf).fEmbedding = false;
0172       logic.push_back(smallerPi ? land : lor);
0173       logic.push_back(isurf);
0174       logic.push_back(lminus); // ')'
0175     }
0176
0177     for (size_t iside = 0; iside < sideCount; iside++) {
0178       getSinCos(getPhi(iside), getPhi(iside + 1));
0179       // Add bottom frames
0180       if (iseg == 0) {
0181         vertices = {{rMin[iseg] * conv * cephi, rMin[iseg] * conv * sephi, z1},
0182                     {rMax[iseg] * conv * cephi, rMax[iseg] * conv * sephi, z1},
0183                     {rMax[iseg] * conv * csphi, rMax[iseg] * conv * ssphi, z1},
0184                     {rMin[iseg] * conv * csphi, rMin[iseg] * conv * ssphi, z1}};
0185         isurf    = builder::CreateLocalSurfaceFromVertices<Real_t>(vertices, logical_id);
0186         if (intersection) builder::GetSurface<Real_t>(isurf).fSkipConvexity = true;
0187         // the bottom surface may be degenerated
0188         if (isurf >= 0) {
0189           builder::GetSurface<Real_t>(isurf).fEmbedding = false;
0190           logic.push_back(land);
0191           logic.push_back(isurf);
0192         }
0193       }
0194
0195       // Add top frames
0196       if (iseg == (nseg - 1)) {
0197         vertices = {{rMin[iseg + 1] * conv * csphi, rMin[iseg + 1] * conv * ssphi, z2},
0198                     {rMax[iseg + 1] * conv * csphi, rMax[iseg + 1] * conv * ssphi, z2},
0199                     {rMax[iseg + 1] * conv * cephi, rMax[iseg + 1] * conv * sephi, z2},
0200                     {rMin[iseg + 1] * conv * cephi, rMin[iseg + 1] * conv * sephi, z2}};
0201         isurf    = builder::CreateLocalSurfaceFromVertices<Real_t>(vertices, logical_id);
0202         if (intersection) builder::GetSurface<Real_t>(isurf).fSkipConvexity = true;
0203         // the top surface may be degenerated
0204         if (isurf >= 0) {
0205           builder::GetSurface<Real_t>(isurf).fEmbedding = false;
0206           logic.push_back(land);
0207           logic.push_back(isurf);
0208         }
0209       }
0210     }
0211
0212     if (realSeg) logic.push_back(lminus); // ')' end section logic
0213   };
0214
0215   for (size_t i = 0; i < nseg; ++i)
0216     convertZsegment(i);
0217
0218   builder::AddLogicToShell<Real_t>(logical_id, logic);
0219   return true;
0220 }
0221
0222 } // namespace conv
0223 } // namespace vgbrep
0224 #endif