From 0cba45dfcfc06906750fcfd073de90f868355c16 Mon Sep 17 00:00:00 2001 From: Dylan Pulver Date: Sun, 30 Aug 2026 15:50:44 +0300 Subject: [PATCH 1/6] Compute bounding boxes for general planes and tori in C++ The C++ Surface::bounding_box overrides cover the axis-aligned planes, the axis-aligned cylinders and the sphere, but not SurfacePlane or the three torus classes, while the Python API computes a box for all four. A cell bounded by a general plane that is axis-aligned to within roundoff therefore reports a different bounding box through openmc.lib.Cell.bounding_box than through openmc.Cell.bounding_box. Add the four missing overrides, mirroring the Python implementations. PlaneMixin.bounding_box tested for axis alignment using the normalized normal but chose per-axis intercepts using the raw coefficients. Those criteria disagree for off-axis coefficients between the tolerance and roughly 1e-6, where the resulting box excluded points that lie inside the half-space. Apply the alignment tolerance consistently so both implementations agree and neither produces such a box. --- include/openmc/constants.h | 5 + include/openmc/surface.h | 4 + openmc/surface.py | 26 ++- src/surface.cpp | 62 +++++++ tests/cpp_unit_tests/CMakeLists.txt | 1 + tests/cpp_unit_tests/test_surface.cpp | 240 ++++++++++++++++++++++++++ 6 files changed, 323 insertions(+), 15 deletions(-) create mode 100644 tests/cpp_unit_tests/test_surface.cpp diff --git a/include/openmc/constants.h b/include/openmc/constants.h index a1d94e5819e..1c13bc071d7 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -58,6 +58,11 @@ constexpr double FP_COINCIDENT {1e-12}; constexpr double TORUS_TOL {1e-10}; constexpr double RADIAL_MESH_TOL {1e-10}; +// Tolerance on the normalized normal of a general plane for treating that +// plane as axis-aligned when computing a bounding box. Matches the value of +// Surface._atol used by PlaneMixin.bounding_box in openmc/surface.py. +constexpr double PLANE_ALIGNMENT_TOL {1e-12}; + // Maximum number of random samples per history constexpr int MAX_SAMPLE {100000}; diff --git a/include/openmc/surface.h b/include/openmc/surface.h index 2d8580345a4..b459d847dc1 100644 --- a/include/openmc/surface.h +++ b/include/openmc/surface.h @@ -171,6 +171,7 @@ class SurfacePlane : public Surface { double distance(Position r, Direction u, bool coincident) const override; Direction normal(Position r) const override; void to_hdf5_inner(hid_t group_id) const override; + BoundingBox bounding_box(bool pos_side) const override; double A_, B_, C_, D_; }; @@ -337,6 +338,7 @@ class SurfaceXTorus : public Surface { double distance(Position r, Direction u, bool coincident) const override; Direction normal(Position r) const override; void to_hdf5_inner(hid_t group_id) const override; + BoundingBox bounding_box(bool pos_side) const override; double x0_, y0_, z0_, A_, B_, C_; }; @@ -354,6 +356,7 @@ class SurfaceYTorus : public Surface { double distance(Position r, Direction u, bool coincident) const override; Direction normal(Position r) const override; void to_hdf5_inner(hid_t group_id) const override; + BoundingBox bounding_box(bool pos_side) const override; double x0_, y0_, z0_, A_, B_, C_; }; @@ -371,6 +374,7 @@ class SurfaceZTorus : public Surface { double distance(Position r, Direction u, bool coincident) const override; Direction normal(Position r) const override; void to_hdf5_inner(hid_t group_id) const override; + BoundingBox bounding_box(bool pos_side) const override; double x0_, y0_, z0_, A_, B_, C_; }; diff --git a/openmc/surface.py b/openmc/surface.py index c2afeb613ac..3c8ef823c93 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -561,22 +561,18 @@ def bounding_box(self, side): nhat = self._get_normal() ll = np.array([-np.inf, -np.inf, -np.inf]) ur = np.array([np.inf, np.inf, np.inf]) - # If the plane is axis aligned, find the proper bounding box - if np.any(np.isclose(np.abs(nhat), 1., rtol=0., atol=self._atol)): + # A plane only bounds a half-space when its normal is parallel to a + # coordinate axis, in which case it bounds it along that axis alone. + aligned = np.isclose(np.abs(nhat), 1., rtol=0., atol=self._atol) + if aligned.any(): + axis = int(np.argmax(aligned)) sign = nhat.sum() - a, b, c, d = self._get_base_coeffs() - vals = [d/val if not np.isclose(val, 0., rtol=0., atol=self._atol) - else np.nan for val in (a, b, c)] - if side == '-': - if sign > 0: - ur = np.array([v if not np.isnan(v) else np.inf for v in vals]) - else: - ll = np.array([v if not np.isnan(v) else -np.inf for v in vals]) - elif side == '+': - if sign > 0: - ll = np.array([v if not np.isnan(v) else -np.inf for v in vals]) - else: - ur = np.array([v if not np.isnan(v) else np.inf for v in vals]) + coeffs = self._get_base_coeffs() + intercept = coeffs[3]/coeffs[axis] + if (side == '+') == (sign > 0): + ll[axis] = intercept + else: + ur[axis] = intercept return BoundingBox(ll, ur) diff --git a/src/surface.cpp b/src/surface.cpp index 81b756deae7..6b09a2c6dec 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -352,6 +352,41 @@ double SurfacePlane::evaluate(Position r) const return A_ * r.x + B_ * r.y + C_ * r.z - D_; } +BoundingBox SurfacePlane::bounding_box(bool pos_side) const +{ + // A general plane bounds a half-space in one direction only when its normal + // is parallel to a coordinate axis; otherwise both half-spaces are unbounded + // along every axis. This mirrors PlaneMixin.bounding_box on the Python side, + // so that a plane whose off-axis coefficients are rotation-matrix roundoff + // (e.g. B = 1 with A = C = 6.1e-17) yields the same box through both APIs. + const array coeffs {A_, B_, C_}; + const double norm = std::sqrt(A_ * A_ + B_ * B_ + C_ * C_); + if (norm == 0.0) + return {}; + + int axis = -1; + double sign = 0.0; + for (int i = 0; i < 3; ++i) { + const double n = coeffs[i] / norm; + sign += n; + if (axis == -1 && std::abs(std::abs(n) - 1.0) <= PLANE_ALIGNMENT_TOL) + axis = i; + } + if (axis == -1) + return {}; + + // The half-space is bounded below when the outward normal points along the + // positive axis direction and we are on the positive side, or vice versa. + BoundingBox bbox; + const double intercept = D_ / coeffs[axis]; + if (pos_side == (sign > 0.0)) { + bbox.min[axis] = intercept; + } else { + bbox.max[axis] = intercept; + } + return bbox; +} + double SurfacePlane::distance(Position r, Direction u, bool coincident) const { const double f = A_ * r.x + B_ * r.y + C_ * r.z - D_; @@ -1034,6 +1069,17 @@ double SurfaceXTorus::evaluate(Position r) const std::pow(std::sqrt(y * y + z * z) - A_, 2) / (C_ * C_) - 1.; } +BoundingBox SurfaceXTorus::bounding_box(bool pos_side) const +{ + // The torus interior is compact: it extends +/-B_ along the axis of + // revolution and +/-(A_ + C_) in the two perpendicular directions. Mirrors + // XTorus.bounding_box on the Python side. + if (pos_side) + return {}; + return {{x0_ - B_, y0_ - A_ - C_, z0_ - A_ - C_}, + {x0_ + B_, y0_ + A_ + C_, z0_ + A_ + C_}}; +} + double SurfaceXTorus::distance(Position r, Direction u, bool coincident) const { double x = r.x - x0_; @@ -1087,6 +1133,14 @@ double SurfaceYTorus::evaluate(Position r) const std::pow(std::sqrt(x * x + z * z) - A_, 2) / (C_ * C_) - 1.; } +BoundingBox SurfaceYTorus::bounding_box(bool pos_side) const +{ + if (pos_side) + return {}; + return {{x0_ - A_ - C_, y0_ - B_, z0_ - A_ - C_}, + {x0_ + A_ + C_, y0_ + B_, z0_ + A_ + C_}}; +} + double SurfaceYTorus::distance(Position r, Direction u, bool coincident) const { double x = r.x - x0_; @@ -1140,6 +1194,14 @@ double SurfaceZTorus::evaluate(Position r) const std::pow(std::sqrt(x * x + y * y) - A_, 2) / (C_ * C_) - 1.; } +BoundingBox SurfaceZTorus::bounding_box(bool pos_side) const +{ + if (pos_side) + return {}; + return {{x0_ - A_ - C_, y0_ - A_ - C_, z0_ - B_}, + {x0_ + A_ + C_, y0_ + A_ + C_, z0_ + B_}}; +} + double SurfaceZTorus::distance(Position r, Direction u, bool coincident) const { double x = r.x - x0_; diff --git a/tests/cpp_unit_tests/CMakeLists.txt b/tests/cpp_unit_tests/CMakeLists.txt index 991f219f528..b6892549bde 100644 --- a/tests/cpp_unit_tests/CMakeLists.txt +++ b/tests/cpp_unit_tests/CMakeLists.txt @@ -9,6 +9,7 @@ set(TEST_NAMES test_photon test_ray test_region + test_surface test_tensor test_geometry # Add additional unit test files here diff --git a/tests/cpp_unit_tests/test_surface.cpp b/tests/cpp_unit_tests/test_surface.cpp new file mode 100644 index 00000000000..d315806ce31 --- /dev/null +++ b/tests/cpp_unit_tests/test_surface.cpp @@ -0,0 +1,240 @@ +#include + +#include +#include + +#include + +#include "openmc/cell.h" +#include "openmc/constants.h" +#include "openmc/surface.h" + +using namespace openmc; + +namespace { + +template +std::unique_ptr make_surface( + pugi::xml_document& doc, int id, const char* type, const char* coeffs) +{ + pugi::xml_node n = doc.append_child("surface"); + n.append_attribute("id") = id; + n.append_attribute("type") = type; + n.append_attribute("coeffs") = coeffs; + return std::make_unique(n); +} + +// Register a surface under the given 1-based index so that Region can find it +template +void add_surface( + pugi::xml_document& doc, int id, const char* type, const char* coeffs) +{ + model::surfaces.push_back(make_surface(doc, id, type, coeffs)); + model::surface_map[id] = id - 1; +} + +// Builds the cell from the model attached to issue #2632 +class Issue2632Fixture { +public: + Issue2632Fixture() + { + // s19 (1), s48 (2), s58 (3), s62 (4), s64 (5), s65 (6), s68 (7), s101 (8) + add_surface(doc_, 1, "y-cylinder", "0.0 0.0 17.7"); + add_surface(doc_, 2, "plane", + "0.7071067811865476 6.123233995736766e-17 0.7071067811865476 11.45"); + add_surface(doc_, 3, "plane", + "0.7071067811865476 6.123233995736766e-17 0.7071067811865476 14.35"); + add_surface(doc_, 4, "plane", + "6.123233995736766e-17 1.0 6.123233995736766e-17 1.5999999999999999"); + add_surface( + doc_, 5, "plane", "6.123233995736766e-17 1.0 6.123233995736766e-17 -1.3"); + add_surface(doc_, 6, "plane", + "-0.7071067811865475 6.123233995736766e-17 0.7071067811865476 1.45"); + add_surface(doc_, 7, "plane", + "-0.7071067811865475 6.123233995736766e-17 0.7071067811865476 -1.45"); + add_surface(doc_, 8, "y-plane", "5.6"); + } + + ~Issue2632Fixture() + { + model::surfaces.clear(); + model::surface_map.clear(); + } + +private: + pugi::xml_document doc_; +}; + +} // anonymous namespace + +TEST_CASE("General plane bounding box") +{ + pugi::xml_document doc; + + SECTION("Exactly axis-aligned planes bound one axis only") + { + // +x normal: the positive half-space starts at x = D/A + auto px = make_surface(doc, 1, "plane", "1.0 0.0 0.0 5.0"); + BoundingBox pos = px->bounding_box(true); + CHECK(pos.min.x == Catch::Approx(5.0)); + CHECK(pos.min.y == -INFTY); + CHECK(pos.min.z == -INFTY); + CHECK(pos.max.x == INFTY); + + BoundingBox neg = px->bounding_box(false); + CHECK(neg.max.x == Catch::Approx(5.0)); + CHECK(neg.min.x == -INFTY); + CHECK(neg.max.y == INFTY); + + // -y normal: the sense of the bound flips with the sign of the normal + auto ny = make_surface(doc, 2, "plane", "0.0 -1.0 0.0 3.0"); + BoundingBox ny_pos = ny->bounding_box(true); + CHECK(ny_pos.max.y == Catch::Approx(-3.0)); + CHECK(ny_pos.min.y == -INFTY); + CHECK(ny_pos.min.x == -INFTY); + + BoundingBox ny_neg = ny->bounding_box(false); + CHECK(ny_neg.min.y == Catch::Approx(-3.0)); + CHECK(ny_neg.max.y == INFTY); + } + + SECTION("Coefficients need not be normalized") + { + // 4z - 10 = 0 is the same plane as z = 2.5 + auto pz = make_surface(doc, 3, "plane", "0.0 0.0 4.0 10.0"); + CHECK(pz->bounding_box(true).min.z == Catch::Approx(2.5)); + CHECK(pz->bounding_box(false).max.z == Catch::Approx(2.5)); + } + + SECTION("Rotation roundoff is still axis aligned (issue #2632)") + { + // The surface s64 from the reported model: a y-plane at -1.3 written with + // cos(pi/2) in the x and z slots. + auto p = make_surface( + doc, 4, "plane", "6.123233995736766e-17 1.0 6.123233995736766e-17 -1.3"); + BoundingBox pos = p->bounding_box(true); + CHECK(pos.min.y == Catch::Approx(-1.3)); + // The two off-axis directions must stay unbounded + CHECK(pos.min.x == -INFTY); + CHECK(pos.min.z == -INFTY); + CHECK(pos.max.x == INFTY); + CHECK(pos.max.z == INFTY); + } + + SECTION("Oblique planes bound nothing") + { + auto p = make_surface( + doc, 5, "plane", "0.7071067811865476 0.0 0.7071067811865476 11.45"); + for (bool side : {false, true}) { + BoundingBox bb = p->bounding_box(side); + CHECK(bb.min.x == -INFTY); + CHECK(bb.min.y == -INFTY); + CHECK(bb.min.z == -INFTY); + CHECK(bb.max.x == INFTY); + CHECK(bb.max.y == INFTY); + CHECK(bb.max.z == INFTY); + } + } + + SECTION("A plane tilted well beyond tolerance bounds nothing") + { + // 1e-5 is far outside PLANE_ALIGNMENT_TOL, so this must not be treated as + // an x-plane, and in particular must not acquire a bound on y. + auto p = make_surface(doc, 6, "plane", "1.0 1e-5 0.0 5.0"); + BoundingBox bb = p->bounding_box(true); + CHECK(bb.min.x == -INFTY); + CHECK(bb.min.y == -INFTY); + } + + SECTION("An almost-aligned plane bounds only the aligned axis") + { + // The off-axis coefficient is above PLANE_ALIGNMENT_TOL while the + // normalized normal is still within it. Only x may be bounded; a spurious + // y bound here would exclude points that lie inside the half-space. + auto p = make_surface(doc, 7, "plane", "1.0 1e-11 0.0 5.0"); + BoundingBox bb = p->bounding_box(true); + CHECK(bb.min.x == Catch::Approx(5.0)); + CHECK(bb.min.y == -INFTY); + CHECK(bb.min.z == -INFTY); + } + + SECTION("A degenerate plane bounds nothing") + { + auto p = make_surface(doc, 8, "plane", "0.0 0.0 0.0 1.0"); + BoundingBox bb = p->bounding_box(true); + CHECK(bb.min.x == -INFTY); + CHECK(bb.max.x == INFTY); + } +} + +TEST_CASE("Torus bounding box") +{ + pugi::xml_document doc; + + // x0 y0 z0 A B C, so the interior spans +/-B along the axis of revolution + // and +/-(A + C) in the perpendicular directions. + SECTION("x-torus") + { + auto t = make_surface( + doc, 1, "x-torus", "1.0 2.0 3.0 5.0 0.5 0.25"); + BoundingBox in = t->bounding_box(false); + CHECK(in.min.x == Catch::Approx(0.5)); + CHECK(in.max.x == Catch::Approx(1.5)); + CHECK(in.min.y == Catch::Approx(-3.25)); + CHECK(in.max.y == Catch::Approx(7.25)); + CHECK(in.min.z == Catch::Approx(-2.25)); + CHECK(in.max.z == Catch::Approx(8.25)); + + // The exterior of a torus is unbounded + BoundingBox out = t->bounding_box(true); + CHECK(out.min.x == -INFTY); + CHECK(out.max.z == INFTY); + } + + SECTION("y-torus") + { + auto t = make_surface( + doc, 2, "y-torus", "1.0 2.0 3.0 5.0 0.5 0.25"); + BoundingBox in = t->bounding_box(false); + CHECK(in.min.y == Catch::Approx(1.5)); + CHECK(in.max.y == Catch::Approx(2.5)); + CHECK(in.min.x == Catch::Approx(-4.25)); + CHECK(in.max.x == Catch::Approx(6.25)); + CHECK(in.min.z == Catch::Approx(-2.25)); + CHECK(in.max.z == Catch::Approx(8.25)); + + CHECK(t->bounding_box(true).max.y == INFTY); + } + + SECTION("z-torus") + { + auto t = make_surface( + doc, 3, "z-torus", "1.0 2.0 3.0 5.0 0.5 0.25"); + BoundingBox in = t->bounding_box(false); + CHECK(in.min.z == Catch::Approx(2.5)); + CHECK(in.max.z == Catch::Approx(3.5)); + CHECK(in.min.x == Catch::Approx(-4.25)); + CHECK(in.max.x == Catch::Approx(6.25)); + CHECK(in.min.y == Catch::Approx(-3.25)); + CHECK(in.max.y == Catch::Approx(7.25)); + + CHECK(t->bounding_box(true).min.z == -INFTY); + } +} + +TEST_CASE("Cell bounding box matches the Python API for issue #2632") +{ + Issue2632Fixture fixture; + + // +s64 -s101 -s19 +s48 (+s58 | +s62 | -s64 | +s65 | -s68) + Region region("5 -8 -1 2 (3 | 4 | -5 | 6 | -7)", 0); + BoundingBox bb = region.bounding_box(0); + + // The values reported by openmc.Cell.bounding_box in the issue + CHECK(bb.min.x == Catch::Approx(-17.7)); + CHECK(bb.min.y == Catch::Approx(-1.3)); + CHECK(bb.min.z == Catch::Approx(-17.7)); + CHECK(bb.max.x == Catch::Approx(17.7)); + CHECK(bb.max.y == Catch::Approx(5.6)); + CHECK(bb.max.z == Catch::Approx(17.7)); +} From 93138e5891e3a31e4139890a27f0d29073847ddb Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 9 Sep 2026 18:15:47 -0500 Subject: [PATCH 2/6] Simplify checks; add test for Python behavior change --- openmc/surface.py | 3 +-- src/surface.cpp | 9 ++++---- tests/unit_tests/test_surface.py | 38 ++++++++++++++++++++++++++++++++ 3 files changed, 43 insertions(+), 7 deletions(-) diff --git a/openmc/surface.py b/openmc/surface.py index 3c8ef823c93..b63ec806867 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -566,10 +566,9 @@ def bounding_box(self, side): aligned = np.isclose(np.abs(nhat), 1., rtol=0., atol=self._atol) if aligned.any(): axis = int(np.argmax(aligned)) - sign = nhat.sum() coeffs = self._get_base_coeffs() intercept = coeffs[3]/coeffs[axis] - if (side == '+') == (sign > 0): + if (side == '+') == (coeffs[axis] > 0): ll[axis] = intercept else: ur[axis] = intercept diff --git a/src/surface.cpp b/src/surface.cpp index 6b09a2c6dec..d6159c7995a 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -365,12 +365,11 @@ BoundingBox SurfacePlane::bounding_box(bool pos_side) const return {}; int axis = -1; - double sign = 0.0; for (int i = 0; i < 3; ++i) { - const double n = coeffs[i] / norm; - sign += n; - if (axis == -1 && std::abs(std::abs(n) - 1.0) <= PLANE_ALIGNMENT_TOL) + if (std::abs(std::abs(coeffs[i] / norm) - 1.0) <= PLANE_ALIGNMENT_TOL) { axis = i; + break; + } } if (axis == -1) return {}; @@ -379,7 +378,7 @@ BoundingBox SurfacePlane::bounding_box(bool pos_side) const // positive axis direction and we are on the positive side, or vice versa. BoundingBox bbox; const double intercept = D_ / coeffs[axis]; - if (pos_side == (sign > 0.0)) { + if (pos_side == (coeffs[axis] > 0.0)) { bbox.min[axis] = intercept; } else { bbox.max[axis] = intercept; diff --git a/tests/unit_tests/test_surface.py b/tests/unit_tests/test_surface.py index f5c0b8f8b62..1e7ff2af3c9 100644 --- a/tests/unit_tests/test_surface.py +++ b/tests/unit_tests/test_surface.py @@ -72,6 +72,44 @@ def test_plane_from_points(): assert s.d == 1.0 +def test_plane_bounding_box(): + # A plane written by a rotation has roundoff in the off-axis slots but is + # still a y-plane and should bound the half-space along y alone (#2632) + eps = math.cos(math.pi/2) + s = openmc.Plane(eps, 1., eps, -1.3) + ll, ur = (+s).bounding_box + assert ll == pytest.approx((-np.inf, -1.3, -np.inf)) + assert np.all(np.isinf(ur)) + ll, ur = (-s).bounding_box + assert ur == pytest.approx((np.inf, -1.3, np.inf)) + assert np.all(np.isinf(ll)) + + # The intercept comes from the raw coefficients, which need not be + # normalized: 4z - 10 = 0 is the plane z = 2.5 + s = openmc.Plane(0., 0., 4., 10.) + assert (+s).bounding_box[0] == pytest.approx((-np.inf, -np.inf, 2.5)) + assert (-s).bounding_box[1] == pytest.approx((np.inf, np.inf, 2.5)) + + # The sense of the bound flips with the sign of the normal + s = openmc.Plane(0., -1., 0., 3.) + assert (+s).bounding_box[1] == pytest.approx((np.inf, -3., np.inf)) + assert (-s).bounding_box[0] == pytest.approx((-np.inf, -3., -np.inf)) + + # An off-axis coefficient that is small enough to leave the plane axis + # aligned must not bound the off-axis directions, or the box would exclude + # points that lie inside the half-space it describes + s = openmc.Plane(1., 1.e-11, 0., 5.) + ll, ur = (+s).bounding_box + assert ll == pytest.approx((5., -np.inf, -np.inf)) + assert np.all(np.isinf(ur)) + p = (5., 1.e6, 0.) + assert p in +s + assert np.all(ll <= p) and np.all(p <= ur) + + # A plane tilted well beyond the alignment tolerance bounds nothing + assert_infinite_bb(openmc.Plane(1., 1.e-5, 0., 5.)) + + def test_xplane(): s = openmc.XPlane(3., boundary_type='reflective') assert s.x0 == 3. From aa69455bfc83d787bd684de4f4860807f42a0b34 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 9 Sep 2026 23:37:42 -0500 Subject: [PATCH 3/6] Stricter alignment check on plane normal --- openmc/surface.py | 7 ++++--- src/surface.cpp | 14 +++++++++++++- tests/cpp_unit_tests/test_surface.cpp | 20 ++++++++++++-------- tests/unit_tests/test_surface.py | 12 ++++-------- 4 files changed, 33 insertions(+), 20 deletions(-) diff --git a/openmc/surface.py b/openmc/surface.py index b63ec806867..441d9486ab2 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -563,9 +563,10 @@ def bounding_box(self, side): ur = np.array([np.inf, np.inf, np.inf]) # A plane only bounds a half-space when its normal is parallel to a # coordinate axis, in which case it bounds it along that axis alone. - aligned = np.isclose(np.abs(nhat), 1., rtol=0., atol=self._atol) - if aligned.any(): - axis = int(np.argmax(aligned)) + axis = int(np.argmax(np.abs(nhat))) + on_axis = np.isclose(abs(nhat[axis]), 1., rtol=0., atol=self._atol) + off_axis = np.delete(nhat, axis) + if on_axis and np.all(np.isclose(off_axis, 0., rtol=0., atol=self._atol)): coeffs = self._get_base_coeffs() intercept = coeffs[3]/coeffs[axis] if (side == '+') == (coeffs[axis] > 0): diff --git a/src/surface.cpp b/src/surface.cpp index d6159c7995a..d568d9f0a1a 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -364,9 +364,21 @@ BoundingBox SurfacePlane::bounding_box(bool pos_side) const if (norm == 0.0) return {}; + const array normal { + coeffs[0] / norm, coeffs[1] / norm, coeffs[2] / norm}; int axis = -1; for (int i = 0; i < 3; ++i) { - if (std::abs(std::abs(coeffs[i] / norm) - 1.0) <= PLANE_ALIGNMENT_TOL) { + if (std::abs(std::abs(normal[i]) - 1.0) > PLANE_ALIGNMENT_TOL) + continue; + + bool aligned = true; + for (int j = 0; j < 3; ++j) { + if (j != i && std::abs(normal[j]) > PLANE_ALIGNMENT_TOL) { + aligned = false; + break; + } + } + if (aligned) { axis = i; break; } diff --git a/tests/cpp_unit_tests/test_surface.cpp b/tests/cpp_unit_tests/test_surface.cpp index d315806ce31..5904ac2bb16 100644 --- a/tests/cpp_unit_tests/test_surface.cpp +++ b/tests/cpp_unit_tests/test_surface.cpp @@ -146,16 +146,20 @@ TEST_CASE("General plane bounding box") CHECK(bb.min.y == -INFTY); } - SECTION("An almost-aligned plane bounds only the aligned axis") + SECTION("An off-axis coefficient above tolerance bounds nothing") { - // The off-axis coefficient is above PLANE_ALIGNMENT_TOL while the - // normalized normal is still within it. Only x may be bounded; a spurious - // y bound here would exclude points that lie inside the half-space. + // Although the x component of the normalized normal is within + // PLANE_ALIGNMENT_TOL of one, the y component is above the tolerance. auto p = make_surface(doc, 7, "plane", "1.0 1e-11 0.0 5.0"); - BoundingBox bb = p->bounding_box(true); - CHECK(bb.min.x == Catch::Approx(5.0)); - CHECK(bb.min.y == -INFTY); - CHECK(bb.min.z == -INFTY); + for (bool side : {false, true}) { + BoundingBox bb = p->bounding_box(side); + CHECK(bb.min.x == -INFTY); + CHECK(bb.min.y == -INFTY); + CHECK(bb.min.z == -INFTY); + CHECK(bb.max.x == INFTY); + CHECK(bb.max.y == INFTY); + CHECK(bb.max.z == INFTY); + } } SECTION("A degenerate plane bounds nothing") diff --git a/tests/unit_tests/test_surface.py b/tests/unit_tests/test_surface.py index 1e7ff2af3c9..dc6843f8af1 100644 --- a/tests/unit_tests/test_surface.py +++ b/tests/unit_tests/test_surface.py @@ -95,16 +95,12 @@ def test_plane_bounding_box(): assert (+s).bounding_box[1] == pytest.approx((np.inf, -3., np.inf)) assert (-s).bounding_box[0] == pytest.approx((-np.inf, -3., -np.inf)) - # An off-axis coefficient that is small enough to leave the plane axis - # aligned must not bound the off-axis directions, or the box would exclude - # points that lie inside the half-space it describes + # An off-axis normalized coefficient above the alignment tolerance makes + # the half-space unbounded along every axis s = openmc.Plane(1., 1.e-11, 0., 5.) - ll, ur = (+s).bounding_box - assert ll == pytest.approx((5., -np.inf, -np.inf)) - assert np.all(np.isinf(ur)) - p = (5., 1.e6, 0.) + assert_infinite_bb(s) + p = (4., 2.e11, 0.) assert p in +s - assert np.all(ll <= p) and np.all(p <= ur) # A plane tilted well beyond the alignment tolerance bounds nothing assert_infinite_bb(openmc.Plane(1., 1.e-5, 0., 5.)) From 9fee3d67ab41a6321ee47927660efe0f0a6194f2 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 9 Sep 2026 23:44:28 -0500 Subject: [PATCH 4/6] Use normal() in SurfacePlane::bounding_box --- src/surface.cpp | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/src/surface.cpp b/src/surface.cpp index d568d9f0a1a..06f0f88b3c4 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -359,21 +359,19 @@ BoundingBox SurfacePlane::bounding_box(bool pos_side) const // along every axis. This mirrors PlaneMixin.bounding_box on the Python side, // so that a plane whose off-axis coefficients are rotation-matrix roundoff // (e.g. B = 1 with A = C = 6.1e-17) yields the same box through both APIs. - const array coeffs {A_, B_, C_}; - const double norm = std::sqrt(A_ * A_ + B_ * B_ + C_ * C_); + const Direction n = normal({}); + const double norm = n.norm(); if (norm == 0.0) return {}; - const array normal { - coeffs[0] / norm, coeffs[1] / norm, coeffs[2] / norm}; int axis = -1; for (int i = 0; i < 3; ++i) { - if (std::abs(std::abs(normal[i]) - 1.0) > PLANE_ALIGNMENT_TOL) + if (std::abs(std::abs(n[i] / norm) - 1.0) > PLANE_ALIGNMENT_TOL) continue; bool aligned = true; for (int j = 0; j < 3; ++j) { - if (j != i && std::abs(normal[j]) > PLANE_ALIGNMENT_TOL) { + if (j != i && std::abs(n[j] / norm) > PLANE_ALIGNMENT_TOL) { aligned = false; break; } @@ -389,8 +387,8 @@ BoundingBox SurfacePlane::bounding_box(bool pos_side) const // The half-space is bounded below when the outward normal points along the // positive axis direction and we are on the positive side, or vice versa. BoundingBox bbox; - const double intercept = D_ / coeffs[axis]; - if (pos_side == (coeffs[axis] > 0.0)) { + const double intercept = D_ / n[axis]; + if (pos_side == (n[axis] > 0.0)) { bbox.min[axis] = intercept; } else { bbox.max[axis] = intercept; From f28d07c3ef73124c62d25a62312adb4d0f4a4094 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 10 Sep 2026 00:13:45 -0500 Subject: [PATCH 5/6] Add test based on #2632 MWE --- tests/unit_tests/test_surface.py | 43 ++++++++++++++++++++++++-------- 1 file changed, 32 insertions(+), 11 deletions(-) diff --git a/tests/unit_tests/test_surface.py b/tests/unit_tests/test_surface.py index dc6843f8af1..887ac319bfc 100644 --- a/tests/unit_tests/test_surface.py +++ b/tests/unit_tests/test_surface.py @@ -4,6 +4,7 @@ import numpy as np import math import openmc +import openmc.lib import pytest @@ -106,6 +107,26 @@ def test_plane_bounding_box(): assert_infinite_bb(openmc.Plane(1., 1.e-5, 0., 5.)) +def test_plane_bounding_box_lib(mpi_intracomm): + """Check Python and C++ bounding consistency.""" + openmc.reset_auto_ids() + + cyl = openmc.YCylinder(r=17.7) + y_min = openmc.Plane(6.123233995736766e-17, 1., 6.123233995736766e-17, -1.3) + y_max = openmc.YPlane(5.6, boundary_type='vacuum') + cell = openmc.Cell(region=+y_min & -y_max & -cyl) + model = openmc.Model(geometry=openmc.Geometry([cell])) + model.settings.batches = 1 + model.settings.particles = 100 + + python_box = cell.bounding_box + with openmc.lib.TemporarySession(model, intracomm=mpi_intracomm): + lib_box = openmc.lib.cells[cell.id].bounding_box + + np.testing.assert_array_equal(python_box.lower_left, lib_box.lower_left) + np.testing.assert_array_equal(python_box.upper_right, lib_box.upper_right) + + def test_xplane(): s = openmc.XPlane(3., boundary_type='reflective') assert s.x0 == 3. @@ -229,12 +250,12 @@ def test_cylinder(): assert s.dy == -1 assert s.dz == 1 assert s.r == 2 - + # Check radius must be positive with pytest.raises(ValueError): openmc.Cylinder(x0=x0, y0=y0, z0=z0, dx=dx, dy=dy, dz=dz, r=0.0) with pytest.raises(ValueError): - openmc.Cylinder(x0=x0, y0=y0, z0=z0, dx=dx, dy=dy, dz=dz, r=-1.0) + openmc.Cylinder(x0=x0, y0=y0, z0=z0, dx=dx, dy=dy, dz=dz, r=-1.0) # Check bounding box assert_infinite_bb(s) @@ -284,7 +305,7 @@ def test_xcylinder(): assert s.y0 == y assert s.z0 == z assert s.r == r - + # Check radius must be positive with pytest.raises(ValueError): openmc.XCylinder(y0=y, z0=z, r=0.0) @@ -336,7 +357,7 @@ def test_ycylinder(): assert s.x0 == x assert s.z0 == z assert s.r == r - + # Check radius must be positive with pytest.raises(ValueError): openmc.YCylinder(x0=x, z0=z, r=0.0) @@ -379,7 +400,7 @@ def test_zcylinder(): assert s.x0 == x assert s.y0 == y assert s.r == r - + # Check radius must be positive with pytest.raises(ValueError): openmc.ZCylinder(x0=x, y0=y, r=0.0) @@ -423,7 +444,7 @@ def test_sphere(): assert s.y0 == y assert s.z0 == z assert s.r == r - + # Check radius must be positive with pytest.raises(ValueError): openmc.Sphere(x0=x, y0=y, z0=z, r=0.0) @@ -468,7 +489,7 @@ def cone_common(apex, r2, cls): assert s.y0 == y assert s.z0 == z assert s.r2 == r2 - + # Check radius must be positive with pytest.raises(ValueError): cls(x0=x, y0=y, z0=z, r2=0.0) @@ -512,12 +533,12 @@ def test_cone(): assert s.dy == -1 assert s.dz == 1 assert s.r2 == 4 - + # Check radius must be positive with pytest.raises(ValueError): openmc.Cone(x0=x0, y0=y0, z0=z0, dx=dx, dy=dy, dz=dz, r2=0.0) with pytest.raises(ValueError): - openmc.Cone(x0=x0, y0=y0, z0=z0, dx=dx, dy=dy, dz=dz, r2=-1.0) + openmc.Cone(x0=x0, y0=y0, z0=z0, dx=dx, dy=dy, dz=dz, r2=-1.0) # Check bounding box assert_infinite_bb(s) @@ -698,9 +719,9 @@ def torus_common(center, R, r1, r2, cls): assert s.a == R assert s.b == r1 assert s.c == r2 - + # Check radius must be positive - params = [(0.0, r1, r2), (R, 0.0, r2), (R, r1, 0.0), + params = [(0.0, r1, r2), (R, 0.0, r2), (R, r1, 0.0), (-1.0, r1, r2), (R, -1.0, r2), (R, r1, -1.0)] for a,b,c in params: with pytest.raises(ValueError): From 01ea0507cdb364c0724ca9fcf100c4de675ec5e0 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 10 Sep 2026 00:32:06 -0500 Subject: [PATCH 6/6] Simplify C++ tests --- tests/cpp_unit_tests/test_surface.cpp | 58 --------------------------- 1 file changed, 58 deletions(-) diff --git a/tests/cpp_unit_tests/test_surface.cpp b/tests/cpp_unit_tests/test_surface.cpp index 5904ac2bb16..1529cf21031 100644 --- a/tests/cpp_unit_tests/test_surface.cpp +++ b/tests/cpp_unit_tests/test_surface.cpp @@ -24,47 +24,6 @@ std::unique_ptr make_surface( return std::make_unique(n); } -// Register a surface under the given 1-based index so that Region can find it -template -void add_surface( - pugi::xml_document& doc, int id, const char* type, const char* coeffs) -{ - model::surfaces.push_back(make_surface(doc, id, type, coeffs)); - model::surface_map[id] = id - 1; -} - -// Builds the cell from the model attached to issue #2632 -class Issue2632Fixture { -public: - Issue2632Fixture() - { - // s19 (1), s48 (2), s58 (3), s62 (4), s64 (5), s65 (6), s68 (7), s101 (8) - add_surface(doc_, 1, "y-cylinder", "0.0 0.0 17.7"); - add_surface(doc_, 2, "plane", - "0.7071067811865476 6.123233995736766e-17 0.7071067811865476 11.45"); - add_surface(doc_, 3, "plane", - "0.7071067811865476 6.123233995736766e-17 0.7071067811865476 14.35"); - add_surface(doc_, 4, "plane", - "6.123233995736766e-17 1.0 6.123233995736766e-17 1.5999999999999999"); - add_surface( - doc_, 5, "plane", "6.123233995736766e-17 1.0 6.123233995736766e-17 -1.3"); - add_surface(doc_, 6, "plane", - "-0.7071067811865475 6.123233995736766e-17 0.7071067811865476 1.45"); - add_surface(doc_, 7, "plane", - "-0.7071067811865475 6.123233995736766e-17 0.7071067811865476 -1.45"); - add_surface(doc_, 8, "y-plane", "5.6"); - } - - ~Issue2632Fixture() - { - model::surfaces.clear(); - model::surface_map.clear(); - } - -private: - pugi::xml_document doc_; -}; - } // anonymous namespace TEST_CASE("General plane bounding box") @@ -225,20 +184,3 @@ TEST_CASE("Torus bounding box") CHECK(t->bounding_box(true).min.z == -INFTY); } } - -TEST_CASE("Cell bounding box matches the Python API for issue #2632") -{ - Issue2632Fixture fixture; - - // +s64 -s101 -s19 +s48 (+s58 | +s62 | -s64 | +s65 | -s68) - Region region("5 -8 -1 2 (3 | 4 | -5 | 6 | -7)", 0); - BoundingBox bb = region.bounding_box(0); - - // The values reported by openmc.Cell.bounding_box in the issue - CHECK(bb.min.x == Catch::Approx(-17.7)); - CHECK(bb.min.y == Catch::Approx(-1.3)); - CHECK(bb.min.z == Catch::Approx(-17.7)); - CHECK(bb.max.x == Catch::Approx(17.7)); - CHECK(bb.max.y == Catch::Approx(5.6)); - CHECK(bb.max.z == Catch::Approx(17.7)); -}