From 2870da61e63b7dfa4d9c425d94f2736a7cdf8a54 Mon Sep 17 00:00:00 2001 From: Ethan Peterson Date: Fri, 21 Feb 2020 20:49:06 -0500 Subject: [PATCH] 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__()}