From 2870da61e63b7dfa4d9c425d94f2736a7cdf8a54 Mon Sep 17 00:00:00 2001 From: Ethan Peterson Date: Fri, 21 Feb 2020 20:49:06 -0500 Subject: [PATCH 1/8] Refactoring improvements --- openmc/surface.py | 333 ++++++++++++++++++++++++++++++++++------------ 1 file changed, 246 insertions(+), 87 deletions(-) diff --git a/openmc/surface.py b/openmc/surface.py index 92c828d65f..0c1513a527 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -24,6 +24,8 @@ _WARNING_KWARGS = """\ will not accept positional parameters for superclass arguments.\ """ +_SURFACE_CLASSES = {} + class Surface(IDManagerMixin, metaclass=ABCMeta): """An implicit surface with an associated boundary condition. @@ -77,6 +79,10 @@ class Surface(IDManagerMixin, metaclass=ABCMeta): # Value - coefficient value self._coefficients = {} + def __init_subclass__(cls, **kwargs): + super().__init_subclass__(**kwargs) + _SURFACE_CLASSES[cls._type] = cls + def __neg__(self): return Halfspace(self, '-') @@ -361,6 +367,36 @@ class PlaneMixin(metaclass=ABCMeta): self._periodic_surface = periodic_surface periodic_surface._periodic_surface = self + def _get_base_coeffs(self): + return (self.a, self.b, self.c, self.d) + + def _get_normal(self): + a, b, c = self._get_base_coeffs()[0:3] + return np.array((a, b, c)) / np.sqrt(a*a + b*b + c*c) + + def bounding_box(self, side): + # Compute the bounding box based on the normal vector to the plane + nhat = self._get_normal() + lb = np.array([-np.inf, -np.inf, -np.inf]) + ub = 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)): + sign = np.dot(np.array([1., 1., 1.]), nhat) + a, b, c, d = self._get_base_coeffs() + vals = tuple(d/i if round(i) != 0 else np.nan for i in (a, b, c)) + if side == '-': + if sign > 0: + ub = np.array([v if ~np.isnan(v) else np.inf for v in vals]) + else: + lb = np.array([v if ~np.isnan(v) else -np.inf for v in vals]) + elif side == '+': + if sign > 0: + lb = np.array([v if ~np.isnan(v) else -np.inf for v in vals]) + else: + ub = np.array([v if ~np.isnan(v) else np.inf for v in vals]) + + return (lb, ub) + def evaluate(self, point): """Evaluate the surface equation at a given point. @@ -509,6 +545,12 @@ class Plane(PlaneMixin, Surface): FutureWarning) setattr(self, k.lower(), val) + @classmethod + def __subclasshook__(cls, c): + if cls is Plane and c in (XPlane, YPlane, ZPlane): + return True + return NotImplemented + @property def a(self): return self.coefficients['a'] @@ -545,9 +587,6 @@ class Plane(PlaneMixin, Surface): check_type('D coefficient', d, Real) self._coefficients['d'] = d - def _get_base_coeffs(self): - return (self.a, self.b, self.c, self.d) - @classmethod def from_points(cls, p1, p2, p3, **kwargs): """Return a plane given three points that pass through it. @@ -638,29 +677,31 @@ class XPlane(PlaneMixin, Surface): def x0(self): return self.coefficients['x0'] + @property + def a(self): + return 1. + + @property + def b(self): + return 0. + + @property + def c(self): + return 0. + + @property + def d(self): + return self.x0 + @x0.setter def x0(self, x0): check_type('x0 coefficient', x0, Real) self._coefficients['x0'] = x0 - def _get_base_coeffs(self): - return (1., 0., 0., self.x0) - - def bounding_box(self, side): - 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])) - def evaluate(self, point): return point[0] - self.x0 -Plane.register(XPlane) - - class YPlane(PlaneMixin, Surface): """A plane perpendicular to the y axis of the form :math:`y - y_0 = 0` @@ -718,29 +759,31 @@ class YPlane(PlaneMixin, Surface): def y0(self): return self.coefficients['y0'] + @property + def a(self): + return 0. + + @property + def b(self): + return 1. + + @property + def c(self): + return 0. + + @property + def d(self): + return self.y0 + @y0.setter def y0(self, y0): check_type('y0 coefficient', y0, Real) self._coefficients['y0'] = y0 - def _get_base_coeffs(self): - return (0., 1., 0., self.y0) - - def bounding_box(self, side): - 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])) - def evaluate(self, point): return point[1] - self.y0 -Plane.register(YPlane) - - class ZPlane(PlaneMixin, Surface): """A plane perpendicular to the z axis of the form :math:`z - z_0 = 0` @@ -798,29 +841,31 @@ class ZPlane(PlaneMixin, Surface): def z0(self): return self.coefficients['z0'] + @property + def a(self): + return 0. + + @property + def b(self): + return 0. + + @property + def c(self): + return 1. + + @property + def d(self): + return self.z0 + @z0.setter def z0(self, z0): check_type('z0 coefficient', z0, Real) self._coefficients['z0'] = z0 - def _get_base_coeffs(self): - return (0., 0., 1., self.z0) - - def bounding_box(self, side): - 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])) - def evaluate(self, point): return point[2] - self.z0 -Plane.register(ZPlane) - - class QuadricMixin(metaclass=ABCMeta): """A Mixin class implementing common functionality for quadric surfaces""" @@ -884,7 +929,7 @@ class QuadricMixin(metaclass=ABCMeta): """ x = np.asarray(point) A, b, c = self.get_Abc() - return np.matmul(x.T, np.matmul(A, x)) + np.matmul(b.T, x) + c + return x.T @ A @ x + b.T @ x + c def translate(self, vector, inplace=False): """Translate surface in given direction @@ -907,17 +952,20 @@ class QuadricMixin(metaclass=ABCMeta): surf = self if inplace else self.clone() - if set(('x0', 'y0', 'z0')).intersection(set(surf._coeff_keys)): + if any((isinstance(self, cls) for cls in (Cylinder, Sphere, Cone))): for vi, xi in zip(vector, ('x0', 'y0', 'z0')): - val = getattr(surf, xi, None) - if val is not None: + val = getattr(surf, xi) + try: setattr(surf, xi, val + vi) - else: + except AttributeError: + # That attribute is read only i.e x0 for XCylinder + pass + + elif isinstance(self, Quadric): A, bvec, cnst = self.get_Abc() - g, h, j = bvec - 2*np.matmul(vector.T, A) - k = cnst + np.matmul(vector.T, np.matmul(A, vector)) \ - - np.matmul(bvec.T, vector) + g, h, j = bvec - 2*vector.T @ A + k = cnst + vector.T @ A @ vector - bvec.T @ vector for key, val in zip(('g', 'h', 'j', 'k'), (g, h, j, k)): setattr(surf, key, val) @@ -993,16 +1041,17 @@ class Cylinder(QuadricMixin, Surface): _coeff_keys = ('x0', 'y0', 'z0', 'r', 'dx', 'dy', 'dz') def __init__(self, x0=0., y0=0., z0=0., r=1., dx=0., dy=0., dz=1., **kwargs): - raise NotImplementedError('There is no C++ implementation for general ' - 'Cylinders yet, please use ' - 'openmc.model.funcs.cylinder_from_points to ' - 'return a Quadric instance instead for now') - super().__init__(**kwargs) for key, val in zip(self._coeff_keys, (x0, y0, z0, r, dx, dy, dz)): setattr(self, key, val) + @classmethod + def __subclasshook__(cls, c): + if cls is Cylinder and c in (XCylinder, YCylinder, ZCylinder): + return True + return NotImplemented + @property def x0(self): return self.coefficients['x0'] @@ -1118,10 +1167,6 @@ class Cylinder(QuadricMixin, Surface): radius r. """ - raise NotImplementedError('There is no C++ implementation for general ' - 'Cylinders yet, please use ' - 'openmc.model.funcs.cylinder_from_points to ' - 'return a Quadric instance instead for now') # Convert to numpy arrays p1 = np.asarray(p1) p2 = np.asarray(p2) @@ -1130,6 +1175,30 @@ class Cylinder(QuadricMixin, Surface): return cls(x0=x0, y0=y0, z0=z0, r=r, dx=dx, dy=dy, dz=dz, **kwargs) + def to_xml_element(self): + """Return XML representation of the surface + + Returns + ------- + element : xml.etree.ElementTree.Element + XML element containing source data + + """ + # This method overrides Surface.to_xml_element to generate a Quadric + # since the C++ layer doesn't support Cylinders right now + element = ET.Element("surface") + element.set("id", str(self._id)) + + if len(self._name) > 0: + element.set("name", str(self._name)) + + element.set("type", 'quadric') + if self.boundary_type != 'transmission': + element.set("boundary", self.boundary_type) + element.set("coeffs", ' '.join([str(c) for c in self._get_base_coeffs()])) + + return element + class XCylinder(QuadricMixin, Surface): """An infinite cylinder whose length is parallel to the x-axis of the form @@ -1202,6 +1271,22 @@ class XCylinder(QuadricMixin, Surface): def r(self): return self.coefficients['r'] + @property + def x0(self): + return 0. + + @property + def dx(self): + return 1. + + @property + def dy(self): + return 0. + + @property + def dz(self): + return 0. + @y0.setter def y0(self, y0): check_type('y0 coefficient', y0, Real) @@ -1240,9 +1325,6 @@ class XCylinder(QuadricMixin, Surface): return y*y + z*z - self.r**2 -Cylinder.register(XCylinder) - - class YCylinder(QuadricMixin, Surface): """An infinite cylinder whose length is parallel to the y-axis of the form :math:`(x - x_0)^2 + (z - z_0)^2 = r^2`. @@ -1314,6 +1396,22 @@ class YCylinder(QuadricMixin, Surface): def r(self): return self.coefficients['r'] + @property + def y0(self): + return 0. + + @property + def dx(self): + return 0. + + @property + def dy(self): + return 1. + + @property + def dz(self): + return 0. + @x0.setter def x0(self, x0): check_type('x0 coefficient', x0, Real) @@ -1352,9 +1450,6 @@ class YCylinder(QuadricMixin, Surface): return x*x + z*z - self.r**2 -Cylinder.register(YCylinder) - - class ZCylinder(QuadricMixin, Surface): """An infinite cylinder whose length is parallel to the z-axis of the form :math:`(x - x_0)^2 + (y - y_0)^2 = r^2`. @@ -1426,6 +1521,22 @@ class ZCylinder(QuadricMixin, Surface): def r(self): return self.coefficients['r'] + @property + def z0(self): + return 0. + + @property + def dx(self): + return 0. + + @property + def dy(self): + return 0. + + @property + def dz(self): + return 1. + @x0.setter def x0(self, x0): check_type('x0 coefficient', x0, Real) @@ -1464,9 +1575,6 @@ class ZCylinder(QuadricMixin, Surface): return x*x + y*y - self.r**2 -Cylinder.register(ZCylinder) - - class Sphere(QuadricMixin, Surface): """A sphere of the form :math:`(x - x_0)^2 + (y - y_0)^2 + (z - z_0)^2 = r^2`. @@ -1656,9 +1764,6 @@ class Cone(QuadricMixin, Surface): _coeff_keys = ('x0', 'y0', 'z0', 'r2', 'dx', 'dy', 'dz') def __init__(self, x0=0., y0=0., z0=0., r2=1., dx=0., dy=0., dz=1., **kwargs): - raise NotImplementedError('There is no C++ implementation for general ' - 'Cones yet, this functionality should be ' - 'added soon.') R2 = kwargs.pop('R2', None) if R2 is not None: warn(_WARNING_UPPER.format(type(self).__name__, 'r2', 'R2'), @@ -1669,6 +1774,12 @@ class Cone(QuadricMixin, Surface): for key, val in zip(self._coeff_keys, (x0, y0, z0, r2, dx, dy, dz)): setattr(self, key, val) + @classmethod + def __subclasshook__(cls, c): + if cls is Cone and c in (XCone, YCone, ZCone): + return True + return NotImplemented + @property def x0(self): return self.coefficients['x0'] @@ -1761,11 +1872,35 @@ class Cone(QuadricMixin, Surface): f = 2*dx*dz g = -2*(dx*dx*x0 + dx*dy*y0 + dx*dz*z0 - cos2) h = -2*(dy*dy*y0 + dx*dy*x0 + dy*dz*z0 - cos2) - j = -2*(dz*dz*y0 + dx*dz*x0 + dy*dz*y0 - cos2) + j = -2*(dz*dz*z0 + dx*dz*x0 + dy*dz*y0 - cos2) k = (dx*x0 + dy*y0 + dz*z0)**2 - cos2*(x0*x0 + y0*y0 + z0*z0) return (a, b, c, d, e, f, g, h, j, k) + def to_xml_element(self): + """Return XML representation of the surface + + Returns + ------- + element : xml.etree.ElementTree.Element + XML element containing source data + + """ + # This method overrides Surface.to_xml_element to generate a Quadric + # since the C++ layer doesn't support Cylinders right now + element = ET.Element("surface") + element.set("id", str(self._id)) + + if len(self._name) > 0: + element.set("name", str(self._name)) + + element.set("type", 'quadric') + if self.boundary_type != 'transmission': + element.set("boundary", self.boundary_type) + element.set("coeffs", ' '.join([str(c) for c in self._get_base_coeffs()])) + + return element + class XCone(QuadricMixin, Surface): """A cone parallel to the x-axis of the form :math:`(y - y_0)^2 + (z - z_0)^2 = @@ -1845,6 +1980,18 @@ class XCone(QuadricMixin, Surface): def r2(self): return self.coefficients['r2'] + @property + def dx(self): + return 1. + + @property + def dy(self): + return 0. + + @property + def dz(self): + return 0. + @x0.setter def x0(self, x0): check_type('x0 coefficient', x0, Real) @@ -1883,9 +2030,6 @@ class XCone(QuadricMixin, Surface): return y*y + z*z - self.r2*x*x -Cone.register(XCone) - - class YCone(QuadricMixin, Surface): """A cone parallel to the y-axis of the form :math:`(x - x_0)^2 + (z - z_0)^2 = r^2 (y - y_0)^2`. @@ -1964,6 +2108,18 @@ class YCone(QuadricMixin, Surface): def r2(self): return self.coefficients['r2'] + @property + def dx(self): + return 0. + + @property + def dy(self): + return 1. + + @property + def dz(self): + return 0. + @x0.setter def x0(self, x0): check_type('x0 coefficient', x0, Real) @@ -2002,9 +2158,6 @@ class YCone(QuadricMixin, Surface): return x*x + z*z - self.r2*y*y -Cone.register(YCone) - - class ZCone(QuadricMixin, Surface): """A cone parallel to the x-axis of the form :math:`(x - x_0)^2 + (y - y_0)^2 = r^2 (z - z_0)^2`. @@ -2083,6 +2236,18 @@ class ZCone(QuadricMixin, Surface): def r2(self): return self.coefficients['r2'] + @property + def dx(self): + return 0. + + @property + def dy(self): + return 0. + + @property + def dz(self): + return 1. + @x0.setter def x0(self, x0): check_type('x0 coefficient', x0, Real) @@ -2121,9 +2286,6 @@ class ZCone(QuadricMixin, Surface): return x*x + y*y - self.r2*z*z -Cone.register(ZCone) - - class Quadric(QuadricMixin, Surface): """A surface of the form :math:`Ax^2 + By^2 + Cz^2 + Dxy + Eyz + Fxz + Gx + Hy + Jz + K = 0`. @@ -2451,6 +2613,3 @@ class Halfspace(Region): # Return translated surface return type(self)(memo[key], self.side) - - -_SURFACE_CLASSES = {cls._type: cls for cls in Surface.__subclasses__()} From 089b6d04ef5609c5f0c205f773472d0f3afbc722 Mon Sep 17 00:00:00 2001 From: Ethan Peterson Date: Sat, 22 Feb 2020 08:03:35 -0500 Subject: [PATCH 2/8] removed __init_subclass__ since it's not in python 3.5 --- openmc/surface.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/openmc/surface.py b/openmc/surface.py index 0c1513a527..62ea83456c 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -79,10 +79,6 @@ class Surface(IDManagerMixin, metaclass=ABCMeta): # Value - coefficient value self._coefficients = {} - def __init_subclass__(cls, **kwargs): - super().__init_subclass__(**kwargs) - _SURFACE_CLASSES[cls._type] = cls - def __neg__(self): return Halfspace(self, '-') @@ -2613,3 +2609,5 @@ class Halfspace(Region): # Return translated surface return type(self)(memo[key], self.side) + +_SURFACE_CLASSES = {cls._type: cls for cls in Surface.__subclasses__()} From ac5bb4047d4dc4de9fe71877c05b3d26c123b386 Mon Sep 17 00:00:00 2001 From: Ethan Peterson Date: Sun, 23 Feb 2020 10:42:27 -0500 Subject: [PATCH 3/8] removed extraneous definition --- openmc/surface.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/openmc/surface.py b/openmc/surface.py index 62ea83456c..a4c8b8c88d 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -24,8 +24,6 @@ _WARNING_KWARGS = """\ will not accept positional parameters for superclass arguments.\ """ -_SURFACE_CLASSES = {} - class Surface(IDManagerMixin, metaclass=ABCMeta): """An implicit surface with an associated boundary condition. From 378856d434a1778d79bc3d3d34077914b3cbf88b Mon Sep 17 00:00:00 2001 From: Ethan Peterson Date: Mon, 24 Feb 2020 11:27:29 -0500 Subject: [PATCH 4/8] Apply suggestions from code review Co-Authored-By: Paul Romano --- openmc/surface.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/openmc/surface.py b/openmc/surface.py index a4c8b8c88d..a1696b08c4 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -366,7 +366,7 @@ class PlaneMixin(metaclass=ABCMeta): def _get_normal(self): a, b, c = self._get_base_coeffs()[0:3] - return np.array((a, b, c)) / np.sqrt(a*a + b*b + c*c) + return np.array((a, b, c)) / math.sqrt(a*a + b*b + c*c) def bounding_box(self, side): # Compute the bounding box based on the normal vector to the plane @@ -375,12 +375,12 @@ class PlaneMixin(metaclass=ABCMeta): ub = 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)): - sign = np.dot(np.array([1., 1., 1.]), nhat) + sign = nhat.sum() a, b, c, d = self._get_base_coeffs() - vals = tuple(d/i if round(i) != 0 else np.nan for i in (a, b, c)) + vals = [d/val if round(val) != 0 else np.nan for val in (a, b, c)] if side == '-': if sign > 0: - ub = np.array([v if ~np.isnan(v) else np.inf for v in vals]) + ub = np.array([v if not np.isnan(v) else np.inf for v in vals]) else: lb = np.array([v if ~np.isnan(v) else -np.inf for v in vals]) elif side == '+': @@ -946,7 +946,7 @@ class QuadricMixin(metaclass=ABCMeta): surf = self if inplace else self.clone() - if any((isinstance(self, cls) for cls in (Cylinder, Sphere, Cone))): + if hasattr(self, 'x0'): for vi, xi in zip(vector, ('x0', 'y0', 'z0')): val = getattr(surf, xi) try: @@ -955,7 +955,7 @@ class QuadricMixin(metaclass=ABCMeta): # That attribute is read only i.e x0 for XCylinder pass - elif isinstance(self, Quadric): + else: A, bvec, cnst = self.get_Abc() g, h, j = bvec - 2*vector.T @ A From 3d3a27f841b03c765196a84c3a5dfd0e9b76cc70 Mon Sep 17 00:00:00 2001 From: Ethan Peterson Date: Mon, 24 Feb 2020 11:55:42 -0500 Subject: [PATCH 5/8] applied suggestions from code review --- docs/source/usersguide/install.rst | 2 +- openmc/model/funcs.py | 47 +----------------- openmc/surface.py | 79 +++++++++++++++++------------- setup.py | 3 +- tests/unit_tests/test_surface.py | 27 +++++----- 5 files changed, 63 insertions(+), 95 deletions(-) diff --git a/docs/source/usersguide/install.rst b/docs/source/usersguide/install.rst index 66440a3c37..ec10fea4fa 100644 --- a/docs/source/usersguide/install.rst +++ b/docs/source/usersguide/install.rst @@ -404,7 +404,7 @@ to install the Python package in :ref:`"editable" mode `. Prerequisites ------------- -The Python API works with Python 3.4+. In addition to Python itself, the API +The Python API works with Python 3.5+. In addition to Python itself, the API relies on a number of third-party packages. All prerequisites can be installed using Conda_ (recommended), pip_, or through the package manager in most Linux distributions. To run simulations in parallel using MPI, it is recommended to diff --git a/openmc/model/funcs.py b/openmc/model/funcs.py index 3ee5052402..9999c383ff 100644 --- a/openmc/model/funcs.py +++ b/openmc/model/funcs.py @@ -373,52 +373,7 @@ def get_hexagonal_prism(*args, **kwargs): return hexagonal_prism(*args, **kwargs) -def cylinder_from_points(p1, p2, r, **kwargs): - """Return cylinder defined by two points passing through its center. - - Parameters - ---------- - p1, p2 : 3-tuples - Coordinates of two points that pass through the center of the cylinder - r : float - Radius of the cylinder - kwargs : dict - Keyword arguments passed to the :class:`openmc.Quadric` constructor - - Returns - ------- - openmc.Quadric - Quadric surface representing the cylinder. - - """ - # Get x, y, z coordinates of two points - x1, y1, z1 = p1 - x2, y2, z2 = p2 - - # Define intermediate terms - dx = x2 - x1 - dy = y2 - y1 - dz = z2 - z1 - cx = y1*z2 - y2*z1 - cy = x2*z1 - x1*z2 - cz = x1*y2 - x2*y1 - - # Given p=(x,y,z), p1=(x1, y1, z1), p2=(x2, y2, z2), the equation for the - # cylinder can be derived as r = |(p - p1) ⨯ (p - p2)| / |p2 - p1|. - # Expanding out all terms and grouping according to what Quadric expects - # gives the following coefficients. - kwargs['a'] = dy*dy + dz*dz - kwargs['b'] = dx*dx + dz*dz - kwargs['c'] = dx*dx + dy*dy - kwargs['d'] = -2*dx*dy - kwargs['e'] = -2*dy*dz - kwargs['f'] = -2*dx*dz - kwargs['g'] = 2*(cy*dz - cz*dy) - kwargs['h'] = 2*(cz*dx - cx*dz) - kwargs['j'] = 2*(cx*dy - cy*dx) - kwargs['k'] = cx*cx + cy*cy + cz*cz - (dx*dx + dy*dy + dz*dz)*r*r - - return Quadric(**kwargs) +cylinder_from_points = Cylinder.from_points def subdivide(surfaces): diff --git a/openmc/surface.py b/openmc/surface.py index a1696b08c4..4d0434a4b0 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -3,13 +3,14 @@ from collections import OrderedDict from copy import deepcopy from numbers import Real from xml.etree import ElementTree as ET -from warnings import warn +from warnings import warn, catch_warnings, simplefilter +import math import numpy as np from openmc.checkvalue import check_type, check_value from openmc.region import Region, Intersection, Union -from openmc.mixin import IDManagerMixin +from openmc.mixin import IDManagerMixin, IDWarning _BOUNDARY_TYPES = ['transmission', 'vacuum', 'reflective', 'periodic', 'white'] @@ -369,10 +370,32 @@ class PlaneMixin(metaclass=ABCMeta): return np.array((a, b, c)) / math.sqrt(a*a + b*b + c*c) def bounding_box(self, side): + """Determine an axis-aligned bounding box. + + An axis-aligned bounding box for Plane 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.ndarray + Lower-left coordinates of the axis-aligned bounding box for the + desired half-space + numpy.ndarray + Upper-right coordinates of the axis-aligned bounding box for the + desired half-space + + """ # Compute the bounding box based on the normal vector to the plane nhat = self._get_normal() - lb = np.array([-np.inf, -np.inf, -np.inf]) - ub = np.array([np.inf, np.inf, np.inf]) + 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)): sign = nhat.sum() @@ -380,16 +403,16 @@ class PlaneMixin(metaclass=ABCMeta): vals = [d/val if round(val) != 0 else np.nan for val in (a, b, c)] if side == '-': if sign > 0: - ub = np.array([v if not np.isnan(v) else np.inf for v in vals]) + ur = np.array([v if not np.isnan(v) else np.inf for v in vals]) else: - lb = np.array([v if ~np.isnan(v) else -np.inf for v in vals]) + ll = np.array([v if not np.isnan(v) else -np.inf for v in vals]) elif side == '+': if sign > 0: - lb = np.array([v if ~np.isnan(v) else -np.inf for v in vals]) + ll = np.array([v if not np.isnan(v) else -np.inf for v in vals]) else: - ub = np.array([v if ~np.isnan(v) else np.inf for v in vals]) + ur = np.array([v if not np.isnan(v) else np.inf for v in vals]) - return (lb, ub) + return (ll, ur) def evaluate(self, point): """Evaluate the surface equation at a given point. @@ -1180,18 +1203,12 @@ class Cylinder(QuadricMixin, Surface): """ # This method overrides Surface.to_xml_element to generate a Quadric # since the C++ layer doesn't support Cylinders right now - element = ET.Element("surface") - element.set("id", str(self._id)) - - if len(self._name) > 0: - element.set("name", str(self._name)) - - element.set("type", 'quadric') - if self.boundary_type != 'transmission': - element.set("boundary", self.boundary_type) - element.set("coeffs", ' '.join([str(c) for c in self._get_base_coeffs()])) - - return element + with catch_warnings(): + simplefilter('ignore', IDWarning) + kwargs = {'boundary_type': self.boundary_type, 'name': self.name, + 'surface_id': self.id} + quad_rep = Quadric(*self._get_base_coeffs(), **kwargs) + return quad_rep.to_xml_element() class XCylinder(QuadricMixin, Surface): @@ -1881,19 +1898,13 @@ class Cone(QuadricMixin, Surface): """ # This method overrides Surface.to_xml_element to generate a Quadric - # since the C++ layer doesn't support Cylinders right now - element = ET.Element("surface") - element.set("id", str(self._id)) - - if len(self._name) > 0: - element.set("name", str(self._name)) - - element.set("type", 'quadric') - if self.boundary_type != 'transmission': - element.set("boundary", self.boundary_type) - element.set("coeffs", ' '.join([str(c) for c in self._get_base_coeffs()])) - - return element + # since the C++ layer doesn't support Cones right now + with catch_warnings(): + simplefilter('ignore', IDWarning) + kwargs = {'boundary_type': self.boundary_type, 'name': self.name, + 'surface_id': self.id} + quad_rep = Quadric(*self._get_base_coeffs(), **kwargs) + return quad_rep.to_xml_element() class XCone(QuadricMixin, Surface): diff --git a/setup.py b/setup.py index 99db831484..c193645c35 100755 --- a/setup.py +++ b/setup.py @@ -56,14 +56,13 @@ kwargs = { 'Topic :: Scientific/Engineering' 'Programming Language :: C++', 'Programming Language :: Python :: 3', - 'Programming Language :: Python :: 3.4', 'Programming Language :: Python :: 3.5', 'Programming Language :: Python :: 3.6', 'Programming Language :: Python :: 3.7', ], # Dependencies - 'python_requires': '>=3.4', + 'python_requires': '>=3.5', 'install_requires': [ 'numpy>=1.9', 'h5py', 'scipy', 'ipython', 'matplotlib', 'pandas', 'lxml', 'uncertainties' diff --git a/tests/unit_tests/test_surface.py b/tests/unit_tests/test_surface.py index 703010a059..86e01b6af3 100644 --- a/tests/unit_tests/test_surface.py +++ b/tests/unit_tests/test_surface.py @@ -370,23 +370,26 @@ def test_cylinder_from_points_axis(): # (x - 3)^2 + (y - 4)^2 = 2^2 # x^2 + y^2 - 6x - 8y + 21 = 0 s = openmc.model.cylinder_from_points((3., 4., 0.), (3., 4., 1.), 2.) - assert (s.a, s.b, s.c) == pytest.approx((1., 1., 0.)) - assert (s.d, s.e, s.f) == pytest.approx((0., 0., 0.)) - assert (s.g, s.h, s.j) == pytest.approx((-6., -8., 0.)) - assert s.k == pytest.approx(21.) + a, b, c, d, e, f, g, h, j, k = s._get_base_coeffs() + assert (a, b, c) == pytest.approx((1., 1., 0.)) + assert (d, e, f) == pytest.approx((0., 0., 0.)) + assert (g, h, j) == pytest.approx((-6., -8., 0.)) + assert k == pytest.approx(21.) # (y + 7)^2 + (z - 1)^2 = 3^2 # y^2 + z^2 + 14y - 2z + 41 = 0 s = openmc.model.cylinder_from_points((0., -7, 1.), (1., -7., 1.), 3.) - assert (s.a, s.b, s.c) == pytest.approx((0., 1., 1.)) - assert (s.d, s.e, s.f) == pytest.approx((0., 0., 0.)) - assert (s.g, s.h, s.j) == pytest.approx((0., 14., -2.)) - assert s.k == 41. + a, b, c, d, e, f, g, h, j, k = s._get_base_coeffs() + assert (a, b, c) == pytest.approx((0., 1., 1.)) + assert (d, e, f) == pytest.approx((0., 0., 0.)) + assert (g, h, j) == pytest.approx((0., 14., -2.)) + assert k == 41. # (x - 2)^2 + (z - 5)^2 = 4^2 # x^2 + z^2 - 4x - 10z + 13 = 0 s = openmc.model.cylinder_from_points((2., 0., 5.), (2., 1., 5.), 4.) - assert (s.a, s.b, s.c) == pytest.approx((1., 0., 1.)) - assert (s.d, s.e, s.f) == pytest.approx((0., 0., 0.)) - assert (s.g, s.h, s.j) == pytest.approx((-4., 0., -10.)) - assert s.k == pytest.approx(13.) + a, b, c, d, e, f, g, h, j, k = s._get_base_coeffs() + assert (a, b, c) == pytest.approx((1., 0., 1.)) + assert (d, e, f) == pytest.approx((0., 0., 0.)) + assert (g, h, j) == pytest.approx((-4., 0., -10.)) + assert k == pytest.approx(13.) From 1e5a8e56ff8b2195fdeb09161f945ac128f3913d Mon Sep 17 00:00:00 2001 From: Ethan Peterson Date: Mon, 24 Feb 2020 12:06:12 -0500 Subject: [PATCH 6/8] removed use of round in favor of isclose --- openmc/surface.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/openmc/surface.py b/openmc/surface.py index 4d0434a4b0..cc92a810b3 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -400,7 +400,8 @@ class PlaneMixin(metaclass=ABCMeta): if np.any(np.isclose(np.abs(nhat), 1., rtol=0., atol=self._atol)): sign = nhat.sum() a, b, c, d = self._get_base_coeffs() - vals = [d/val if round(val) != 0 else np.nan for val in (a, b, c)] + 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]) From e3af035095bc7783b127d97df956f06edef21f70 Mon Sep 17 00:00:00 2001 From: Ethan Peterson Date: Mon, 24 Feb 2020 15:30:28 -0500 Subject: [PATCH 7/8] finished unit tests and fixed wrong math --- openmc/surface.py | 44 ++++++----- tests/unit_tests/test_surface.py | 121 +++++++++++++++++++++++++++++++ 2 files changed, 146 insertions(+), 19 deletions(-) diff --git a/openmc/surface.py b/openmc/surface.py index cc92a810b3..e4c11692a1 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -887,6 +887,15 @@ class ZPlane(PlaneMixin, Surface): class QuadricMixin(metaclass=ABCMeta): """A Mixin class implementing common functionality for quadric surfaces""" + @property + def _origin(self): + return np.array((self.x0, self.y0, self.z0)) + + @property + def _axis(self): + axis = np.array((self.dx, self.dy, self.dz)) + return axis / np.linalg.norm(axis) + def get_Abc(self, coeffs=None): """Compute matrix, vector, and scalar coefficients for this surface or for a specified set of coefficients. @@ -1135,8 +1144,8 @@ class Cylinder(QuadricMixin, Surface): def _get_base_coeffs(self): # Get x, y, z coordinates of two points - x1, y1, z1 = self.x0, self.y0, self.z0 - x2, y2, z2 = x1 + self.dx, y1 + self.dy, z1 + self.dz + x1, y1, z1 = self._origin + x2, y2, z2 = self._origin + self._axis r = self.r # Define intermediate terms @@ -1868,24 +1877,21 @@ class Cone(QuadricMixin, Surface): # The argument r2 for cones is actually tan^2(theta) so that # cos^2(theta) = 1 / (1 + r2) - x0, y0, z0, r2 = self.x0, self.y0, self.z0, self.r2 - dx, dy, dz = self.dx, self.dy, self.dz - dnorm = dx*dx + dy*dy + dz*dz - dx /= dnorm - dy /= dnorm - dz /= dnorm - cos2 = 1 / (1 + r2) + x0, y0, z0 = self._origin + dx, dy, dz = self._axis + cos2 = 1 / (1 + self.r2) - a = dx*dx - cos2 - b = dy*dy - cos2 - c = dz*dz - cos2 - d = 2*dx*dy - e = 2*dy*dz - f = 2*dx*dz - g = -2*(dx*dx*x0 + dx*dy*y0 + dx*dz*z0 - cos2) - h = -2*(dy*dy*y0 + dx*dy*x0 + dy*dz*z0 - cos2) - j = -2*(dz*dz*z0 + dx*dz*x0 + dy*dz*y0 - cos2) - k = (dx*x0 + dy*y0 + dz*z0)**2 - cos2*(x0*x0 + y0*y0 + z0*z0) + a = cos2 - dx*dx + b = cos2 - dy*dy + c = cos2 - dz*dz + d = -2*dx*dy + e = -2*dy*dz + f = -2*dx*dz + g = 2*(dx*(dy*y0 + dz*z0) - a*x0) + h = 2*(dy*(dx*x0 + dz*z0) - b*y0) + j = 2*(dz*(dx*x0 + dy*y0) - c*z0) + k = a*x0*x0 + b*y0*y0 + c*z0*z0 - 2*(dx*dy*x0*y0 + dy*dz*y0*z0 + + dx*dz*x0*z0) return (a, b, c, d, e, f, g, h, j, k) diff --git a/tests/unit_tests/test_surface.py b/tests/unit_tests/test_surface.py index 86e01b6af3..b6c7398778 100644 --- a/tests/unit_tests/test_surface.py +++ b/tests/unit_tests/test_surface.py @@ -2,6 +2,7 @@ from functools import partial from random import uniform, seed import numpy as np +import math import openmc import pytest @@ -139,6 +140,63 @@ def test_zplane(): repr(s) +def test_cylinder(): + x0, y0, z0, r = 2, 3, 4, 2 + dx, dy, dz = 1, -1, 1 + s = openmc.Cylinder(x0=x0, y0=y0, z0=z0, dx=dx, dy=dy, dz=dz, r=r) + assert s.x0 == 2 + assert s.y0 == 3 + assert s.z0 == 4 + assert s.dx == 1 + assert s.dy == -1 + assert s.dz == 1 + assert s.r == 2 + + # Check bounding box + assert_infinite_bb(s) + + # evaluate method + # |(p - p1) ⨯ (p - p2)|^2 / |p2 - p1|^2 - r^2 + # point inside + p = s._origin + 5*s._axis + p1 = s._origin + p2 = p1 + s._axis + c1 = np.linalg.norm(np.cross(p - p1, p - p2)) / np.linalg.norm(p2 - p1) + val = c1*c1 - s.r*s.r + assert s.evaluate(p) < 0 + assert s.evaluate(p) == pytest.approx(val) + # point outside + p = np.array((4., 0., 2.5)) + p1 = s._origin + p2 = p1 + s._axis + c1 = np.linalg.norm(np.cross(p - p1, p - p2)) / np.linalg.norm(p2 - p1) + val = c1*c1 - s.r*s.r + assert s.evaluate(p) > 0 + assert s.evaluate(p) == pytest.approx(val) + # point on cylinder + perp = np.array((1, -2, 1))*(1 / s._axis) + p = s._origin + s.r*perp / np.linalg.norm(perp) + p1 = s._origin + p2 = p1 + s._axis + c1 = np.linalg.norm(np.cross(p - p1, p - p2)) / np.linalg.norm(p2 - p1) + val = c1*c1 - s.r*s.r + assert 0 == pytest.approx(s.evaluate(p)) + assert val == pytest.approx(s.evaluate(p)) + + # translate method + st = s.translate((1.0, 1.0, 1.0)) + assert st.x0 == s.x0 + 1 + assert st.y0 == s.y0 + 1 + assert st.z0 == s.z0 + 1 + assert st.dx == s.dx + assert st.dy == s.dy + assert st.dz == s.dz + assert st.r == s.r + + # Make sure repr works + repr(s) + + def test_xcylinder(): y, z, r = 3, 5, 2 s = openmc.XCylinder(y0=y, z0=z, r=r) @@ -284,6 +342,69 @@ def cone_common(apex, r2, cls): repr(s) +def test_cone(): + x0, y0, z0, r2 = 2, 3, 4, 4 + dx, dy, dz = 1, -1, 1 + s = openmc.Cone(x0=x0, y0=y0, z0=z0, dx=dx, dy=dy, dz=dz, r2=r2) + assert s.x0 == 2 + assert s.y0 == 3 + assert s.z0 == 4 + assert s.dx == 1 + assert s.dy == -1 + assert s.dz == 1 + assert s.r2 == 4 + + # Check bounding box + assert_infinite_bb(s) + + # evaluate method + # cos^2(theta) * ((p - p1))**2 - (d @ (p - p1))^2 + # The argument r2 for cones is actually tan^2(theta) so that + # cos^2(theta) = 1 / (1 + r2) + # + # This makes the evaluation equation shown below where p is the evaluation + # point (x, y, z) p1 is the apex (origin) of the cone and r2 is related to + # the aperature of the cone as described above + # (p - p1) @ (p - p1) / (1 + r2) - (d @ (p - p1))^2 + # point inside + p1 = s._origin + d = s._axis + p = p1 + 5*d + val = np.sum((p - p1)**2) / (1 + s.r2) - np.sum((d @ (p - p1))**2) + assert s.evaluate(p) < 0 + assert s.evaluate(p) == pytest.approx(val) + # point outside + p1 = s._origin + d = s._axis + perp = np.array((1, -2, 1))*(1 / d) + p = p1 + 3.2*perp + val = np.sum((p - p1)**2) / (1 + s.r2) - np.sum((d @ (p - p1))**2) + assert s.evaluate(p) > 0 + assert s.evaluate(p) == pytest.approx(val) + # point on cone + p1 = s._origin + d = s._axis + perp = np.array((1, -2, 1))*(1 / d) + perp /= np.linalg.norm(perp) + p = p1 + 3.2*d + 3.2*math.sqrt(s.r2)*perp + val = np.sum((p - p1)**2) / (1 + s.r2) - np.sum((d @ (p - p1))**2) + assert 0 == pytest.approx(s.evaluate(p)) + assert val == pytest.approx(s.evaluate(p)) + + # translate method + st = s.translate((1.0, 1.0, 1.0)) + assert st.x0 == s.x0 + 1 + assert st.y0 == s.y0 + 1 + assert st.z0 == s.z0 + 1 + assert st.dx == s.dx + assert st.dy == s.dy + assert st.dz == s.dz + assert st.r2 == s.r2 + + # Make sure repr works + repr(s) + + def test_xcone(): apex = (10, 0, 0) r2 = 4 From e6cf49275ab18b6cef7845d570b0fab408b2f906 Mon Sep 17 00:00:00 2001 From: Ethan Peterson Date: Tue, 25 Feb 2020 08:40:58 -0500 Subject: [PATCH 8/8] applied suggestions from code review --- openmc/surface.py | 2 +- tests/unit_tests/test_surface.py | 67 +++++++++++--------------------- 2 files changed, 23 insertions(+), 46 deletions(-) diff --git a/openmc/surface.py b/openmc/surface.py index e4c11692a1..567968b2c8 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -366,7 +366,7 @@ class PlaneMixin(metaclass=ABCMeta): return (self.a, self.b, self.c, self.d) def _get_normal(self): - a, b, c = self._get_base_coeffs()[0:3] + a, b, c = self._get_base_coeffs()[:3] return np.array((a, b, c)) / math.sqrt(a*a + b*b + c*c) def bounding_box(self, side): diff --git a/tests/unit_tests/test_surface.py b/tests/unit_tests/test_surface.py index b6c7398778..749ae96212 100644 --- a/tests/unit_tests/test_surface.py +++ b/tests/unit_tests/test_surface.py @@ -157,31 +157,19 @@ def test_cylinder(): # evaluate method # |(p - p1) ⨯ (p - p2)|^2 / |p2 - p1|^2 - r^2 - # point inside - p = s._origin + 5*s._axis p1 = s._origin p2 = p1 + s._axis - c1 = np.linalg.norm(np.cross(p - p1, p - p2)) / np.linalg.norm(p2 - p1) - val = c1*c1 - s.r*s.r - assert s.evaluate(p) < 0 - assert s.evaluate(p) == pytest.approx(val) - # point outside - p = np.array((4., 0., 2.5)) - p1 = s._origin - p2 = p1 + s._axis - c1 = np.linalg.norm(np.cross(p - p1, p - p2)) / np.linalg.norm(p2 - p1) - val = c1*c1 - s.r*s.r - assert s.evaluate(p) > 0 - assert s.evaluate(p) == pytest.approx(val) - # point on cylinder perp = np.array((1, -2, 1))*(1 / s._axis) - p = s._origin + s.r*perp / np.linalg.norm(perp) - p1 = s._origin - p2 = p1 + s._axis - c1 = np.linalg.norm(np.cross(p - p1, p - p2)) / np.linalg.norm(p2 - p1) - val = c1*c1 - s.r*s.r - assert 0 == pytest.approx(s.evaluate(p)) - assert val == pytest.approx(s.evaluate(p)) + divisor = np.linalg.norm(p2 - p1) + pin = p1 + 5*s._axis # point inside cylinder + pout = np.array((4., 0., 2.5)) # point outside the cylinder + pon = p1 + s.r*perp / np.linalg.norm(perp) # point on cylinder + for p, fn in zip((pin, pout, pon), (np.less, np.greater, np.isclose)): + c1 = np.linalg.norm(np.cross(p - p1, p - p2)) / divisor + val = c1*c1 - s.r*s.r + p_eval = s.evaluate(p) + assert fn(p_eval, 0.) + assert p_eval == pytest.approx(val) # translate method st = s.translate((1.0, 1.0, 1.0)) @@ -369,27 +357,16 @@ def test_cone(): # point inside p1 = s._origin d = s._axis - p = p1 + 5*d - val = np.sum((p - p1)**2) / (1 + s.r2) - np.sum((d @ (p - p1))**2) - assert s.evaluate(p) < 0 - assert s.evaluate(p) == pytest.approx(val) - # point outside - p1 = s._origin - d = s._axis - perp = np.array((1, -2, 1))*(1 / d) - p = p1 + 3.2*perp - val = np.sum((p - p1)**2) / (1 + s.r2) - np.sum((d @ (p - p1))**2) - assert s.evaluate(p) > 0 - assert s.evaluate(p) == pytest.approx(val) - # point on cone - p1 = s._origin - d = s._axis perp = np.array((1, -2, 1))*(1 / d) perp /= np.linalg.norm(perp) - p = p1 + 3.2*d + 3.2*math.sqrt(s.r2)*perp - val = np.sum((p - p1)**2) / (1 + s.r2) - np.sum((d @ (p - p1))**2) - assert 0 == pytest.approx(s.evaluate(p)) - assert val == pytest.approx(s.evaluate(p)) + pin = p1 + 5*d # point inside cone + pout = p1 + 3.2*perp # point outside cone + pon = p1 + 3.2*d + 3.2*math.sqrt(s.r2)*perp # point on cone + for p, fn in zip((pin, pout, pon), (np.less, np.greater, np.isclose)): + val = np.sum((p - p1)**2) / (1 + s.r2) - np.sum((d @ (p - p1))**2) + p_eval = s.evaluate(p) + assert fn(p_eval, 0.) + assert p_eval == pytest.approx(val) # translate method st = s.translate((1.0, 1.0, 1.0)) @@ -460,7 +437,7 @@ def test_cylinder_from_points(): p1 = np.array([xi(), xi(), xi()]) p2 = np.array([xi(), xi(), xi()]) r = uniform(1.0, 100.0) - s = openmc.model.cylinder_from_points(p1, p2, r) + s = openmc.Cylinder.from_points(p1, p2, r) # Points p1 and p2 need to be inside cylinder assert p1 in -s @@ -490,7 +467,7 @@ def test_cylinder_from_points_axis(): # (x - 3)^2 + (y - 4)^2 = 2^2 # x^2 + y^2 - 6x - 8y + 21 = 0 - s = openmc.model.cylinder_from_points((3., 4., 0.), (3., 4., 1.), 2.) + s = openmc.Cylinder.from_points((3., 4., 0.), (3., 4., 1.), 2.) a, b, c, d, e, f, g, h, j, k = s._get_base_coeffs() assert (a, b, c) == pytest.approx((1., 1., 0.)) assert (d, e, f) == pytest.approx((0., 0., 0.)) @@ -499,7 +476,7 @@ def test_cylinder_from_points_axis(): # (y + 7)^2 + (z - 1)^2 = 3^2 # y^2 + z^2 + 14y - 2z + 41 = 0 - s = openmc.model.cylinder_from_points((0., -7, 1.), (1., -7., 1.), 3.) + s = openmc.Cylinder.from_points((0., -7, 1.), (1., -7., 1.), 3.) a, b, c, d, e, f, g, h, j, k = s._get_base_coeffs() assert (a, b, c) == pytest.approx((0., 1., 1.)) assert (d, e, f) == pytest.approx((0., 0., 0.)) @@ -508,7 +485,7 @@ def test_cylinder_from_points_axis(): # (x - 2)^2 + (z - 5)^2 = 4^2 # x^2 + z^2 - 4x - 10z + 13 = 0 - s = openmc.model.cylinder_from_points((2., 0., 5.), (2., 1., 5.), 4.) + s = openmc.Cylinder.from_points((2., 0., 5.), (2., 1., 5.), 4.) a, b, c, d, e, f, g, h, j, k = s._get_base_coeffs() assert (a, b, c) == pytest.approx((1., 0., 1.)) assert (d, e, f) == pytest.approx((0., 0., 0.))