surfaces/conv/PolyconeConverter.h

0001 #ifndef VECGEOM_SURFACE_POLYCONECONVERTER_H_
0002 #define VECGEOM_SURFACE_POLYCONECONVERTER_H_
0003
0004 #include <VecGeom/surfaces/conv/ConvHelper.h>
0005
0006 #include <VecGeom/volumes/Polycone.h>
0007
0008 namespace vgbrep {
0009 namespace conv {
0010
0011 /// @brief Converter for Polycone
0012 /// @tparam Real_t Precision type
0013 /// @param polycone Polycone solid to be converted
0014 /// @param logical_id Id of the logical volume
0015 /// @return Conversion success
0016 template <typename Real_t>
0017 bool CreatePolyconeSurfaces(vecgeom::UnplacedPolycone const &polycone, int logical_id, bool intersection = false)
0018 {
0019   using RingMask_t = RingMask<Real_t>;
0020   using ZPhiMask_t = ZPhiMask<Real_t>;
0021   using Vector3D   = vecgeom::Vector3D<vecgeom::Precision>;
0022
0023   LogicExpressionCPU
0024       logic; // top & bottom & [rmin] & rmax & (dphi < 180) ? sphi * ephi : sphi | ephi  auto rmin1 = cone.GetRmin1();
0025
0026   auto nSect = polycone.GetNSections();
0027   auto sphi  = polycone.GetStartPhi();
0028   auto dphi  = polycone.GetDeltaPhi();
0029   auto ephi  = polycone.GetEndPhi();
0030   VECGEOM_ASSERT(dphi > vecgeom::kTolerance);
0031
0032   bool fullCirc  = ApproxEqual(dphi, vecgeom::kTwoPi);
0033   bool smallerPi = dphi < (vecgeom::kPi - vecgeom::kTolerance);
0034
0035   int isurf;
0036   vecgeom::Precision surfdata[2];
0037
0038   // We need angles in degrees for transformations
0039   auto sphid = vecgeom::kRadToDeg * sphi;
0040   auto ephid = vecgeom::kRadToDeg * ephi;
0041
0042   vecgeom::Transformation3D transformation;
0043   std::vector<Vector3D> vert; // Stores 3D coordinates of the four corners that create a quadrilateral.
0044   auto csphi = std::cos(sphi);
0045   auto ssphi = std::sin(sphi);
0046   auto cephi = std::cos(ephi);
0047   auto sephi = std::sin(ephi);
0048
0049   // generate z, rmin and rmax arrays, 2*nSect because each section has a bottom and a top value
0050   std::unique_ptr<vecgeom::Precision[]> zArray(new vecgeom::Precision[2 * nSect]);
0051   std::unique_ptr<vecgeom::Precision[]> rMinArray(new vecgeom::Precision[2 * nSect]);
0052   std::unique_ptr<vecgeom::Precision[]> rMaxArray(new vecgeom::Precision[2 * nSect]);
0053   // fill rArrays per section
0054   for (int i = 0; i < nSect; ++i) {
0055     rMinArray[2 * i]     = polycone.GetRmin1AtSection(i);
0056     rMinArray[2 * i + 1] = polycone.GetRmin2AtSection(i);
0057     rMaxArray[2 * i]     = polycone.GetRmax1AtSection(i);
0058     rMaxArray[2 * i + 1] = polycone.GetRmax2AtSection(i);
0059   }
0060   // fill zArray per section
0061   zArray[0] = polycone.GetZAtPlane(0);
0062   for (int i = 1; i < nSect; ++i) {
0063     zArray[2 * i - 1] = polycone.GetZAtPlane(i);
0064     zArray[2 * i]     = polycone.GetZAtPlane(i);
0065   }
0066   zArray[2 * nSect - 1] = polycone.GetZAtPlane(nSect);
0067
0068   // loop over sections, at each section the full cone-like solid is constructed
0069   for (int i = 0; i < nSect; i++) {
0070     logic.push_back(lplus); // '(' parenthesis around section
0071
0072     auto z1    = polycone.GetZAtPlane(i);     // note: indices for GetZAtPlane is the
0073     auto z2    = polycone.GetZAtPlane(i + 1); // number of distinct z points of the polycone
0074     auto rmin1 = polycone.GetRmin1AtSection(i);
0075     auto rmax1 = polycone.GetRmax1AtSection(i);
0076     auto rmin2 = polycone.GetRmin2AtSection(i);
0077     auto rmax2 = polycone.GetRmax2AtSection(i);
0078
0079     VECGEOM_ASSERT(rmax1 - rmin1 > -vecgeom::kTolerance);
0080     VECGEOM_ASSERT(rmax2 - rmin2 > -vecgeom::kTolerance);
0081
0082     // if rmax == rmin then for safety it is set to rmax = rmin + ConeTolerance already in the solid model
0083     // For using mixed precision, the compiled tolerance must be replaced by the mixed tolerance
0084     if (rmax1 - rmin1 == vecgeom::kConeTolerance) rmax1 = rmin1 + vecgeom::kToleranceCone<Real_t>;
0085     if (rmax2 - rmin2 == vecgeom::kConeTolerance) rmax2 = rmin2 + vecgeom::kToleranceCone<Real_t>;
0086
0087     // construct top and bottom surfaces:
0088     // if there are previous or following sections, the top and bottom surface consist
0089     // of up to two rings (inner and outer ring). All of those are real surfaces
0090
0091     // lets start with the adjusted bottom rings:
0092     // case where there is a previous section
0093     if (i > 0) {
0094       auto rmin_prev = polycone.GetRmin2AtSection(i - 1);
0095       auto rmax_prev = polycone.GetRmax2AtSection(i - 1);
0096
0097       bool has_inner_ring = rmin1 < rmin_prev - vecgeom::kTolerance;
0098       bool has_outer_ring = rmax1 > rmax_prev + vecgeom::kTolerance;
0099
0100       if (has_inner_ring) {
0101         isurf = builder::CreateLocalSurface<Real_t>(
0102             builder::CreateUnplacedSurface<Real_t>(SurfaceType::kPlanar),
0103             builder::CreateFrame<Real_t>(FrameType::kRing, RingMask_t{rmin1, rmin_prev, fullCirc, sphi, ephi}),
0104             vecgeom::Transformation3DMP<Precision>(0., 0., z1, 0., 180., -sphid - ephid));
0105         builder::GetSurface<Real_t>(isurf).fEmbedding = false;
0106         builder::AddSurfaceToShell<Real_t>(logical_id, isurf);
0107         if (intersection) builder::GetSurface<Real_t>(isurf).fSkipConvexity = true;
0108         logic.push_back(isurf);
0109         logic.push_back(land);
0110       }
0111       if (has_outer_ring) {
0112         isurf = builder::CreateLocalSurface<Real_t>(
0113             builder::CreateUnplacedSurface<Real_t>(SurfaceType::kPlanar),
0114             builder::CreateFrame<Real_t>(FrameType::kRing, RingMask_t{rmax_prev, rmax1, fullCirc, sphi, ephi}),
0115             vecgeom::Transformation3DMP<Precision>(0., 0., z1, 0., 180., -sphid - ephid));
0116         builder::GetSurface<Real_t>(isurf).fEmbedding = false;
0117         if (intersection) builder::GetSurface<Real_t>(isurf).fSkipConvexity = true;
0118         builder::AddSurfaceToShell<Real_t>(logical_id, isurf);
0119         logic.push_back(isurf);
0120         logic.push_back(land);
0121       }
0122
0123       // add virtual surface at z1 if there is no real one
0124       if (!has_inner_ring && !has_outer_ring) {
0125         isurf = builder::CreateLocalSurface<Real_t>(builder::CreateUnplacedSurface<Real_t>(SurfaceType::kPlanar),
0126                                                     Frame{FrameType::kNoFrame},
0127                                                     vecgeom::Transformation3DMP<Precision>(0., 0., z1, 0., 180., 0.));
0128         builder::AddSurfaceToShell<Real_t>(logical_id, isurf);
0129         if (intersection) builder::GetSurface<Real_t>(isurf).fSkipConvexity = true;
0130         logic.push_back(isurf);
0131         logic.push_back(land);
0132       }
0133     }
0134     // adjusted top rings in case there is a following section
0135     if (i < nSect - 1) {
0136       auto rmin_next = polycone.GetRmin1AtSection(i + 1);
0137       auto rmax_next = polycone.GetRmax1AtSection(i + 1);
0138
0139       bool has_inner_ring = rmin2 < rmin_next - vecgeom::kTolerance;
0140       bool has_outer_ring = rmax2 > rmax_next + vecgeom::kTolerance;
0141
0142       if (has_inner_ring) {
0143         isurf = builder::CreateLocalSurface<Real_t>(
0144             builder::CreateUnplacedSurface<Real_t>(SurfaceType::kPlanar),
0145             builder::CreateFrame<Real_t>(FrameType::kRing, RingMask_t{rmin2, rmin_next, fullCirc, sphi, ephi}),
0146             vecgeom::Transformation3DMP<Precision>(0., 0., z2, 0., 0., 0.));
0147         builder::AddSurfaceToShell<Real_t>(logical_id, isurf);
0148         if (intersection) builder::GetSurface<Real_t>(isurf).fSkipConvexity = true;
0149         builder::GetSurface<Real_t>(isurf).fEmbedding = false;
0150         logic.push_back(isurf);
0151         logic.push_back(land);
0152       }
0153       if (has_outer_ring) {
0154         isurf = builder::CreateLocalSurface<Real_t>(
0155             builder::CreateUnplacedSurface<Real_t>(SurfaceType::kPlanar),
0156             builder::CreateFrame<Real_t>(FrameType::kRing, RingMask_t{rmax_next, rmax2, fullCirc, sphi, ephi}),
0157             vecgeom::Transformation3DMP<Precision>(0., 0., z2, 0., 0., 0.));
0158         builder::AddSurfaceToShell<Real_t>(logical_id, isurf);
0159         if (intersection) builder::GetSurface<Real_t>(isurf).fSkipConvexity = true;
0160         builder::GetSurface<Real_t>(isurf).fEmbedding = false;
0161         logic.push_back(isurf);
0162         logic.push_back(land);
0163       }
0164
0165       // add virtual surface at z2 if there is no real one
0166       if (!has_inner_ring && !has_outer_ring) {
0167         isurf = builder::CreateLocalSurface<Real_t>(builder::CreateUnplacedSurface<Real_t>(SurfaceType::kPlanar),
0168                                                     Frame{FrameType::kNoFrame},
0169                                                     vecgeom::Transformation3DMP<Precision>(0., 0., z2, 0., 0., 0.));
0170         builder::AddSurfaceToShell<Real_t>(logical_id, isurf);
0171         if (intersection) builder::GetSurface<Real_t>(isurf).fSkipConvexity = true;
0172         logic.push_back(isurf);
0173         logic.push_back(land);
0174       }
0175     }
0176
0177     // first section has the full initial ring as bottom surface
0178     if (i == 0) {
0179       // real surface at z1
0180       if (rmax1 - rmin1 > vecgeom::kTolerance) {
0181         isurf = builder::CreateLocalSurface<Real_t>(
0182             builder::CreateUnplacedSurface<Real_t>(SurfaceType::kPlanar),
0183             builder::CreateFrame<Real_t>(FrameType::kRing, RingMask_t{rmin1, rmax1, fullCirc, sphi, ephi}),
0184             vecgeom::Transformation3DMP<Precision>(0., 0., z1, 0., 180., -sphid - ephid));
0185         builder::AddSurfaceToShell<Real_t>(logical_id, isurf);
0186         builder::GetSurface<Real_t>(isurf).fEmbedding = false;
0187         if (intersection) builder::GetSurface<Real_t>(isurf).fSkipConvexity = true;
0188         logic.push_back(isurf);
0189         logic.push_back(land);
0190       }
0191     }
0192
0193     // last section has the full initial ring as top surface
0194     if (i == nSect - 1) {
0195       // real surface at z2
0196       if (rmax1 - rmin1 > vecgeom::kTolerance) {
0197         isurf = builder::CreateLocalSurface<Real_t>(
0198             builder::CreateUnplacedSurface<Real_t>(SurfaceType::kPlanar),
0199             builder::CreateFrame<Real_t>(FrameType::kRing, RingMask_t{rmin2, rmax2, fullCirc, sphi, ephi}),
0200             vecgeom::Transformation3DMP<Precision>(0., 0., z2, 0., 0., 0.));
0201         builder::AddSurfaceToShell<Real_t>(logical_id, isurf);
0202         builder::GetSurface<Real_t>(isurf).fEmbedding = false;
0203         if (intersection) builder::GetSurface<Real_t>(isurf).fSkipConvexity = true;
0204         logic.push_back(isurf);
0205         logic.push_back(land);
0206       }
0207     }
0208
0209     // Cones:
0210     const Real_t z_shift = 0.5 * (z1 + z2);
0211     // inner cone
0212     if (rmin2 > vecgeom::kTolerance || rmin1 > vecgeom::kTolerance) {
0213       surfdata[0] = 0.5 * (rmin1 + rmin2);
0214       surfdata[1] = (rmin2 - rmin1) / (z2 - z1);
0215       auto stype  = std::abs(surfdata[1]) > vecgeom::kTolerance ? SurfaceType::kConical : SurfaceType::kCylindrical;
0216       isurf       = builder::CreateLocalSurface<Real_t>(
0217           builder::CreateUnplacedSurface<Real_t>(stype, surfdata, /*flipped=*/true),
0218           builder::CreateFrame<Real_t>(FrameType::kZPhi,
0219                                        ZPhiMask_t{z1 - z_shift, z2 - z_shift, fullCirc, rmin1, rmin2, sphi, ephi}),
0220           vecgeom::Transformation3DMP<Precision>(0., 0., z_shift, 0., 0., 0.));
0221       builder::AddSurfaceToShell<Real_t>(logical_id, isurf);
0222       // Only a convex inner cone is embedding, since this is practically a boolean union.
0223       auto is_convex = IsConvexConcave<vecgeom::Precision>(2 * i, zArray, rMinArray, /*convex_check=*/1, 2 * nSect);
0224       if (!(is_convex)) builder::GetSurface<Real_t>(isurf).fEmbedding = false;
0225       if (intersection)
0226         builder::GetSurface<Real_t>(isurf).fSkipConvexity =
0227             true; // this is for the convexity check of making booleans treated as normal surfaces and is not related to
0228                   // the convexity check above
0229       logic.push_back(isurf);
0230       logic.push_back(land);
0231     }
0232     // outer cone
0233     surfdata[0] = 0.5 * (rmax1 + rmax2);
0234     surfdata[1] = (rmax2 - rmax1) / (z2 - z1);
0235     auto stype  = std::abs(surfdata[1]) > vecgeom::kTolerance ? SurfaceType::kConical : SurfaceType::kCylindrical;
0236     isurf       = builder::CreateLocalSurface<Real_t>(
0237         builder::CreateUnplacedSurface<Real_t>(stype, surfdata),
0238         builder::CreateFrame<Real_t>(FrameType::kZPhi,
0239                                      ZPhiMask_t{z1 - z_shift, z2 - z_shift, fullCirc, rmax1, rmax2, sphi, ephi}),
0240         vecgeom::Transformation3DMP<Precision>(0., 0., z_shift, 0., 0., 0.));
0241     builder::AddSurfaceToShell<Real_t>(logical_id, isurf);
0242     // Only a concave outer cone is embedding, since this is practically a boolean union.
0243     auto is_concave = IsConvexConcave<vecgeom::Precision>(2 * i, zArray, rMaxArray, /*convex_check=*/0, 2 * nSect);
0244     if (!(is_concave)) builder::GetSurface<Real_t>(isurf).fEmbedding = false;
0245     if (intersection) builder::GetSurface<Real_t>(isurf).fSkipConvexity = true;
0246     logic.push_back(isurf);
0247
0248     if (!fullCirc) {
0249       // corners of the sphi and ephi faces
0250       std::vector<Vector3D> corners = {{rmin1 * csphi, rmin1 * ssphi, z1}, {rmax1 * csphi, rmax1 * ssphi, z1},
0251                                        {rmax2 * csphi, rmax2 * ssphi, z2}, {rmin2 * csphi, rmin2 * ssphi, z2},
0252                                        {rmax1 * cephi, rmax1 * sephi, z1}, {rmin1 * cephi, rmin1 * sephi, z1},
0253                                        {rmin2 * cephi, rmin2 * sephi, z2}, {rmax2 * cephi, rmax2 * sephi, z2}};
0254
0255       // plane cap at sphi
0256       vert  = {corners[0], corners[1], corners[2], corners[3]};
0257       isurf = builder::CreateLocalSurfaceFromVertices<Real_t>(vert, logical_id);
0258       if (intersection) builder::GetSurface<Real_t>(isurf).fSkipConvexity = true;
0259       VECGEOM_ASSERT(isurf >= 0);
0260       builder::GetSurface<Real_t>(isurf).fEmbedding = false;
0261       logic.push_back(land);
0262       logic.push_back(lplus); // '('
0263       logic.push_back(isurf);
0264
0265       // plane cap at sphi+dphi
0266       vert  = {corners[4], corners[5], corners[6], corners[7]};
0267       isurf = builder::CreateLocalSurfaceFromVertices<Real_t>(vert, logical_id);
0268       if (intersection) builder::GetSurface<Real_t>(isurf).fSkipConvexity = true;
0269       VECGEOM_ASSERT(isurf >= 0);
0270       builder::GetSurface<Real_t>(isurf).fEmbedding = false;
0271       logic.push_back(smallerPi ? land : lor);
0272       logic.push_back(isurf);
0273       logic.push_back(lminus); // ')'
0274     }
0275
0276     logic.push_back(lminus);                 // ')' // close parenthesis after section
0277     if (i < nSect - 1) logic.push_back(lor); // union with next section if any
0278   }                                          // end loop over sections
0279   builder::AddLogicToShell<Real_t>(logical_id, logic);
0280
0281   return true;
0282 }
0283
0284 } // namespace conv
0285 } // namespace vgbrep
0286 #endif