implemented BoundingBox class

This commit is contained in:
Jonathan Shimwell 2023-04-14 13:11:01 +01:00
parent 93de45d044
commit 6eb26b80c6
7 changed files with 188 additions and 40 deletions

View file

@ -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 *

44
openmc/bounding_box.py Normal file
View file

@ -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]))

View file

@ -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):

View file

@ -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):

View file

@ -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):

View file

@ -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):

View file

@ -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])))