From 6eb26b80c6a8fe4613b2183df61d6c3c49856366 Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Fri, 14 Apr 2023 13:11:01 +0100 Subject: [PATCH] implemented BoundingBox class --- openmc/__init__.py | 1 + openmc/bounding_box.py | 44 +++++++++++ openmc/cell.py | 5 +- openmc/lib/cell.py | 4 +- openmc/region.py | 5 +- openmc/surface.py | 102 +++++++++++++++++--------- tests/unit_tests/test_bounding_box.py | 67 +++++++++++++++++ 7 files changed, 188 insertions(+), 40 deletions(-) create mode 100644 openmc/bounding_box.py create mode 100644 tests/unit_tests/test_bounding_box.py diff --git a/openmc/__init__.py b/openmc/__init__.py index ac9e02e168..5bb150fb07 100644 --- a/openmc/__init__.py +++ b/openmc/__init__.py @@ -1,4 +1,5 @@ from openmc.arithmetic import * +from openmc.bounding_box import * from openmc.cell import * from openmc.checkvalue import * from openmc.mesh import * diff --git a/openmc/bounding_box.py b/openmc/bounding_box.py new file mode 100644 index 0000000000..0818937975 --- /dev/null +++ b/openmc/bounding_box.py @@ -0,0 +1,44 @@ +import numpy as np +from .checkvalue import check_type, check_length + + +class BoundingBox(tuple): + """Axis-aligned bounding box. + + Parameters + ---------- + corners : 2-tuple of numpy.array + Lower-left and upper-right coordinates of an axis-aligned bounding box + of the domain. + + Attributes + ---------- + volume: float + The volume of the bounding box in cm3 + center: numpy.array + x, y, z coordinates of the center of the bounding box in cm. + """ + + def __init__(self, corners): + + check_type("corners", corners, tuple) + check_type("corners", corners[0], np.ndarray) + check_type("corners", corners[1], np.ndarray) + check_length("corners", corners, 2, 2) + check_length("corners", corners[0], 3, 3) + check_length("corners", corners[1], 3, 3) + self.corners = corners + + @property + def center(self): + """The center x, y, z coordinates of the bounding box""" + # x_center = (self.corners[0][0] - self.corners[1][0]) / 2 + # y_center = (self.corners[0][1] - self.corners[1][1]) / 2 + # z_center = (self.corners[0][2] - self.corners[1][2]) / 2 + # return np.array([x_center, y_center, z_center]) + return (self.corners[0] + self.corners[1]) / 2 + + @property + def volume(self): + """The volume of the bounding box""" + return np.abs(np.prod(self.corners[1] - self.corners[0])) diff --git a/openmc/cell.py b/openmc/cell.py index a8e1178f41..db8ffd48be 100644 --- a/openmc/cell.py +++ b/openmc/cell.py @@ -13,6 +13,7 @@ from ._xml import get_text from .mixin import IDManagerMixin from .region import Region, Complement from .surface import Halfspace +from .bounding_box import BoundingBox class Cell(IDManagerMixin): @@ -249,8 +250,8 @@ class Cell(IDManagerMixin): if self.region is not None: return self.region.bounding_box else: - return (np.array([-np.inf, -np.inf, -np.inf]), - np.array([np.inf, np.inf, np.inf])) + return BoundingBox((np.array([-np.inf, -np.inf, -np.inf]), + np.array([np.inf, np.inf, np.inf]))) @property def num_instances(self): diff --git a/openmc/lib/cell.py b/openmc/lib/cell.py index b53a6be7a7..c88d6da519 100644 --- a/openmc/lib/cell.py +++ b/openmc/lib/cell.py @@ -11,6 +11,8 @@ from . import _dll from .core import _FortranObjectWithID from .error import _error_handler from .material import Material +from ..bounding_box import BoundingBox + __all__ = ['Cell', 'cells'] @@ -289,7 +291,7 @@ class Cell(_FortranObjectWithID): llc[llc == -inf] = -np.inf urc[urc == -inf] = -np.inf - return llc, urc + return BoundingBox((llc, urc)) class _CellMapping(Mapping): diff --git a/openmc/region.py b/openmc/region.py index 5d0e680204..d8935a0705 100644 --- a/openmc/region.py +++ b/openmc/region.py @@ -6,6 +6,7 @@ from copy import deepcopy import numpy as np from .checkvalue import check_type +from .bounding_box import BoundingBox class Region(ABC): @@ -418,7 +419,7 @@ class Intersection(Region, MutableSequence): lower_left_n, upper_right_n = n.bounding_box lower_left[:] = np.maximum(lower_left, lower_left_n) upper_right[:] = np.minimum(upper_right, upper_right_n) - return lower_left, upper_right + return BoundingBox((lower_left, upper_right)) class Union(Region, MutableSequence): @@ -506,7 +507,7 @@ class Union(Region, MutableSequence): lower_left_n, upper_right_n = n.bounding_box lower_left[:] = np.minimum(lower_left, lower_left_n) upper_right[:] = np.maximum(upper_right, upper_right_n) - return lower_left, upper_right + return BoundingBox((lower_left, upper_right)) class Complement(Region): diff --git a/openmc/surface.py b/openmc/surface.py index b8ea5479d5..98ced13c69 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -12,6 +12,7 @@ import numpy as np from .checkvalue import check_type, check_value, check_length from .mixin import IDManagerMixin, IDWarning from .region import Region, Intersection, Union +from .bounding_box import BoundingBox _BOUNDARY_TYPES = ['transmission', 'vacuum', 'reflective', 'periodic', 'white'] @@ -233,8 +234,10 @@ class Surface(IDManagerMixin, ABC): desired half-space """ - return (np.array([-np.inf, -np.inf, -np.inf]), - np.array([np.inf, np.inf, np.inf])) + return BoundingBox(( + np.array([-np.inf, -np.inf, -np.inf]), + np.array([np.inf, np.inf, np.inf]) + )) def clone(self, memo=None): """Create a copy of this surface with a new unique ID. @@ -539,7 +542,7 @@ class PlaneMixin: else: ur = np.array([v if not np.isnan(v) else np.inf for v in vals]) - return (ll, ur) + return BoundingBox((ll, ur)) def evaluate(self, point): """Evaluate the surface equation at a given point. @@ -1196,8 +1199,10 @@ class Cylinder(QuadricMixin, Surface): return (np.array(ll), np.array(ur)) elif side == '+': - return (np.array([-np.inf, -np.inf, -np.inf]), - np.array([np.inf, np.inf, np.inf])) + return BoundingBox(( + np.array([-np.inf, -np.inf, -np.inf]), + np.array([np.inf, np.inf, np.inf]) + )) def _get_base_coeffs(self): # Get x, y, z coordinates of two points @@ -1359,11 +1364,15 @@ class XCylinder(QuadricMixin, Surface): def bounding_box(self, side): if side == '-': - return (np.array([-np.inf, self.y0 - self.r, self.z0 - self.r]), - np.array([np.inf, self.y0 + self.r, self.z0 + self.r])) + return BoundingBox(( + np.array([-np.inf, self.y0 - self.r, self.z0 - self.r]), + np.array([np.inf, self.y0 + self.r, self.z0 + self.r]) + )) elif side == '+': - return (np.array([-np.inf, -np.inf, -np.inf]), - np.array([np.inf, np.inf, np.inf])) + return BoundingBox(( + np.array([-np.inf, -np.inf, -np.inf]), + np.array([np.inf, np.inf, np.inf]) + )) def evaluate(self, point): y = point[1] - self.y0 @@ -1450,11 +1459,15 @@ class YCylinder(QuadricMixin, Surface): def bounding_box(self, side): if side == '-': - return (np.array([self.x0 - self.r, -np.inf, self.z0 - self.r]), - np.array([self.x0 + self.r, np.inf, self.z0 + self.r])) + return BoundingBox(( + np.array([self.x0 - self.r, -np.inf, self.z0 - self.r]), + np.array([self.x0 + self.r, np.inf, self.z0 + self.r]) + )) elif side == '+': - return (np.array([-np.inf, -np.inf, -np.inf]), - np.array([np.inf, np.inf, np.inf])) + return BoundingBox(( + np.array([-np.inf, -np.inf, -np.inf]), + np.array([np.inf, np.inf, np.inf]) + )) def evaluate(self, point): x = point[0] - self.x0 @@ -1541,11 +1554,15 @@ class ZCylinder(QuadricMixin, Surface): def bounding_box(self, side): if side == '-': - return (np.array([self.x0 - self.r, self.y0 - self.r, -np.inf]), - np.array([self.x0 + self.r, self.y0 + self.r, np.inf])) + return BoundingBox(( + np.array([self.x0 - self.r, self.y0 - self.r, -np.inf]), + np.array([self.x0 + self.r, self.y0 + self.r, np.inf]) + )) elif side == '+': - return (np.array([-np.inf, -np.inf, -np.inf]), - np.array([np.inf, np.inf, np.inf])) + return BoundingBox(( + np.array([-np.inf, -np.inf, -np.inf]), + np.array([np.inf, np.inf, np.inf]) + )) def evaluate(self, point): x = point[0] - self.x0 @@ -1631,13 +1648,15 @@ class Sphere(QuadricMixin, Surface): def bounding_box(self, side): if side == '-': - return (np.array([self.x0 - self.r, self.y0 - self.r, - self.z0 - self.r]), - np.array([self.x0 + self.r, self.y0 + self.r, - self.z0 + self.r])) + return BoundingBox(( + np.array([self.x0 - self.r, self.y0 - self.r, self.z0 - self.r]), + np.array([self.x0 + self.r, self.y0 + self.r, self.z0 + self.r]) + )) elif side == '+': - return (np.array([-np.inf, -np.inf, -np.inf]), - np.array([np.inf, np.inf, np.inf])) + return BoundingBox(( + np.array([-np.inf, -np.inf, -np.inf]), + np.array([np.inf, np.inf, np.inf]) + )) def evaluate(self, point): x = point[0] - self.x0 @@ -2253,11 +2272,15 @@ class XTorus(TorusMixin, Surface): x0, y0, z0 = self.x0, self.y0, self.z0 a, b, c = self.a, self.b, self.c if side == '-': - return (np.array([x0 - b, y0 - a - c, z0 - a - c]), - np.array([x0 + b, y0 + a + c, z0 + a + c])) + return BoundingBox(( + np.array([x0 - b, y0 - a - c, z0 - a - c]), + np.array([x0 + b, y0 + a + c, z0 + a + c]) + )) elif side == '+': - return (np.array([-np.inf, -np.inf, -np.inf]), - np.array([np.inf, np.inf, np.inf])) + return BoundingBox(( + np.array([-np.inf, -np.inf, -np.inf]), + np.array([np.inf, np.inf, np.inf]) + )) class YTorus(TorusMixin, Surface): r"""A torus of the form :math:`(y - y_0)^2/B^2 + (\sqrt{(x - x_0)^2 + (z - @@ -2324,11 +2347,16 @@ class YTorus(TorusMixin, Surface): x0, y0, z0 = self.x0, self.y0, self.z0 a, b, c = self.a, self.b, self.c if side == '-': - return (np.array([x0 - a - c, y0 - b, z0 - a - c]), - np.array([x0 + a + c, y0 + b, z0 + a + c])) + return BoundingBox(( + np.array([x0 - a - c, y0 - b, z0 - a - c]), + np.array([x0 + a + c, y0 + b, z0 + a + c]) + )) elif side == '+': - return (np.array([-np.inf, -np.inf, -np.inf]), - np.array([np.inf, np.inf, np.inf])) + return BoundingBox(( + np.array([-np.inf, -np.inf, -np.inf]), + np.array([np.inf, np.inf, np.inf]) + )) + class ZTorus(TorusMixin, Surface): r"""A torus of the form :math:`(z - z_0)^2/B^2 + (\sqrt{(x - x_0)^2 + (y - @@ -2395,11 +2423,15 @@ class ZTorus(TorusMixin, Surface): x0, y0, z0 = self.x0, self.y0, self.z0 a, b, c = self.a, self.b, self.c if side == '-': - return (np.array([x0 - a - c, y0 - a - c, z0 - b]), - np.array([x0 + a + c, y0 + a + c, z0 + b])) + return BoundingBox(( + np.array([x0 - a - c, y0 - a - c, z0 - b]), + np.array([x0 + a + c, y0 + a + c, z0 + b]) + )) elif side == '+': - return (np.array([-np.inf, -np.inf, -np.inf]), - np.array([np.inf, np.inf, np.inf])) + return BoundingBox(( + np.array([-np.inf, -np.inf, -np.inf]), + np.array([np.inf, np.inf, np.inf]) + )) class Halfspace(Region): diff --git a/tests/unit_tests/test_bounding_box.py b/tests/unit_tests/test_bounding_box.py new file mode 100644 index 0000000000..72c1baf486 --- /dev/null +++ b/tests/unit_tests/test_bounding_box.py @@ -0,0 +1,67 @@ +import numpy as np +import openmc +import pytest + + +test_bb_1 = openmc.BoundingBox( + (np.array([-10.0, -20.0, -30.0]), np.array([1.0, 2.0, 3.0])) +) +test_bb_2 = openmc.BoundingBox( + (np.array([1.0, 2.0, 3.0]), np.array([11.0, 22.0, 33.0])) +) +test_bb_3 = openmc.BoundingBox( + (np.array([-10.0, -20.0, -30.0]), np.array([-1.0, -2.0, -3.0])) +) + + +@pytest.mark.parametrize( + "bb, expected", + [ + (test_bb_1, 7986), # 11 * 22 * 33 + (test_bb_2, 6000), # 10 * 20 * 30 + (test_bb_3, 4374), # 9 * 18 * 27 + ], +) +def test_bounding_box_volume(bb, expected): + assert bb.volume == expected + + +@pytest.mark.parametrize( + "bb, expected", + [ + (test_bb_1, np.array([-4.5, -9.0, -13.5])), + (test_bb_2, np.array([6.0, 12.0, 18.0])), + (test_bb_3, np.array([-5.5, -11.0, -16.5])), + ], +) +def test_bounding_box_center(bb, expected): + assert np.array_equiv(expected, bb.center) + + +def test_bounding_box_input_checking(): + # checks that a list is not accepted instead of np.array + with pytest.raises(TypeError): + openmc.BoundingBox(([-10, -20, -30], [1, 2, 3])) + # checks that a list is not accepted instead of tuple + with pytest.raises(TypeError): + openmc.BoundingBox([np.array([-10, -20, -30]), np.array([1, 2, 3])]) + # checks that a tuple with one entry is not accepted + with pytest.raises(TypeError): + openmc.BoundingBox((np.array([-10, -20, -30]))) + # checks that a tuple with three entry is not accepted + with pytest.raises(ValueError): + openmc.BoundingBox( + (np.array([-1, -2, -3]), np.array([-1, -2, -3]), np.array([-1, -2, -3])) + ) + # checks that a numpy array with two entries is not accepted + with pytest.raises(ValueError): + openmc.BoundingBox((np.array([-10, -30]), np.array([1, 2, 3]))) + # checks that a numpy array with two entries is not accepted + with pytest.raises(ValueError): + openmc.BoundingBox((np.array([-10, -20, -30]), np.array([1, 3]))) + # checks that a numpy array with four entries is not accepted + with pytest.raises(ValueError): + openmc.BoundingBox((np.array([-10, -20, -3, -4]), np.array([1, 2, 3]))) + # checks that a numpy array with four entries is not accepted + with pytest.raises(ValueError): + openmc.BoundingBox((np.array([-10, -20, -4]), np.array([1, 2, 3, 4])))