From 0cba45dfcfc06906750fcfd073de90f868355c16 Mon Sep 17 00:00:00 2001 From: Dylan Pulver Date: Sun, 30 Aug 2026 15:50:44 +0300 Subject: [PATCH] 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)); +}