Compute bounding boxes for general planes and tori in C++

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};
 

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_;
 };

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)):
-            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])
+        # 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.
+        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):
+                ll[axis] = intercept
+            else:
+                ur[axis] = intercept
 
         return BoundingBox(ll, ur)
 

src/surface.cpp

@@ -352,6 +352,50 @@ 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 Direction n = normal({});
+  const double norm = n.norm();
+  if (norm == 0.0)
+    return {};
+
+  int axis = -1;
+  for (int i = 0; i < 3; ++i) {
+    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(n[j] / norm) > PLANE_ALIGNMENT_TOL) {
+        aligned = false;
+        break;
+      }
+    }
+    if (aligned) {
+      axis = i;
+      break;
+    }
+  }
+  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_ / n[axis];
+  if (pos_side == (n[axis] > 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 +1078,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 +1142,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 +1203,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_;

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

tests/cpp_unit_tests/test_surface.cpp

@@ -0,0 +1,186 @@
+#include <memory>
+
+#include <catch2/catch_approx.hpp>
+#include <catch2/catch_test_macros.hpp>
+
+#include <pugixml.hpp>
+
+#include "openmc/cell.h"
+#include "openmc/constants.h"
+#include "openmc/surface.h"
+
+using namespace openmc;
+
+namespace {
+
+template<typename T>
+std::unique_ptr<T> 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<T>(n);
+}
+
+} // 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<SurfacePlane>(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<SurfacePlane>(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<SurfacePlane>(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<SurfacePlane>(
+      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<SurfacePlane>(
+      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<SurfacePlane>(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 off-axis coefficient above tolerance bounds nothing")
+  {
+    // 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<SurfacePlane>(doc, 7, "plane", "1.0 1e-11 0.0 5.0");
+    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")
+  {
+    auto p = make_surface<SurfacePlane>(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<SurfaceXTorus>(
+      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<SurfaceYTorus>(
+      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<SurfaceZTorus>(
+      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);
+  }
+}