diff --git a/Core/include/Acts/Geometry/GridPortalLink.hpp b/Core/include/Acts/Geometry/GridPortalLink.hpp index 7ba790b4d82..d94ba5fc749 100644 --- a/Core/include/Acts/Geometry/GridPortalLink.hpp +++ b/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 @@ -67,6 +67,9 @@ class GridPortalLink : public PortalLinkBase { } else if (dynamic_cast(surface.get()) != nullptr && direction != AxisR && direction != AxisPhi) { throw std::invalid_argument{"Invalid binning direction"}; + } else if (dynamic_cast(surface.get()) != nullptr && + direction != AxisX && direction != AxisY) { + throw std::invalid_argument{"Invalid binning direction"}; } return std::make_unique>( @@ -90,6 +93,8 @@ class GridPortalLink : public PortalLinkBase { direction = AxisDirection::AxisRPhi; } else if (dynamic_cast(surface.get()) != nullptr) { direction = AxisDirection::AxisR; + } else if (dynamic_cast(surface.get()) != nullptr) { + direction = AxisDirection::AxisX; } return std::make_unique>( @@ -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 extendTo2dImpl( const std::shared_ptr& 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 extendTo2dImpl( + const std::shared_ptr& 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(m_surface.get())) { + checkConsistency(*plane); + + if (direction == AxisX) { + m_projection = &projection; + } else if (direction == AxisDirection::AxisY) { + m_projection = &projection; + } 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(m_surface)) { return extendTo2dImpl(disc, other); + } else if (auto plane = + std::dynamic_pointer_cast(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 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)); } diff --git a/Core/src/Geometry/CompositePortalLink.cpp b/Core/src/Geometry/CompositePortalLink.cpp index 19d9ea3c02a..461c0cd72a5 100644 --- a/Core/src/Geometry/CompositePortalLink.cpp +++ b/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 +#include #include #include #include @@ -48,7 +51,9 @@ std::shared_ptr mergedSurface(const Surface& a, return merged; } else if (const auto* planeA = dynamic_cast(&a); planeA != nullptr) { - throw std::logic_error{"Plane surfaces not implemented yet"}; + const auto& planeB = dynamic_cast(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 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 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(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"}; } diff --git a/Core/src/Geometry/GridPortalLink.cpp b/Core/src/Geometry/GridPortalLink.cpp index a0e429708d2..6a0917ba2ce 100644 --- a/Core/src/Geometry/GridPortalLink.cpp +++ b/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 #include namespace Acts { @@ -66,8 +69,19 @@ std::unique_ptr GridPortalLink::make( } else if (const auto* plane = dynamic_cast(surface.get()); plane != nullptr) { - throw std::invalid_argument{"Plane surface is not implemented yet"}; - + const auto& bounds = dynamic_cast(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(&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::extendTo2dImpl( } } +std::unique_ptr GridPortalLink::extendTo2dImpl( + const std::shared_ptr& surface, const IAxis* other) const { + assert(dim() == 1); + + const auto* bounds = dynamic_cast(&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 { + if (dirX) { + return GridPortalLink::make(surface, axis0, axis1); + } else { + return GridPortalLink::make(surface, axis1, axis0); + } + }); + }); + + fillGrid1dTo2d(fillDir, *this, *grid); + return grid; +} + } // namespace Acts diff --git a/Core/src/Geometry/GridPortalLinkMerging.cpp b/Core/src/Geometry/GridPortalLinkMerging.cpp index ae49bedfa46..c013e7a19b0 100644 --- a/Core/src/Geometry/GridPortalLinkMerging.cpp +++ b/Core/src/Geometry/GridPortalLinkMerging.cpp @@ -257,7 +257,8 @@ std::unique_ptr mergeGridPortals( dynamic_cast(*surfaceB), direction, true, logger); } else if (const auto* planeA = dynamic_cast(surfaceA); planeA != nullptr) { - throw std::logic_error{"Plane surfaces not implemented yet"}; + std::tie(mergedSurface, reversed) = planeA->mergedWith( + dynamic_cast(*surfaceB), direction, logger); } else { throw std::invalid_argument{"Unsupported surface type"}; } @@ -442,6 +443,7 @@ std::unique_ptr mergeGridPortals(const GridPortalLink* a, const auto* cylinder = dynamic_cast(&a->surface()); const auto* disc = dynamic_cast(&a->surface()); + const auto* plane = dynamic_cast(&a->surface()); if (a->dim() == b->dim()) { ACTS_VERBOSE("Grid both have same dimension: " << a->dim()); @@ -454,8 +456,11 @@ std::unique_ptr mergeGridPortals(const GridPortalLink* a, return mergeGridPortals(a, b, disc, &dynamic_cast(b->surface()), AxisR, AxisPhi, direction, logger); + } else if (plane != nullptr) { + return mergeGridPortals(a, b, plane, + &dynamic_cast(b->surface()), + AxisX, AxisY, direction, logger); } else { - // @TODO: Support PlaneSurface ACTS_VERBOSE("Surface type is not supported here, falling back"); return nullptr; } diff --git a/Core/src/Geometry/PortalLinkBase.cpp b/Core/src/Geometry/PortalLinkBase.cpp index 16826172a2d..58560103633 100644 --- a/Core/src/Geometry/PortalLinkBase.cpp +++ b/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(&discB->bounds()), "DiscSurface bounds must be RadialBounds"); + } else if (const auto* planeA = dynamic_cast(&surfaceA); + planeA != nullptr) { + const auto* planeB = dynamic_cast(&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(&planeA->bounds()) && + dynamic_cast(&planeB->bounds()), + "PlaneSurface bounds must be RectangleBounds"); } else { throw std::logic_error{"Surface type is not supported"}; } diff --git a/Tests/UnitTests/Core/Geometry/PortalLinkTests.cpp b/Tests/UnitTests/Core/Geometry/PortalLinkTests.cpp index 90b39c63d4d..14f6ee23d33 100644 --- a/Tests/UnitTests/Core/Geometry/PortalLinkTests.cpp +++ b/Tests/UnitTests/Core/Geometry/PortalLinkTests.cpp @@ -15,12 +15,14 @@ #include "Acts/Definitions/Algebra.hpp" #include "Acts/Definitions/Units.hpp" #include "Acts/Geometry/CompositePortalLink.hpp" +#include "Acts/Geometry/CuboidVolumeBounds.hpp" #include "Acts/Geometry/CylinderVolumeBounds.hpp" #include "Acts/Geometry/GridPortalLink.hpp" #include "Acts/Geometry/TrackingVolume.hpp" #include "Acts/Geometry/TrivialPortalLink.hpp" #include "Acts/Surfaces/CylinderSurface.hpp" #include "Acts/Surfaces/RadialBounds.hpp" +#include "Acts/Surfaces/RectangleBounds.hpp" #include "Acts/Surfaces/SurfaceMergingException.hpp" #include "Acts/Tests/CommonHelpers/FloatComparisons.hpp" #include "Acts/Utilities/AxisDefinitions.hpp" @@ -465,7 +467,124 @@ BOOST_AUTO_TEST_CASE(Disc) { } BOOST_AUTO_TEST_CASE(Plane) { - // @TODO: Add plane tests + using enum AxisType; + using enum AxisBoundaryType; + BOOST_TEST_CONTEXT("1D") { + // Initialize surfaces + auto rBounds = std::make_shared(30_mm, 100_mm); + auto planeX = + Surface::makeShared(Transform3::Identity(), rBounds); + auto planeY = + Surface::makeShared(Transform3::Identity(), rBounds); + + // Volume for bin testing + auto vol = std::make_shared( + Transform3::Identity(), + std::make_shared(30_mm, 40_mm, 100_mm)); + + // Initialize grids + auto gridX = GridPortalLink::make(planeX, AxisDirection::AxisX, + Axis{AxisBound, -30_mm, 30_mm, 3}); + auto gridY = GridPortalLink::make(planeY, AxisDirection::AxisY, + Axis{AxisBound, -100_mm, 100_mm, 3}); + + // Check grid X + BOOST_REQUIRE_NE(gridX, nullptr); + BOOST_CHECK_EQUAL(gridX->grid().axes().size(), 1); + const auto& axisX = *gridX->grid().axes().front(); + Axis axisXExpected{AxisBound, -30_mm, 30_mm, 3}; + BOOST_CHECK_EQUAL(axisX, axisXExpected); + + // Check grid Y + BOOST_REQUIRE_NE(gridY, nullptr); + BOOST_CHECK_EQUAL(gridY->grid().axes().size(), 1); + const auto& axisY = *gridY->grid().axes().front(); + Axis axisYExpected{AxisBound, -100_mm, 100_mm, 3}; + BOOST_CHECK_EQUAL(axisY, axisYExpected); + + // Check gridX/gridX bin content + auto checkAllBins = [&](const auto& grid) { + visitBins(grid, [&](const TrackingVolume* content) { + BOOST_CHECK_EQUAL(content, vol.get()); + }); + }; + + gridX->setVolume(vol.get()); + checkAllBins(*gridX); + + gridY->setVolume(vol.get()); + checkAllBins(*gridY); + + // Test making 2D grids from the 1D ones + + // Test grid X + auto gridX2d = gridX->extendTo2d(nullptr); + BOOST_REQUIRE(gridX2d); + BOOST_CHECK_EQUAL(gridX2d->grid().axes().size(), 2); + const auto* axisX1 = gridX2d->grid().axes().front(); + BOOST_CHECK_EQUAL(*axisX1, axisXExpected); + const auto* axisX2 = gridX2d->grid().axes().back(); + BOOST_CHECK_EQUAL(axisX2->getMin(), -100_mm); + BOOST_CHECK_EQUAL(axisX2->getMax(), 100_mm); + BOOST_CHECK_EQUAL(axisX2->getNBins(), 1); + BOOST_CHECK_EQUAL(axisX2->getType(), AxisType::Equidistant); + BOOST_CHECK_EQUAL(axisX2->getBoundaryType(), AxisBoundaryType::Bound); + + // Test grid X explicit + Axis axisYExplicit{AxisClosed, -100_mm, 100_mm, 3}; + auto gridX2dExplicit = gridX->extendTo2d(&axisYExplicit); + BOOST_REQUIRE(gridX2dExplicit); + BOOST_CHECK_EQUAL(gridX2dExplicit->grid().axes().size(), 2); + axisX1 = gridX2dExplicit->grid().axes().front(); + axisX2 = gridX2dExplicit->grid().axes().back(); + BOOST_CHECK_EQUAL(*axisX1, axisXExpected); + BOOST_CHECK_EQUAL(*axisX2, axisYExplicit); + + checkAllBins( + dynamic_cast< + GridPortalLinkT&>( + *gridX2dExplicit)); + + // Test grid Y + auto gridY2d = gridY->extendTo2d(nullptr); + BOOST_REQUIRE(gridY2d); + BOOST_CHECK_EQUAL(gridY2d->grid().axes().size(), 2); + const auto* axisY1 = gridY2d->grid().axes().front(); + BOOST_CHECK_EQUAL(axisY1->getMin(), -30_mm); + BOOST_CHECK_EQUAL(axisY1->getMax(), 30_mm); + BOOST_CHECK_EQUAL(axisY1->getNBins(), 1); + BOOST_CHECK_EQUAL(axisY1->getType(), AxisType::Equidistant); + BOOST_CHECK_EQUAL(axisY1->getBoundaryType(), AxisBoundaryType::Bound); + const auto* axisY2 = gridY2d->grid().axes().back(); + BOOST_CHECK_EQUAL(*axisY2, axisYExpected); + + // Test grid Y explicit + Axis axisXExplicit{AxisClosed, -30_mm, 30_mm, 3}; + auto gridY2dExplicit = gridY->extendTo2d(&axisXExplicit); + BOOST_REQUIRE(gridY2dExplicit); + BOOST_CHECK_EQUAL(gridY2dExplicit->grid().axes().size(), 2); + axisY1 = gridY2dExplicit->grid().axes().front(); + axisY2 = gridY2dExplicit->grid().axes().back(); + BOOST_CHECK_EQUAL(*axisY1, axisXExplicit); + BOOST_CHECK_EQUAL(*axisY2, axisYExpected); + + checkAllBins( + dynamic_cast< + GridPortalLinkT&>( + *gridY2dExplicit)); + } + BOOST_TEST_CONTEXT("2D") { + auto rBounds = std::make_shared(30_mm, 100_mm); + + auto plane = + Surface::makeShared(Transform3::Identity(), rBounds); + + Axis xAxis{AxisBound, -30_mm, 30_mm, 5}; + Axis yAxis{AxisBound, -100_mm, 100_mm, 5}; + + auto grid = GridPortalLink::make(plane, xAxis, yAxis); + BOOST_REQUIRE_NE(grid, nullptr); + } } BOOST_AUTO_TEST_CASE(FromTrivial) { @@ -567,6 +686,45 @@ BOOST_AUTO_TEST_CASE(FromTrivial) { gridPhi->resolveVolume(gctx, Vector2{90_mm, 10_degree}).value(), vol.get()); } + BOOST_TEST_CONTEXT("Plane") { + auto rBounds = std::make_shared(30_mm, 100_mm); + + auto plane = + Surface::makeShared(Transform3::Identity(), rBounds); + + auto vol = std::make_shared( + Transform3::Identity(), + std::make_shared(30_mm, 40_mm, 100_mm)); + + auto trivial = std::make_unique(plane, *vol); + BOOST_REQUIRE(trivial); + + // Doesn't matter which position + BOOST_CHECK_EQUAL( + trivial->resolveVolume(gctx, Vector2{10_mm, 20_mm}).value(), vol.get()); + + auto gridX = trivial->makeGrid(AxisDirection::AxisX); + BOOST_REQUIRE(gridX); + + BOOST_CHECK_EQUAL(gridX->grid().axes().size(), 1); + BOOST_CHECK_EQUAL(gridX->surface().bounds(), plane->bounds()); + Axis axisXExpected{AxisBound, -30_mm, 30_mm, 1}; + BOOST_CHECK_EQUAL(*gridX->grid().axes().front(), axisXExpected); + + BOOST_CHECK_EQUAL(gridX->resolveVolume(gctx, Vector2{20_mm, 10_mm}).value(), + vol.get()); + + auto gridY = trivial->makeGrid(AxisDirection::AxisY); + BOOST_REQUIRE(gridY); + + BOOST_CHECK_EQUAL(gridY->grid().axes().size(), 1); + BOOST_CHECK_EQUAL(gridY->surface().bounds(), plane->bounds()); + Axis axisYExpected{AxisBound, -100_mm, 100_mm, 1}; + BOOST_CHECK_EQUAL(*gridY->grid().axes().front(), axisYExpected); + + BOOST_CHECK_EQUAL(gridY->resolveVolume(gctx, Vector2{15_mm, 20_mm}).value(), + vol.get()); + } } BOOST_AUTO_TEST_SUITE_END() // GridConstruction @@ -2048,6 +2206,391 @@ BOOST_AUTO_TEST_CASE(PhiDirection) { BOOST_AUTO_TEST_SUITE_END() // MergeCrossDisc +BOOST_AUTO_TEST_SUITE(Merging1dPlane) + +BOOST_AUTO_TEST_CASE(ColinearMerge) { + auto rBounds1 = std::make_shared(30_mm, 100_mm); + + auto plane1 = + Surface::makeShared(Transform3::Identity(), rBounds1); + + auto gridX1 = GridPortalLink::make(plane1, AxisDirection::AxisX, + Axis{AxisBound, -30_mm, 30_mm, 6}); + auto gridY1 = GridPortalLink::make(plane1, AxisDirection::AxisY, + Axis{AxisBound, -100_mm, 100_mm, 10}); + + auto rBounds2 = std::make_shared(80_mm, 100_mm); + Translation3 offsetX{110_mm, 0., 0.}; + auto planeX2 = Surface::makeShared( + Transform3::Identity() * offsetX, rBounds2); + + auto rBounds3 = std::make_shared(30_mm, 20_mm); + Translation3 offsetY{0, 120_mm, 0.}; + auto planeY2 = Surface::makeShared( + Transform3::Identity() * offsetY, rBounds3); + + auto gridX2 = GridPortalLink::make(planeX2, AxisDirection::AxisX, + Axis{AxisBound, -80_mm, 80_mm, 16}); + auto gridY2 = GridPortalLink::make(planeY2, AxisDirection::AxisY, + Axis{AxisBound, -20_mm, 20_mm, 2}); + + auto mergedPtrX = + GridPortalLink::merge(*gridX1, *gridX2, AxisDirection::AxisX, *logger); + BOOST_REQUIRE(mergedPtrX); + const auto* mergedX = dynamic_cast(mergedPtrX.get()); + BOOST_REQUIRE_NE(mergedX, nullptr); + + BOOST_CHECK_EQUAL(mergedX->grid().axes().size(), 1); + Axis axisXExpected{AxisBound, -110_mm, 110_mm, 22}; + BOOST_CHECK_EQUAL(*mergedX->grid().axes().front(), axisXExpected); + + auto mergedPtrY = + GridPortalLink::merge(*gridY1, *gridY2, AxisDirection::AxisY, *logger); + BOOST_REQUIRE(mergedPtrY); + const auto* mergedY = dynamic_cast(mergedPtrY.get()); + BOOST_REQUIRE_NE(mergedY, nullptr); + + BOOST_CHECK_EQUAL(mergedY->grid().axes().size(), 1); + Axis axisYExpected{AxisBound, -120_mm, 120_mm, 12}; + BOOST_CHECK_EQUAL(*mergedY->grid().axes().front(), axisYExpected); +} + +BOOST_AUTO_TEST_CASE(ParallelMerge) { + auto rBounds = std::make_shared(30_mm, 100_mm); + auto plane1 = + Surface::makeShared(Transform3::Identity(), rBounds); + + auto grid1X = GridPortalLink::make(plane1, AxisDirection::AxisX, + Axis{AxisBound, -30_mm, 30_mm, 6}); + + auto grid1Y = GridPortalLink::make(plane1, AxisDirection::AxisY, + Axis{AxisBound, -100_mm, 100_mm, 5}); + + Translation3 offsetX{60_mm, 0, 0.}; + auto plane2 = Surface::makeShared( + Transform3::Identity() * offsetX, rBounds); + auto grid2 = GridPortalLink::make(plane2, AxisDirection::AxisY, + Axis{AxisBound, -100_mm, 100_mm, 5}); + + Translation3 offsetY{0, 200_mm, 0.}; + auto plane3 = Surface::makeShared( + Transform3::Identity() * offsetY, rBounds); + auto grid3 = GridPortalLink::make(plane3, AxisDirection::AxisX, + Axis{AxisBound, -30_mm, 30_mm, 6}); + + auto mergedPtrX = + GridPortalLink::merge(*grid1Y, *grid2, AxisDirection::AxisX, *logger); + BOOST_REQUIRE(mergedPtrX); + const auto* mergedX = dynamic_cast(mergedPtrX.get()); + BOOST_REQUIRE_NE(mergedX, nullptr); + + BOOST_CHECK_EQUAL(mergedX->grid().axes().size(), 2); + const auto& axisX1 = *mergedX->grid().axes().front(); + const auto& axisX2 = *mergedX->grid().axes().back(); + Axis axisX1Expected{AxisBound, -60_mm, 60_mm, 2}; + BOOST_CHECK_EQUAL(axisX1, axisX1Expected); + Axis axisX2Expected{AxisBound, -100_mm, 100_mm, 5}; + BOOST_CHECK_EQUAL(axisX2, axisX2Expected); + + auto mergedPtrY = + GridPortalLink::merge(*grid1X, *grid3, AxisDirection::AxisY, *logger); + BOOST_REQUIRE(mergedPtrY); + const auto* mergedY = dynamic_cast(mergedPtrY.get()); + BOOST_REQUIRE_NE(mergedY, nullptr); + + BOOST_CHECK_EQUAL(mergedY->grid().axes().size(), 2); + const auto& axisY1 = *mergedY->grid().axes().front(); + const auto& axisY2 = *mergedY->grid().axes().back(); + Axis axisY1Expected{AxisBound, -30_mm, 30_mm, 6}; + BOOST_CHECK_EQUAL(axisY1, axisY1Expected); + Axis axisY2Expected{AxisBound, -200_mm, 200_mm, 2}; + BOOST_CHECK_EQUAL(axisY2, axisY2Expected); +} + +BOOST_AUTO_TEST_CASE(BinFilling) { + // Volumes for bin content checking + auto vol1 = std::make_shared( + Transform3::Identity(), + std::make_shared(30_mm, 40_mm, 100_mm)); + auto vol2 = std::make_shared( + Transform3::Identity(), + std::make_shared(30_mm, 40_mm, 100_mm)); + + auto rBounds = std::make_shared(30_mm, 100_mm); + auto plane1 = + Surface::makeShared(Transform3::Identity(), rBounds); + + auto grid1X = GridPortalLink::make(plane1, AxisDirection::AxisX, + Axis{AxisBound, -30_mm, 30_mm, 2}); + grid1X->setVolume(vol1.get()); + + auto grid1Y = GridPortalLink::make(plane1, AxisDirection::AxisY, + Axis{AxisBound, -100_mm, 100_mm, 2}); + grid1Y->setVolume(vol1.get()); + + Translation3 offsetX{60_mm, 0., 0.}; + auto plane2 = Surface::makeShared( + Transform3::Identity() * offsetX, rBounds); + auto grid2 = GridPortalLink::make(plane2, AxisDirection::AxisX, + Axis{AxisBound, -30_mm, 30_mm, 2}); + grid2->setVolume(vol2.get()); + + Translation3 offsetY{0., 200_mm, 0.}; + auto plane3 = Surface::makeShared( + Transform3::Identity() * offsetY, rBounds); + auto grid3 = GridPortalLink::make(plane3, AxisDirection::AxisY, + Axis{AxisBound, -100_mm, 100_mm, 2}); + grid3->setVolume(vol2.get()); + + auto mergedPtrX = + GridPortalLink::merge(*grid1X, *grid2, AxisDirection::AxisX, *logger); + + auto mergedPtrY = + GridPortalLink::merge(*grid1Y, *grid3, AxisDirection::AxisY, *logger); + + using merged_type = + GridPortalLinkT>; + + const auto* mergedX = dynamic_cast(mergedPtrX.get()); + BOOST_REQUIRE(mergedX); + + grid1X->printContents(std::cout); + grid2->printContents(std::cout); + mergedX->printContents(std::cout); + + BOOST_CHECK_EQUAL(mergedX->grid().atLocalBins({1}), vol1.get()); + BOOST_CHECK_EQUAL(mergedX->grid().atLocalBins({2}), vol1.get()); + BOOST_CHECK_EQUAL(mergedX->grid().atLocalBins({3}), vol2.get()); + BOOST_CHECK_EQUAL(mergedX->grid().atLocalBins({4}), vol2.get()); + + const auto* mergedY = dynamic_cast(mergedPtrX.get()); + BOOST_REQUIRE(mergedY); + + BOOST_CHECK_EQUAL(mergedY->grid().atLocalBins({1}), vol1.get()); + BOOST_CHECK_EQUAL(mergedY->grid().atLocalBins({2}), vol1.get()); + BOOST_CHECK_EQUAL(mergedY->grid().atLocalBins({3}), vol2.get()); + BOOST_CHECK_EQUAL(mergedY->grid().atLocalBins({4}), vol2.get()); + + grid1X->printContents(std::cout); + grid2->printContents(std::cout); + grid3->printContents(std::cout); + mergedX->printContents(std::cout); + mergedY->printContents(std::cout); +} + +BOOST_AUTO_TEST_SUITE_END() // Merging1dPlane + +BOOST_AUTO_TEST_SUITE(Merging2dPlane) + +BOOST_AUTO_TEST_CASE(XYDirection) { + // Basic, because the parallel 1D case already tests this to some degree + auto rBounds = std::make_shared(30_mm, 100_mm); + auto plane1 = + Surface::makeShared(Transform3::Identity(), rBounds); + + auto grid1 = GridPortalLink::make(plane1, Axis{AxisBound, -30_mm, 30_mm, 10}, + Axis{AxisBound, -100_mm, 100_mm, 3}); + + Translation3 offsetX{60_mm, 0., 0.}; + auto plane2 = Surface::makeShared( + Transform3::Identity() * offsetX, rBounds); + auto grid2 = GridPortalLink::make(plane2, Axis{AxisBound, -30_mm, 30_mm, 10}, + Axis{AxisBound, -100_mm, 100_mm, 3}); + + Translation3 offsetY{0., 200_mm, 0.}; + auto plane3 = Surface::makeShared( + Transform3::Identity() * offsetY, rBounds); + auto grid3 = GridPortalLink::make(plane3, Axis{AxisBound, -30_mm, 30_mm, 10}, + Axis{AxisBound, -100_mm, 100_mm, 3}); + + auto mergedPtrX = + GridPortalLink::merge(*grid1, *grid2, AxisDirection::AxisX, *logger); + + auto mergedPtrY = + GridPortalLink::merge(*grid1, *grid3, AxisDirection::AxisY, *logger); + + BOOST_REQUIRE(mergedPtrX); + const auto* mergedX = dynamic_cast(mergedPtrX.get()); + BOOST_REQUIRE_NE(mergedX, nullptr); + BOOST_CHECK_EQUAL(mergedX->grid().axes().size(), 2); + const auto& axisX1 = *mergedX->grid().axes().front(); + const auto& axisX2 = *mergedX->grid().axes().back(); + + BOOST_CHECK_EQUAL(axisX1.getMin(), -60_mm); + BOOST_CHECK_EQUAL(axisX1.getMax(), 60_mm); + BOOST_CHECK_EQUAL(axisX1.getNBins(), 20); + BOOST_CHECK_EQUAL(axisX1.getType(), AxisType::Equidistant); + BOOST_CHECK_EQUAL(axisX1.getBoundaryType(), AxisBoundaryType::Bound); + BOOST_CHECK_EQUAL(axisX2.getMin(), -100_mm); + BOOST_CHECK_EQUAL(axisX2.getMax(), 100_mm); + BOOST_CHECK_EQUAL(axisX2.getNBins(), 3); + BOOST_CHECK_EQUAL(axisX2.getType(), AxisType::Equidistant); + + BOOST_REQUIRE(mergedPtrY); + const auto* mergedY = dynamic_cast(mergedPtrY.get()); + BOOST_REQUIRE_NE(mergedY, nullptr); + BOOST_CHECK_EQUAL(mergedY->grid().axes().size(), 2); + const auto& axisY1 = *mergedY->grid().axes().front(); + const auto& axisY2 = *mergedY->grid().axes().back(); + + BOOST_CHECK_EQUAL(axisY1.getMin(), -30_mm); + BOOST_CHECK_EQUAL(axisY1.getMax(), 30_mm); + BOOST_CHECK_EQUAL(axisY1.getNBins(), 10); + BOOST_CHECK_EQUAL(axisY1.getType(), AxisType::Equidistant); + BOOST_CHECK_EQUAL(axisY1.getBoundaryType(), AxisBoundaryType::Bound); + BOOST_CHECK_EQUAL(axisY2.getMin(), -200_mm); + BOOST_CHECK_EQUAL(axisY2.getMax(), 200_mm); + BOOST_CHECK_EQUAL(axisY2.getNBins(), 6); + BOOST_CHECK_EQUAL(axisY2.getType(), AxisType::Equidistant); +} + +BOOST_AUTO_TEST_CASE(BinFilling) { + // Volumes for bin content checking + // Volume shape/transform is irrelevant, only used for pointer identity + auto vol1 = std::make_shared( + Transform3::Identity(), + std::make_shared(30_mm, 40_mm, 100_mm)); + + auto vol2 = std::make_shared( + Transform3::Identity(), + std::make_shared(30_mm, 40_mm, 100_mm)); + + auto fillCheckerBoard = [&](auto& grid) { + auto loc = grid.numLocalBins(); + for (std::size_t i = 1; i <= loc[0]; ++i) { + for (std::size_t j = 1; j <= loc[1]; ++j) { + grid.atLocalBins({i, j}) = (i + j) % 2 == 0 ? vol1.get() : vol2.get(); + } + } + }; + + auto checkCheckerBoard = [&](const auto& grid) { + auto loc = grid.numLocalBins(); + for (std::size_t i = 1; i <= loc[0]; ++i) { + for (std::size_t j = 1; j <= loc[1]; ++j) { + const auto* vol = grid.atLocalBins({i, j}); + if (vol != ((i + j) % 2 == 0 ? vol1.get() : vol2.get())) { + BOOST_ERROR("Is not a checkerboard pattern"); + return; + } + } + } + }; + + // Basic, because the parallel 1D case already tests this to some degree + auto rBounds = std::make_shared(30_mm, 100_mm); + auto plane1 = + Surface::makeShared(Transform3::Identity(), rBounds); + + auto grid1 = GridPortalLink::make(plane1, Axis{AxisBound, -30_mm, 30_mm, 2}, + Axis{AxisBound, -100_mm, 100_mm, 2}); + + Translation3 offsetX{60_mm, 0., 0.}; + auto plane2 = Surface::makeShared( + Transform3::Identity() * offsetX, rBounds); + auto grid2 = GridPortalLink::make(plane2, Axis{AxisBound, -30_mm, 30_mm, 2}, + Axis{AxisBound, -100_mm, 100_mm, 2}); + + Translation3 offsetY{0., 200_mm, 0.}; + auto plane3 = Surface::makeShared( + Transform3::Identity() * offsetY, rBounds); + auto grid3 = GridPortalLink::make(plane3, Axis{AxisBound, -30_mm, 30_mm, 2}, + Axis{AxisBound, -100_mm, 100_mm, 2}); + + fillCheckerBoard(grid1->grid()); + checkCheckerBoard(grid1->grid()); + + fillCheckerBoard(grid2->grid()); + checkCheckerBoard(grid2->grid()); + + fillCheckerBoard(grid3->grid()); + checkCheckerBoard(grid3->grid()); + + auto mergedPtrX = + GridPortalLink::merge(*grid1, *grid2, AxisDirection::AxisX, *logger); + auto mergedPtrY = + GridPortalLink::merge(*grid1, *grid3, AxisDirection::AxisY, *logger); + + using merged_type = + GridPortalLinkT, + Axis>; + + const auto* mergedX = dynamic_cast(mergedPtrX.get()); + BOOST_REQUIRE(mergedX); + checkCheckerBoard(mergedX->grid()); + + const auto* mergedY = dynamic_cast(mergedPtrY.get()); + BOOST_REQUIRE(mergedY); + checkCheckerBoard(mergedY->grid()); + + // Fill a / b + grid1->setVolume(vol1.get()); + grid2->setVolume(vol2.get()); + grid3->setVolume(vol2.get()); + + mergedPtrX = + GridPortalLink::merge(*grid1, *grid2, AxisDirection::AxisX, *logger); + mergedPtrY = + GridPortalLink::merge(*grid1, *grid3, AxisDirection::AxisY, *logger); + + mergedX = dynamic_cast(mergedPtrX.get()); + BOOST_REQUIRE(mergedX); + + mergedY = dynamic_cast(mergedPtrY.get()); + BOOST_REQUIRE(mergedY); + + const auto* v1 = vol1.get(); + const auto* v2 = vol2.get(); + + std::vector> locationsX = { + {{-45_mm, 0_mm}, v1}, + {{-15_mm, 0_mm}, v1}, + {{15_mm, 0_mm}, v2}, + {{45_mm, 0_mm}, v2}, + }; + std::vector> locationsY = { + {{0_mm, -150_mm}, v1}, + {{0_mm, -50_mm}, v1}, + {{0_mm, 50_mm}, v2}, + {{0_mm, 150_mm}, v2}, + }; + + for (const auto& [loc, vol] : locationsX) { + BOOST_TEST_CONTEXT(loc.transpose()) + BOOST_CHECK_EQUAL(mergedX->resolveVolume(gctx, loc).value(), vol); + } + for (const auto& [loc, vol] : locationsY) { + BOOST_TEST_CONTEXT(loc.transpose()) + BOOST_CHECK_EQUAL(mergedY->resolveVolume(gctx, loc).value(), vol); + } + + std::vector> contentsX = { + {v1, v1}, + {v1, v1}, + {v2, v2}, + {v2, v2}, + }; + std::vector> contentsY = { + {v1, v1, v2, v2}, + {v1, v1, v2, v2}, + }; + + for (std::size_t i = 0; i < 4; ++i) { + for (std::size_t j = 0; j < 2; ++j) { + BOOST_CHECK_EQUAL(mergedX->grid().atLocalBins({i + 1, j + 1}), + contentsX.at(i).at(j)); + } + } + for (std::size_t i = 0; i < 2; ++i) { + for (std::size_t j = 0; j < 4; ++j) { + BOOST_CHECK_EQUAL(mergedY->grid().atLocalBins({i + 1, j + 1}), + contentsY.at(i).at(j)); + } + } +} + +BOOST_AUTO_TEST_SUITE_END() // Merging2dPlane + BOOST_AUTO_TEST_SUITE_END() // GridMerging BOOST_AUTO_TEST_CASE(CompositeConstruction) { @@ -2298,6 +2841,152 @@ BOOST_DATA_TEST_CASE(TrivialTrivial, grid123->resolveVolume(gctx, Vector2{40_mm, 190_mm}).value(), vol3.get()); } + BOOST_TEST_CONTEXT("PlaneXDirection") { + auto vol1 = std::make_shared( + Transform3::Identity(), + std::make_shared(30_mm, 40_mm, 100_mm)); + + auto vol2 = std::make_shared( + Transform3::Identity() * Translation3{Vector3::UnitX() * 60}, + std::make_shared(30_mm, 40_mm, 100_mm)); + + auto vol3 = std::make_shared( + Transform3::Identity() * Translation3{Vector3::UnitX() * 120}, + std::make_shared(30_mm, 40_mm, 100_mm)); + + auto rBounds = std::make_shared(30_mm, 40_mm); + auto plane1 = + Surface::makeShared(Transform3::Identity(), rBounds); + + auto plane2 = Surface::makeShared( + Transform3::Identity() * Translation3{Vector3::UnitX() * 60}, rBounds); + + auto plane3 = Surface::makeShared( + Transform3::Identity() * Translation3{Vector3::UnitX() * 120}, rBounds); + + auto trivial1 = std::make_unique(plane1, *vol1); + BOOST_REQUIRE(trivial1); + auto trivial2 = std::make_unique(plane2, *vol2); + BOOST_REQUIRE(trivial2); + auto trivial3 = std::make_unique(plane3, *vol3); + BOOST_REQUIRE(trivial3); + + auto grid1 = trivial1->makeGrid(AxisDirection::AxisX); + auto compGridTrivial = PortalLinkBase::merge( + std::move(grid1), std::make_unique(*trivial2), + AxisDirection::AxisX, *logger); + BOOST_REQUIRE(compGridTrivial); + BOOST_CHECK_EQUAL(dynamic_cast(*compGridTrivial) + .makeGrid(gctx, *logger), + nullptr); + + auto composite = + PortalLinkBase::merge(std::move(trivial1), std::move(trivial2), + AxisDirection::AxisX, *logger); + BOOST_REQUIRE(composite); + + auto grid12 = + dynamic_cast(*composite).makeGrid(gctx, *logger); + BOOST_REQUIRE(grid12); + + BOOST_CHECK_EQUAL(grid12->resolveVolume(gctx, Vector2{-30_mm, 0}).value(), + vol1.get()); + + BOOST_CHECK_EQUAL(grid12->resolveVolume(gctx, Vector2{30_mm, 0}).value(), + vol2.get()); + + composite = PortalLinkBase::merge(std::move(composite), std::move(trivial3), + AxisDirection::AxisX, *logger); + BOOST_REQUIRE(composite); + + auto grid123 = + dynamic_cast(*composite).makeGrid(gctx, *logger); + BOOST_REQUIRE(grid123); + + BOOST_CHECK_EQUAL( + grid123->resolveVolume(gctx, Vector2{-80_mm, 0_mm}).value(), + vol1.get()); + + BOOST_CHECK_EQUAL(grid123->resolveVolume(gctx, Vector2{0_mm, 0_mm}).value(), + vol2.get()); + + BOOST_CHECK_EQUAL( + grid123->resolveVolume(gctx, Vector2{80_mm, 0_mm}).value(), vol3.get()); + } + BOOST_TEST_CONTEXT("PlaneYDirection") { + auto vol1 = std::make_shared( + Transform3::Identity(), + std::make_shared(30_mm, 40_mm, 100_mm)); + + auto vol2 = std::make_shared( + Transform3::Identity() * Translation3{Vector3::UnitY() * 80}, + std::make_shared(30_mm, 40_mm, 100_mm)); + + auto vol3 = std::make_shared( + Transform3::Identity() * Translation3{Vector3::UnitY() * 160}, + std::make_shared(30_mm, 40_mm, 100_mm)); + + auto rBounds = std::make_shared(30_mm, 40_mm); + auto plane1 = + Surface::makeShared(Transform3::Identity(), rBounds); + + auto plane2 = Surface::makeShared( + Transform3::Identity() * Translation3{Vector3::UnitY() * 80}, rBounds); + + auto plane3 = Surface::makeShared( + Transform3::Identity() * Translation3{Vector3::UnitY() * 160}, rBounds); + + auto trivial1 = std::make_unique(plane1, *vol1); + BOOST_REQUIRE(trivial1); + auto trivial2 = std::make_unique(plane2, *vol2); + BOOST_REQUIRE(trivial2); + auto trivial3 = std::make_unique(plane3, *vol3); + BOOST_REQUIRE(trivial3); + + auto grid1 = trivial1->makeGrid(AxisDirection::AxisY); + auto compGridTrivial = PortalLinkBase::merge( + std::move(grid1), std::make_unique(*trivial2), + AxisDirection::AxisY, *logger); + BOOST_REQUIRE(compGridTrivial); + BOOST_CHECK_EQUAL(dynamic_cast(*compGridTrivial) + .makeGrid(gctx, *logger), + nullptr); + + auto composite = + PortalLinkBase::merge(std::move(trivial1), std::move(trivial2), + AxisDirection::AxisY, *logger); + BOOST_REQUIRE(composite); + + auto grid12 = + dynamic_cast(*composite).makeGrid(gctx, *logger); + BOOST_REQUIRE(grid12); + + BOOST_CHECK_EQUAL( + grid12->resolveVolume(gctx, Vector2{0_mm, -40_mm}).value(), vol1.get()); + + BOOST_CHECK_EQUAL(grid12->resolveVolume(gctx, Vector2{0_mm, 40_mm}).value(), + vol2.get()); + + composite = PortalLinkBase::merge(std::move(composite), std::move(trivial3), + AxisDirection::AxisY, *logger); + BOOST_REQUIRE(composite); + + auto grid123 = + dynamic_cast(*composite).makeGrid(gctx, *logger); + BOOST_REQUIRE(grid123); + + BOOST_CHECK_EQUAL( + grid123->resolveVolume(gctx, Vector2{0_mm, -110_mm}).value(), + vol1.get()); + + BOOST_CHECK_EQUAL( + grid123->resolveVolume(gctx, Vector2{0_mm, -10_mm}).value(), + vol2.get()); + + BOOST_CHECK_EQUAL( + grid123->resolveVolume(gctx, Vector2{0_mm, 110_mm}).value(), + vol3.get()); + } } BOOST_AUTO_TEST_CASE(TrivialGridR) { @@ -2384,6 +3073,48 @@ BOOST_AUTO_TEST_CASE(TrivialGridPhi) { } } +BOOST_AUTO_TEST_CASE(TrivialGridXY) { + auto vol1 = std::make_shared( + Transform3::Identity(), + std::make_shared(30_mm, 30_mm, 30_mm)); + + auto vol2 = std::make_shared( + Transform3::Identity() * Translation3{Vector3::UnitX() * 60}, + std::make_shared(30_mm, 30_mm, 30_mm)); + + auto rBounds = std::make_shared(30_mm, 30_mm); + auto plane1 = + Surface::makeShared(Transform3::Identity(), rBounds); + + auto plane2 = Surface::makeShared( + Transform3::Identity() * Translation3{Vector3::UnitX() * 60}, rBounds); + + auto gridX = GridPortalLink::make(plane1, AxisDirection::AxisX, + Axis{AxisBound, -30_mm, 30_mm, 2}); + gridX->setVolume(vol1.get()); + + auto gridY = GridPortalLink::make(plane1, AxisDirection::AxisY, + Axis{AxisBound, -30_mm, 30_mm, 2}); + gridY->setVolume(vol1.get()); + + auto trivial = std::make_unique(plane2, *vol2); + BOOST_REQUIRE(trivial); + + BOOST_TEST_CONTEXT("Colinear") { + auto merged = PortalLinkBase::merge(copy(trivial), copy(gridX), + AxisDirection::AxisX, *logger); + BOOST_REQUIRE(merged); + BOOST_CHECK_NE(dynamic_cast(merged.get()), nullptr); + } + + BOOST_TEST_CONTEXT("Orthogonal") { + auto merged = PortalLinkBase::merge(copy(gridY), copy(trivial), + AxisDirection::AxisX, *logger); + BOOST_REQUIRE(merged); + BOOST_CHECK_NE(dynamic_cast(merged.get()), nullptr); + } +} + BOOST_AUTO_TEST_CASE(CompositeOther) { auto vol1 = std::make_shared( Transform3::Identity(),