surfaces/conv/Builder.h

0001 #ifndef VECGEOM_SURFACE_BUILDER_H
0002 #define VECGEOM_SURFACE_BUILDER_H
0003
0004 #include <VecGeom/surfaces/base/CpuTypes.h>
0005
0006 #include <initializer_list>
0007
0008 namespace vgbrep {
0009
0010 namespace builder {
0011
0012 template <typename Real_t>
0013 UnplacedSurface<Real_t> CreateUnplacedSurface(SurfaceType type, vecgeom::Precision *data = nullptr, bool flip = false)
0014 {
0015   auto &cpudata = CPUsurfData<Real_t>::Instance();
0016   switch (type) {
0017   case SurfaceType::kPlanar:
0018     return UnplacedSurface<Real_t>(type);
0019   case SurfaceType::kCylindrical:
0020   case SurfaceType::kSpherical:
0021     return UnplacedSurface<Real_t>(type, -1, data[0], 0., flip);
0022   case SurfaceType::kElliptical:
0023     cpudata.fEllipData.push_back({data[0], data[1], data[2]});
0024     return UnplacedSurface<Real_t>(type, cpudata.fEllipData.size() - 1);
0025   case SurfaceType::kConical:
0026     return UnplacedSurface<Real_t>(type, -1, data[0], data[1], flip);
0027   case SurfaceType::kTorus:
0028     cpudata.fTorusData.push_back({data[0], data[1], data[2], data[3], flip});
0029     return UnplacedSurface<Real_t>(type, cpudata.fTorusData.size() - 1);
0030   case SurfaceType::kArb4:
0031     cpudata.fArb4Data.push_back(
0032         {data[0], data[1], data[2], data[3], data[4], data[5], data[6], data[7], data[8], data[9]});
0033     return UnplacedSurface<Real_t>(type, cpudata.fArb4Data.size() - 1);
0034   default:
0035     return UnplacedSurface<Real_t>(type);
0036   };
0037 }
0038
0039 // Creators for different types of frames.
0040 template <typename Real_t>
0041 Frame CreateFrame(FrameType type, WindowMask<Real_t> const &mask)
0042 {
0043   auto &cpudata = CPUsurfData<Real_t>::Instance();
0044   int id        = cpudata.fWindowMasks.size();
0045   cpudata.fWindowMasks.push_back(mask);
0046   return Frame(type, id);
0047 }
0048
0049 template <typename Real_t>
0050 Frame CreateFrame(FrameType type, RingMask<Real_t> const &mask)
0051 {
0052   auto &cpudata = CPUsurfData<Real_t>::Instance();
0053   int id        = cpudata.fRingMasks.size();
0054   cpudata.fRingMasks.push_back(mask);
0055   return Frame(type, id);
0056 }
0057
0058 template <typename Real_t>
0059 Frame CreateFrame(FrameType type, ZPhiMask<Real_t> const &mask)
0060 {
0061   auto &cpudata = CPUsurfData<Real_t>::Instance();
0062   int id        = cpudata.fZPhiMasks.size();
0063   cpudata.fZPhiMasks.push_back(mask);
0064   return Frame(type, id);
0065 }
0066
0067 template <typename Real_t>
0068 Frame CreateFrame(FrameType type, QuadrilateralMask<Real_t> const &mask)
0069 {
0070   auto &cpudata = CPUsurfData<Real_t>::Instance();
0071   int id        = cpudata.fQuadMasks.size();
0072   cpudata.fQuadMasks.push_back(mask);
0073   return Frame(type, id);
0074 }
0075
0076 template <typename Real_t>
0077 Frame CreateFrame(FrameType type, TriangleMask<Real_t> const &mask)
0078 {
0079   auto &cpudata = CPUsurfData<Real_t>::Instance();
0080   int id        = cpudata.fTriangleMasks.size();
0081   cpudata.fTriangleMasks.push_back(mask);
0082   return Frame(type, id);
0083 }
0084
0085 template <typename Real_t>
0086 int CreateLocalSurface(UnplacedSurface<Real_t> const &unplaced, Frame const &frame, TransformationMP<Real_t> trans,
0087                        bool never_check = false)
0088 {
0089   auto &cpudata = CPUsurfData<Real_t>::Instance();
0090   int id        = cpudata.fLocalSurfaces.size();
0091   cpudata.fLocalSurfaces.push_back({unplaced, frame, trans, /*NavIndex=*/0, never_check});
0092   return id;
0093 }
0094
0095 template <typename Real_t>
0096 FramedSurface<Real_t, TransformationMP<Real_t>> &GetSurface(int isurf)
0097 {
0098   auto &cpudata = CPUsurfData<Real_t>::Instance();
0099   VECGEOM_ASSERT(size_t(isurf) < cpudata.fLocalSurfaces.size());
0100   return cpudata.fLocalSurfaces[isurf];
0101 }
0102
0103 template <typename Real_t>
0104 int AddSurfaceToShell(int logical_id, int isurf)
0105 {
0106   auto &cpudata = CPUsurfData<Real_t>::Instance();
0107   if (cpudata.fShells.size() == 0) {
0108     std::cout << "BrepHelper::AddSurfaceToShell: need to call SetNvolumes first\n";
0109     return -1;
0110   }
0111   VECGEOM_VALIDATE(logical_id < (int)cpudata.fShells.size(), << "surface shell id exceeding number of volumes");
0112   int id = cpudata.fShells[logical_id].fSurfaces.size();
0113   cpudata.fShells[logical_id].fSurfaces.push_back(isurf);
0114   cpudata.fSceneShells[logical_id].fSurfaces.push_back(isurf);
0115   cpudata.fLocalSurfaces[isurf].fSurfIndex = id;
0116   return id;
0117 }
0118
0119 template <typename Real_t>
0120 void AddLogicToShell(int logical_id, LogicExpressionCPU &logic)
0121 {
0122   // Add solid logic to existing shell logic
0123   auto &cpudata  = CPUsurfData<Real_t>::Instance();
0124   auto &crtlogic = cpudata.fShells[logical_id].fLogic;
0125   crtlogic.insert(crtlogic.end(), logic.begin(), logic.end());
0126 }
0127
0128 /// @brief Creates a quadrilateral frame, triangular frame of no frame based on a vector of vertices (using only XY coordinates)
0129 /// @tparam Real_t Precision type
0130 /// @tparam Container Container type
0131 /// @param points Vector of points
0132 /// @return Created frame. If this has the type kNoFrame, the user must abort framed surface creation
0133 template <typename Real_t, typename Container>
0134 Frame CreateFrameFromVertices(Container &points, TransformationMP<Real_t> &trans)
0135 {
0136   if (points.size() == 3)
0137     return CreateFrame<Real_t>(FrameType::kTriangle, TriangleMask<Real_t>{points[0].x(), points[0].y(), points[1].x(),
0138                                                                           points[1].y(), points[2].x(), points[2].y()});
0139   if (points.size() == 4) {
0140     // Check for rectangular frame
0141     // The 0->1 vector is aligned with the local Ox
0142     bool rectangle = ApproxEqualVector(points[1] - points[0], points[2] - points[3]);
0143     rectangle &=
0144         ApproxEqual(static_cast<Real_t>((points[1] - points[0]).Dot(points[3] - points[0])), static_cast<Real_t>(0.));
0145     if (rectangle) {
0146       auto dx = 0.5 * (points[1] - points[0]).Mag();
0147       auto dy = 0.5 * (points[3] - points[0]).Mag();
0148       return CreateFrame<Real_t>(FrameType::kWindow, WindowMask<Real_t>{dx, dy});
0149     } else {
0150       return CreateFrame<Real_t>(FrameType::kQuadrilateral,
0151                                  QuadrilateralMask<Real_t>{points[0].x(), points[0].y(), points[1].x(), points[1].y(),
0152                                                            points[2].x(), points[2].y(), points[3].x(), points[3].y()});
0153     }
0154   }
0155   return Frame(FrameType::kNoFrame);
0156 }
0157
0158 /// @brief Compute transformation (rotation + translation) for a surface defined by a set of co-planar points.
0159 /// @details The points must be ordered such that the cross product of any two consecutive segments has the same
0160 /// direction as the normal. All points must be different and not all colinear.
0161 /// @param points Input container holding the ordered Vector3D<Real_t> points, defined in the solid frame.
0162 //  The container will return the points in the local surface reference frame
0163 /// @return Transformation moving a surface from the (XOY) plane to the final position. The container will hold the
0164 /// input points transformed with the inverse transformation, lying in the (XOY) plane
0165 template <typename Real_t, typename Container>
0166 vecgeom::Transformation3DMP<Real_t> TransformationFromPlanarPoints(Container &points)
0167 {
0168   using Vector3 = vecgeom::Vector3D<Real_t>;
0169
0170   int npoints = points.size();
0171   VECGEOM_VALIDATE(npoints > 2, << "TransformationFromPlanarPoints takes at least three points");
0172   int istart            = 0;
0173   Real_t cross_mag2_max = 0.;
0174   Vector3 normal;
0175   Vector3 center;
0176   for (int i = 0; i < npoints; ++i) {
0177     center += points[i];
0178     auto a          = points[(i + 1) % npoints] - points[i];
0179     auto b          = points[(i + 2) % npoints] - points[(i + 1) % npoints];
0180     auto a_cross_b  = a.Cross(b);
0181     auto cross_mag2 = a_cross_b.Mag2();
0182     if (cross_mag2 > cross_mag2_max + cross_mag2_max * vecgeom::kToleranceDistSquared<Real_t>) {
0183       normal         = a_cross_b.Unit();
0184       istart         = i;
0185       cross_mag2_max = cross_mag2;
0186     }
0187   }
0188   VECGEOM_ASSERT(cross_mag2_max > vecgeom::kToleranceDistSquared<Real_t> &&
0189                  "TransformationFromPlanarPoints: degenerated polygon");
0190   center *= 1. / npoints;
0191
0192   Vector3 zref = normal;
0193   Vector3 xref = (points[(istart + 1) % npoints] - points[istart]).Unit();
0194   Vector3 yref = zref.Cross(xref);
0195   vecgeom::Transformation3DMP<Real_t> transformation(center[0], center[1], center[2], xref[0], yref[0], zref[0],
0196                                                      xref[1], yref[1], zref[1], xref[2], yref[2], zref[2]);
0197   // Convert points to the local frame
0198   for (int i = 0; i < npoints; ++i) {
0199     Vector3 local = transformation.Transform(points[i]);
0200     points[i]     = local;
0201   }
0202   return transformation;
0203 }
0204
0205 /// @brief Create local surface starting from a vector of maximum four vertices
0206 /// @tparam Real_t Precision type
0207 /// @tparam Container Container type
0208 /// @param points Vertices vector
0209 /// @return Index of created local surface
0210 template <typename Real_t, typename Container>
0211 int CreateLocalSurfaceFromVertices(Container &points, int logical_id)
0212 {
0213   // helper function to find non-coplanar surfaces (Arb4)
0214   using Vector3      = vecgeom::Vector3D<Real_t>;
0215   auto PointsOnPlane = [](Container &points) {
0216     // this method of checking the co-planarity was taken from the UnplacedTrapezoid
0217     Vector3 normalVector = (points[2] - points[0]).Cross(points[3] - points[1]);
0218     normalVector.Normalize();
0219     Vector3 centr = 0.25 * (points[0] + points[1] + points[2] + points[3]);
0220
0221     Real_t dist_scale = Max(Max((points[0] - centr).Length(), (points[1] - centr).Length()),
0222                             Max((points[2] - centr).Length(), (points[3] - centr).Length()));
0223     // check co-planarity
0224     Real_t resid1 = normalVector.Dot(points[0] - centr);
0225     Real_t resid2 = normalVector.Dot(points[1] - centr);
0226     Real_t resid3 = normalVector.Dot(points[2] - centr);
0227     Real_t resid4 = normalVector.Dot(points[3] - centr);
0228     Real_t resid  = Max(Max(fabs(resid1), fabs(resid2)), Max(fabs(resid3), fabs(resid4)));
0229
0230     // the residue should be small compared to the length scale of the quadrilateral
0231     return resid / dist_scale < 100 * vecgeom::kTolerance;
0232   };
0233
0234   // copy container because the content may get changed due to degenerated vertices
0235   Container vertices(points);
0236   // Remove duplicated vertices
0237   size_t i = 0;
0238   for (i = 0; i < vertices.size(); ++i) {
0239     for (size_t j = 0; j < i; ++j) {
0240       if (ApproxEqualVector(vertices[i], vertices[j])) {
0241         // remove duplicate vertex
0242         vertices.erase(vertices.begin() + i--);
0243         break;
0244       }
0245     }
0246   }
0247   if (vertices.size() < 3) return -1;
0248
0249   int isurf = 0;
0250   // Create transformation
0251   auto transformation = TransformationFromPlanarPoints<Real_t>(vertices);
0252   auto frame          = CreateFrameFromVertices<Real_t>(vertices, transformation);
0253   if (vertices.size() != 4 || PointsOnPlane(vertices)) {
0254     isurf =
0255         builder::CreateLocalSurface<Real_t>(CreateUnplacedSurface<Real_t>(SurfaceType::kPlanar), frame, transformation);
0256   } else { // creating Arb4 surface
0257     vecgeom::Precision surfdata[10];
0258     surfdata[0] = points[0].x();
0259     surfdata[1] = points[0].y();
0260     surfdata[2] = points[0].z();
0261     surfdata[3] = points[1].x();
0262     surfdata[4] = points[1].y();
0263     surfdata[5] = points[2].x();
0264     surfdata[6] = points[2].y();
0265     surfdata[7] = points[2].z();
0266     surfdata[8] = points[3].x();
0267     surfdata[9] = points[3].y();
0268     // we are using the frame above although it is never used for the Arb4
0269     vecgeom::Transformation3DMP<Real_t> identity;
0270     isurf = builder::CreateLocalSurface<Real_t>(CreateUnplacedSurface<Real_t>(SurfaceType::kArb4, surfdata), frame,
0271                                                 /*identity transformation*/ identity, /*never_check=*/1);
0272   }
0273   AddSurfaceToShell<Real_t>(logical_id, isurf);
0274   return isurf;
0275 }
0276
0277 } // namespace builder
0278
0279 } // namespace vgbrep
0280
0281 #endif