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)):
+        # 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)
 

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<double, 3> 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_;

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