crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-27 21:55:41 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Ability to determine bounding boxes for Regions

Browse source
This commit is contained in:
Paul Romano 2015-12-01 07:22:20 -06:00
parent 4e31f4314c
commit 8506f32c4f
2 changed files with 319 additions and 0 deletions
Download patch file Download diff file Expand all files Collapse all files

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):
ll = np.array([-np.inf, -np.inf, -np.inf])
ur = np.array([np.inf, np.inf, np.inf])
for n in self.nodes:
ll_n, ur_n = n.bounding_box
ll[:] = np.maximum(ll, ll_n)
ur[:] = np.minimum(ur, ur_n)
return ll, ur
@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):
ll = np.array([np.inf, np.inf, np.inf])
ur = np.array([-np.inf, -np.inf, -np.inf])
for n in self.nodes:
ll_n, ur_n = n.bounding_box
ll[:] = np.minimum(ll, ll_n)
ur[:] = np.maximum(ur, ur_n)
return ll, ur
@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

278
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. 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
"""
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. 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
"""
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. 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
"""
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,37 @@ 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. 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
"""
if side == '-':
return (np.array([-np.inf, self.y0 - self.r, self.z0 - self.r]),
np.array([np.inf, self.y0 + self.r, self.y0 + 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 +696,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 +723,37 @@ 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. 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
"""
if side == '-':
return (np.array([self.x0 - self.r, -np.inf, self.z0 - self.r]),
np.array([self.x0 + self.r, np.inf, self.y0 + 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 +794,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 +821,37 @@ 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. 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
"""
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 +896,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 +949,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. 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
"""
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 +1027,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 +1252,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 +1399,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 +1429,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)