feat: Adding planar portal link (#4044)

Extending the `PortalLinkBase`, `GridPortalLink`, and `CompositePortalLink` classes to support `PlaneSurface`. 

<!-- This is an auto-generated comment: release notes by coderabbit.ai -->
## Summary by CodeRabbit

- **New Features**
	- Added support for `PlaneSurface` in grid portal link functionality.
	- Enhanced merging capabilities for plane surfaces.
	- Improved error handling and consistency checks for surface grid operations.

- **Bug Fixes**
	- Resolved issues with surface merging and grid creation for plane surfaces.
	- Added robust validation for surface bounds and axis directions.

- **Tests**
	- Introduced comprehensive test cases for plane surface grid merging.
	- Added validation for 1D and 2D grid portal link scenarios.
<!-- end of auto-generated comment: release notes by coderabbit.ai -->

Core/include/Acts/Geometry/GridPortalLink.hpp

@@ -11,13 +11,13 @@
 #include "Acts/Definitions/Tolerance.hpp"
 #include "Acts/Geometry/PortalLinkBase.hpp"
 #include "Acts/Geometry/TrackingVolume.hpp"
-#include "Acts/Surfaces/BoundaryTolerance.hpp"
 #include "Acts/Surfaces/CylinderSurface.hpp"
 #include "Acts/Surfaces/DiscSurface.hpp"
 #include "Acts/Surfaces/PlaneSurface.hpp"
 #include "Acts/Utilities/AxisDefinitions.hpp"
 #include "Acts/Utilities/Grid.hpp"
 #include "Acts/Utilities/Logger.hpp"
+#include "Acts/Utilities/ThrowAssert.hpp"
 
 #include <iosfwd>
 
@@ -67,6 +67,9 @@ class GridPortalLink : public PortalLinkBase {
     } else if (dynamic_cast<const DiscSurface*>(surface.get()) != nullptr &&
                direction != AxisR && direction != AxisPhi) {
       throw std::invalid_argument{"Invalid binning direction"};
+    } else if (dynamic_cast<const PlaneSurface*>(surface.get()) != nullptr &&
+               direction != AxisX && direction != AxisY) {
+      throw std::invalid_argument{"Invalid binning direction"};
     }
 
     return std::make_unique<GridPortalLinkT<axis_t>>(
@@ -90,6 +93,8 @@ class GridPortalLink : public PortalLinkBase {
       direction = AxisDirection::AxisRPhi;
     } else if (dynamic_cast<const DiscSurface*>(surface.get()) != nullptr) {
       direction = AxisDirection::AxisR;
+    } else if (dynamic_cast<const PlaneSurface*>(surface.get()) != nullptr) {
+      direction = AxisDirection::AxisX;
     }
 
     return std::make_unique<GridPortalLinkT<axis_1_t, axis_2_t>>(
@@ -353,6 +358,10 @@ class GridPortalLink : public PortalLinkBase {
   /// @param disc The disc surface
   void checkConsistency(const DiscSurface& disc) const;
 
+  /// Helper function to check consistency for grid on a plane surface
+  /// @param plane The plane surface
+  void checkConsistency(const PlaneSurface& plane) const;
+
   /// Expand a 1D grid to a 2D one for a cylinder surface
   /// @param surface The cylinder surface
   /// @param other The axis to use for the missing direction,
@@ -370,6 +379,14 @@ class GridPortalLink : public PortalLinkBase {
   std::unique_ptr<GridPortalLink> extendTo2dImpl(
       const std::shared_ptr<DiscSurface>& surface, const IAxis* other) const;
 
+  /// Expand a 1D grid to a 2D one for a plane surface
+  /// @param surface The plane surface
+  /// @param other The axis to use for the missing direction,
+  ///              can be null for auto determination
+  /// @return A unique pointer to the 2D grid portal link
+  std::unique_ptr<GridPortalLink> extendTo2dImpl(
+      const std::shared_ptr<PlaneSurface>& surface, const IAxis* other) const;
+
   /// Helper enum to declare which local direction to fill
   enum class FillDirection {
     loc0,
@@ -451,6 +468,17 @@ class GridPortalLinkT : public GridPortalLink {
       } else {
         throw std::invalid_argument{"Invalid binning direction"};
       }
+    } else if (const auto* plane =
+                   dynamic_cast<const PlaneSurface*>(m_surface.get())) {
+      checkConsistency(*plane);
+
+      if (direction == AxisX) {
+        m_projection = &projection<PlaneSurface, AxisX>;
+      } else if (direction == AxisDirection::AxisY) {
+        m_projection = &projection<PlaneSurface, AxisY>;
+      } else {
+        throw std::invalid_argument{"Invalid binning direction"};
+      }
 
     } else {
       throw std::logic_error{"Surface type is not supported"};
@@ -490,6 +518,9 @@ class GridPortalLinkT : public GridPortalLink {
       } else if (auto disc =
                      std::dynamic_pointer_cast<DiscSurface>(m_surface)) {
         return extendTo2dImpl(disc, other);
+      } else if (auto plane =
+                     std::dynamic_pointer_cast<PlaneSurface>(m_surface)) {
+        return extendTo2dImpl(plane, other);
       } else {
         throw std::logic_error{
             "Surface type is not supported (this should not happen)"};
@@ -539,7 +570,8 @@ class GridPortalLinkT : public GridPortalLink {
   Result<const TrackingVolume*> resolveVolume(
       const GeometryContext& /*gctx*/, const Vector2& position,
       double /*tolerance*/ = s_onSurfaceTolerance) const override {
-    assert(surface().insideBounds(position, BoundaryTolerance::None()));
+    throw_assert(surface().insideBounds(position, BoundaryTolerance::None()),
+                 "Checking volume outside of bounds");
     return m_grid.atPosition(m_projection(position));
   }
 

Core/src/Geometry/CompositePortalLink.cpp

@@ -15,10 +15,13 @@
 #include "Acts/Surfaces/DiscSurface.hpp"
 #include "Acts/Surfaces/PlaneSurface.hpp"
 #include "Acts/Surfaces/RadialBounds.hpp"
+#include "Acts/Surfaces/RectangleBounds.hpp"
 #include "Acts/Surfaces/RegularSurface.hpp"
 #include "Acts/Utilities/Axis.hpp"
+#include "Acts/Utilities/AxisDefinitions.hpp"
 
 #include <algorithm>
+#include <cstddef>
 #include <iostream>
 #include <iterator>
 #include <stdexcept>
@@ -48,7 +51,9 @@ std::shared_ptr<RegularSurface> mergedSurface(const Surface& a,
     return merged;
   } else if (const auto* planeA = dynamic_cast<const PlaneSurface*>(&a);
              planeA != nullptr) {
-    throw std::logic_error{"Plane surfaces not implemented yet"};
+    const auto& planeB = dynamic_cast<const PlaneSurface&>(b);
+    auto [merged, reversed] = planeA->mergedWith(planeB, direction);
+    return merged;
   } else {
     throw std::invalid_argument{"Unsupported surface type"};
   }
@@ -289,7 +294,55 @@ std::unique_ptr<GridPortalLink> CompositePortalLink::makeGrid(
 
     return grid;
   } else if (surface().type() == Surface::SurfaceType::Plane) {
-    throw std::runtime_error{"Plane surfaces not implemented yet"};
+    ACTS_VERBOSE("Combining composite into plane grid");
+
+    if (m_direction != AxisDirection::AxisX &&
+        m_direction != AxisDirection::AxisY) {
+      ACTS_ERROR("Plane grid only supports binning in x/y direction");
+      throw std::runtime_error{"Unsupported binning direction"};
+    }
+
+    bool dirX = m_direction == AxisDirection::AxisX;
+
+    std::vector<double> edges;
+    edges.reserve(m_children.size() + 1);
+
+    const Transform3& groupTransform = m_surface->transform(gctx);
+    Transform3 itransform = groupTransform.inverse();
+
+    std::size_t sortingDir = dirX ? eX : eY;
+    std::ranges::sort(trivialLinks, [&itransform, &gctx, sortingDir](
+                                        const auto& a, const auto& b) {
+      return (itransform * a->surface().transform(gctx))
+                 .translation()[sortingDir] <
+             (itransform * b->surface().transform(gctx))
+                 .translation()[sortingDir];
+    });
+
+    for (const auto& [i, child] : enumerate(trivialLinks)) {
+      const auto& bounds =
+          dynamic_cast<const RectangleBounds&>(child->surface().bounds());
+      Transform3 ltransform = itransform * child->surface().transform(gctx);
+      double half = dirX ? bounds.halfLengthX() : bounds.halfLengthY();
+      double min = ltransform.translation()[sortingDir] - half;
+      double max = ltransform.translation()[sortingDir] + half;
+      if (i == 0) {
+        edges.push_back(min);
+      }
+      edges.push_back(max);
+    }
+
+    ACTS_VERBOSE("~> Determined bin edges to be " << printEdges(edges));
+
+    Axis axis{AxisBound, edges};
+
+    auto grid = GridPortalLink::make(m_surface, m_direction, std::move(axis));
+    for (const auto& [i, child] : enumerate(trivialLinks)) {
+      grid->atLocalBins({i + 1}) = &child->volume();
+    }
+
+    return grid;
+
   } else {
     throw std::invalid_argument{"Unsupported surface type"};
   }

Core/src/Geometry/GridPortalLink.cpp

@@ -10,7 +10,10 @@
 
 #include "Acts/Surfaces/PlaneSurface.hpp"
 #include "Acts/Surfaces/RadialBounds.hpp"
+#include "Acts/Surfaces/RectangleBounds.hpp"
+#include "Acts/Utilities/AxisDefinitions.hpp"
 
+#include <iostream>
 #include <numbers>
 
 namespace Acts {
@@ -66,8 +69,19 @@ std::unique_ptr<GridPortalLink> GridPortalLink::make(
   } else if (const auto* plane =
                  dynamic_cast<const PlaneSurface*>(surface.get());
              plane != nullptr) {
-    throw std::invalid_argument{"Plane surface is not implemented yet"};
-
+    const auto& bounds = dynamic_cast<const RectangleBounds&>(plane->bounds());
+    if (direction != AxisDirection::AxisX &&
+        direction != AxisDirection::AxisY) {
+      throw std::invalid_argument{"Invalid binning direction"};
+    }
+    double min = (direction == AxisDirection::AxisX)
+                     ? bounds.get(RectangleBounds::eMinX)
+                     : bounds.get(RectangleBounds::eMinY);
+    double max = (direction == AxisDirection::AxisX)
+                     ? bounds.get(RectangleBounds::eMaxX)
+                     : bounds.get(RectangleBounds::eMaxY);
+    grid =
+        GridPortalLink::make(surface, direction, Axis{AxisBound, min, max, 1});
   } else {
     throw std::invalid_argument{"Surface type is not supported"};
   }
@@ -201,6 +215,39 @@ void GridPortalLink::checkConsistency(const DiscSurface& disc) const {
   }
 }
 
+void GridPortalLink::checkConsistency(const PlaneSurface& plane) const {
+  constexpr auto tolerance = s_onSurfaceTolerance;
+  auto same = [](auto a, auto b) { return std::abs(a - b) < tolerance; };
+
+  const auto* bounds = dynamic_cast<const RectangleBounds*>(&plane.bounds());
+  if (bounds == nullptr) {
+    throw std::invalid_argument(
+        "GridPortalLink: PlaneBounds: invalid bounds type.");
+  }
+  auto check = [&bounds, same](const IAxis& axis, AxisDirection dir) {
+    double min = (dir == AxisDirection::AxisX)
+                     ? bounds->get(RectangleBounds::eMinX)
+                     : bounds->get(RectangleBounds::eMinY);
+    double max = (dir == AxisDirection::AxisX)
+                     ? bounds->get(RectangleBounds::eMaxX)
+                     : bounds->get(RectangleBounds::eMaxY);
+    if (!same(axis.getMin(), min) || !same(axis.getMax(), max)) {
+      throw std::invalid_argument(
+          "GridPortalLink: PlaneBounds: invalid setup.");
+    }
+  };
+
+  if (dim() == 1) {
+    const IAxis& axisLoc0 = *grid().axes().front();
+    check(axisLoc0, direction());
+  } else {  // DIM == 2
+    const auto& axisLoc0 = *grid().axes().front();
+    const auto& axisLoc1 = *grid().axes().back();
+    check(axisLoc0, AxisDirection::AxisX);
+    check(axisLoc1, AxisDirection::AxisY);
+  }
+}
+
 void GridPortalLink::printContents(std::ostream& os) const {
   std::size_t dim = grid().axes().size();
   os << "----- GRID " << dim << "d -----" << std::endl;
@@ -419,4 +466,41 @@ std::unique_ptr<GridPortalLink> GridPortalLink::extendTo2dImpl(
   }
 }
 
+std::unique_ptr<GridPortalLink> GridPortalLink::extendTo2dImpl(
+    const std::shared_ptr<PlaneSurface>& surface, const IAxis* other) const {
+  assert(dim() == 1);
+
+  const auto* bounds = dynamic_cast<const RectangleBounds*>(&surface->bounds());
+  if (bounds == nullptr) {
+    throw std::invalid_argument(
+        "GridPortalLink: RectangleBounds: invalid bounds type.");
+  }
+
+  bool dirX = direction() == AxisDirection::AxisX;
+  const auto& thisAxis = *grid().axes().front();
+
+  double minExtend = dirX ? bounds->get(RectangleBounds::BoundValues::eMinY)
+                          : bounds->get(RectangleBounds::BoundValues::eMinX);
+  double maxExtend = dirX ? bounds->get(RectangleBounds::BoundValues::eMaxY)
+                          : bounds->get(RectangleBounds::BoundValues::eMaxX);
+
+  FillDirection fillDir = dirX ? FillDirection::loc1 : FillDirection::loc0;
+
+  auto grid = thisAxis.visit([&](const auto& axis0) {
+    Axis axisExtend{AxisBound, minExtend, maxExtend, 1};
+    const IAxis* axis = other != nullptr ? other : &axisExtend;
+    return axis->visit(
+        [&](const auto& axis1) -> std::unique_ptr<GridPortalLink> {
+          if (dirX) {
+            return GridPortalLink::make(surface, axis0, axis1);
+          } else {
+            return GridPortalLink::make(surface, axis1, axis0);
+          }
+        });
+  });
+
+  fillGrid1dTo2d(fillDir, *this, *grid);
+  return grid;
+}
+
 }  // namespace Acts

Core/src/Geometry/GridPortalLinkMerging.cpp

@@ -257,7 +257,8 @@ std::unique_ptr<PortalLinkBase> mergeGridPortals(
         dynamic_cast<const DiscSurface&>(*surfaceB), direction, true, logger);
   } else if (const auto* planeA = dynamic_cast<const PlaneSurface*>(surfaceA);
              planeA != nullptr) {
-    throw std::logic_error{"Plane surfaces not implemented yet"};
+    std::tie(mergedSurface, reversed) = planeA->mergedWith(
+        dynamic_cast<const PlaneSurface&>(*surfaceB), direction, logger);
   } else {
     throw std::invalid_argument{"Unsupported surface type"};
   }
@@ -442,6 +443,7 @@ std::unique_ptr<PortalLinkBase> mergeGridPortals(const GridPortalLink* a,
 
   const auto* cylinder = dynamic_cast<const CylinderSurface*>(&a->surface());
   const auto* disc = dynamic_cast<const DiscSurface*>(&a->surface());
+  const auto* plane = dynamic_cast<const PlaneSurface*>(&a->surface());
 
   if (a->dim() == b->dim()) {
     ACTS_VERBOSE("Grid both have same dimension: " << a->dim());
@@ -454,8 +456,11 @@ std::unique_ptr<PortalLinkBase> mergeGridPortals(const GridPortalLink* a,
       return mergeGridPortals(a, b, disc,
                               &dynamic_cast<const DiscSurface&>(b->surface()),
                               AxisR, AxisPhi, direction, logger);
+    } else if (plane != nullptr) {
+      return mergeGridPortals(a, b, plane,
+                              &dynamic_cast<const PlaneSurface&>(b->surface()),
+                              AxisX, AxisY, direction, logger);
     } else {
-      // @TODO: Support PlaneSurface
       ACTS_VERBOSE("Surface type is not supported here, falling back");
       return nullptr;
     }

Core/src/Geometry/PortalLinkBase.cpp

@@ -14,6 +14,7 @@
 #include "Acts/Surfaces/CylinderSurface.hpp"
 #include "Acts/Surfaces/DiscSurface.hpp"
 #include "Acts/Surfaces/RadialBounds.hpp"
+#include "Acts/Surfaces/RectangleBounds.hpp"
 #include "Acts/Surfaces/RegularSurface.hpp"
 #include "Acts/Utilities/ThrowAssert.hpp"
 
@@ -56,6 +57,18 @@ void PortalLinkBase::checkMergePreconditions(const PortalLinkBase& a,
                      dynamic_cast<const RadialBounds*>(&discB->bounds()),
                  "DiscSurface bounds must be RadialBounds");
 
+  } else if (const auto* planeA = dynamic_cast<const PlaneSurface*>(&surfaceA);
+             planeA != nullptr) {
+    const auto* planeB = dynamic_cast<const PlaneSurface*>(&surfaceB);
+    throw_assert(planeB != nullptr,
+                 "Cannot merge PlaneSurface with non-PlaneSurface");
+    throw_assert(
+        direction == AxisDirection::AxisX || direction == AxisDirection::AxisY,
+        "Invalid binning direction: " + axisDirectionName(direction));
+
+    throw_assert(dynamic_cast<const RectangleBounds*>(&planeA->bounds()) &&
+                     dynamic_cast<const RectangleBounds*>(&planeB->bounds()),
+                 "PlaneSurface bounds must be RectangleBounds");
   } else {
     throw std::logic_error{"Surface type is not supported"};
   }