crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 14:15:42 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Merge pull request #516 from paulromano/bounding-box

Browse source 
Ability to determine bounding boxes for Regions
This commit is contained in:
Will Boyd 2015-12-03 09:12:01 -05:00
commit 43d13428aa
5 changed files with 329 additions and 3 deletions
Download patch file Download diff file Expand all files Collapse all files

4
examples/python/boxes/build-xml.py
View file

@ -1,3 +1,5 @@
import numpy as np
import openmc
###############################################################################
@ -115,7 +117,7 @@ settings_file = openmc.SettingsFile()
settings_file.batches = batches
settings_file.inactive = inactive
settings_file.particles = particles
settings_file.set_source_space('point', [0., 0., 0.])
settings_file.set_source_space('box', np.concatenate(outer_cube.bounding_box))
settings_file.export_to_xml()
###############################################################################

4
examples/python/reflective/build-xml.py
View file

@ -1,3 +1,5 @@
import numpy as np
import openmc
###############################################################################
@ -82,5 +84,5 @@ settings_file = openmc.SettingsFile()
settings_file.batches = batches
settings_file.inactive = inactive
settings_file.particles = particles
settings_file.set_source_space('box', [-1, -1, -1, 1, 1, 1])
settings_file.set_source_space('box', np.concatenate(cell.region.bounding_box))
settings_file.export_to_xml()

2
examples/xml/boxes/settings.xml
View file

@ -10,7 +10,7 @@
<!-- Starting source -->
<source>
<space type="point" parameters="0. 0. 0." />
<space type="box" parameters="-10. -10. -10. 10. 10. 10." />
</source>
</settings>

41
openmc/region.py
View file

@ -1,6 +1,8 @@
from abc import ABCMeta, abstractmethod
from collections import Iterable
import numpy as np
from openmc.checkvalue import check_type
@ -218,6 +220,8 @@ class Intersection(Region):
----------
nodes : tuple of Region
Regions to take the intersection of
bounding_box : tuple of numpy.array
Lower-left and upper-right coordinates of an axis-aligned bounding box
"""
@ -231,6 +235,16 @@ class Intersection(Region):
def nodes(self):
return self._nodes
@property
def bounding_box(self):
lower_left = np.array([-np.inf, -np.inf, -np.inf])
upper_right = np.array([np.inf, np.inf, np.inf])
for n in self.nodes:
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
@nodes.setter
def nodes(self, nodes):
check_type('nodes', nodes, Iterable, Region)
@ -257,6 +271,8 @@ class Union(Region):
----------
nodes : tuple of Region
Regions to take the union of
bounding_box : tuple of numpy.array
Lower-left and upper-right coordinates of an axis-aligned bounding box
"""
@ -270,6 +286,16 @@ class Union(Region):
def nodes(self):
return self._nodes
@property
def bounding_box(self):
lower_left = np.array([np.inf, np.inf, np.inf])
upper_right = np.array([-np.inf, -np.inf, -np.inf])
for n in self.nodes:
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
@nodes.setter
def nodes(self, nodes):
check_type('nodes', nodes, Iterable, Region)
@ -300,6 +326,8 @@ class Complement(Region):
----------
node : Region
Regions to take the complement of
bounding_box : tuple of numpy.array
Lower-left and upper-right coordinates of an axis-aligned bounding box
"""
@ -317,3 +345,16 @@ class Complement(Region):
def node(self, node):
check_type('node', node, Region)
self._node = node
@property
def bounding_box(self):
# Use De Morgan's laws to distribute the complement operator so that it
# only applies to surface half-spaces, thus allowing us to calculate the
# bounding box in the usual recursive manner.
if isinstance(self.node, Union):
temp_region = Intersection(*[~n for n in self.node.nodes])
elif isinstance(self.node, Intersection):
temp_region = Union(*[~n for n in self.node.nodes])
else:
temp_region = ~n
return temp_region.bounding_box

281
openmc/surface.py
View file

@ -3,6 +3,8 @@ from numbers import Real, Integral
from xml.etree import ElementTree as ET
import sys
import numpy as np
from openmc.checkvalue import check_type, check_value, check_greater_than
from openmc.region import Region
@ -136,6 +138,33 @@ class Surface(object):
check_value('boundary type', boundary_type, _BC_TYPES)
self._boundary_type = boundary_type
def bounding_box(self, side):
"""Determine an axis-aligned bounding box.
An axis-aligned bounding box for surface half-spaces is represented by
its lower-left and upper-right coordinates. If the half-space is
unbounded in a particular direction, numpy.inf is used to represent
infinity.
Parameters
----------
side : {'+', '-'}
Indicates the negative or positive half-space
Returns
-------
numpy.array
Lower-left coordinates of the axis-aligned bounding box for the
desired half-space
numpy.array
Upper-right coordinates of the axis-aligned bounding box for the
desired half-space
"""
return (np.array([-np.inf, -np.inf, -np.inf]),
np.array([np.inf, np.inf, np.inf]))
def create_xml_subelement(self):
element = ET.Element("surface")
element.set("id", str(self._id))
@ -194,6 +223,10 @@ class Plane(Surface):
self._type = 'plane'
self._coeff_keys = ['A', 'B', 'C', 'D']
self._coeffs['A'] = 1.
self._coeffs['B'] = 0.
self._coeffs['C'] = 0.
self._coeffs['D'] = 0.
if A is not None:
self.a = A
@ -276,6 +309,7 @@ class XPlane(Plane):
self._type = 'x-plane'
self._coeff_keys = ['x0']
self._coeffs['x0'] = 0.
if x0 is not None:
self.x0 = x0
@ -289,6 +323,37 @@ class XPlane(Plane):
check_type('x0 coefficient', x0, Real)
self._coeffs['x0'] = x0
def bounding_box(self, side):
"""Determine an axis-aligned bounding box.
An axis-aligned bounding box for surface half-spaces is represented by
its lower-left and upper-right coordinates. For the x-plane surface, the
half-spaces are unbounded in their y- and z- directions. To represent
infinity, numpy.inf is used.
Parameters
----------
side : {'+', '-'}
Indicates the negative or positive half-space
Returns
-------
numpy.array
Lower-left coordinates of the axis-aligned bounding box for the
desired half-space
numpy.array
Upper-right coordinates of the axis-aligned bounding box for the
desired half-space
"""
if side == '-':
return (np.array([-np.inf, -np.inf, -np.inf]),
np.array([self.x0, np.inf, np.inf]))
elif side == '+':
return (np.array([self.x0, -np.inf, -np.inf]),
np.array([np.inf, np.inf, np.inf]))
class YPlane(Plane):
"""A plane perpendicular to the y axis, i.e. a surface of the form :math:`y -
@ -322,6 +387,7 @@ class YPlane(Plane):
self._type = 'y-plane'
self._coeff_keys = ['y0']
self._coeffs['y0'] = 0.
if y0 is not None:
self.y0 = y0
@ -335,6 +401,37 @@ class YPlane(Plane):
check_type('y0 coefficient', y0, Real)
self._coeffs['y0'] = y0
def bounding_box(self, side):
"""Determine an axis-aligned bounding box.
An axis-aligned bounding box for surface half-spaces is represented by
its lower-left and upper-right coordinates. For the y-plane surface, the
half-spaces are unbounded in their x- and z- directions. To represent
infinity, numpy.inf is used.
Parameters
----------
side : {'+', '-'}
Indicates the negative or positive half-space
Returns
-------
numpy.array
Lower-left coordinates of the axis-aligned bounding box for the
desired half-space
numpy.array
Upper-right coordinates of the axis-aligned bounding box for the
desired half-space
"""
if side == '-':
return (np.array([-np.inf, -np.inf, -np.inf]),
np.array([np.inf, self.y0, np.inf]))
elif side == '+':
return (np.array([-np.inf, self.y0, -np.inf]),
np.array([np.inf, np.inf, np.inf]))
class ZPlane(Plane):
"""A plane perpendicular to the z axis, i.e. a surface of the form :math:`z -
@ -368,6 +465,7 @@ class ZPlane(Plane):
self._type = 'z-plane'
self._coeff_keys = ['z0']
self._coeffs['z0'] = 0.
if z0 is not None:
self.z0 = z0
@ -381,6 +479,37 @@ class ZPlane(Plane):
check_type('z0 coefficient', z0, Real)
self._coeffs['z0'] = z0
def bounding_box(self, side):
"""Determine an axis-aligned bounding box.
An axis-aligned bounding box for surface half-spaces is represented by
its lower-left and upper-right coordinates. For the z-plane surface, the
half-spaces are unbounded in their x- and y- directions. To represent
infinity, numpy.inf is used.
Parameters
----------
side : {'+', '-'}
Indicates the negative or positive half-space
Returns
-------
numpy.array
Lower-left coordinates of the axis-aligned bounding box for the
desired half-space
numpy.array
Upper-right coordinates of the axis-aligned bounding box for the
desired half-space
"""
if side == '-':
return (np.array([-np.inf, -np.inf, -np.inf]),
np.array([np.inf, np.inf, self.z0]))
elif side == '+':
return (np.array([-np.inf, -np.inf, self.z0]),
np.array([np.inf, np.inf, np.inf]))
class Cylinder(Surface):
"""A cylinder whose length is parallel to the x-, y-, or z-axis.
@ -415,6 +544,7 @@ class Cylinder(Surface):
super(Cylinder, self).__init__(surface_id, boundary_type, name=name)
self._coeff_keys = ['R']
self._coeffs['R'] = 1.
if R is not None:
self.r = R
@ -468,6 +598,8 @@ class XCylinder(Cylinder):
self._type = 'x-cylinder'
self._coeff_keys = ['y0', 'z0', 'R']
self._coeffs['y0'] = 0.
self._coeffs['z0'] = 0.
if y0 is not None:
self.y0 = y0
@ -493,6 +625,38 @@ class XCylinder(Cylinder):
check_type('z0 coefficient', z0, Real)
self._coeffs['z0'] = z0
def bounding_box(self, side):
"""Determine an axis-aligned bounding box.
An axis-aligned bounding box for surface half-spaces is represented by
its lower-left and upper-right coordinates. For the x-cylinder surface,
the negative half-space is unbounded in the x- direction and the
positive half-space is unbounded in all directions. To represent
infinity, numpy.inf is used.
Parameters
----------
side : {'+', '-'}
Indicates the negative or positive half-space
Returns
-------
numpy.array
Lower-left coordinates of the axis-aligned bounding box for the
desired half-space
numpy.array
Upper-right coordinates of the axis-aligned bounding box for the
desired half-space
"""
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]))
elif side == '+':
return (np.array([-np.inf, -np.inf, -np.inf]),
np.array([np.inf, np.inf, np.inf]))
class YCylinder(Cylinder):
"""An infinite cylinder whose length is parallel to the y-axis. This is a
@ -533,6 +697,8 @@ class YCylinder(Cylinder):
self._type = 'y-cylinder'
self._coeff_keys = ['x0', 'z0', 'R']
self._coeffs['x0'] = 0.
self._coeffs['z0'] = 0.
if x0 is not None:
self.x0 = x0
@ -558,6 +724,38 @@ class YCylinder(Cylinder):
check_type('z0 coefficient', z0, Real)
self._coeffs['z0'] = z0
def bounding_box(self, side):
"""Determine an axis-aligned bounding box.
An axis-aligned bounding box for surface half-spaces is represented by
its lower-left and upper-right coordinates. For the y-cylinder surface,
the negative half-space is unbounded in the y- direction and the
positive half-space is unbounded in all directions. To represent
infinity, numpy.inf is used.
Parameters
----------
side : {'+', '-'}
Indicates the negative or positive half-space
Returns
-------
numpy.array
Lower-left coordinates of the axis-aligned bounding box for the
desired half-space
numpy.array
Upper-right coordinates of the axis-aligned bounding box for the
desired half-space
"""
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]))
elif side == '+':
return (np.array([-np.inf, -np.inf, -np.inf]),
np.array([np.inf, np.inf, np.inf]))
class ZCylinder(Cylinder):
"""An infinite cylinder whose length is parallel to the z-axis. This is a
@ -598,6 +796,8 @@ class ZCylinder(Cylinder):
self._type = 'z-cylinder'
self._coeff_keys = ['x0', 'y0', 'R']
self._coeffs['x0'] = 0.
self._coeffs['y0'] = 0.
if x0 is not None:
self.x0 = x0
@ -623,6 +823,38 @@ class ZCylinder(Cylinder):
check_type('y0 coefficient', y0, Real)
self._coeffs['y0'] = y0
def bounding_box(self, side):
"""Determine an axis-aligned bounding box.
An axis-aligned bounding box for surface half-spaces is represented by
its lower-left and upper-right coordinates. For the z-cylinder surface,
the negative half-space is unbounded in the z- direction and the
positive half-space is unbounded in all directions. To represent
infinity, numpy.inf is used.
Parameters
----------
side : {'+', '-'}
Indicates the negative or positive half-space
Returns
-------
numpy.array
Lower-left coordinates of the axis-aligned bounding box for the
desired half-space
numpy.array
Upper-right coordinates of the axis-aligned bounding box for the
desired half-space
"""
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]))
elif side == '+':
return (np.array([-np.inf, -np.inf, -np.inf]),
np.array([np.inf, np.inf, np.inf]))
class Sphere(Surface):
"""A sphere of the form :math:`(x - x_0)^2 + (y - y_0)^2 + (z - z_0)^2 = R^2`.
@ -667,6 +899,10 @@ class Sphere(Surface):
self._type = 'sphere'
self._coeff_keys = ['x0', 'y0', 'z0', 'R']
self._coeffs['x0'] = 0.
self._coeffs['y0'] = 0.
self._coeffs['z0'] = 0.
self._coeffs['R'] = 1.
if x0 is not None:
self.x0 = x0
@ -716,6 +952,39 @@ class Sphere(Surface):
check_type('R coefficient', R, Real)
self._coeffs['R'] = R
def bounding_box(self, side):
"""Determine an axis-aligned bounding box.
An axis-aligned bounding box for surface half-spaces is represented by
its lower-left and upper-right coordinates. The positive half-space of a
sphere is unbounded in all directions. To represent infinity, numpy.inf
is used.
Parameters
----------
side : {'+', '-'}
Indicates the negative or positive half-space
Returns
-------
numpy.array
Lower-left coordinates of the axis-aligned bounding box for the
desired half-space
numpy.array
Upper-right coordinates of the axis-aligned bounding box for the
desired half-space
"""
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]))
elif side == '+':
return (np.array([-np.inf, -np.inf, -np.inf]),
np.array([np.inf, np.inf, np.inf]))
class Cone(Surface):
"""A conical surface parallel to the x-, y-, or z-axis.
@ -761,6 +1030,10 @@ class Cone(Surface):
super(Cone, self).__init__(surface_id, boundary_type, name=name)
self._coeff_keys = ['x0', 'y0', 'z0', 'R2']
self._coeffs['x0'] = 0.
self._coeffs['y0'] = 0.
self._coeffs['z0'] = 0.
self._coeffs['R2'] = 1.
if x0 is not None:
self.x0 = x0
@ -982,6 +1255,8 @@ class Quadric(Surface):
self._type = 'quadric'
self._coeff_keys = ['a', 'b', 'c', 'd', 'e', 'f', 'g', 'h', 'j', 'k']
for key in self._coeff_keys:
self._coeffs[key] = 0.
if a is not None:
self.a = a
@ -1127,6 +1402,8 @@ class Halfspace(Region):
Surface which divides Euclidean space.
side : {'+', '-'}
Indicates whether the positive or negative half-space is used.
bounding_box : tuple of numpy.array
Lower-left and upper-right coordinates of an axis-aligned bounding box
"""
@ -1155,6 +1432,10 @@ class Halfspace(Region):
check_value('side', side, ('+', '-'))
self._side = side
@property
def bounding_box(self):
return self.surface.bounding_box(self.side)
def __str__(self):
return '-' + str(self.surface.id) if self.side == '-' \
else str(self.surface.id)