From 6d6988b2c8cb2ce611234e1ce546580f1cbacb82 Mon Sep 17 00:00:00 2001 From: Ethan Peterson Date: Mon, 10 Feb 2020 14:15:28 -0500 Subject: [PATCH 01/33] refactoring surface classes --- openmc/surface.py | 715 +++++++++++++++++++++++++++------------------- 1 file changed, 423 insertions(+), 292 deletions(-) diff --git a/openmc/surface.py b/openmc/surface.py index 55e3f459d3..d64ca242fd 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -1,5 +1,6 @@ from abc import ABCMeta, abstractmethod from collections import OrderedDict +from abc.collections import Iterable from copy import deepcopy from numbers import Real, Integral from xml.etree import ElementTree as ET @@ -7,7 +8,7 @@ from warnings import warn import numpy as np -from openmc.checkvalue import check_type, check_value +from openmc.checkvalue import check_type, check_value, check_length from openmc.region import Region, Intersection, Union from openmc.mixin import IDManagerMixin @@ -187,6 +188,17 @@ class Surface(IDManagerMixin, metaclass=ABCMeta): def translate(self, vector): pass + @classmethod + def get_subclasses(cls): + """Recursively find all subclasses of this class""" + return set(cls.__subclasses__()).union([s for c in cls.__subclasses__() + for s in get_subclasses(c)]) + + @classmethod + def get_subclass_map(cls): + """Generate mapping of class _type attributes to classes""" + return {c._type: c for c in cls.get_subclasses()} + def to_xml_element(self): """Return XML representation of the surface @@ -228,21 +240,7 @@ class Surface(IDManagerMixin, metaclass=ABCMeta): # Determine appropriate class surf_type = elem.get('type') - surface_classes = { - 'plane': Plane, - 'x-plane': XPlane, - 'y-plane': YPlane, - 'z-plane': ZPlane, - 'x-cylinder': XCylinder, - 'y-cylinder': YCylinder, - 'z-cylinder': ZCylinder, - 'sphere': Sphere, - 'x-cone': XCone, - 'y-cone': YCone, - 'z-cone': ZCone, - 'quadric': Quadric, - } - cls = surface_classes[surf_type] + cls = _SURFACE_CLASSES[surf_type] # Determine ID, boundary type, coefficients kwargs = {} @@ -268,62 +266,215 @@ class Surface(IDManagerMixin, metaclass=ABCMeta): Instance of surface subclass """ + surface_id = int(group.name.split('/')[-1].lstrip('surface ')) name = group['name'][()].decode() if 'name' in group else '' surf_type = group['type'][()].decode() bc = group['boundary_type'][()].decode() coeffs = group['coefficients'][...] - # Create the Surface based on its type - if surf_type == 'x-plane': - x0 = coeffs[0] - surface = XPlane(x0, bc, name, surface_id) + cls = _SURFACE_CLASSES[surf_type] - elif surf_type == 'y-plane': - y0 = coeffs[0] - surface = YPlane(y0, bc, name, surface_id) - - elif surf_type == 'z-plane': - z0 = coeffs[0] - surface = ZPlane(z0, bc, name, surface_id) - - elif surf_type == 'plane': - A, B, C, D = coeffs - surface = Plane(A, B, C, D, bc, name, surface_id) - - elif surf_type == 'x-cylinder': - y0, z0, r = coeffs - surface = XCylinder(y0, z0, r, bc, name, surface_id) - - elif surf_type == 'y-cylinder': - x0, z0, r = coeffs - surface = YCylinder(x0, z0, r, bc, name, surface_id) - - elif surf_type == 'z-cylinder': - x0, y0, r = coeffs - surface = ZCylinder(x0, y0, r, bc, name, surface_id) - - elif surf_type == 'sphere': - x0, y0, z0, r = coeffs - surface = Sphere(x0, y0, z0, r, bc, name, surface_id) - - elif surf_type in ['x-cone', 'y-cone', 'z-cone']: - x0, y0, z0, r2 = coeffs - if surf_type == 'x-cone': - surface = XCone(x0, y0, z0, r2, bc, name, surface_id) - elif surf_type == 'y-cone': - surface = YCone(x0, y0, z0, r2, bc, name, surface_id) - elif surf_type == 'z-cone': - surface = ZCone(x0, y0, z0, r2, bc, name, surface_id) - - elif surf_type == 'quadric': - a, b, c, d, e, f, g, h, j, k = coeffs - surface = Quadric(a, b, c, d, e, f, g, h, j, k, bc, name, surface_id) - - return surface + return cls(*coeffs, bc, name, surface_id) -class Plane(Surface): +_SURFACE_CLASSES = Surface.get_subclass_map() + + +class PlaneMeta(metaclass=ABCMeta): + """A Plane Meta class for all operations on order 1 surfaces""" + + def __new__(cls, *args, **kwargs): + """Simplify this plane if possible to an XPlane, YPlane, or ZPlane""" + if cls is Plane: + (a,b,c,d,bt, name, surfidlen(args) + if np.all(np.isclose((self.b, self.c), 0., atol=atol)): + x0 = self.d / self.a + return XPlane(x0=x0, boundary_type=self.boundary_type, + name=self.name, surface_id=self.id) + + elif np.all(np.isclose((self.a, self.c), 0., atol=atol)): + y0 = self.d / self.b + return YPlane(y0=y0, boundary_type=self.boundary_type, + name=self.name, surface_id=self.id) + + elif np.all(np.isclose((self.a, self.b), 0., atol=atol)): + z0 = self.d / self.c + return ZPlane(z0=z0, boundary_type=self.boundary_type, + name=self.name, surface_id=self.id) + + def __init__(self): + pass + + @property + def periodic_surface(self): + return self._periodic_surface + + @periodic_surface.setter + def periodic_surface(self, periodic_surface): + check_type('periodic surface', periodic_surface, Plane) + self._periodic_surface = periodic_surface + periodic_surface._periodic_surface = self + + @abstractmethod + def _get_base_coeffs(self): + pass + + @abstractmethod + def _update_from_base_coeffs(self): + pass + + def _neg_bounds(self): + return (np.array([-np.inf, -np.inf, -np.inf]), + np.array([np.inf, np.inf, np.inf])) + + def _pos_bounds(self): + return (np.array([-np.inf, -np.inf, -np.inf]), + np.array([np.inf, np.inf, np.inf])) + + def evaluate(self, point): + """Evaluate the surface equation at a given point. + + Parameters + ---------- + point : 3-tuple of float + The Cartesian coordinates, :math:`(x',y',z')`, at which the surface + equation should be evaluated. + + Returns + ------- + float + :math:`Ax' + By' + Cz' - D` + + """ + + x, y, z = point + a, b, c, d = self._get_base_coeffs() + return a*x + b*y + c*z - d + + def translate(self, vector, clone=False): + """Translate surface in given direction + + Parameters + ---------- + vector : iterable of float + Direction in which surface should be translated + clone : boolean + Whether or not to return a new instance of a Plane or to modify the + coefficients of this plane. + + Returns + ------- + openmc.Plane + Translated surface + + """ + vx, vy, vz = vector + a, b, c, d = self._get_base_coeffs() + d = d + a*vx + b*vy + c*vz + + if clone: + surf = self.clone() + else: + surf = self + + surf._update_from_base_coeffs(a, b, c, d) + return surf + + def rotate(self, rotation, frame='lab', clone=False): + """Rotate surface by given Tait-Bryan angles + + Parameters + ---------- + rotation : iterable of float + Intrinsic Tait-Bryan angles in degrees used to rotate the surface + frame : str, one of 'lab' or 'body' + clone : boolean + Whether or not to return a new instance of a Plane or to modify the + coefficients of this plane. + + Returns + ------- + openmc.Plane or None + Rotated surface + + """ + + check_type('surface rotation', rotation, Iterable, Real) + check_length('surface rotation', rotation, 3) + + # Calculate rotation matrix Rmat from angles phi, theta, psi + phi, theta, psi = rotation*(-np.pi/180.) + c3, s3 = np.cos(phi), np.sin(phi) + c2, s2 = np.cos(theta), np.sin(theta) + c1, s1 = np.cos(psi), np.sin(psi) + Rmat = np.array([[c1*c2, c1*s2*s3 - c3*s1, s1*s3 + c1*c3*s2], + [c2*s1, c1*c3 + s1*s2*s3, c3*s1*s2 - c1*s3], + [-s2, c2*s3, c2*c3]]) + + a, b, c, d = self._get_base_coeffs() + bvec = np.array([a, b, c]) + + # Compute new rotated coefficients a, b, c + a, b, c = np.dot(bvec.T, Rmat.T) + + if clone: + surf = self.clone() + else: + surf = self + + surf._update_from_base_coeffs(a, b, c, d) + return surf + + 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.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 + + """ + + if side == '-': + return self._neg_bounds() + elif side == '+': + return self._pos_bounds() + + def to_xml_element(self): + """Return XML representation of the surface + + Returns + ------- + element : xml.etree.ElementTree.Element + XML element containing source data + + """ + element = super().to_xml_element() + + # Add periodic surface pair information + if self.boundary_type == 'periodic': + if self.periodic_surface is not None: + element.set("periodic_surface_id", + str(self.periodic_surface.id)) + return element + + +class Plane(PlaneMeta, Surface): """An arbitrary plane of the form :math:`Ax + By + Cz = D`. Parameters @@ -406,10 +557,6 @@ class Plane(Surface): def d(self): return self.coefficients['d'] - @property - def periodic_surface(self): - return self._periodic_surface - @a.setter def a(self, a): check_type('A coefficient', a, Real) @@ -430,68 +577,15 @@ class Plane(Surface): check_type('D coefficient', d, Real) self._coefficients['d'] = d - @periodic_surface.setter - def periodic_surface(self, periodic_surface): - check_type('periodic surface', periodic_surface, Plane) - 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 evaluate(self, point): - """Evaluate the surface equation at a given point. - - Parameters - ---------- - point : 3-tuple of float - The Cartesian coordinates, :math:`(x',y',z')`, at which the surface - equation should be evaluated. - - Returns - ------- - float - :math:`Ax' + By' + Cz' - D` - - """ - - x, y, z = point - return self.a*x + self.b*y + self.c*z - self.d - - def translate(self, vector): - """Translate surface in given direction - - Parameters - ---------- - vector : iterable of float - Direction in which surface should be translated - - Returns - ------- - openmc.Plane - Translated surface - - """ - vx, vy, vz = vector - d = self.d + self.a*vx + self.b*vy + self.c*vz - if d == self.d: - return self - else: - return type(self)(a=self.a, b=self.b, c=self.c, d=d) - - def to_xml_element(self): - """Return XML representation of the surface - - Returns - ------- - element : xml.etree.ElementTree.Element - XML element containing source data - - """ - element = super().to_xml_element() - - # Add periodic surface pair information - if self.boundary_type == 'periodic': - if self.periodic_surface is not None: - element.set("periodic_surface_id", str(self.periodic_surface.id)) - return element + def _update_from_base_coeffs(self, coeffs): + a, b, c, d = coeffs + self.a = a + self.b = b + self.c = c + self.d = d @classmethod def from_points(cls, p1, p2, p3, **kwargs): @@ -526,7 +620,7 @@ class Plane(Surface): return cls(a=a, b=b, c=c, d=d, **kwargs) -class XPlane(Plane): +class XPlane(PlaneMeta, Surface): """A plane perpendicular to the x axis of the form :math:`x - x_0 = 0` Parameters @@ -582,76 +676,25 @@ class XPlane(Plane): check_type('x0 coefficient', x0, Real) self._coefficients['x0'] = x0 - def bounding_box(self, side): - """Determine an axis-aligned bounding box. + def _get_base_coeffs(self): + return (1., 0., 0., self.x0) - 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. + def _update_from_base_coeffs(self, coeffs): + a, b, c, d = coeffs + self.x0 = d - Parameters - ---------- - side : {'+', '-'} - Indicates the negative or positive half-space + def _neg_bounds(self): + """Return the lower and upper bounds of the negative half space""" + return (np.array([-np.inf, -np.inf, -np.inf]), + np.array([self.x0, np.inf, np.inf])) - 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 - - """ - - 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): - """Evaluate the surface equation at a given point. - - Parameters - ---------- - point : 3-tuple of float - The Cartesian coordinates, :math:`(x',y',z')`, at which the surface - equation should be evaluated. - - Returns - ------- - float - :math:`x' - x_0` - - """ - return point[0] - self.x0 - - def translate(self, vector): - """Translate surface in given direction - - Parameters - ---------- - vector : iterable of float - Direction in which surface should be translated - - Returns - ------- - openmc.XPlane - Translated surface - - """ - vx = vector[0] - if vx == 0: - return self - else: - return type(self)(x0=self.x0 + vx) + def _pos_bounds(self): + """Return the lower and upper bounds of the positive half space""" + return (np.array([self.x0, -np.inf, -np.inf]), + np.array([np.inf, np.inf, np.inf])) -class YPlane(Plane): +class YPlane(PlaneMeta, Surface): """A plane perpendicular to the y axis of the form :math:`y - y_0 = 0` Parameters @@ -707,76 +750,26 @@ class YPlane(Plane): check_type('y0 coefficient', y0, Real) self._coefficients['y0'] = y0 - def bounding_box(self, side): - """Determine an axis-aligned bounding box. + def _get_base_coeffs(self): + """Return generalized coefficients for a plane""" + return (0., 1., 0., self.y0) - 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. + def _update_from_base_coeffs(self, coeffs): + a, b, c, d = coeffs + self.y0 = d - Parameters - ---------- - side : {'+', '-'} - Indicates the negative or positive half-space + def _neg_bounds(self): + """Return the lower and upper bounds of the negative half space""" + return (np.array([-np.inf, -np.inf, -np.inf]), + np.array([np.inf, self.y0, np.inf])) - 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 - - """ - - 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): - """Evaluate the surface equation at a given point. - - Parameters - ---------- - point : 3-tuple of float - The Cartesian coordinates, :math:`(x',y',z')`, at which the surface - equation should be evaluated. - - Returns - ------- - float - :math:`y' - y_0` - - """ - return point[1] - self.y0 - - def translate(self, vector): - """Translate surface in given direction - - Parameters - ---------- - vector : iterable of float - Direction in which surface should be translated - - Returns - ------- - openmc.YPlane - Translated surface - - """ - vy = vector[1] - if vy == 0.0: - return self - else: - return type(self)(y0=self.y0 + vy) + def _pos_bounds(self): + """Return the lower and upper bounds of the positive half space""" + return (np.array([-np.inf, self.y0, -np.inf]), + np.array([np.inf, np.inf, np.inf])) -class ZPlane(Plane): +class ZPlane(PlaneMeta, Surface): """A plane perpendicular to the z axis of the form :math:`z - z_0 = 0` Parameters @@ -832,36 +825,166 @@ class ZPlane(Plane): check_type('z0 coefficient', z0, Real) self._coefficients['z0'] = z0 - def bounding_box(self, side): - """Determine an axis-aligned bounding box. + def _get_base_coeffs(self): + """Return generalized coefficients for a plane""" + return (0., 0., 1., self.z0) - 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. + def _update_from_base_coeffs(self, coeffs): + a, b, c, d = coeffs + self.z0 = d - Parameters - ---------- - side : {'+', '-'} - Indicates the negative or positive half-space + def _neg_bounds(self): + """Return the lower and upper bounds of the negative half space""" + return (np.array([-np.inf, -np.inf, -np.inf]), + np.array([np.inf, np.inf, self.z0])) - 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 + def _pos_bounds(self): + """Return the lower and upper bounds of the positive half space""" + return (np.array([-np.inf, -np.inf, self.z0]), + np.array([np.inf, np.inf, np.inf])) - """ +class QuadricMeta(metaclass=ABCMeta): + """A surface of the form :math:`Ax^2 + By^2 + Cz^2 + Dxy + Eyz + Fxz + Gx + Hy + + Jz + K = 0`. - 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])) + Parameters + ---------- + a, b, c, d, e, f, g, h, j, k : float, optional + coefficients for the surface. All default to 0. + boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic', 'white'}, optional + Boundary condition that defines the behavior for particles hitting the + surface. Defaults to transmissive boundary condition where particles + freely pass through the surface. + name : str, optional + Name of the surface. If not specified, the name will be the empty string. + surface_id : int, optional + Unique identifier for the surface. If not specified, an identifier will + automatically be assigned. + + Attributes + ---------- + a, b, c, d, e, f, g, h, j, k : float + coefficients for the surface + boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic', 'white'} + Boundary condition that defines the behavior for particles hitting the + surface. + coefficients : dict + Dictionary of surface coefficients + id : int + Unique identifier for the surface + name : str + Name of the surface + type : str + Type of the surface + + """ + + _type = 'quadric' + _coeff_keys = ('a', 'b', 'c', 'd', 'e', 'f', 'g', 'h', 'j', 'k') + + def __init__(self, a=0., b=0., c=0., d=0., e=0., f=0., g=0., h=0., j=0., + k=0., boundary_type='transmission', name='', surface_id=None): + super().__init__(surface_id, boundary_type, name=name) + self.a = a + self.b = b + self.c = c + self.d = d + self.e = e + self.f = f + self.g = g + self.h = h + self.j = j + self.k = k + + @property + def a(self): + return self.coefficients['a'] + + @property + def b(self): + return self.coefficients['b'] + + @property + def c(self): + return self.coefficients['c'] + + @property + def d(self): + return self.coefficients['d'] + + @property + def e(self): + return self.coefficients['e'] + + @property + def f(self): + return self.coefficients['f'] + + @property + def g(self): + return self.coefficients['g'] + + @property + def h(self): + return self.coefficients['h'] + + @property + def j(self): + return self.coefficients['j'] + + @property + def k(self): + return self.coefficients['k'] + + @a.setter + def a(self, a): + check_type('a coefficient', a, Real) + self._coefficients['a'] = a + + @b.setter + def b(self, b): + check_type('b coefficient', b, Real) + self._coefficients['b'] = b + + @c.setter + def c(self, c): + check_type('c coefficient', c, Real) + self._coefficients['c'] = c + + @d.setter + def d(self, d): + check_type('d coefficient', d, Real) + self._coefficients['d'] = d + + @e.setter + def e(self, e): + check_type('e coefficient', e, Real) + self._coefficients['e'] = e + + @f.setter + def f(self, f): + check_type('f coefficient', f, Real) + self._coefficients['f'] = f + + @g.setter + def g(self, g): + check_type('g coefficient', g, Real) + self._coefficients['g'] = g + + @h.setter + def h(self, h): + check_type('h coefficient', h, Real) + self._coefficients['h'] = h + + @j.setter + def j(self, j): + check_type('j coefficient', j, Real) + self._coefficients['j'] = j + + @k.setter + def k(self, k): + check_type('k coefficient', k, Real) + self._coefficients['k'] = k def evaluate(self, point): """Evaluate the surface equation at a given point. @@ -875,10 +998,14 @@ class ZPlane(Plane): Returns ------- float - :math:`z' - z_0` + :math:`Ax'^2 + By'^2 + Cz'^2 + Dx'y' + Ey'z' + Fx'z' + Gx' + Hy' + + Jz' + K = 0` """ - return point[2] - self.z0 + x, y, z = point + return x*(self.a*x + self.d*y + self.g) + \ + y*(self.b*y + self.e*z + self.h) + \ + z*(self.c*z + self.f*x + self.j) + self.k def translate(self, vector): """Translate surface in given direction @@ -890,15 +1017,19 @@ class ZPlane(Plane): Returns ------- - openmc.ZPlane + openmc.Quadric Translated surface """ - vz = vector[2] - if vz == 0.0: - return self - else: - return type(self)(z0=self.z0 + vz) + vx, vy, vz = vector + a, b, c, d, e, f, g, h, j, k = (getattr(self, key) for key in + self._coeff_keys) + k = (k + vx*vx + vy*vy + vz*vz + d*vx*vy + e*vy*vz + f*vx*vz + - g*vx - h*vy - j*vz) + g = g - 2*a*vx - d*vy - f*vz + h = h - 2*b*vy - d*vx - e*vz + j = j - 2*c*vz - e*vy - f*vx + return type(self)(a=a, b=b, c=c, d=d, e=e, f=f, g=g, h=h, j=j, k=k) class Cylinder(Surface): From 5af686b4fb49b6a9b24cca787a33f3091d26e5e2 Mon Sep 17 00:00:00 2001 From: Ethan Peterson Date: Mon, 10 Feb 2020 17:36:14 -0500 Subject: [PATCH 02/33] working on quadric refactoring --- openmc/surface.py | 382 ++++++++++++++++++---------------------------- 1 file changed, 151 insertions(+), 231 deletions(-) diff --git a/openmc/surface.py b/openmc/surface.py index d64ca242fd..de1ec4cc2d 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -1,6 +1,6 @@ from abc import ABCMeta, abstractmethod from collections import OrderedDict -from abc.collections import Iterable +from collections.abc import Iterable from copy import deepcopy from numbers import Real, Integral from xml.etree import ElementTree as ET @@ -283,28 +283,30 @@ _SURFACE_CLASSES = Surface.get_subclass_map() class PlaneMeta(metaclass=ABCMeta): """A Plane Meta class for all operations on order 1 surfaces""" +# def __new__(cls, *args, **kwargs): +# for key in cls._coeff_keys, kwargs.pop(key) to get arguments for class +# """Simplify this plane if possible to an XPlane, YPlane, or ZPlane""" +# pass +# +# if cls is Plane: +# if np.all(np.isclose((self.b, self.c), 0., atol=atol)): +# x0 = self.d / self.a +# return XPlane(x0=x0, boundary_type=self.boundary_type, +# name=self.name, surface_id=self.id) +# +# elif np.all(np.isclose((self.a, self.c), 0., atol=atol)): +# y0 = self.d / self.b +# return YPlane(y0=y0, boundary_type=self.boundary_type, +# name=self.name, surface_id=self.id) +# +# elif np.all(np.isclose((self.a, self.b), 0., atol=atol)): +# z0 = self.d / self.c +# return ZPlane(z0=z0, boundary_type=self.boundary_type, +# name=self.name, surface_id=self.id) - def __new__(cls, *args, **kwargs): - """Simplify this plane if possible to an XPlane, YPlane, or ZPlane""" - if cls is Plane: - (a,b,c,d,bt, name, surfidlen(args) - if np.all(np.isclose((self.b, self.c), 0., atol=atol)): - x0 = self.d / self.a - return XPlane(x0=x0, boundary_type=self.boundary_type, - name=self.name, surface_id=self.id) - - elif np.all(np.isclose((self.a, self.c), 0., atol=atol)): - y0 = self.d / self.b - return YPlane(y0=y0, boundary_type=self.boundary_type, - name=self.name, surface_id=self.id) - - elif np.all(np.isclose((self.a, self.b), 0., atol=atol)): - z0 = self.d / self.c - return ZPlane(z0=z0, boundary_type=self.boundary_type, - name=self.name, surface_id=self.id) - - def __init__(self): - pass + def __init__(self, **kwargs): + super().__init__(**kwargs) + self._periodic_surface = None @property def periodic_surface(self): @@ -381,51 +383,6 @@ class PlaneMeta(metaclass=ABCMeta): surf._update_from_base_coeffs(a, b, c, d) return surf - def rotate(self, rotation, frame='lab', clone=False): - """Rotate surface by given Tait-Bryan angles - - Parameters - ---------- - rotation : iterable of float - Intrinsic Tait-Bryan angles in degrees used to rotate the surface - frame : str, one of 'lab' or 'body' - clone : boolean - Whether or not to return a new instance of a Plane or to modify the - coefficients of this plane. - - Returns - ------- - openmc.Plane or None - Rotated surface - - """ - - check_type('surface rotation', rotation, Iterable, Real) - check_length('surface rotation', rotation, 3) - - # Calculate rotation matrix Rmat from angles phi, theta, psi - phi, theta, psi = rotation*(-np.pi/180.) - c3, s3 = np.cos(phi), np.sin(phi) - c2, s2 = np.cos(theta), np.sin(theta) - c1, s1 = np.cos(psi), np.sin(psi) - Rmat = np.array([[c1*c2, c1*s2*s3 - c3*s1, s1*s3 + c1*c3*s2], - [c2*s1, c1*c3 + s1*s2*s3, c3*s1*s2 - c1*s3], - [-s2, c2*s3, c2*c3]]) - - a, b, c, d = self._get_base_coeffs() - bvec = np.array([a, b, c]) - - # Compute new rotated coefficients a, b, c - a, b, c = np.dot(bvec.T, Rmat.T) - - if clone: - surf = self.clone() - else: - surf = self - - surf._update_from_base_coeffs(a, b, c, d) - return surf - def bounding_box(self, side): """Determine an axis-aligned bounding box. @@ -527,15 +484,18 @@ class Plane(PlaneMeta, Surface): _type = 'plane' _coeff_keys = ('a', 'b', 'c', 'd') - def __init__(self, a=1., b=0., c=0., d=0., boundary_type='transmission', - name='', surface_id=None, **kwargs): - super().__init__(surface_id, boundary_type, name=name) - self._periodic_surface = None + def __init__(self, a=1., b=0., c=0., d=0., **kwargs): + # work around until capital letter kwargs are deprecated + oldkwargs = deepcopy(kwargs) + for k in 'ABCD': + kwargs.pop(k, None) + + super().__init__(**kwargs) self.a = a self.b = b self.c = c self.d = d - for k, v in kwargs.items(): + for k, v in oldkwargs.items(): if k in 'ABCD': warn(_WARNING_UPPER.format(type(self).__name__, k.lower(), k), FutureWarning) @@ -662,9 +622,8 @@ class XPlane(PlaneMeta, Surface): _type = 'x-plane' _coeff_keys = ('x0',) - def __init__(self, x0=0., boundary_type='transmission', - name='', surface_id=None): - super().__init__(surface_id=surface_id, boundary_type=boundary_type, name=name) + def __init__(self, x0=0., **kwargs): + super().__init__(**kwargs) self.x0 = x0 @property @@ -736,9 +695,8 @@ class YPlane(PlaneMeta, Surface): _type = 'y-plane' _coeff_keys = ('y0',) - def __init__(self, y0=0., boundary_type='transmission', - name='', surface_id=None): - super().__init__(surface_id=surface_id, boundary_type=boundary_type, name=name) + def __init__(self, y0=0., **kwargs): + super().__init__(**kwargs) self.y0 = y0 @property @@ -811,9 +769,8 @@ class ZPlane(PlaneMeta, Surface): _type = 'z-plane' _coeff_keys = ('z0',) - def __init__(self, z0=0., boundary_type='transmission', - name='', surface_id=None): - super().__init__(surface_id=surface_id, boundary_type=boundary_type, name=name) + def __init__(self, z0=0., **kwargs): + super().__init__(**kwargs) self.z0 = z0 @property @@ -843,148 +800,28 @@ class ZPlane(PlaneMeta, Surface): return (np.array([-np.inf, -np.inf, self.z0]), np.array([np.inf, np.inf, np.inf])) + class QuadricMeta(metaclass=ABCMeta): - """A surface of the form :math:`Ax^2 + By^2 + Cz^2 + Dxy + Eyz + Fxz + Gx + Hy + - Jz + K = 0`. + """A Meta class implementing common functionality for quadric surfaces""" - Parameters - ---------- - a, b, c, d, e, f, g, h, j, k : float, optional - coefficients for the surface. All default to 0. - boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic', 'white'}, optional - Boundary condition that defines the behavior for particles hitting the - surface. Defaults to transmissive boundary condition where particles - freely pass through the surface. - name : str, optional - Name of the surface. If not specified, the name will be the empty string. - surface_id : int, optional - Unique identifier for the surface. If not specified, an identifier will - automatically be assigned. + def __init__(self, **kwargs): + super().__init__(**kwargs) - Attributes - ---------- - a, b, c, d, e, f, g, h, j, k : float - coefficients for the surface - boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic', 'white'} - Boundary condition that defines the behavior for particles hitting the - surface. - coefficients : dict - Dictionary of surface coefficients - id : int - Unique identifier for the surface - name : str - Name of the surface - type : str - Type of the surface + @abstractmethod + def _get_base_coeffs(self): + pass - """ + @abstractmethod + def _update_from_base_coeffs(self): + pass - _type = 'quadric' - _coeff_keys = ('a', 'b', 'c', 'd', 'e', 'f', 'g', 'h', 'j', 'k') + def _neg_bounds(self): + return (np.array([-np.inf, -np.inf, -np.inf]), + np.array([np.inf, np.inf, np.inf])) - def __init__(self, a=0., b=0., c=0., d=0., e=0., f=0., g=0., h=0., j=0., - k=0., boundary_type='transmission', name='', surface_id=None): - super().__init__(surface_id, boundary_type, name=name) - self.a = a - self.b = b - self.c = c - self.d = d - self.e = e - self.f = f - self.g = g - self.h = h - self.j = j - self.k = k - - @property - def a(self): - return self.coefficients['a'] - - @property - def b(self): - return self.coefficients['b'] - - @property - def c(self): - return self.coefficients['c'] - - @property - def d(self): - return self.coefficients['d'] - - @property - def e(self): - return self.coefficients['e'] - - @property - def f(self): - return self.coefficients['f'] - - @property - def g(self): - return self.coefficients['g'] - - @property - def h(self): - return self.coefficients['h'] - - @property - def j(self): - return self.coefficients['j'] - - @property - def k(self): - return self.coefficients['k'] - - @a.setter - def a(self, a): - check_type('a coefficient', a, Real) - self._coefficients['a'] = a - - @b.setter - def b(self, b): - check_type('b coefficient', b, Real) - self._coefficients['b'] = b - - @c.setter - def c(self, c): - check_type('c coefficient', c, Real) - self._coefficients['c'] = c - - @d.setter - def d(self, d): - check_type('d coefficient', d, Real) - self._coefficients['d'] = d - - @e.setter - def e(self, e): - check_type('e coefficient', e, Real) - self._coefficients['e'] = e - - @f.setter - def f(self, f): - check_type('f coefficient', f, Real) - self._coefficients['f'] = f - - @g.setter - def g(self, g): - check_type('g coefficient', g, Real) - self._coefficients['g'] = g - - @h.setter - def h(self, h): - check_type('h coefficient', h, Real) - self._coefficients['h'] = h - - @j.setter - def j(self, j): - check_type('j coefficient', j, Real) - self._coefficients['j'] = j - - @k.setter - def k(self, k): - check_type('k coefficient', k, Real) - self._coefficients['k'] = k + def _pos_bounds(self): + return (np.array([-np.inf, -np.inf, -np.inf]), + np.array([np.inf, np.inf, np.inf])) def evaluate(self, point): """Evaluate the surface equation at a given point. @@ -1003,11 +840,10 @@ class QuadricMeta(metaclass=ABCMeta): """ x, y, z = point - return x*(self.a*x + self.d*y + self.g) + \ - y*(self.b*y + self.e*z + self.h) + \ - z*(self.c*z + self.f*x + self.j) + self.k + a, b, c, d, e, f, g, h, j, k = self._get_base_coeffs() + return x*(a*x + d*y + g) + y*(b*y + e*z + h) + z*(c*z + f*x + j) + k - def translate(self, vector): + def translate(self, vector, clone=False): """Translate surface in given direction Parameters @@ -1022,18 +858,23 @@ class QuadricMeta(metaclass=ABCMeta): """ vx, vy, vz = vector - a, b, c, d, e, f, g, h, j, k = (getattr(self, key) for key in - self._coeff_keys) + a, b, c, d, e, f, g, h, j, k = self._get_base_coeffs() k = (k + vx*vx + vy*vy + vz*vz + d*vx*vy + e*vy*vz + f*vx*vz - g*vx - h*vy - j*vz) g = g - 2*a*vx - d*vy - f*vz h = h - 2*b*vy - d*vx - e*vz j = j - 2*c*vz - e*vy - f*vx - return type(self)(a=a, b=b, c=c, d=d, e=e, f=f, g=g, h=h, j=j, k=k) + if clone: + surf = self.clone() + else: + surf = self + + surf._update_from_base_coeffs((a, b, c, d, e, f, g, h, j, k)) + return surf -class Cylinder(Surface): - """A cylinder whose length is parallel to the x-, y-, or z-axis. +class Cylinder(QuadricMeta, Surface): + """A cylinder Parameters ---------- @@ -1067,22 +908,90 @@ class Cylinder(Surface): Type of the surface """ - def __init__(self, r=1., boundary_type='transmission', - name='', surface_id=None): - super().__init__(surface_id, boundary_type, name=name) + _type = 'cylinder' + _coeff_keys = ('x0', 'y0', 'z0', 'r', 'u', 'v','w') + def __init__(self, x0=0., y0=0., z0=0. r=1., u=0., v=0., w=1., **kwargs): + super().__init__(**kwargs) + self.x0 = x0 + self.y0 = y0 + self.z0 = z0 self.r = r + self.u = u + self.v = v + self.w = w + + @property + def x0(self): + return self.coefficients['x0'] + + @property + def y0(self): + return self.coefficients['y0'] + + @property + def z0(self): + return self.coefficients['z0'] @property def r(self): return self.coefficients['r'] + @property + def u(self): + return self.coefficients['u'] + + @property + def v(self): + return self.coefficients['v'] + + @property + def w(self): + return self.coefficients['w'] + + @x0.setter + def x0(self, x0): + check_type('x0 coefficient', x0, Real) + self._coefficients['x0'] = x0 + + @y0.setter + def y0(self, y0): + check_type('y0 coefficient', y0, Real) + self._coefficients['y0'] = y0 + + @z0.setter + def z0(self, z0): + check_type('z0 coefficient', z0, Real) + self._coefficients['z0'] = z0 + @r.setter def r(self, r): check_type('r coefficient', r, Real) self._coefficients['r'] = r + @u.setter + def u(self, u): + check_type('u coefficient', u, Real) + self._coefficients['u'] = u -class XCylinder(Cylinder): + @v.setter + def v(self, v): + check_type('v coefficient', v, Real) + self._coefficients['v'] = v + + @w.setter + def w(self, w): + check_type('w coefficient', w, Real) + self._coefficients['w'] = w + + def _get_base_coeffs(self): + """Return generalized coefficients for a cylinder""" + return (0., 0., 1., self.z0) + + def _update_from_base_coeffs(self, coeffs): + a, b, c, d, e, f, g, h, j, k = coeffs + + +class XCylinder(QuadricMeta, Surface): """An infinite cylinder whose length is parallel to the x-axis of the form :math:`(y - y_0)^2 + (z - z_0)^2 = r^2`. @@ -1131,12 +1040,17 @@ class XCylinder(Cylinder): def __init__(self, y0=0., z0=0., r=1., boundary_type='transmission', name='', surface_id=None, *, R=None): if R is not None: - warn(_WARNING_UPPER.format(type(self).__name__, 'r', 'R'), FutureWarning) + warn(_WARNING_UPPER.format(type(self).__name__, 'r', 'R'), + FutureWarning) r = R super().__init__(r, boundary_type, name, surface_id) self.y0 = y0 self.z0 = z0 + @property + def x0(self): + return self.coefficients['x0'] + @property def y0(self): return self.coefficients['y0'] @@ -1145,6 +1059,11 @@ class XCylinder(Cylinder): def z0(self): return self.coefficients['z0'] + @x0.setter + def x0(self, x0): + check_type('x0 coefficient', x0, Real) + self._coefficients['x0'] = x0 + @y0.setter def y0(self, y0): check_type('y0 coefficient', y0, Real) @@ -1155,6 +1074,7 @@ class XCylinder(Cylinder): check_type('z0 coefficient', z0, Real) self._coefficients['z0'] = z0 + def bounding_box(self, side): """Determine an axis-aligned bounding box. From 64ee9b177f286c5948ad490bdc1d4d8ebf305501 Mon Sep 17 00:00:00 2001 From: Ethan Peterson Date: Mon, 10 Feb 2020 20:13:49 -0500 Subject: [PATCH 03/33] finished skeleton framework for refactoring surfaces --- openmc/surface.py | 650 +++++++++++++++++++++++----------------------- 1 file changed, 323 insertions(+), 327 deletions(-) diff --git a/openmc/surface.py b/openmc/surface.py index de1ec4cc2d..6124494fe0 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -272,10 +272,11 @@ class Surface(IDManagerMixin, metaclass=ABCMeta): surf_type = group['type'][()].decode() bc = group['boundary_type'][()].decode() coeffs = group['coefficients'][...] + kwargs = {'boundary_type': bc, 'name': name, 'surface_id': surface_id} cls = _SURFACE_CLASSES[surf_type] - return cls(*coeffs, bc, name, surface_id) + return cls(*coeffs, **kwargs) _SURFACE_CLASSES = Surface.get_subclass_map() @@ -850,10 +851,12 @@ class QuadricMeta(metaclass=ABCMeta): ---------- vector : iterable of float Direction in which surface should be translated + clone : bool + Whether to return a clone of the Surface or the Surface itself Returns ------- - openmc.Quadric + openmc.QuadricMeta Translated surface """ @@ -864,12 +867,14 @@ class QuadricMeta(metaclass=ABCMeta): g = g - 2*a*vx - d*vy - f*vz h = h - 2*b*vy - d*vx - e*vz j = j - 2*c*vz - e*vy - f*vx + if clone: surf = self.clone() else: surf = self surf._update_from_base_coeffs((a, b, c, d, e, f, g, h, j, k)) + return surf @@ -878,8 +883,23 @@ class Cylinder(QuadricMeta, Surface): Parameters ---------- + x0 : float, optional + x-coordinate for the origin of the Cylinder. Defaults to 0 + y0 : float, optional + y-coordinate for the origin of the Cylinder. Defaults to 0 + z0 : float, optional + z-coordinate for the origin of the Cylinder. Defaults to 0 r : float, optional Radius of the cylinder. Defaults to 1. + u : float, optional + x-component of the vector representing the axis of the cylinder. + Defaults to 0. + v : float, optional + y-component of the vector representing the axis of the cylinder. + Defaults to 0. + w : float, optional + z-component of the vector representing the axis of the cylinder. + Defaults to 1. boundary_type : {'transmission, 'vacuum', 'reflective', 'white'}, optional Boundary condition that defines the behavior for particles hitting the surface. Defaults to transmissive boundary condition where particles @@ -893,8 +913,20 @@ class Cylinder(QuadricMeta, Surface): Attributes ---------- + x0 : float + x-coordinate for the origin of the Cylinder + y0 : float + y-coordinate for the origin of the Cylinder + z0 : float + z-coordinate for the origin of the Cylinder r : float Radius of the cylinder + u : float + x-component of the vector representing the axis of the cylinder + v : float + y-component of the vector representing the axis of the cylinder + w : float + z-component of the vector representing the axis of the cylinder boundary_type : {'transmission, 'vacuum', 'reflective', 'white'} Boundary condition that defines the behavior for particles hitting the surface. @@ -998,11 +1030,11 @@ class XCylinder(QuadricMeta, Surface): Parameters ---------- y0 : float, optional - y-coordinate of the center of the cylinder. Defaults to 0. + y-coordinate for the origin of the Cylinder. Defaults to 0 z0 : float, optional - z-coordinate of the center of the cylinder. Defaults to 0. + z-coordinate for the origin of the Cylinder. Defaults to 0 r : float, optional - Radius of the cylinder. Defaults to 0. + Radius of the cylinder. Defaults to 1. boundary_type : {'transmission, 'vacuum', 'reflective', 'white'}, optional Boundary condition that defines the behavior for particles hitting the surface. Defaults to transmissive boundary condition where particles @@ -1017,9 +1049,11 @@ class XCylinder(QuadricMeta, Surface): Attributes ---------- y0 : float - y-coordinate of the center of the cylinder + y-coordinate for the origin of the Cylinder z0 : float - z-coordinate of the center of the cylinder + z-coordinate for the origin of the Cylinder + r : float + Radius of the cylinder boundary_type : {'transmission, 'vacuum', 'reflective', 'white'} Boundary condition that defines the behavior for particles hitting the surface. @@ -1037,15 +1071,16 @@ class XCylinder(QuadricMeta, Surface): _type = 'x-cylinder' _coeff_keys = ('y0', 'z0', 'r') - def __init__(self, y0=0., z0=0., r=1., boundary_type='transmission', - name='', surface_id=None, *, R=None): + def __init__(self, y0=0., z0=0., r=1., **kwargs): + R = kwargs.pop('R', None) if R is not None: warn(_WARNING_UPPER.format(type(self).__name__, 'r', 'R'), FutureWarning) r = R - super().__init__(r, boundary_type, name, surface_id) + super().__init__(**kwargs) self.y0 = y0 self.z0 = z0 + self.r = r @property def x0(self): @@ -1074,6 +1109,12 @@ class XCylinder(QuadricMeta, Surface): check_type('z0 coefficient', z0, Real) self._coefficients['z0'] = z0 + def _get_base_coeffs(self): + """Return generalized coefficients for a cylinder""" + return (0., 0., 1., self.z0) + + def _update_from_base_coeffs(self, coeffs): + a, b, c, d, e, f, g, h, j, k = coeffs def bounding_box(self, side): """Determine an axis-aligned bounding box. @@ -1107,58 +1148,17 @@ class XCylinder(QuadricMeta, Surface): return (np.array([-np.inf, -np.inf, -np.inf]), np.array([np.inf, np.inf, np.inf])) - def evaluate(self, point): - """Evaluate the surface equation at a given point. - Parameters - ---------- - point : 3-tuple of float - The Cartesian coordinates, :math:`(x',y',z')`, at which the surface - equation should be evaluated. - - Returns - ------- - float - :math:`(y' - y_0)^2 + (z' - z_0)^2 - r^2` - - """ - y = point[1] - self.y0 - z = point[2] - self.z0 - return y**2 + z**2 - self.r**2 - - def translate(self, vector): - """Translate surface in given direction - - Parameters - ---------- - vector : iterable of float - Direction in which surface should be translated - - Returns - ------- - openmc.XCylinder - Translated surface - - """ - vx, vy, vz = vector - if vy == 0.0 and vz == 0.0: - return self - else: - y0 = self.y0 + vy - z0 = self.z0 + vz - return type(self)(y0=y0, z0=z0, r=self.r) - - -class YCylinder(Cylinder): +class YCylinder(QuadricMeta): """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`. Parameters ---------- x0 : float, optional - x-coordinate of the center of the cylinder. Defaults to 0. + x-coordinate for the origin of the Cylinder. Defaults to 0 z0 : float, optional - z-coordinate of the center of the cylinder. Defaults to 0. + z-coordinate for the origin of the Cylinder. Defaults to 0 r : float, optional Radius of the cylinder. Defaults to 1. boundary_type : {'transmission, 'vacuum', 'reflective', 'white'}, optional @@ -1175,9 +1175,11 @@ class YCylinder(Cylinder): Attributes ---------- x0 : float - x-coordinate of the center of the cylinder + x-coordinate for the origin of the Cylinder z0 : float - z-coordinate of the center of the cylinder + z-coordinate for the origin of the Cylinder + r : float + Radius of the cylinder boundary_type : {'transmission, 'vacuum', 'reflective', 'white'} Boundary condition that defines the behavior for particles hitting the surface. @@ -1195,14 +1197,15 @@ class YCylinder(Cylinder): _type = 'y-cylinder' _coeff_keys = ('x0', 'z0', 'r') - def __init__(self, x0=0., z0=0., r=1., boundary_type='transmission', - name='', surface_id=None, *, R=None): + def __init__(self, x0=0., z0=0., r=1., **kwargs): + R = kwargs.pop('R', None) if R is not None: warn(_WARNING_UPPER.format(type(self).__name__, 'r', 'R'), FutureWarning) r = R - super().__init__(r, boundary_type, name, surface_id) + super().__init__(**kwargs) self.x0 = x0 self.z0 = z0 + self.r = r @property def x0(self): @@ -1222,6 +1225,13 @@ class YCylinder(Cylinder): check_type('z0 coefficient', z0, Real) self._coefficients['z0'] = z0 + def _get_base_coeffs(self): + """Return generalized coefficients for a cylinder""" + return (0., 0., 1., self.z0) + + def _update_from_base_coeffs(self, coeffs): + a, b, c, d, e, f, g, h, j, k = coeffs + def bounding_box(self, side): """Determine an axis-aligned bounding box. @@ -1254,77 +1264,38 @@ class YCylinder(Cylinder): return (np.array([-np.inf, -np.inf, -np.inf]), np.array([np.inf, np.inf, np.inf])) - def evaluate(self, point): - """Evaluate the surface equation at a given point. - Parameters - ---------- - point : 3-tuple of float - The Cartesian coordinates, :math:`(x',y',z')`, at which the surface - equation should be evaluated. - - Returns - ------- - float - :math:`(x' - x_0)^2 + (z' - z_0)^2 - r^2` - - """ - x = point[0] - self.x0 - z = point[2] - self.z0 - return x**2 + z**2 - self.r**2 - - def translate(self, vector): - """Translate surface in given direction - - Parameters - ---------- - vector : iterable of float - Direction in which surface should be translated - - Returns - ------- - openmc.YCylinder - Translated surface - - """ - vx, vy, vz = vector - if vx == 0.0 and vz == 0.0: - return self - else: - x0 = self.x0 + vx - z0 = self.z0 + vz - return type(self)(x0=x0, z0=z0, r=self.r) - - -class ZCylinder(Cylinder): +class ZCylinder(QuadricMeta, 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`. Parameters ---------- - surface_id : int, optional - Unique identifier for the surface. If not specified, an identifier will - automatically be assigned. + x0 : float, optional + x-coordinate for the origin of the Cylinder. Defaults to 0 + y0 : float, optional + y-coordinate for the origin of the Cylinder. Defaults to 0 + r : float, optional + Radius of the cylinder. Defaults to 1. boundary_type : {'transmission, 'vacuum', 'reflective', 'white'}, optional Boundary condition that defines the behavior for particles hitting the surface. Defaults to transmissive boundary condition where particles freely pass through the surface. - x0 : float, optional - x-coordinate of the center of the cylinder. Defaults to 0. - y0 : float, optional - y-coordinate of the center of the cylinder. Defaults to 0. - r : float, optional - Radius of the cylinder. Defaults to 1. name : str, optional Name of the cylinder. If not specified, the name will be the empty string. + surface_id : int, optional + Unique identifier for the surface. If not specified, an identifier will + automatically be assigned. Attributes ---------- x0 : float - x-coordinate of the center of the cylinder + x-coordinate for the origin of the Cylinder y0 : float - y-coordinate of the center of the cylinder + y-coordinate for the origin of the Cylinder + r : float + Radius of the cylinder boundary_type : {'transmission, 'vacuum', 'reflective', 'white'} Boundary condition that defines the behavior for particles hitting the surface. @@ -1342,14 +1313,15 @@ class ZCylinder(Cylinder): _type = 'z-cylinder' _coeff_keys = ('x0', 'y0', 'r') - def __init__(self, x0=0., y0=0., r=1., boundary_type='transmission', - name='', surface_id=None, *, R=None): + def __init__(self, x0=0., y0=0., r=1., **kwargs): + R = kwargs.pop('R', None) if R is not None: warn(_WARNING_UPPER.format(type(self).__name__, 'r', 'R'), FutureWarning) r = R - super().__init__(r, boundary_type, name, surface_id) + super().__init__(**kwargs) self.x0 = x0 self.y0 = y0 + self.r = r @property def x0(self): @@ -1369,6 +1341,13 @@ class ZCylinder(Cylinder): check_type('y0 coefficient', y0, Real) self._coefficients['y0'] = y0 + def _get_base_coeffs(self): + """Return generalized coefficients for a cylinder""" + return (0., 0., 1., self.z0) + + def _update_from_base_coeffs(self, coeffs): + a, b, c, d, e, f, g, h, j, k = coeffs + def bounding_box(self, side): """Determine an axis-aligned bounding box. @@ -1401,49 +1380,8 @@ class ZCylinder(Cylinder): return (np.array([-np.inf, -np.inf, -np.inf]), np.array([np.inf, np.inf, np.inf])) - def evaluate(self, point): - """Evaluate the surface equation at a given point. - Parameters - ---------- - point : 3-tuple of float - The Cartesian coordinates, :math:`(x',y',z')`, at which the surface - equation should be evaluated. - - Returns - ------- - float - :math:`(x' - x_0)^2 + (y' - y_0)^2 - r^2` - - """ - x = point[0] - self.x0 - y = point[1] - self.y0 - return x**2 + y**2 - self.r**2 - - def translate(self, vector): - """Translate surface in given direction - - Parameters - ---------- - vector : iterable of float - Direction in which surface should be translated - - Returns - ------- - openmc.ZCylinder - Translated surface - - """ - vx, vy, vz = vector - if vx == 0.0 and vy == 0.0: - return self - else: - x0 = self.x0 + vx - y0 = self.y0 + vy - return type(self)(x0=x0, y0=y0, r=self.r) - - -class Sphere(Surface): +class Sphere(QuadricMeta, Surface): """A sphere of the form :math:`(x - x_0)^2 + (y - y_0)^2 + (z - z_0)^2 = r^2`. Parameters @@ -1493,12 +1431,12 @@ class Sphere(Surface): _type = 'sphere' _coeff_keys = ('x0', 'y0', 'z0', 'r') - def __init__(self, x0=0., y0=0., z0=0., r=1., boundary_type='transmission', - name='', surface_id=None, *, R=None): + def __init__(self, x0=0., y0=0., z0=0., r=1., **kwargs): + R = kwargs.pop('R', None) if R is not None: warn(_WARNING_UPPER.format(type(self).__name__, 'r', 'R'), FutureWarning) r = R - super().__init__(surface_id, boundary_type, name=name) + super().__init__(**kwargs) self.x0 = x0 self.y0 = y0 self.z0 = z0 @@ -1540,6 +1478,13 @@ class Sphere(Surface): check_type('r coefficient', r, Real) self._coefficients['r'] = r + def _get_base_coeffs(self): + """Return generalized coefficients for a cylinder""" + return (0., 0., 1., self.z0) + + def _update_from_base_coeffs(self, coeffs): + a, b, c, d, e, f, g, h, j, k = coeffs + def bounding_box(self, side): """Determine an axis-aligned bounding box. @@ -1573,51 +1518,8 @@ class Sphere(Surface): return (np.array([-np.inf, -np.inf, -np.inf]), np.array([np.inf, np.inf, np.inf])) - def evaluate(self, point): - """Evaluate the surface equation at a given point. - Parameters - ---------- - point : 3-tuple of float - The Cartesian coordinates, :math:`(x',y',z')`, at which the surface - equation should be evaluated. - - Returns - ------- - float - :math:`(x' - x_0)^2 + (y' - y_0)^2 + (z' - z_0)^2 - r^2` - - """ - x = point[0] - self.x0 - y = point[1] - self.y0 - z = point[2] - self.z0 - return x**2 + y**2 + z**2 - self.r**2 - - def translate(self, vector): - """Translate surface in given direction - - Parameters - ---------- - vector : iterable of float - Direction in which surface should be translated - - Returns - ------- - openmc.Sphere - Translated surface - - """ - vx, vy, vz = vector - if vx == 0.0 and vy == 0.0 and vz == 0.0: - return self - else: - x0 = self.x0 + vx - y0 = self.y0 + vy - z0 = self.z0 + vz - return type(self)(x0=x0, y0=y0, z0=z0, r=self.r) - - -class Cone(Surface): +class Cone(QuadricMeta, Surface): """A conical surface parallel to the x-, y-, or z-axis. Parameters @@ -1630,6 +1532,15 @@ class Cone(Surface): z-coordinate of the apex. Defaults to 0. r2 : float, optional Parameter related to the aperature. Defaults to 1. + u : float, optional + x-component of the vector representing the axis of the cone. + Defaults to 0. + v : float, optional + y-component of the vector representing the axis of the cone. + Defaults to 0. + w : float, optional + z-component of the vector representing the axis of the cone. + Defaults to 1. surface_id : int, optional Unique identifier for the surface. If not specified, an identifier will automatically be assigned. @@ -1650,6 +1561,12 @@ class Cone(Surface): z-coordinate of the apex r2 : float Parameter related to the aperature + u : float + x-component of the vector representing the axis of the cone. + v : float + y-component of the vector representing the axis of the cone. + w : float + z-component of the vector representing the axis of the cone. boundary_type : {'transmission, 'vacuum', 'reflective', 'white'} Boundary condition that defines the behavior for particles hitting the surface. @@ -1664,18 +1581,23 @@ class Cone(Surface): """ - _coeff_keys = ('x0', 'y0', 'z0', 'r2') + _type = 'cylinder' + _coeff_keys = ('x0', 'y0', 'z0', 'r2', 'u', 'v', 'w') - def __init__(self, x0=0., y0=0., z0=0., r2=1., boundary_type='transmission', - name='', surface_id=None, *, R2=None): + def __init__(self, x0=0., y0=0., z0=0., r2=1., u=0., v=0., w=1., **kwargs): + R2 = kwargs.pop('R2', None) if R2 is not None: - warn(_WARNING_UPPER.format(type(self).__name__, 'r2', 'R2'), FutureWarning) + warn(_WARNING_UPPER.format(type(self).__name__, 'r2', 'R2'), + FutureWarning) r2 = R2 - super().__init__(surface_id, boundary_type, name=name) + super().__init__(**kwargs) self.x0 = x0 self.y0 = y0 self.z0 = z0 self.r2 = r2 + self.u = u + self.v = v + self.w = w @property def x0(self): @@ -1693,6 +1615,18 @@ class Cone(Surface): def r2(self): return self.coefficients['r2'] + @property + def u(self): + return self.coefficients['u'] + + @property + def v(self): + return self.coefficients['v'] + + @property + def w(self): + return self.coefficients['w'] + @x0.setter def x0(self, x0): check_type('x0 coefficient', x0, Real) @@ -1713,31 +1647,30 @@ class Cone(Surface): check_type('r^2 coefficient', r2, Real) self._coefficients['r2'] = r2 - def translate(self, vector): - """Translate surface in given direction + @u.setter + def u(self, u): + check_type('u coefficient', u, Real) + self._coefficients['u'] = u - Parameters - ---------- - vector : iterable of float - Direction in which surface should be translated + @v.setter + def v(self, v): + check_type('v coefficient', v, Real) + self._coefficients['v'] = v - Returns - ------- - openmc.Cone - Translated surface + @w.setter + def w(self, w): + check_type('w coefficient', w, Real) + self._coefficients['w'] = w - """ - vx, vy, vz = vector - if vx == 0.0 and vy == 0.0 and vz == 0.0: - return self - else: - x0 = self.x0 + vx - y0 = self.y0 + vy - z0 = self.z0 + vz - return type(self)(x0=x0, y0=y0, z0=z0, r2=self.r2) + def _get_base_coeffs(self): + """Return generalized coefficients for a cylinder""" + return (0., 0., 1., self.z0) + + def _update_from_base_coeffs(self, coeffs): + a, b, c, d, e, f, g, h, j, k = coeffs -class XCone(Cone): +class XCone(QuadricMeta, Surface): """A cone parallel to the x-axis of the form :math:`(y - y_0)^2 + (z - z_0)^2 = r^2 (x - x_0)^2`. @@ -1786,29 +1719,65 @@ class XCone(Cone): """ _type = 'x-cone' + _coeff_keys = ('x0', 'y0', 'z0', 'r2') - def evaluate(self, point): - """Evaluate the surface equation at a given point. + def __init__(self, x0=0., y0=0., z0=0., r2=1., **kwargs): + R2 = kwargs.pop('R2', None) + if R2 is not None: + warn(_WARNING_UPPER.format(type(self).__name__, 'r2', 'R2'), + FutureWarning) + r2 = R2 + super().__init__(**kwargs) + self.x0 = x0 + self.y0 = y0 + self.z0 = z0 + self.r2 = r2 - Parameters - ---------- - point : 3-tuple of float - The Cartesian coordinates, :math:`(x',y',z')`, at which the surface - equation should be evaluated. + @property + def x0(self): + return self.coefficients['x0'] - Returns - ------- - float - :math:`(y' - y_0)^2 + (z' - z_0)^2 - r^2(x' - x_0)^2` + @property + def y0(self): + return self.coefficients['y0'] - """ - x = point[0] - self.x0 - y = point[1] - self.y0 - z = point[2] - self.z0 - return y**2 + z**2 - self.r2*x**2 + @property + def z0(self): + return self.coefficients['z0'] + + @property + def r2(self): + return self.coefficients['r2'] + + @x0.setter + def x0(self, x0): + check_type('x0 coefficient', x0, Real) + self._coefficients['x0'] = x0 + + @y0.setter + def y0(self, y0): + check_type('y0 coefficient', y0, Real) + self._coefficients['y0'] = y0 + + @z0.setter + def z0(self, z0): + check_type('z0 coefficient', z0, Real) + self._coefficients['z0'] = z0 + + @r2.setter + def r2(self, r2): + check_type('r^2 coefficient', r2, Real) + self._coefficients['r2'] = r2 + + def _get_base_coeffs(self): + """Return generalized coefficients for a cylinder""" + return (0., 0., 1., self.z0) + + def _update_from_base_coeffs(self, coeffs): + a, b, c, d, e, f, g, h, j, k = coeffs -class YCone(Cone): +class YCone(QuadricMeta, 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`. @@ -1857,29 +1826,65 @@ class YCone(Cone): """ _type = 'y-cone' + _coeff_keys = ('x0', 'y0', 'z0', 'r2') - def evaluate(self, point): - """Evaluate the surface equation at a given point. + def __init__(self, x0=0., y0=0., z0=0., r2=1., **kwargs): + R2 = kwargs.pop('R2', None) + if R2 is not None: + warn(_WARNING_UPPER.format(type(self).__name__, 'r2', 'R2'), + FutureWarning) + r2 = R2 + super().__init__(**kwargs) + self.x0 = x0 + self.y0 = y0 + self.z0 = z0 + self.r2 = r2 - Parameters - ---------- - point : 3-tuple of float - The Cartesian coordinates, :math:`(x',y',z')`, at which the surface - equation should be evaluated. + @property + def x0(self): + return self.coefficients['x0'] - Returns - ------- - float - :math:`(x' - x_0)^2 + (z' - z_0)^2 - r^2(y' - y_0)^2` + @property + def y0(self): + return self.coefficients['y0'] - """ - x = point[0] - self.x0 - y = point[1] - self.y0 - z = point[2] - self.z0 - return x**2 + z**2 - self.r2*y**2 + @property + def z0(self): + return self.coefficients['z0'] + + @property + def r2(self): + return self.coefficients['r2'] + + @x0.setter + def x0(self, x0): + check_type('x0 coefficient', x0, Real) + self._coefficients['x0'] = x0 + + @y0.setter + def y0(self, y0): + check_type('y0 coefficient', y0, Real) + self._coefficients['y0'] = y0 + + @z0.setter + def z0(self, z0): + check_type('z0 coefficient', z0, Real) + self._coefficients['z0'] = z0 + + @r2.setter + def r2(self, r2): + check_type('r^2 coefficient', r2, Real) + self._coefficients['r2'] = r2 + + def _get_base_coeffs(self): + """Return generalized coefficients for a cylinder""" + return (0., 0., 1., self.z0) + + def _update_from_base_coeffs(self, coeffs): + a, b, c, d, e, f, g, h, j, k = coeffs -class ZCone(Cone): +class ZCone(QuadricMeta, 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`. @@ -1928,29 +1933,65 @@ class ZCone(Cone): """ _type = 'z-cone' + _coeff_keys = ('x0', 'y0', 'z0', 'r2') - def evaluate(self, point): - """Evaluate the surface equation at a given point. + def __init__(self, x0=0., y0=0., z0=0., r2=1., **kwargs): + R2 = kwargs.pop('R2', None) + if R2 is not None: + warn(_WARNING_UPPER.format(type(self).__name__, 'r2', 'R2'), + FutureWarning) + r2 = R2 + super().__init__(**kwargs) + self.x0 = x0 + self.y0 = y0 + self.z0 = z0 + self.r2 = r2 - Parameters - ---------- - point : 3-tuple of float - The Cartesian coordinates, :math:`(x',y',z')`, at which the surface - equation should be evaluated. + @property + def x0(self): + return self.coefficients['x0'] - Returns - ------- - float - :math:`(x' - x_0)^2 + (y' - y_0)^2 - r^2(z' - z_0)^2` + @property + def y0(self): + return self.coefficients['y0'] - """ - x = point[0] - self.x0 - y = point[1] - self.y0 - z = point[2] - self.z0 - return x**2 + y**2 - self.r2*z**2 + @property + def z0(self): + return self.coefficients['z0'] + + @property + def r2(self): + return self.coefficients['r2'] + + @x0.setter + def x0(self, x0): + check_type('x0 coefficient', x0, Real) + self._coefficients['x0'] = x0 + + @y0.setter + def y0(self, y0): + check_type('y0 coefficient', y0, Real) + self._coefficients['y0'] = y0 + + @z0.setter + def z0(self, z0): + check_type('z0 coefficient', z0, Real) + self._coefficients['z0'] = z0 + + @r2.setter + def r2(self, r2): + check_type('r^2 coefficient', r2, Real) + self._coefficients['r2'] = r2 + + def _get_base_coeffs(self): + """Return generalized coefficients for a cylinder""" + return (0., 0., 1., self.z0) + + def _update_from_base_coeffs(self, coeffs): + a, b, c, d, e, f, g, h, j, k = coeffs -class Quadric(Surface): +class Quadric(QuadricMeta, Surface): """A surface of the form :math:`Ax^2 + By^2 + Cz^2 + Dxy + Eyz + Fxz + Gx + Hy + Jz + K = 0`. @@ -1990,8 +2031,8 @@ class Quadric(Surface): _coeff_keys = ('a', 'b', 'c', 'd', 'e', 'f', 'g', 'h', 'j', 'k') def __init__(self, a=0., b=0., c=0., d=0., e=0., f=0., g=0., h=0., j=0., - k=0., boundary_type='transmission', name='', surface_id=None): - super().__init__(surface_id, boundary_type, name=name) + k=0., **kwargs): + super().__init__(**kwargs) self.a = a self.b = b self.c = c @@ -2093,51 +2134,6 @@ class Quadric(Surface): check_type('k coefficient', k, Real) self._coefficients['k'] = k - def evaluate(self, point): - """Evaluate the surface equation at a given point. - - Parameters - ---------- - point : 3-tuple of float - The Cartesian coordinates, :math:`(x',y',z')`, at which the surface - equation should be evaluated. - - Returns - ------- - float - :math:`Ax'^2 + By'^2 + Cz'^2 + Dx'y' + Ey'z' + Fx'z' + Gx' + Hy' + - Jz' + K = 0` - - """ - x, y, z = point - return x*(self.a*x + self.d*y + self.g) + \ - y*(self.b*y + self.e*z + self.h) + \ - z*(self.c*z + self.f*x + self.j) + self.k - - def translate(self, vector): - """Translate surface in given direction - - Parameters - ---------- - vector : iterable of float - Direction in which surface should be translated - - Returns - ------- - openmc.Quadric - Translated surface - - """ - vx, vy, vz = vector - a, b, c, d, e, f, g, h, j, k = (getattr(self, key) for key in - self._coeff_keys) - k = (k + vx*vx + vy*vy + vz*vz + d*vx*vy + e*vy*vz + f*vx*vz - - g*vx - h*vy - j*vz) - g = g - 2*a*vx - d*vy - f*vz - h = h - 2*b*vy - d*vx - e*vz - j = j - 2*c*vz - e*vy - f*vx - return type(self)(a=a, b=b, c=c, d=d, e=e, f=f, g=g, h=h, j=j, k=k) - class Halfspace(Region): """A positive or negative half-space region. From 9740f91ffc51f58f6ccb321ca9bd0cf65fe75138 Mon Sep 17 00:00:00 2001 From: Ethan Peterson Date: Mon, 10 Feb 2020 20:59:34 -0500 Subject: [PATCH 04/33] working out some _get_base_coeff methods --- openmc/surface.py | 74 +++++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 69 insertions(+), 5 deletions(-) diff --git a/openmc/surface.py b/openmc/surface.py index 6124494fe0..ddee08aad8 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -641,7 +641,7 @@ class XPlane(PlaneMeta, Surface): def _update_from_base_coeffs(self, coeffs): a, b, c, d = coeffs - self.x0 = d + self.x0 = d / a def _neg_bounds(self): """Return the lower and upper bounds of the negative half space""" @@ -715,7 +715,7 @@ class YPlane(PlaneMeta, Surface): def _update_from_base_coeffs(self, coeffs): a, b, c, d = coeffs - self.y0 = d + self.y0 = d / b def _neg_bounds(self): """Return the lower and upper bounds of the negative half space""" @@ -789,7 +789,7 @@ class ZPlane(PlaneMeta, Surface): def _update_from_base_coeffs(self, coeffs): a, b, c, d = coeffs - self.z0 = d + self.z0 = d / c def _neg_bounds(self): """Return the lower and upper bounds of the negative half space""" @@ -879,7 +879,9 @@ class QuadricMeta(metaclass=ABCMeta): class Cylinder(QuadricMeta, Surface): - """A cylinder + """A cylinder with radius r, centered on the point (x0, y0, z0) with an + axis specified by the line through points (x0, y0, z0) and (x0+u, y0+v, + z0+w) Parameters ---------- @@ -1017,11 +1019,73 @@ class Cylinder(QuadricMeta, Surface): def _get_base_coeffs(self): """Return generalized coefficients for a cylinder""" - return (0., 0., 1., self.z0) + x0, y0, z0 = self.x0, self.y0, self.z0 + x1, y1, z1 = self.x0 + self.u, self.y0 + self.v, self.z0 + self.w + dx, dy, dz = x1 - x0, y1 - y0, z1 - z0 + + # Set coefficients for Quadric surface that represents a cylinder of + # radius r whose axis is the line defined by p1 and p2 + a = dy**2 + dz**2 + b = dx**2 + dz**2 + c = dx**2 + dy**2 + d = -2*dx*dy + e = -2*dy*dz + f = -2*dx*dz + g = -2*((z1 - z0)*(x0*z1 - x1*z0) + (y1 - y0)*(x0*y1-x1*y0)) + h = 2*((x1 - x0)*(x0*y1 - x1*y0) - (z1 - z0)*(y0*z1 - y1*z0)) + j = 2*((x1 - x0)*(x0*z1 - x1*z0) + (y1 - y0)*(y0*z1 - y1*z0)) + k = (y0*z1 - y1*z0)**2 + (x0*z1 - x1*z0)**2 + (x0*y1 - x1*y0)**2 \ + - r**2*(dx**2 + dy**2 + dz**2) + + return (a, b, c, d, e, f, g, h, j, k) def _update_from_base_coeffs(self, coeffs): a, b, c, d, e, f, g, h, j, k = coeffs + @classmethod + def from_points(cls, p1, p2, r=1., **kwargs): + """Return a cylinder given points that define the axis and a radius. + + Parameters + ---------- + p1, p2 : 3-tuples + Points that pass through the plane + r : float, optional + Radius of the cylinder. Defaults to 1. + kwargs : dict + Keyword arguments passed to the :class:`Quadric` constructor + + Returns + ------- + Cylinder + Cylinder that has an axis through the points p1 and p2 + + """ + # Convert to numpy arrays + p1 = np.asarray(p1) + p2 = np.asarray(p2) + x1, y1, z1 = p1 + x2, y2, z2 = p2 + dx, dy, dz = p2 - p1 + + # Set coefficients for Quadric surface that represents a cylinder of + # radius r whose axis is the line defined by p1 and p2 + a = dy**2 + dz**2 + b = dx**2 + dz**2 + c = dx**2 + dy**2 + d = -2*dx*dy + e = -2*dy*dz + f = -2*dx*dz + g = -2*((z2 - z1)*(x1*z2 - x2*z1) + (y2 - y1)*(x1*y2-x2*y1)) + h = 2*((x2 - x1)*(x1*y2 - x2*y1) - (z2 - z1)*(y1*z2 - y2*z1)) + j = 2*((x2 - x1)*(x1*z2 - x2*z1) + (y2 - y1)*(y1*z2 - y2*z1)) + k = (y1*z2 - y2*z1)**2 + (x1*z2 - x2*z1)**2 + (x1*y2 - x2*y1)**2 \ + - r**2*(dx**2 + dy**2 + dz**2) + + cyl = cls(x0=x1, y0=y1, z0=z1, r=r, **kwargs) + cyl._update_from_base_coeffs((a, b, c, d, e, f, g, h, j, k)) + return cyl + class XCylinder(QuadricMeta, Surface): """An infinite cylinder whose length is parallel to the x-axis of the form From 542643f91aebbc54c01093ce8aca7809d38e2a69 Mon Sep 17 00:00:00 2001 From: Ethan Peterson Date: Mon, 10 Feb 2020 21:03:35 -0500 Subject: [PATCH 05/33] fixed typos in translate and __init__ --- openmc/surface.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/openmc/surface.py b/openmc/surface.py index ddee08aad8..aacaf45b4a 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -381,7 +381,7 @@ class PlaneMeta(metaclass=ABCMeta): else: surf = self - surf._update_from_base_coeffs(a, b, c, d) + surf._update_from_base_coeffs((a, b, c, d)) return surf def bounding_box(self, side): @@ -944,7 +944,7 @@ class Cylinder(QuadricMeta, Surface): """ _type = 'cylinder' _coeff_keys = ('x0', 'y0', 'z0', 'r', 'u', 'v','w') - def __init__(self, x0=0., y0=0., z0=0. r=1., u=0., v=0., w=1., **kwargs): + def __init__(self, x0=0., y0=0., z0=0., r=1., u=0., v=0., w=1., **kwargs): super().__init__(**kwargs) self.x0 = x0 self.y0 = y0 From ab465084b60c9fa3170e534d3c362d66b599b369 Mon Sep 17 00:00:00 2001 From: Ethan Peterson Date: Tue, 11 Feb 2020 11:25:51 -0500 Subject: [PATCH 06/33] still working through coefficient transformations, done with (X,Y,Z)Cylinders and Cones, Planes, Quadric --- openmc/surface.py | 626 +++++++++++++++++++++++++++++++++++----------- 1 file changed, 479 insertions(+), 147 deletions(-) diff --git a/openmc/surface.py b/openmc/surface.py index aacaf45b4a..8b05ce8c53 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -282,28 +282,8 @@ class Surface(IDManagerMixin, metaclass=ABCMeta): _SURFACE_CLASSES = Surface.get_subclass_map() -class PlaneMeta(metaclass=ABCMeta): +class PlaneMixin(metaclass=ABCMeta): """A Plane Meta class for all operations on order 1 surfaces""" -# def __new__(cls, *args, **kwargs): -# for key in cls._coeff_keys, kwargs.pop(key) to get arguments for class -# """Simplify this plane if possible to an XPlane, YPlane, or ZPlane""" -# pass -# -# if cls is Plane: -# if np.all(np.isclose((self.b, self.c), 0., atol=atol)): -# x0 = self.d / self.a -# return XPlane(x0=x0, boundary_type=self.boundary_type, -# name=self.name, surface_id=self.id) -# -# elif np.all(np.isclose((self.a, self.c), 0., atol=atol)): -# y0 = self.d / self.b -# return YPlane(y0=y0, boundary_type=self.boundary_type, -# name=self.name, surface_id=self.id) -# -# elif np.all(np.isclose((self.a, self.b), 0., atol=atol)): -# z0 = self.d / self.c -# return ZPlane(z0=z0, boundary_type=self.boundary_type, -# name=self.name, surface_id=self.id) def __init__(self, **kwargs): super().__init__(**kwargs) @@ -321,17 +301,32 @@ class PlaneMeta(metaclass=ABCMeta): @abstractmethod def _get_base_coeffs(self): + """Return coefficients a, b, c, d representing a general plane of the + form :math:`ax + by + cz = d`. + + """ pass @abstractmethod - def _update_from_base_coeffs(self): + def _update_from_base_coeffs(self, coeffs): + """Update the current plane from coefficients representing a general + plane of the form :math:`ax + by + cz = d`. + + Parameters + ---------- + coeffs : tuple + Tuple of general coefficients (a, b, c, d) representing the plane + + """ pass def _neg_bounds(self): + """Return the lower and upper bounds of the negative half space""" return (np.array([-np.inf, -np.inf, -np.inf]), np.array([np.inf, np.inf, np.inf])) def _pos_bounds(self): + """Return the lower and upper bounds of the positive half space""" return (np.array([-np.inf, -np.inf, -np.inf]), np.array([np.inf, np.inf, np.inf])) @@ -432,7 +427,7 @@ class PlaneMeta(metaclass=ABCMeta): return element -class Plane(PlaneMeta, Surface): +class Plane(PlaneMixin, Surface): """An arbitrary plane of the form :math:`Ax + By + Cz = D`. Parameters @@ -539,14 +534,25 @@ class Plane(PlaneMeta, Surface): self._coefficients['d'] = d def _get_base_coeffs(self): + """Return coefficients a, b, c, d representing a general plane of the + form :math:`ax + by + cz = d`. + + """ return (self.a, self.b, self.c, self.d) def _update_from_base_coeffs(self, coeffs): - a, b, c, d = coeffs - self.a = a - self.b = b - self.c = c - self.d = d + """Update the current plane from coefficients representing a general + plane of the form :math:`ax + by + cz = d`. + + Parameters + ---------- + coeffs : tuple + Tuple of general coefficients (a, b, c, d) representing the plane + + """ + + for c, val in zip(self._coeff_keys, coeffs): + setattr(self, c, val) @classmethod def from_points(cls, p1, p2, p3, **kwargs): @@ -581,7 +587,7 @@ class Plane(PlaneMeta, Surface): return cls(a=a, b=b, c=c, d=d, **kwargs) -class XPlane(PlaneMeta, Surface): +class XPlane(PlaneMixin, Surface): """A plane perpendicular to the x axis of the form :math:`x - x_0 = 0` Parameters @@ -637,9 +643,22 @@ class XPlane(PlaneMeta, Surface): self._coefficients['x0'] = x0 def _get_base_coeffs(self): + """Return coefficients a, b, c, d representing a general plane of the + form :math:`ax + by + cz = d`. + + """ return (1., 0., 0., self.x0) def _update_from_base_coeffs(self, coeffs): + """Update the current plane from coefficients representing a general + plane of the form :math:`ax + by + cz = d`. + + Parameters + ---------- + coeffs : tuple + Tuple of general coefficients (a, b, c, d) representing the plane + + """ a, b, c, d = coeffs self.x0 = d / a @@ -654,7 +673,10 @@ class XPlane(PlaneMeta, Surface): np.array([np.inf, np.inf, np.inf])) -class YPlane(PlaneMeta, Surface): +Plane.register(XPlane) + + +class YPlane(PlaneMixin, Surface): """A plane perpendicular to the y axis of the form :math:`y - y_0 = 0` Parameters @@ -710,10 +732,22 @@ class YPlane(PlaneMeta, Surface): self._coefficients['y0'] = y0 def _get_base_coeffs(self): - """Return generalized coefficients for a plane""" + """Return coefficients a, b, c, d representing a general plane of the + form :math:`ax + by + cz = d`. + + """ return (0., 1., 0., self.y0) def _update_from_base_coeffs(self, coeffs): + """Update the current plane from coefficients representing a general + plane of the form :math:`ax + by + cz = d`. + + Parameters + ---------- + coeffs : tuple + Tuple of general coefficients (a, b, c, d) representing the plane + + """ a, b, c, d = coeffs self.y0 = d / b @@ -728,23 +762,26 @@ class YPlane(PlaneMeta, Surface): np.array([np.inf, np.inf, np.inf])) -class ZPlane(PlaneMeta, Surface): +Plane.register(YPlane) + + +class ZPlane(PlaneMixin, Surface): """A plane perpendicular to the z axis of the form :math:`z - z_0 = 0` Parameters ---------- - surface_id : int, optional - Unique identifier for the surface. If not specified, an identifier will - automatically be assigned. + z0 : float, optional + Location of the plane. Defaults to 0. boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic', 'white'}, optional Boundary condition that defines the behavior for particles hitting the surface. Defaults to transmissive boundary condition where particles freely pass through the surface. Only axis-aligned periodicity is supported, i.e., x-planes can only be paired with x-planes. - z0 : float, optional - Location of the plane. Defaults to 0. name : str, optional Name of the plane. If not specified, the name will be the empty string. + surface_id : int, optional + Unique identifier for the surface. If not specified, an identifier will + automatically be assigned. Attributes ---------- @@ -784,10 +821,22 @@ class ZPlane(PlaneMeta, Surface): self._coefficients['z0'] = z0 def _get_base_coeffs(self): - """Return generalized coefficients for a plane""" + """Return coefficients a, b, c, d, e, f, g, h, j, k representing a + general plane form :math:`ax + by + cz = d`. + + """ return (0., 0., 1., self.z0) def _update_from_base_coeffs(self, coeffs): + """Update the current plane from coefficients representing a general + plane of the form :math:`ax + by + cz = d`. + + Parameters + ---------- + coeffs : tuple + Tuple of general coefficients (a, b, c, d) representing the plane + + """ a, b, c, d = coeffs self.z0 = d / c @@ -802,25 +851,46 @@ class ZPlane(PlaneMeta, Surface): np.array([np.inf, np.inf, np.inf])) -class QuadricMeta(metaclass=ABCMeta): - """A Meta class implementing common functionality for quadric surfaces""" +Plane.register(ZPlane) + + +class QuadricMixin(metaclass=ABCMeta): + """A Mixin class implementing common functionality for quadric surfaces""" def __init__(self, **kwargs): super().__init__(**kwargs) @abstractmethod def _get_base_coeffs(self): + """Return coefficients a, b, c, d, e, f, g, h, j, k representing a + general quadric surface of the form: + :math:`ax^2 + by^2 + cz^2 + dxy + eyz + fxz+ gx + hy + jz +k = 0`. + + """ pass @abstractmethod - def _update_from_base_coeffs(self): + def _update_from_base_coeffs(self, coeffs): + """Update the current surface from coefficients representing a general + quadric surface of the form: + :math:`ax^2 + by^2 + cz^2 + dxy + eyz + fxz+ gx + hy + jz +k = 0`. + + Parameters + ---------- + coeffs : tuple + Tuple of general coefficients (a, b, c, d, e, f, g, h, j, k) + representing the quadric surface + + """ pass def _neg_bounds(self): + """Return the lower and upper bounds of the negative half space""" return (np.array([-np.inf, -np.inf, -np.inf]), np.array([np.inf, np.inf, np.inf])) def _pos_bounds(self): + """Return the lower and upper bounds of the positive half space""" return (np.array([-np.inf, -np.inf, -np.inf]), np.array([np.inf, np.inf, np.inf])) @@ -856,7 +926,7 @@ class QuadricMeta(metaclass=ABCMeta): Returns ------- - openmc.QuadricMeta + openmc.QuadricMixin Translated surface """ @@ -878,7 +948,7 @@ class QuadricMeta(metaclass=ABCMeta): return surf -class Cylinder(QuadricMeta, Surface): +class Cylinder(QuadricMixin, Surface): """A cylinder with radius r, centered on the point (x0, y0, z0) with an axis specified by the line through points (x0, y0, z0) and (x0+u, y0+v, z0+w) @@ -893,13 +963,13 @@ class Cylinder(QuadricMeta, Surface): z-coordinate for the origin of the Cylinder. Defaults to 0 r : float, optional Radius of the cylinder. Defaults to 1. - u : float, optional + dx : float, optional x-component of the vector representing the axis of the cylinder. Defaults to 0. - v : float, optional + dy : float, optional y-component of the vector representing the axis of the cylinder. Defaults to 0. - w : float, optional + dz : float, optional z-component of the vector representing the axis of the cylinder. Defaults to 1. boundary_type : {'transmission, 'vacuum', 'reflective', 'white'}, optional @@ -923,11 +993,11 @@ class Cylinder(QuadricMeta, Surface): z-coordinate for the origin of the Cylinder r : float Radius of the cylinder - u : float + dx : float x-component of the vector representing the axis of the cylinder - v : float + dy : float y-component of the vector representing the axis of the cylinder - w : float + dz : float z-component of the vector representing the axis of the cylinder boundary_type : {'transmission, 'vacuum', 'reflective', 'white'} Boundary condition that defines the behavior for particles hitting the @@ -943,16 +1013,16 @@ class Cylinder(QuadricMeta, Surface): """ _type = 'cylinder' - _coeff_keys = ('x0', 'y0', 'z0', 'r', 'u', 'v','w') - def __init__(self, x0=0., y0=0., z0=0., r=1., u=0., v=0., w=1., **kwargs): + _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): super().__init__(**kwargs) self.x0 = x0 self.y0 = y0 self.z0 = z0 self.r = r - self.u = u - self.v = v - self.w = w + self.dx = dx + self.dy = dy + self.dz = dz @property def x0(self): @@ -971,16 +1041,16 @@ class Cylinder(QuadricMeta, Surface): return self.coefficients['r'] @property - def u(self): - return self.coefficients['u'] + def dx(self): + return self.coefficients['dx'] @property - def v(self): - return self.coefficients['v'] + def dy(self): + return self.coefficients['dy'] @property - def w(self): - return self.coefficients['w'] + def dz(self): + return self.coefficients['dz'] @x0.setter def x0(self, x0): @@ -1002,29 +1072,31 @@ class Cylinder(QuadricMeta, Surface): check_type('r coefficient', r, Real) self._coefficients['r'] = r - @u.setter - def u(self, u): - check_type('u coefficient', u, Real) - self._coefficients['u'] = u + @dx.setter + def dx(self, dx): + check_type('dx coefficient', dx, Real) + self._coefficients['dx'] = dx - @v.setter - def v(self, v): - check_type('v coefficient', v, Real) - self._coefficients['v'] = v + @dy.setter + def dy(self, dy): + check_type('dy coefficient', dy, Real) + self._coefficients['dy'] = dy - @w.setter - def w(self, w): - check_type('w coefficient', w, Real) - self._coefficients['w'] = w + @dz.setter + def dz(self, dz): + check_type('dz coefficient', dz, Real) + self._coefficients['dz'] = dz def _get_base_coeffs(self): - """Return generalized coefficients for a cylinder""" - x0, y0, z0 = self.x0, self.y0, self.z0 - x1, y1, z1 = self.x0 + self.u, self.y0 + self.v, self.z0 + self.w - dx, dy, dz = x1 - x0, y1 - y0, z1 - z0 + """Return coefficients a, b, c, d, e, f, g, h, j, k representing a + general quadric surface of the form: + :math:`ax^2 + by^2 + cz^2 + dxy + eyz + fxz+ gx + hy + jz +k = 0`. + + """ + x0, y0, z0, r = self.x0, self.y0, self.z0, self.r + dx, dy, dz = self.dx, self.dy, self.dz + x1, y1, z1 = x0 + dx, y0 + dy, z0 + dz - # Set coefficients for Quadric surface that represents a cylinder of - # radius r whose axis is the line defined by p1 and p2 a = dy**2 + dz**2 b = dx**2 + dz**2 c = dx**2 + dy**2 @@ -1040,6 +1112,17 @@ class Cylinder(QuadricMeta, Surface): return (a, b, c, d, e, f, g, h, j, k) def _update_from_base_coeffs(self, coeffs): + """Update the current surface from coefficients representing a general + quadric surface of the form: + :math:`ax^2 + by^2 + cz^2 + dxy + eyz + fxz+ gx + hy + jz +k = 0`. + + Parameters + ---------- + coeffs : tuple + Tuple of general coefficients (a, b, c, d, e, f, g, h, j, k) + representing the quadric surface + + """ a, b, c, d, e, f, g, h, j, k = coeffs @classmethod @@ -1087,7 +1170,7 @@ class Cylinder(QuadricMeta, Surface): return cyl -class XCylinder(QuadricMeta, Surface): +class XCylinder(QuadricMixin, Surface): """An infinite cylinder whose length is parallel to the x-axis of the form :math:`(y - y_0)^2 + (z - z_0)^2 = r^2`. @@ -1146,10 +1229,6 @@ class XCylinder(QuadricMeta, Surface): self.z0 = z0 self.r = r - @property - def x0(self): - return self.coefficients['x0'] - @property def y0(self): return self.coefficients['y0'] @@ -1158,10 +1237,9 @@ class XCylinder(QuadricMeta, Surface): def z0(self): return self.coefficients['z0'] - @x0.setter - def x0(self, x0): - check_type('x0 coefficient', x0, Real) - self._coefficients['x0'] = x0 + @property + def r(self): + return self.coefficients['r'] @y0.setter def y0(self, y0): @@ -1173,12 +1251,43 @@ class XCylinder(QuadricMeta, Surface): check_type('z0 coefficient', z0, Real) self._coefficients['z0'] = z0 + @r.setter + def r(self, r): + check_type('r coefficient', r, Real) + self._coefficients['r'] = r + def _get_base_coeffs(self): - """Return generalized coefficients for a cylinder""" - return (0., 0., 1., self.z0) + """Return coefficients a, b, c, d, e, f, g, h, j, k representing a + general quadric surface of the form: + :math:`ax^2 + by^2 + cz^2 + dxy + eyz + fxz+ gx + hy + jz +k = 0`. + + """ + y0, z0, r = self.y0, self.z0, self.r + + a = d = e = f = g = 0. + b = c = 1. + h, j, k = -2*y0, -2*z0, y0**2 + z0**2 - r**2 + + return (a, b, c, d, e, f, g, h, j, k) def _update_from_base_coeffs(self, coeffs): + """Update the current surface from coefficients representing a general + quadric surface of the form: + :math:`ax^2 + by^2 + cz^2 + dxy + eyz + fxz+ gx + hy + jz +k = 0`. + + Parameters + ---------- + coeffs : tuple + Tuple of general coefficients (a, b, c, d, e, f, g, h, j, k) + representing the quadric surface + + """ a, b, c, d, e, f, g, h, j, k = coeffs + y0, z0 = -h / 2, -j / 2 + r = np.sqrt(y0**2 + z0**2 - k) + + for c, val in zip(self._coeff_keys, (y0, z0, r)): + setattr(self, c, val) def bounding_box(self, side): """Determine an axis-aligned bounding box. @@ -1213,7 +1322,7 @@ class XCylinder(QuadricMeta, Surface): np.array([np.inf, np.inf, np.inf])) -class YCylinder(QuadricMeta): +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`. @@ -1279,6 +1388,10 @@ class YCylinder(QuadricMeta): def z0(self): return self.coefficients['z0'] + @property + def r(self): + return self.coefficients['r'] + @x0.setter def x0(self, x0): check_type('x0 coefficient', x0, Real) @@ -1289,12 +1402,43 @@ class YCylinder(QuadricMeta): check_type('z0 coefficient', z0, Real) self._coefficients['z0'] = z0 + @r.setter + def r(self, r): + check_type('r coefficient', r, Real) + self._coefficients['r'] = r + def _get_base_coeffs(self): - """Return generalized coefficients for a cylinder""" - return (0., 0., 1., self.z0) + """Return coefficients a, b, c, d, e, f, g, h, j, k representing a + general quadric surface of the form: + :math:`ax^2 + by^2 + cz^2 + dxy + eyz + fxz+ gx + hy + jz +k = 0`. + + """ + x0, z0, r = self.x0, self.z0, self.r + + b = d = e = f = h = 0. + a = c = 1. + g, j, k = -2*x0, -2*z0, x0**2 + z0**2 - r**2 + + return (a, b, c, d, e, f, g, h, j, k) def _update_from_base_coeffs(self, coeffs): + """Update the current surface from coefficients representing a general + quadric surface of the form: + :math:`ax^2 + by^2 + cz^2 + dxy + eyz + fxz+ gx + hy + jz +k = 0`. + + Parameters + ---------- + coeffs : tuple + Tuple of general coefficients (a, b, c, d, e, f, g, h, j, k) + representing the quadric surface + + """ a, b, c, d, e, f, g, h, j, k = coeffs + x0, z0 = -g / 2, -j / 2 + r = np.sqrt(x0**2 + z0**2 - k) + + for c, val in zip(self._coeff_keys, (x0, z0, r)): + setattr(self, c, val) def bounding_box(self, side): """Determine an axis-aligned bounding box. @@ -1329,7 +1473,7 @@ class YCylinder(QuadricMeta): np.array([np.inf, np.inf, np.inf])) -class ZCylinder(QuadricMeta, Surface): +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`. @@ -1380,7 +1524,8 @@ class ZCylinder(QuadricMeta, Surface): def __init__(self, x0=0., y0=0., r=1., **kwargs): R = kwargs.pop('R', None) if R is not None: - warn(_WARNING_UPPER.format(type(self).__name__, 'r', 'R'), FutureWarning) + warn(_WARNING_UPPER.format(type(self).__name__, 'r', 'R'), + FutureWarning) r = R super().__init__(**kwargs) self.x0 = x0 @@ -1395,6 +1540,10 @@ class ZCylinder(QuadricMeta, Surface): def y0(self): return self.coefficients['y0'] + @property + def r(self): + return self.coefficients['r'] + @x0.setter def x0(self, x0): check_type('x0 coefficient', x0, Real) @@ -1405,12 +1554,43 @@ class ZCylinder(QuadricMeta, Surface): check_type('y0 coefficient', y0, Real) self._coefficients['y0'] = y0 + @r.setter + def r(self, r): + check_type('r coefficient', r, Real) + self._coefficients['r'] = r + def _get_base_coeffs(self): - """Return generalized coefficients for a cylinder""" - return (0., 0., 1., self.z0) + """Return coefficients a, b, c, d, e, f, g, h, j, k representing a + general quadric surface of the form: + :math:`ax^2 + by^2 + cz^2 + dxy + eyz + fxz+ gx + hy + jz +k = 0`. + + """ + x0, y0, r = self.x0, self.y0, self.r + + c = d = e = f = j = 0. + a = b = 1. + g, h, k = -2*x0, -2*y0, x0**2 + y0**2 - r**2 + + return (a, b, c, d, e, f, g, h, j, k) def _update_from_base_coeffs(self, coeffs): + """Update the current surface from coefficients representing a general + quadric surface of the form: + :math:`ax^2 + by^2 + cz^2 + dxy + eyz + fxz+ gx + hy + jz +k = 0`. + + Parameters + ---------- + coeffs : tuple + Tuple of general coefficients (a, b, c, d, e, f, g, h, j, k) + representing the quadric surface + + """ a, b, c, d, e, f, g, h, j, k = coeffs + x0, y0 = -g / 2, -4 / 2 + r = np.sqrt(x0**2 + y0**2 - k) + + for c, val in zip(self._coeff_keys, (x0, y0, r)): + setattr(self, c, val) def bounding_box(self, side): """Determine an axis-aligned bounding box. @@ -1445,7 +1625,7 @@ class ZCylinder(QuadricMeta, Surface): np.array([np.inf, np.inf, np.inf])) -class Sphere(QuadricMeta, Surface): +class Sphere(QuadricMixin, Surface): """A sphere of the form :math:`(x - x_0)^2 + (y - y_0)^2 + (z - z_0)^2 = r^2`. Parameters @@ -1543,11 +1723,35 @@ class Sphere(QuadricMeta, Surface): self._coefficients['r'] = r def _get_base_coeffs(self): - """Return generalized coefficients for a cylinder""" - return (0., 0., 1., self.z0) + """Return coefficients a, b, c, d, e, f, g, h, j, k representing a + general quadric surface of the form: + :math:`ax^2 + by^2 + cz^2 + dxy + eyz + fxz+ gx + hy + jz +k = 0`. + + """ + x0, y0, z0, r = self.x0, self.y0, self.z0, self.r + a = b = c = 1. + d = e = f = 0. + g, h, j = -2*x0, -2*y0, -2*z0 + k = x0**2 + y0**2 + z0**2 - r**2 + return (a, b, c, d, e, f, g, h, j, k) def _update_from_base_coeffs(self, coeffs): + """Update the current surface from coefficients representing a general + quadric surface of the form: + :math:`ax^2 + by^2 + cz^2 + dxy + eyz + fxz+ gx + hy + jz +k = 0`. + + Parameters + ---------- + coeffs : tuple + Tuple of general coefficients (a, b, c, d, e, f, g, h, j, k) + representing the quadric surface + + """ a, b, c, d, e, f, g, h, j, k = coeffs + x0, y0, z0 = -g / 2, -h / 2, -j / 2 + r = np.sqrt(x0**2 + y0**2 + z0**2 - k) + for c, val in zip(self._coeff_keys, (x0, y0, z0, r)): + setattr(self, c, val) def bounding_box(self, side): """Determine an axis-aligned bounding box. @@ -1583,7 +1787,7 @@ class Sphere(QuadricMeta, Surface): np.array([np.inf, np.inf, np.inf])) -class Cone(QuadricMeta, Surface): +class Cone(QuadricMixin, Surface): """A conical surface parallel to the x-, y-, or z-axis. Parameters @@ -1625,11 +1829,11 @@ class Cone(QuadricMeta, Surface): z-coordinate of the apex r2 : float Parameter related to the aperature - u : float + dx : float x-component of the vector representing the axis of the cone. - v : float + dy : float y-component of the vector representing the axis of the cone. - w : float + dz : float z-component of the vector representing the axis of the cone. boundary_type : {'transmission, 'vacuum', 'reflective', 'white'} Boundary condition that defines the behavior for particles hitting the @@ -1646,9 +1850,9 @@ class Cone(QuadricMeta, Surface): """ _type = 'cylinder' - _coeff_keys = ('x0', 'y0', 'z0', 'r2', 'u', 'v', 'w') + _coeff_keys = ('x0', 'y0', 'z0', 'r2', 'dx', 'dy', 'dz') - def __init__(self, x0=0., y0=0., z0=0., r2=1., u=0., v=0., w=1., **kwargs): + def __init__(self, x0=0., y0=0., z0=0., r2=1., dx=0., dy=0., dz=1., **kwargs): R2 = kwargs.pop('R2', None) if R2 is not None: warn(_WARNING_UPPER.format(type(self).__name__, 'r2', 'R2'), @@ -1659,9 +1863,9 @@ class Cone(QuadricMeta, Surface): self.y0 = y0 self.z0 = z0 self.r2 = r2 - self.u = u - self.v = v - self.w = w + self.dx = dx + self.dy = dy + self.dz = dz @property def x0(self): @@ -1680,16 +1884,16 @@ class Cone(QuadricMeta, Surface): return self.coefficients['r2'] @property - def u(self): - return self.coefficients['u'] + def dx(self): + return self.coefficients['dx'] @property - def v(self): - return self.coefficients['v'] + def dy(self): + return self.coefficients['dy'] @property - def w(self): - return self.coefficients['w'] + def dz(self): + return self.coefficients['dz'] @x0.setter def x0(self, x0): @@ -1711,30 +1915,45 @@ class Cone(QuadricMeta, Surface): check_type('r^2 coefficient', r2, Real) self._coefficients['r2'] = r2 - @u.setter - def u(self, u): - check_type('u coefficient', u, Real) - self._coefficients['u'] = u + @dx.setter + def dx(self, dx): + check_type('dx coefficient', dx, Real) + self._coefficients['dx'] = dx - @v.setter - def v(self, v): - check_type('v coefficient', v, Real) - self._coefficients['v'] = v + @dy.setter + def dy(self, dy): + check_type('dy coefficient', dy, Real) + self._coefficients['dy'] = dy - @w.setter - def w(self, w): - check_type('w coefficient', w, Real) - self._coefficients['w'] = w + @dz.setter + def dz(self, dz): + check_type('dz coefficient', dz, Real) + self._coefficients['dz'] = dz def _get_base_coeffs(self): - """Return generalized coefficients for a cylinder""" + """Return coefficients a, b, c, d, e, f, g, h, j, k representing a + general quadric surface of the form: + :math:`ax^2 + by^2 + cz^2 + dxy + eyz + fxz+ gx + hy + jz +k = 0`. + + """ return (0., 0., 1., self.z0) def _update_from_base_coeffs(self, coeffs): + """Update the current surface from coefficients representing a general + quadric surface of the form: + :math:`ax^2 + by^2 + cz^2 + dxy + eyz + fxz+ gx + hy + jz +k = 0`. + + Parameters + ---------- + coeffs : tuple + Tuple of general coefficients (a, b, c, d, e, f, g, h, j, k) + representing the quadric surface + + """ a, b, c, d, e, f, g, h, j, k = coeffs -class XCone(QuadricMeta, Surface): +class XCone(QuadricMixin, Surface): """A cone parallel to the x-axis of the form :math:`(y - y_0)^2 + (z - z_0)^2 = r^2 (x - x_0)^2`. @@ -1834,14 +2053,46 @@ class XCone(QuadricMeta, Surface): self._coefficients['r2'] = r2 def _get_base_coeffs(self): - """Return generalized coefficients for a cylinder""" - return (0., 0., 1., self.z0) + """Return coefficients a, b, c, d, e, f, g, h, j, k representing a + general quadric surface of the form: + :math:`ax^2 + by^2 + cz^2 + dxy + eyz + fxz+ gx + hy + jz +k = 0`. + + """ + x0, y0, z0, r2 = self.x0, self.y0, self.z0, self.r2 + + a = -r2 + b = c = 1. + d = e = f = 0. + g, h, j = 2*x0*r2, -2*y0, -2*z0 + k = y0**2 + z0**2 - r2*x0**2 + + return (a, b, c, d, e, f, g, h, j, k) def _update_from_base_coeffs(self, coeffs): + """Update the current surface from coefficients representing a general + quadric surface of the form: + :math:`ax^2 + by^2 + cz^2 + dxy + eyz + fxz+ gx + hy + jz +k = 0`. + + Parameters + ---------- + coeffs : tuple + Tuple of general coefficients (a, b, c, d, e, f, g, h, j, k) + representing the quadric surface + + """ a, b, c, d, e, f, g, h, j, k = coeffs + r2 = -a + x0, y0, z0 = g / (2*r2), -h / 2, -j / 2 -class YCone(QuadricMeta, Surface): + for c, val in zip(self._coeff_keys, (x0, y0, z0, r2)): + setattr(self, c, val) + + +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`. @@ -1941,14 +2192,46 @@ class YCone(QuadricMeta, Surface): self._coefficients['r2'] = r2 def _get_base_coeffs(self): - """Return generalized coefficients for a cylinder""" - return (0., 0., 1., self.z0) + """Return coefficients a, b, c, d, e, f, g, h, j, k representing a + general quadric surface of the form: + :math:`ax^2 + by^2 + cz^2 + dxy + eyz + fxz+ gx + hy + jz +k = 0`. + + """ + x0, y0, z0, r2 = self.x0, self.y0, self.z0, self.r2 + + b = -r2 + a = c = 1. + d = e = f = 0. + g, h, j = -2*x0, 2*y0*r2, -2*z0 + k = x0**2 + z0**2 - r2*y0**2 + + return (a, b, c, d, e, f, g, h, j, k) def _update_from_base_coeffs(self, coeffs): + """Update the current surface from coefficients representing a general + quadric surface of the form: + :math:`ax^2 + by^2 + cz^2 + dxy + eyz + fxz+ gx + hy + jz +k = 0`. + + Parameters + ---------- + coeffs : tuple + Tuple of general coefficients (a, b, c, d, e, f, g, h, j, k) + representing the quadric surface + + """ a, b, c, d, e, f, g, h, j, k = coeffs + r2 = -b + x0, y0, z0 = -g / 2, h / (2*r2), -j / 2 -class ZCone(QuadricMeta, Surface): + for c, val in zip(self._coeff_keys, (x0, y0, z0, r2)): + setattr(self, c, val) + + +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`. @@ -2048,14 +2331,46 @@ class ZCone(QuadricMeta, Surface): self._coefficients['r2'] = r2 def _get_base_coeffs(self): - """Return generalized coefficients for a cylinder""" - return (0., 0., 1., self.z0) + """Return coefficients a, b, c, d, e, f, g, h, j, k representing a + general quadric surface of the form: + :math:`ax^2 + by^2 + cz^2 + dxy + eyz + fxz+ gx + hy + jz +k = 0`. + + """ + x0, y0, z0, r2 = self.x0, self.y0, self.z0, self.r2 + + c = -r2 + a = b = 1. + d = e = f = 0. + g, h, j = -2*x0, -2*y0, 2*z0*r2 + k = x0**2 + y0**2 - r2*z0**2 + + return (a, b, c, d, e, f, g, h, j, k) def _update_from_base_coeffs(self, coeffs): + """Update the current surface from coefficients representing a general + quadric surface of the form: + :math:`ax^2 + by^2 + cz^2 + dxy + eyz + fxz+ gx + hy + jz +k = 0`. + + Parameters + ---------- + coeffs : tuple + Tuple of general coefficients (a, b, c, d, e, f, g, h, j, k) + representing the quadric surface + + """ a, b, c, d, e, f, g, h, j, k = coeffs + r2 = -c + x0, y0, z0 = -g / 2, -h / 2, j / (2*r2) -class Quadric(QuadricMeta, Surface): + for c, val in zip(self._coeff_keys, (x0, y0, z0, r2)): + setattr(self, c, val) + + +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`. @@ -2096,17 +2411,11 @@ class Quadric(QuadricMeta, Surface): def __init__(self, a=0., b=0., c=0., d=0., e=0., f=0., g=0., h=0., j=0., k=0., **kwargs): + super().__init__(**kwargs) - self.a = a - self.b = b - self.c = c - self.d = d - self.e = e - self.f = f - self.g = g - self.h = h - self.j = j - self.k = k + + for c, val in zip(self._coeff_keys, (a, b, c, d, e, f, g, h, j, k)): + setattr(self, c, val) @property def a(self): @@ -2198,6 +2507,29 @@ class Quadric(QuadricMeta, Surface): check_type('k coefficient', k, Real) self._coefficients['k'] = k + def _get_base_coeffs(self): + """Return coefficients a, b, c, d, e, f, g, h, j, k representing a + general quadric surface of the form: + :math:`ax^2 + by^2 + cz^2 + dxy + eyz + fxz+ gx + hy + jz +k = 0`. + + """ + return tuple(getattr(self, c) for c in self._coeff_keys) + + def _update_from_base_coeffs(self, coeffs): + """Update the current surface from coefficients representing a general + quadric surface of the form: + :math:`ax^2 + by^2 + cz^2 + dxy + eyz + fxz+ gx + hy + jz +k = 0`. + + Parameters + ---------- + coeffs : tuple + Tuple of general coefficients (a, b, c, d, e, f, g, h, j, k) + representing the quadric surface + + """ + for c, val in zip(self._coeff_keys, coeffs): + setattr(self, c, val) + class Halfspace(Region): """A positive or negative half-space region. From 0ad5dbd446ef4e9dd434fff4b648b6af553fb6cc Mon Sep 17 00:00:00 2001 From: Ethan Peterson Date: Tue, 11 Feb 2020 12:10:17 -0500 Subject: [PATCH 07/33] cleaning up and making code more self-consistent --- openmc/surface.py | 49 +++++++++++++++++++++-------------------------- 1 file changed, 22 insertions(+), 27 deletions(-) diff --git a/openmc/surface.py b/openmc/surface.py index 8b05ce8c53..44178f6bc6 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -1103,9 +1103,9 @@ class Cylinder(QuadricMixin, Surface): d = -2*dx*dy e = -2*dy*dz f = -2*dx*dz - g = -2*((z1 - z0)*(x0*z1 - x1*z0) + (y1 - y0)*(x0*y1-x1*y0)) - h = 2*((x1 - x0)*(x0*y1 - x1*y0) - (z1 - z0)*(y0*z1 - y1*z0)) - j = 2*((x1 - x0)*(x0*z1 - x1*z0) + (y1 - y0)*(y0*z1 - y1*z0)) + g = -2*(dz*(x0*z1 - x1*z0) + dy*(x0*y1-x1*y0)) + h = 2*(dx*(x0*y1 - x1*y0) - dz*(y0*z1 - y1*z0)) + j = 2*(dx*(x0*z1 - x1*z0) + dy*(y0*z1 - y1*z0)) k = (y0*z1 - y1*z0)**2 + (x0*z1 - x1*z0)**2 + (x0*y1 - x1*y0)**2 \ - r**2*(dx**2 + dy**2 + dz**2) @@ -1132,42 +1132,26 @@ class Cylinder(QuadricMixin, Surface): Parameters ---------- p1, p2 : 3-tuples - Points that pass through the plane + Points that pass through the plane, p1 will be used as (x0, y0, z0) r : float, optional Radius of the cylinder. Defaults to 1. kwargs : dict - Keyword arguments passed to the :class:`Quadric` constructor + Keyword arguments passed to the :class:`Cylinder` constructor Returns ------- Cylinder - Cylinder that has an axis through the points p1 and p2 + Cylinder that has an axis through the points p1 and p2, and a + radius r. """ # Convert to numpy arrays p1 = np.asarray(p1) p2 = np.asarray(p2) - x1, y1, z1 = p1 - x2, y2, z2 = p2 + x0, y0, z0 = p1 dx, dy, dz = p2 - p1 - # Set coefficients for Quadric surface that represents a cylinder of - # radius r whose axis is the line defined by p1 and p2 - a = dy**2 + dz**2 - b = dx**2 + dz**2 - c = dx**2 + dy**2 - d = -2*dx*dy - e = -2*dy*dz - f = -2*dx*dz - g = -2*((z2 - z1)*(x1*z2 - x2*z1) + (y2 - y1)*(x1*y2-x2*y1)) - h = 2*((x2 - x1)*(x1*y2 - x2*y1) - (z2 - z1)*(y1*z2 - y2*z1)) - j = 2*((x2 - x1)*(x1*z2 - x2*z1) + (y2 - y1)*(y1*z2 - y2*z1)) - k = (y1*z2 - y2*z1)**2 + (x1*z2 - x2*z1)**2 + (x1*y2 - x2*y1)**2 \ - - r**2*(dx**2 + dy**2 + dz**2) - - cyl = cls(x0=x1, y0=y1, z0=z1, r=r, **kwargs) - cyl._update_from_base_coeffs((a, b, c, d, e, f, g, h, j, k)) - return cyl + return cls(x0=x0, y0=y0, z0=z0, r=r, dx=dx, dy=dy, dz=dz, **kwargs) class XCylinder(QuadricMixin, Surface): @@ -1322,6 +1306,9 @@ class XCylinder(QuadricMixin, Surface): np.array([np.inf, np.inf, np.inf])) +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`. @@ -1373,7 +1360,8 @@ class YCylinder(QuadricMixin, Surface): def __init__(self, x0=0., z0=0., r=1., **kwargs): R = kwargs.pop('R', None) if R is not None: - warn(_WARNING_UPPER.format(type(self).__name__, 'r', 'R'), FutureWarning) + warn(_WARNING_UPPER.format(type(self).__name__, 'r', 'R'), + FutureWarning) r = R super().__init__(**kwargs) self.x0 = x0 @@ -1473,6 +1461,9 @@ class YCylinder(QuadricMixin, Surface): np.array([np.inf, np.inf, np.inf])) +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`. @@ -1625,6 +1616,9 @@ class ZCylinder(QuadricMixin, Surface): np.array([np.inf, np.inf, np.inf])) +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`. @@ -1678,7 +1672,8 @@ class Sphere(QuadricMixin, Surface): def __init__(self, x0=0., y0=0., z0=0., r=1., **kwargs): R = kwargs.pop('R', None) if R is not None: - warn(_WARNING_UPPER.format(type(self).__name__, 'r', 'R'), FutureWarning) + warn(_WARNING_UPPER.format(type(self).__name__, 'r', 'R'), + FutureWarning) r = R super().__init__(**kwargs) self.x0 = x0 From 6b3a23199427e4e71892b78cce495ab012b9f863 Mon Sep 17 00:00:00 2001 From: Ethan Peterson Date: Tue, 11 Feb 2020 12:34:31 -0500 Subject: [PATCH 08/33] fixed issue related to finding all surface subclasses --- openmc/surface.py | 31 ++++++++++++++++--------------- 1 file changed, 16 insertions(+), 15 deletions(-) diff --git a/openmc/surface.py b/openmc/surface.py index 44178f6bc6..24b21cf288 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -188,17 +188,6 @@ class Surface(IDManagerMixin, metaclass=ABCMeta): def translate(self, vector): pass - @classmethod - def get_subclasses(cls): - """Recursively find all subclasses of this class""" - return set(cls.__subclasses__()).union([s for c in cls.__subclasses__() - for s in get_subclasses(c)]) - - @classmethod - def get_subclass_map(cls): - """Generate mapping of class _type attributes to classes""" - return {c._type: c for c in cls.get_subclasses()} - def to_xml_element(self): """Return XML representation of the surface @@ -279,9 +268,6 @@ class Surface(IDManagerMixin, metaclass=ABCMeta): return cls(*coeffs, **kwargs) -_SURFACE_CLASSES = Surface.get_subclass_map() - - class PlaneMixin(metaclass=ABCMeta): """A Plane Meta class for all operations on order 1 surfaces""" @@ -1844,7 +1830,7 @@ class Cone(QuadricMixin, Surface): """ - _type = 'cylinder' + _type = 'cone' _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): @@ -2712,3 +2698,18 @@ class Halfspace(Region): # Return translated surface return type(self)(memo[key], self.side) + +def get_subclasses(cls): + """Recursively find all subclasses of this class""" + return set(cls.__subclasses__()).union([s for c in cls.__subclasses__() + for s in get_subclasses(c)]) + + +def get_subclass_map(cls): + """Generate mapping of class _type attributes to classes""" + return {c._type: c for c in get_subclasses(cls)} + + +_SURFACE_CLASSES = get_subclass_map(Surface) + + From a2f96ffe3438ab635effc28458f47c1ac7a1d406 Mon Sep 17 00:00:00 2001 From: Ethan Peterson Date: Tue, 11 Feb 2020 15:11:49 -0500 Subject: [PATCH 09/33] fixed bug in quadric translate function --- openmc/surface.py | 62 ++++++++++++++++++++++++++++++++++++++--------- 1 file changed, 51 insertions(+), 11 deletions(-) diff --git a/openmc/surface.py b/openmc/surface.py index 24b21cf288..41befbb416 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -20,6 +20,11 @@ _WARNING_UPPER = """\ will not accept the capitalized version.\ """ +_WARNING_KWARGS = """ +"{}(...) accepts keyword arguments only for '{}'. Future versions of OpenMC \ +will not accept positional parameters for superclass arguments.\ +""" + class Surface(IDManagerMixin, metaclass=ABCMeta): """An implicit surface with an associated boundary condition. @@ -336,7 +341,7 @@ class PlaneMixin(metaclass=ABCMeta): a, b, c, d = self._get_base_coeffs() return a*x + b*y + c*z - d - def translate(self, vector, clone=False): + def translate(self, vector, clone=True): """Translate surface in given direction Parameters @@ -345,7 +350,7 @@ class PlaneMixin(metaclass=ABCMeta): Direction in which surface should be translated clone : boolean Whether or not to return a new instance of a Plane or to modify the - coefficients of this plane. + coefficients of this plane. Defaults to True Returns ------- @@ -466,12 +471,21 @@ class Plane(PlaneMixin, Surface): _type = 'plane' _coeff_keys = ('a', 'b', 'c', 'd') - def __init__(self, a=1., b=0., c=0., d=0., **kwargs): + def __init__(self, a=1., b=0., c=0., d=0., *args, **kwargs): # work around until capital letter kwargs are deprecated oldkwargs = deepcopy(kwargs) + # work around for accepting Surface kwargs as positional parameters + # until they are deprecated + argsdict = {k: v for k, v in zip(('boundary_type', 'name', + 'surface_id'), args)} + for k, v in argsdict.items(): + warn(_WARNING_KWARGS.format(type(self).__name__, k), FutureWarning) + for k in 'ABCD': kwargs.pop(k, None) + kwargs.update(argsdict) + super().__init__(**kwargs) self.a = a self.b = b @@ -615,7 +629,15 @@ class XPlane(PlaneMixin, Surface): _type = 'x-plane' _coeff_keys = ('x0',) - def __init__(self, x0=0., **kwargs): + def __init__(self, x0=0., *args, **kwargs): + # work around for accepting Surface kwargs as positional parameters + # until they are deprecated + argsdict = {k: v for k, v in zip(('boundary_type', 'name', + 'surface_id'), args)} + for k, v in argsdict.items(): + warn(_WARNING_KWARGS.format(type(self).__name__, k), FutureWarning) + kwargs.update(argsdict) + super().__init__(**kwargs) self.x0 = x0 @@ -704,7 +726,15 @@ class YPlane(PlaneMixin, Surface): _type = 'y-plane' _coeff_keys = ('y0',) - def __init__(self, y0=0., **kwargs): + def __init__(self, y0=0., *args, **kwargs): + # work around for accepting Surface kwargs as positional parameters + # until they are deprecated + argsdict = {k: v for k, v in zip(('boundary_type', 'name', + 'surface_id'), args)} + for k, v in argsdict.items(): + warn(_WARNING_KWARGS.format(type(self).__name__, k), FutureWarning) + kwargs.update(argsdict) + super().__init__(**kwargs) self.y0 = y0 @@ -793,7 +823,15 @@ class ZPlane(PlaneMixin, Surface): _type = 'z-plane' _coeff_keys = ('z0',) - def __init__(self, z0=0., **kwargs): + def __init__(self, z0=0., *args, **kwargs): + # work around for accepting Surface kwargs as positional parameters + # until they are deprecated + argsdict = {k: v for k, v in zip(('boundary_type', 'name', + 'surface_id'), args)} + for k, v in argsdict.items(): + warn(_WARNING_KWARGS.format(type(self).__name__, k), FutureWarning) + kwargs.update(argsdict) + super().__init__(**kwargs) self.z0 = z0 @@ -900,7 +938,7 @@ class QuadricMixin(metaclass=ABCMeta): a, b, c, d, e, f, g, h, j, k = self._get_base_coeffs() return x*(a*x + d*y + g) + y*(b*y + e*z + h) + z*(c*z + f*x + j) + k - def translate(self, vector, clone=False): + def translate(self, vector, clone=True): """Translate surface in given direction Parameters @@ -908,17 +946,18 @@ class QuadricMixin(metaclass=ABCMeta): vector : iterable of float Direction in which surface should be translated clone : bool - Whether to return a clone of the Surface or the Surface itself + Whether to return a clone of the Surface or the Surface itself. + Defaults to True Returns ------- - openmc.QuadricMixin + openmc.Surface Translated surface """ vx, vy, vz = vector a, b, c, d, e, f, g, h, j, k = self._get_base_coeffs() - k = (k + vx*vx + vy*vy + vz*vz + d*vx*vy + e*vy*vz + f*vx*vz + k = (k + a*vx**2 + b*vy**2 + c*vz**2 + d*vx*vy + e*vy*vz + f*vx*vz - g*vx - h*vy - j*vz) g = g - 2*a*vx - d*vy - f*vz h = h - 2*b*vy - d*vx - e*vz @@ -1237,6 +1276,7 @@ class XCylinder(QuadricMixin, Surface): a = d = e = f = g = 0. b = c = 1. h, j, k = -2*y0, -2*z0, y0**2 + z0**2 - r**2 + print("to base", (a, b, c, d, e, f, g, h, j, k)) return (a, b, c, d, e, f, g, h, j, k) @@ -1563,7 +1603,7 @@ class ZCylinder(QuadricMixin, Surface): """ a, b, c, d, e, f, g, h, j, k = coeffs - x0, y0 = -g / 2, -4 / 2 + x0, y0 = -g / 2, -h / 2 r = np.sqrt(x0**2 + y0**2 - k) for c, val in zip(self._coeff_keys, (x0, y0, r)): From 2a1c48e9878c463df62465cc3a927448272845fe Mon Sep 17 00:00:00 2001 From: Ethan Peterson Date: Tue, 11 Feb 2020 16:15:32 -0500 Subject: [PATCH 10/33] refactoring bounding_box and cleaning up constructors --- openmc/surface.py | 297 ++++++++++++++++++---------------------------- 1 file changed, 113 insertions(+), 184 deletions(-) diff --git a/openmc/surface.py b/openmc/surface.py index 41befbb416..5f5aecfc71 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -487,10 +487,10 @@ class Plane(PlaneMixin, Surface): kwargs.update(argsdict) super().__init__(**kwargs) - self.a = a - self.b = b - self.c = c - self.d = d + + for key, val in zip(self._coeff_keys, (a, b, c, d)): + setattr(self, key, val) + for k, v in oldkwargs.items(): if k in 'ABCD': warn(_WARNING_UPPER.format(type(self).__name__, k.lower(), k), @@ -551,8 +551,8 @@ class Plane(PlaneMixin, Surface): """ - for c, val in zip(self._coeff_keys, coeffs): - setattr(self, c, val) + for key, val in zip(self._coeff_keys, coeffs): + setattr(self, key, val) @classmethod def from_points(cls, p1, p2, p3, **kwargs): @@ -918,6 +918,34 @@ class QuadricMixin(metaclass=ABCMeta): return (np.array([-np.inf, -np.inf, -np.inf]), np.array([np.inf, np.inf, np.inf])) + 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. To represent infinity, + numpy.inf is used. + + 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 + + """ + + if side == '-': + return self._neg_bounds() + elif side == '+': + return self._pos_bounds() + def evaluate(self, point): """Evaluate the surface equation at a given point. @@ -1041,13 +1069,9 @@ 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): super().__init__(**kwargs) - self.x0 = x0 - self.y0 = y0 - self.z0 = z0 - self.r = r - self.dx = dx - self.dy = dy - self.dz = dz + + for key, val in zip(self._coeff_keys, (x0, y0, z0, r, dx, dy, dz)): + setattr(self, key, val) @property def x0(self): @@ -1234,9 +1258,9 @@ class XCylinder(QuadricMixin, Surface): FutureWarning) r = R super().__init__(**kwargs) - self.y0 = y0 - self.z0 = z0 - self.r = r + + for key, val in zip(self._coeff_keys, (y0, z0, r)): + setattr(self, key, val) @property def y0(self): @@ -1296,40 +1320,18 @@ class XCylinder(QuadricMixin, Surface): y0, z0 = -h / 2, -j / 2 r = np.sqrt(y0**2 + z0**2 - k) - for c, val in zip(self._coeff_keys, (y0, z0, r)): - setattr(self, c, val) + for key, val in zip(self._coeff_keys, (y0, z0, r)): + setattr(self, key, val) - def bounding_box(self, side): - """Determine an axis-aligned bounding box. + def _neg_bounds(self): + """Return the lower and upper bounds of the negative half space""" + 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])) - 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.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 - - """ - - 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])) + def _pos_bounds(self): + """Return the lower and upper bounds of the positive half space""" + return (np.array([-np.inf, -np.inf, -np.inf]), + np.array([np.inf, np.inf, np.inf])) Cylinder.register(XCylinder) @@ -1390,9 +1392,9 @@ class YCylinder(QuadricMixin, Surface): FutureWarning) r = R super().__init__(**kwargs) - self.x0 = x0 - self.z0 = z0 - self.r = r + + for key, val in zip(self._coeff_keys, (x0, z0, r)): + setattr(self, key, val) @property def x0(self): @@ -1451,40 +1453,18 @@ class YCylinder(QuadricMixin, Surface): x0, z0 = -g / 2, -j / 2 r = np.sqrt(x0**2 + z0**2 - k) - for c, val in zip(self._coeff_keys, (x0, z0, r)): - setattr(self, c, val) + for key, val in zip(self._coeff_keys, (x0, z0, r)): + setattr(self, key, val) - def bounding_box(self, side): - """Determine an axis-aligned bounding box. + def _neg_bounds(self): + """Return the lower and upper bounds of the negative half space""" + 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])) - 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.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 - - """ - - 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])) + def _pos_bounds(self): + """Return the lower and upper bounds of the positive half space""" + return (np.array([-np.inf, -np.inf, -np.inf]), + np.array([np.inf, np.inf, np.inf])) Cylinder.register(YCylinder) @@ -1545,9 +1525,9 @@ class ZCylinder(QuadricMixin, Surface): FutureWarning) r = R super().__init__(**kwargs) - self.x0 = x0 - self.y0 = y0 - self.r = r + + for key, val in zip(self._coeff_keys, (x0, y0, r)): + setattr(self, key, val) @property def x0(self): @@ -1606,40 +1586,18 @@ class ZCylinder(QuadricMixin, Surface): x0, y0 = -g / 2, -h / 2 r = np.sqrt(x0**2 + y0**2 - k) - for c, val in zip(self._coeff_keys, (x0, y0, r)): - setattr(self, c, val) + for key, val in zip(self._coeff_keys, (x0, y0, r)): + setattr(self, key, val) - def bounding_box(self, side): - """Determine an axis-aligned bounding box. + def _neg_bounds(self): + """Return the lower and upper bounds of the negative half space""" + 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])) - 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.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 - - """ - - 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])) + def _pos_bounds(self): + """Return the lower and upper bounds of the positive half space""" + return (np.array([-np.inf, -np.inf, -np.inf]), + np.array([np.inf, np.inf, np.inf])) Cylinder.register(ZCylinder) @@ -1702,10 +1660,9 @@ class Sphere(QuadricMixin, Surface): FutureWarning) r = R super().__init__(**kwargs) - self.x0 = x0 - self.y0 = y0 - self.z0 = z0 - self.r = r + + for key, val in zip(self._coeff_keys, (x0, y0, z0, r)): + setattr(self, key, val) @property def x0(self): @@ -1771,41 +1728,20 @@ class Sphere(QuadricMixin, Surface): a, b, c, d, e, f, g, h, j, k = coeffs x0, y0, z0 = -g / 2, -h / 2, -j / 2 r = np.sqrt(x0**2 + y0**2 + z0**2 - k) - for c, val in zip(self._coeff_keys, (x0, y0, z0, r)): - setattr(self, c, val) + for key, val in zip(self._coeff_keys, (x0, y0, z0, r)): + setattr(self, key, val) - def bounding_box(self, side): - """Determine an axis-aligned bounding box. + def _neg_bounds(self): + """Return the lower and upper bounds of the negative half space""" + 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])) - 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.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 - - """ - - 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])) + def _pos_bounds(self): + """Return the lower and upper bounds of the positive half space""" + return (np.array([-np.inf, -np.inf, -np.inf]), + np.array([np.inf, np.inf, np.inf])) class Cone(QuadricMixin, Surface): @@ -1880,13 +1816,9 @@ class Cone(QuadricMixin, Surface): FutureWarning) r2 = R2 super().__init__(**kwargs) - self.x0 = x0 - self.y0 = y0 - self.z0 = z0 - self.r2 = r2 - self.dx = dx - self.dy = dy - self.dz = dz + + for key, val in zip(self._coeff_keys, (x0, y0, z0, r2, dx, dy, dz)): + setattr(self, key, val) @property def x0(self): @@ -2032,10 +1964,9 @@ class XCone(QuadricMixin, Surface): FutureWarning) r2 = R2 super().__init__(**kwargs) - self.x0 = x0 - self.y0 = y0 - self.z0 = z0 - self.r2 = r2 + + for key, val in zip(self._coeff_keys, (x0, y0, z0, r2)): + setattr(self, key, val) @property def x0(self): @@ -2106,8 +2037,8 @@ class XCone(QuadricMixin, Surface): r2 = -a x0, y0, z0 = g / (2*r2), -h / 2, -j / 2 - for c, val in zip(self._coeff_keys, (x0, y0, z0, r2)): - setattr(self, c, val) + for key, val in zip(self._coeff_keys, (x0, y0, z0, r2)): + setattr(self, key, val) Cone.register(XCone) @@ -2171,10 +2102,9 @@ class YCone(QuadricMixin, Surface): FutureWarning) r2 = R2 super().__init__(**kwargs) - self.x0 = x0 - self.y0 = y0 - self.z0 = z0 - self.r2 = r2 + + for key, val in zip(self._coeff_keys, (x0, y0, z0, r2)): + setattr(self, key, val) @property def x0(self): @@ -2245,8 +2175,8 @@ class YCone(QuadricMixin, Surface): r2 = -b x0, y0, z0 = -g / 2, h / (2*r2), -j / 2 - for c, val in zip(self._coeff_keys, (x0, y0, z0, r2)): - setattr(self, c, val) + for key, val in zip(self._coeff_keys, (x0, y0, z0, r2)): + setattr(self, key, val) Cone.register(YCone) @@ -2310,10 +2240,9 @@ class ZCone(QuadricMixin, Surface): FutureWarning) r2 = R2 super().__init__(**kwargs) - self.x0 = x0 - self.y0 = y0 - self.z0 = z0 - self.r2 = r2 + + for key, val in zip(self._coeff_keys, (x0, y0, z0, r2)): + setattr(self, key, val) @property def x0(self): @@ -2384,8 +2313,8 @@ class ZCone(QuadricMixin, Surface): r2 = -c x0, y0, z0 = -g / 2, -h / 2, j / (2*r2) - for c, val in zip(self._coeff_keys, (x0, y0, z0, r2)): - setattr(self, c, val) + for key, val in zip(self._coeff_keys, (x0, y0, z0, r2)): + setattr(self, key, val) Cone.register(ZCone) @@ -2435,8 +2364,8 @@ class Quadric(QuadricMixin, Surface): super().__init__(**kwargs) - for c, val in zip(self._coeff_keys, (a, b, c, d, e, f, g, h, j, k)): - setattr(self, c, val) + for key, val in zip(self._coeff_keys, (a, b, c, d, e, f, g, h, j, k)): + setattr(self, key, val) @property def a(self): @@ -2548,8 +2477,8 @@ class Quadric(QuadricMixin, Surface): representing the quadric surface """ - for c, val in zip(self._coeff_keys, coeffs): - setattr(self, c, val) + for key, val in zip(self._coeff_keys, coeffs): + setattr(self, key, val) class Halfspace(Region): From 29255732caec4ec2ba4a1652b5e29d097dd66080 Mon Sep 17 00:00:00 2001 From: Ethan Peterson Date: Tue, 11 Feb 2020 16:43:49 -0500 Subject: [PATCH 11/33] changed subclass finding methods to private --- openmc/surface.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/openmc/surface.py b/openmc/surface.py index 5f5aecfc71..ca1501b887 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -25,6 +25,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. @@ -275,7 +277,6 @@ class Surface(IDManagerMixin, metaclass=ABCMeta): class PlaneMixin(metaclass=ABCMeta): """A Plane Meta class for all operations on order 1 surfaces""" - def __init__(self, **kwargs): super().__init__(**kwargs) self._periodic_surface = None @@ -880,7 +881,6 @@ Plane.register(ZPlane) class QuadricMixin(metaclass=ABCMeta): """A Mixin class implementing common functionality for quadric surfaces""" - def __init__(self, **kwargs): super().__init__(**kwargs) @@ -1300,7 +1300,6 @@ class XCylinder(QuadricMixin, Surface): a = d = e = f = g = 0. b = c = 1. h, j, k = -2*y0, -2*z0, y0**2 + z0**2 - r**2 - print("to base", (a, b, c, d, e, f, g, h, j, k)) return (a, b, c, d, e, f, g, h, j, k) @@ -2668,17 +2667,18 @@ class Halfspace(Region): # Return translated surface return type(self)(memo[key], self.side) -def get_subclasses(cls): + +def _get_subclasses(cls): """Recursively find all subclasses of this class""" return set(cls.__subclasses__()).union([s for c in cls.__subclasses__() - for s in get_subclasses(c)]) + for s in _get_subclasses(c)]) -def get_subclass_map(cls): +def _get_subclass_map(cls): """Generate mapping of class _type attributes to classes""" - return {c._type: c for c in get_subclasses(cls)} + return {c._type: c for c in _get_subclasses(cls)} -_SURFACE_CLASSES = get_subclass_map(Surface) +_SURFACE_CLASSES = _get_subclass_map(Surface) From caeb4d3ce5b14936a8dc56921929063af03f46f5 Mon Sep 17 00:00:00 2001 From: Ethan Peterson Date: Tue, 11 Feb 2020 16:49:38 -0500 Subject: [PATCH 12/33] removed empty dictionary definition --- openmc/surface.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/openmc/surface.py b/openmc/surface.py index ca1501b887..5dc25103f3 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -25,8 +25,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 994bd221badf2fd72e1a262f0a636060b01d5e67 Mon Sep 17 00:00:00 2001 From: Ethan Peterson Date: Tue, 11 Feb 2020 21:51:34 -0500 Subject: [PATCH 13/33] fixed general cylinder _to_base_coeffs method --- openmc/surface.py | 20 ++++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/openmc/surface.py b/openmc/surface.py index 5dc25103f3..badbd630c2 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -1137,12 +1137,11 @@ class Cylinder(QuadricMixin, Surface): def _get_base_coeffs(self): """Return coefficients a, b, c, d, e, f, g, h, j, k representing a general quadric surface of the form: - :math:`ax^2 + by^2 + cz^2 + dxy + eyz + fxz+ gx + hy + jz +k = 0`. + :math:`ax^2 + by^2 + cz^2 + dxy + eyz + fxz+ gx + hy + jz + k = 0`. """ x0, y0, z0, r = self.x0, self.y0, self.z0, self.r dx, dy, dz = self.dx, self.dy, self.dz - x1, y1, z1 = x0 + dx, y0 + dy, z0 + dz a = dy**2 + dz**2 b = dx**2 + dz**2 @@ -1150,18 +1149,19 @@ class Cylinder(QuadricMixin, Surface): d = -2*dx*dy e = -2*dy*dz f = -2*dx*dz - g = -2*(dz*(x0*z1 - x1*z0) + dy*(x0*y1-x1*y0)) - h = 2*(dx*(x0*y1 - x1*y0) - dz*(y0*z1 - y1*z0)) - j = 2*(dx*(x0*z1 - x1*z0) + dy*(y0*z1 - y1*z0)) - k = (y0*z1 - y1*z0)**2 + (x0*z1 - x1*z0)**2 + (x0*y1 - x1*y0)**2 \ - - r**2*(dx**2 + dy**2 + dz**2) + g = 2*(dx*(z0*dz + y0*dy) - x0*(dy**2 + dz**2)) + h = 2*(dy*(z0*dz + x0*dx) - y0*(dx**2 + dz**2)) + j = 2*(dz*(y0*dy + x0*dx) - z0*(dx**2 + dy**2)) + k = x0**2*(dy**2 + dz**2) + y0**2*(dx**2 + dz**2) \ + + z0**2*(dx**2 + dy**2) - 2*dy*dz*y0*z0 \ + - 2*dx*dz*x0*z0 - 2*dx*dy*x0*y0 - r**2*(dx**2 + dy**2 + dz**2) return (a, b, c, d, e, f, g, h, j, k) def _update_from_base_coeffs(self, coeffs): """Update the current surface from coefficients representing a general quadric surface of the form: - :math:`ax^2 + by^2 + cz^2 + dxy + eyz + fxz+ gx + hy + jz +k = 0`. + :math:`ax^2 + by^2 + cz^2 + dxy + eyz + fxz+ gx + hy + jz + k = 0`. Parameters ---------- @@ -1170,6 +1170,10 @@ class Cylinder(QuadricMixin, Surface): representing the quadric surface """ + pre_coeffs = self._get_base_coeffs + # infer transformation based on changes between two matrices + # perform same transformation on x0,y0,z0,dx,dy,dz if A matrix + # unchanged then it was a translation if not it was a rotation a, b, c, d, e, f, g, h, j, k = coeffs @classmethod From 777803c029d45499a4b2f40de8f0208d9830230f Mon Sep 17 00:00:00 2001 From: Ethan Peterson Date: Wed, 12 Feb 2020 14:40:12 -0500 Subject: [PATCH 14/33] implemented helper functions --- openmc/surface.py | 114 ++++++++++++++++++++++++++++++++++++++++------ 1 file changed, 100 insertions(+), 14 deletions(-) diff --git a/openmc/surface.py b/openmc/surface.py index badbd630c2..29ec06ae37 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -66,6 +66,7 @@ class Surface(IDManagerMixin, metaclass=ABCMeta): next_id = 1 used_ids = set() + _atol = 1.e-12 def __init__(self, surface_id=None, boundary_type='transmission', name=''): self.id = surface_id @@ -185,6 +186,30 @@ class Surface(IDManagerMixin, metaclass=ABCMeta): return memo[self] + def is_equal(self, coeffs): + """Determine if this Surface is the same as another to within a certain + tolerance + + Parameters + ---------- + coeffs : tuple + Tuple of surface coefficients to compare to the current Surface + object + + """ + def normalize(coeffs): + """Normalize coefficients by first nonzero value""" + coeffs = np.asarray(coeffs) + nonzeros = ~np.isclose(coeffs, 0., rtol=0., atol=self._atol) + norm_factor = coeffs[nonzeros][0] + return tuple([c/norm_factor for c in coeffs]) + + # get normalized sets of coefficients + coeffs = normalize(coeffs) + bcoeffs = normalize(self._get_base_coeffs()) + + return np.all(np.isclose(bcoeffs, coeffs, rtol=0., atol=self._atol)) + @abstractmethod def evaluate(self, point): pass @@ -480,11 +505,11 @@ class Plane(PlaneMixin, Surface): for k, v in argsdict.items(): warn(_WARNING_KWARGS.format(type(self).__name__, k), FutureWarning) + kwargs.update(argsdict) + for k in 'ABCD': kwargs.pop(k, None) - kwargs.update(argsdict) - super().__init__(**kwargs) for key, val in zip(self._coeff_keys, (a, b, c, d)): @@ -493,7 +518,7 @@ class Plane(PlaneMixin, Surface): for k, v in oldkwargs.items(): if k in 'ABCD': warn(_WARNING_UPPER.format(type(self).__name__, k.lower(), k), - FutureWarning) + FutureWarning) setattr(self, k.lower(), v) @property @@ -916,6 +941,38 @@ class QuadricMixin(metaclass=ABCMeta): return (np.array([-np.inf, -np.inf, -np.inf]), np.array([np.inf, np.inf, np.inf])) + def _A(self, coeffs=None): + """Helper function to generate the quadric matrix""" + if coeffs is None: + a, b, c, d, e, f, g, h, j, k = self._get_base_coeffs() + else: + a, b, c, d, e, f, g, h, j, k = coeffs + + return np.asarray([[a, d/2, f/2], [d/2, b, e/2], [f/2, e/2, c]]) + + def _b(self, coeffs=None): + """Helper function to generate the linear coefficients""" + if coeffs is None: + a, b, c, d, e, f, g, h, j, k = self._get_base_coeffs() + else: + a, b, c, d, e, f, g, h, j, k = coeffs + + return np.asarray([g, h, j]) + + def _c(self, coeffs=None): + """Helper function to generate the constant coefficient""" + if coeffs is None: + return self._get_base_coeffs()[-1] + return coeffs[-1] + + def eigh(self, coeffs=None): + """Wrapper method for returning eigenvalues and eigenvectors of this + quadric surface which is used for transformations. + + """ + A = self._A(coeffs=coeffs) + return np.linalg.eigh(A) + def bounding_box(self, side): """Determine an axis-aligned bounding box. @@ -981,13 +1038,22 @@ class QuadricMixin(metaclass=ABCMeta): Translated surface """ - vx, vy, vz = vector + vector = np.asarray(vector) + A = self._A() + bvec = self._b() + const = self._c() + a, b, c, d, e, f, g, h, j, k = self._get_base_coeffs() - k = (k + a*vx**2 + b*vy**2 + c*vz**2 + d*vx*vy + e*vy*vz + f*vx*vz - - g*vx - h*vy - j*vz) - g = g - 2*a*vx - d*vy - f*vz - h = h - 2*b*vy - d*vx - e*vz - j = j - 2*c*vz - e*vy - f*vx + g, h, j = bvec - 2*(vector.T @ A) + k = const + vector.T @ A @ vector - bvec.T @ vector + + #vx, vy, vz = vector + #a, b, c, d, e, f, g, h, j, k = self._get_base_coeffs() + #k = (k + a*vx**2 + b*vy**2 + c*vz**2 + d*vx*vy + e*vy*vz + f*vx*vz + # - g*vx - h*vy - j*vz) + #g = g - 2*a*vx - d*vy - f*vz + #h = h - 2*b*vy - d*vx - e*vz + #j = j - 2*c*vz - e*vy - f*vx if clone: surf = self.clone() @@ -1065,7 +1131,9 @@ class Cylinder(QuadricMixin, Surface): """ _type = 'cylinder' _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): + super().__init__(**kwargs) for key, val in zip(self._coeff_keys, (x0, y0, z0, r, dx, dy, dz)): @@ -1170,11 +1238,29 @@ class Cylinder(QuadricMixin, Surface): representing the quadric surface """ - pre_coeffs = self._get_base_coeffs - # infer transformation based on changes between two matrices - # perform same transformation on x0,y0,z0,dx,dy,dz if A matrix - # unchanged then it was a translation if not it was a rotation - a, b, c, d, e, f, g, h, j, k = coeffs + if self.is_equal(coeffs): + print("surfaces are identical") + return + + if np.isclose(self.Amat(), self.Amat(coeffs), rtol=0., atol=self._atol): + # transformation was translation + vec = np.asarray([self.x0, self.y0, self.z0]) + A = self.Amat() + bvec = self.bvec() + 2*np.dot(vec.T, A) + + # get eigenvalues and eigenvectors of quadric matrix + w, v = self.eigh() + w, u = self.eigh(coeffs=coeffs) + Rmat = np.dot(u.T, v) + xyz0 = np.asarray([self.x0, self.y0, self.z0]) + dxdydz0 = np.asarray([self.dx, self.dy, self.dz]) + # principal axis for cylinder is first eigenvector which has + # an eigenvalue of 0 + x0, y0, z0 = np.dot(Rmat, xyz0) + dx, dy, dz = np.dot(Rmat, dxdydz0) + r = self.r + for key, val in zip(self._coeff_keys, (x0, y0, z0, r, dx, dy, dz)): + setattr(self, key, val) @classmethod def from_points(cls, p1, p2, r=1., **kwargs): From 90d41aeabefd5df975fea5be1128206d084294eb Mon Sep 17 00:00:00 2001 From: Ethan Peterson Date: Wed, 12 Feb 2020 21:21:47 -0500 Subject: [PATCH 15/33] finished general cone transformation --- openmc/surface.py | 25 ++++++++++++++++++++----- 1 file changed, 20 insertions(+), 5 deletions(-) diff --git a/openmc/surface.py b/openmc/surface.py index 29ec06ae37..71112427ec 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -1239,7 +1239,6 @@ class Cylinder(QuadricMixin, Surface): """ if self.is_equal(coeffs): - print("surfaces are identical") return if np.isclose(self.Amat(), self.Amat(coeffs), rtol=0., atol=self._atol): @@ -1844,13 +1843,13 @@ class Cone(QuadricMixin, Surface): z-coordinate of the apex. Defaults to 0. r2 : float, optional Parameter related to the aperature. Defaults to 1. - u : float, optional + dx : float, optional x-component of the vector representing the axis of the cone. Defaults to 0. - v : float, optional + dy : float, optional y-component of the vector representing the axis of the cone. Defaults to 0. - w : float, optional + dz : float, optional z-component of the vector representing the axis of the cone. Defaults to 1. surface_id : int, optional @@ -1976,7 +1975,23 @@ class Cone(QuadricMixin, Surface): :math:`ax^2 + by^2 + cz^2 + dxy + eyz + fxz+ gx + hy + jz +k = 0`. """ - return (0., 0., 1., self.z0) + x0, y0, z0, r2 = self.x0, self.y0, self.z0, self.r2 + dx, dy, dz = self.dx, self.dy, self.dz + cos2 = 1 / (1 + r2) + c1 = (dx**2 + dy**2 + dz**2)*cos2 + + a = dx**2 - c1 + b = dy**2 - c1 + c = dz**2 - c1 + d = 2*dx*dy + e = 2*dy*dz + f = 2*dx*dz + g = -2*(dx**2*x0 + dx*dy*y0 + dx*dz*z0 - c1) + h = -2*(dy**2*y0 + dx*dy*x0 + dy*dz*z0 - c1) + j = -2*(dz**2*y0 + dx*dz*x0 + dy*dz*y0 - c1) + k = (dx*x0 + dy*y0 + dz*z0)**2 - c1*(x0**2 + y0**2 + z0**2) + + return (a, b, c, d, e, f, g, h, j, k) def _update_from_base_coeffs(self, coeffs): """Update the current surface from coefficients representing a general From 2919b2791aeb6f96c74f34135d02b9f786f05635 Mon Sep 17 00:00:00 2001 From: Ethan Peterson Date: Thu, 13 Feb 2020 12:06:01 -0500 Subject: [PATCH 16/33] finished implementation --- openmc/surface.py | 223 ++++++++++++++++++++++++++++------------------ 1 file changed, 137 insertions(+), 86 deletions(-) diff --git a/openmc/surface.py b/openmc/surface.py index 71112427ec..ba6ec15f79 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -1,8 +1,7 @@ from abc import ABCMeta, abstractmethod from collections import OrderedDict -from collections.abc import Iterable from copy import deepcopy -from numbers import Real, Integral +from numbers import Real from xml.etree import ElementTree as ET from warnings import warn @@ -186,29 +185,38 @@ class Surface(IDManagerMixin, metaclass=ABCMeta): return memo[self] - def is_equal(self, coeffs): - """Determine if this Surface is the same as another to within a certain - tolerance + def normalize(self, coeffs=None): + """Normalize coefficients by first nonzero value + Parameters + ---------- + coeffs : tuple, optional + Tuple of surface coefficients to normalize. If none are supplied, + the coefficients from the current surface will be used. + + Returns + ------- + tuple of normalized coefficients + + """ + coeffs = np.asarray(coeffs) + nonzeros = ~np.isclose(coeffs, 0., rtol=0., atol=self._atol) + norm_factor = coeffs[nonzeros][0] + return tuple([c/norm_factor for c in coeffs]) + + def is_equal(self, other): + """Determine if this Surface is equivalent to another Parameters ---------- - coeffs : tuple - Tuple of surface coefficients to compare to the current Surface - object + other : instance of openmc.Surface + Instance of openmc.Surface that should be compared to the current + surface """ - def normalize(coeffs): - """Normalize coefficients by first nonzero value""" - coeffs = np.asarray(coeffs) - nonzeros = ~np.isclose(coeffs, 0., rtol=0., atol=self._atol) - norm_factor = coeffs[nonzeros][0] - return tuple([c/norm_factor for c in coeffs]) + coeffs1 = normalize(self._get_base_coeffs()) + coeffs2 = normalize(other._get_base_coeffs()) - # get normalized sets of coefficients - coeffs = normalize(coeffs) - bcoeffs = normalize(self._get_base_coeffs()) - - return np.all(np.isclose(bcoeffs, coeffs, rtol=0., atol=self._atol)) + return np.all(np.isclose(coeffs1, coeffs2, rtol=0., atol=self._atol)) @abstractmethod def evaluate(self, point): @@ -299,7 +307,7 @@ class Surface(IDManagerMixin, metaclass=ABCMeta): class PlaneMixin(metaclass=ABCMeta): - """A Plane Meta class for all operations on order 1 surfaces""" + """A Plane mixin class for all operations on order 1 surfaces""" def __init__(self, **kwargs): super().__init__(**kwargs) self._periodic_surface = None @@ -941,37 +949,47 @@ class QuadricMixin(metaclass=ABCMeta): return (np.array([-np.inf, -np.inf, -np.inf]), np.array([np.inf, np.inf, np.inf])) - def _A(self, coeffs=None): - """Helper function to generate the quadric matrix""" + def get_Abc(self, coeffs=None): + """Compute matrix, vector, and scalar coefficients for this surface or + for a specified set of coefficients. + + Parameters + ---------- + coeffs : tuple, optional + Tuple of coefficients from which to compute the quadric elements. + If none are supplied the coefficients of this surface will be used. + """ if coeffs is None: a, b, c, d, e, f, g, h, j, k = self._get_base_coeffs() else: a, b, c, d, e, f, g, h, j, k = coeffs - return np.asarray([[a, d/2, f/2], [d/2, b, e/2], [f/2, e/2, c]]) + A = np.asarray([[a, d/2, f/2], [d/2, b, e/2], [f/2, e/2, c]]) + bvec = np.asarray([g, h, j]) - def _b(self, coeffs=None): - """Helper function to generate the linear coefficients""" - if coeffs is None: - a, b, c, d, e, f, g, h, j, k = self._get_base_coeffs() - else: - a, b, c, d, e, f, g, h, j, k = coeffs - - return np.asarray([g, h, j]) - - def _c(self, coeffs=None): - """Helper function to generate the constant coefficient""" - if coeffs is None: - return self._get_base_coeffs()[-1] - return coeffs[-1] + return A, bvec, k def eigh(self, coeffs=None): """Wrapper method for returning eigenvalues and eigenvectors of this quadric surface which is used for transformations. + Parameters + ---------- + coeffs : tuple, optional + Tuple of coefficients from which to compute the quadric elements. + If none are supplied the coefficients of this surface will be used. + + Returns + ------- + w, v : tuple of numpy arrays with shapes (3,) and (3,3) respectively + Returns the eigenvalues and eigenvectors of the quadric matrix A + that represents the supplied coefficients. The vector w contains + the eigenvalues in ascending order and the matrix v contains the + eigenvectors such that v[:,i] is the eigenvector corresponding to + the eigenvalue w[i]. + """ - A = self._A(coeffs=coeffs) - return np.linalg.eigh(A) + return np.linalg.eigh(self.get_Abc(coeffs=coeffs)[0]) def bounding_box(self, side): """Determine an axis-aligned bounding box. @@ -1017,9 +1035,9 @@ class QuadricMixin(metaclass=ABCMeta): Jz' + K = 0` """ - x, y, z = point - a, b, c, d, e, f, g, h, j, k = self._get_base_coeffs() - return x*(a*x + d*y + g) + y*(b*y + e*z + h) + z*(c*z + f*x + j) + k + 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 def translate(self, vector, clone=True): """Translate surface in given direction @@ -1039,21 +1057,12 @@ class QuadricMixin(metaclass=ABCMeta): """ vector = np.asarray(vector) - A = self._A() - bvec = self._b() - const = self._c() + A, bvec, cnst = self.get_Abc() - a, b, c, d, e, f, g, h, j, k = self._get_base_coeffs() - g, h, j = bvec - 2*(vector.T @ A) - k = const + vector.T @ A @ vector - bvec.T @ vector - - #vx, vy, vz = vector - #a, b, c, d, e, f, g, h, j, k = self._get_base_coeffs() - #k = (k + a*vx**2 + b*vy**2 + c*vz**2 + d*vx*vy + e*vy*vz + f*vx*vz - # - g*vx - h*vy - j*vz) - #g = g - 2*a*vx - d*vy - f*vz - #h = h - 2*b*vy - d*vx - e*vz - #j = j - 2*c*vz - e*vy - f*vx + a, b, c, d, e, f = self._get_base_coeffs()[0:6] + 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) if clone: surf = self.clone() @@ -1210,19 +1219,19 @@ class Cylinder(QuadricMixin, Surface): """ x0, y0, z0, r = self.x0, self.y0, self.z0, self.r dx, dy, dz = self.dx, self.dy, self.dz + dx2, dy2, dz2 = dx**2, dy**2, dz**2 - a = dy**2 + dz**2 - b = dx**2 + dz**2 - c = dx**2 + dy**2 + a = dy2 + dz2 + b = dx2 + dz2 + c = dx2 + dy2 d = -2*dx*dy e = -2*dy*dz f = -2*dx*dz - g = 2*(dx*(z0*dz + y0*dy) - x0*(dy**2 + dz**2)) - h = 2*(dy*(z0*dz + x0*dx) - y0*(dx**2 + dz**2)) - j = 2*(dz*(y0*dy + x0*dx) - z0*(dx**2 + dy**2)) - k = x0**2*(dy**2 + dz**2) + y0**2*(dx**2 + dz**2) \ - + z0**2*(dx**2 + dy**2) - 2*dy*dz*y0*z0 \ - - 2*dx*dz*x0*z0 - 2*dx*dy*x0*y0 - r**2*(dx**2 + dy**2 + dz**2) + g = 2*(dx*(z0*dz + y0*dy) - x0*(dy2 + dz2)) + h = 2*(dy*(z0*dz + x0*dx) - y0*(dx2 + dz2)) + j = 2*(dz*(y0*dy + x0*dx) - z0*(dx2 + dy2)) + k = x0**2*(dy2 + dz2) + y0**2*(dx2 + dz2) + z0**2*(dx2 + dy2) \ + + e*y0*z0 + f*x0*z0 + d*x0*y0 - r**2*(dx2 + dy2 + dz2) return (a, b, c, d, e, f, g, h, j, k) @@ -1238,28 +1247,38 @@ class Cylinder(QuadricMixin, Surface): representing the quadric surface """ - if self.is_equal(coeffs): - return + raise NotImplementedError('_update_from_base_coeffs for {} is ' + 'undefined'.format(self.__cls__.__name__)) - if np.isclose(self.Amat(), self.Amat(coeffs), rtol=0., atol=self._atol): - # transformation was translation - vec = np.asarray([self.x0, self.y0, self.z0]) - A = self.Amat() - bvec = self.bvec() + 2*np.dot(vec.T, A) + def translate(self, vector, clone=True): + """Translate surface in given direction + + Parameters + ---------- + vector : iterable of float + Direction in which surface should be translated + clone : bool + Whether to return a clone of the Surface or the Surface itself. + Defaults to True + + Returns + ------- + openmc.Surface + Translated surface + + """ + + if clone: + surf = self.clone() + else: + surf = self + + surf.x0 += vector[0] + surf.y0 += vector[1] + surf.z0 += vector[2] + + return surf - # get eigenvalues and eigenvectors of quadric matrix - w, v = self.eigh() - w, u = self.eigh(coeffs=coeffs) - Rmat = np.dot(u.T, v) - xyz0 = np.asarray([self.x0, self.y0, self.z0]) - dxdydz0 = np.asarray([self.dx, self.dy, self.dz]) - # principal axis for cylinder is first eigenvector which has - # an eigenvalue of 0 - x0, y0, z0 = np.dot(Rmat, xyz0) - dx, dy, dz = np.dot(Rmat, dxdydz0) - r = self.r - for key, val in zip(self._coeff_keys, (x0, y0, z0, r, dx, dy, dz)): - setattr(self, key, val) @classmethod def from_points(cls, p1, p2, r=1., **kwargs): @@ -1797,6 +1816,7 @@ class Sphere(QuadricMixin, Surface): d = e = f = 0. g, h, j = -2*x0, -2*y0, -2*z0 k = x0**2 + y0**2 + z0**2 - r**2 + return (a, b, c, d, e, f, g, h, j, k) def _update_from_base_coeffs(self, coeffs): @@ -1814,6 +1834,7 @@ class Sphere(QuadricMixin, Surface): a, b, c, d, e, f, g, h, j, k = coeffs x0, y0, z0 = -g / 2, -h / 2, -j / 2 r = np.sqrt(x0**2 + y0**2 + z0**2 - k) + for key, val in zip(self._coeff_keys, (x0, y0, z0, r)): setattr(self, key, val) @@ -2005,7 +2026,37 @@ class Cone(QuadricMixin, Surface): representing the quadric surface """ - a, b, c, d, e, f, g, h, j, k = coeffs + raise NotImplementedError('_update_from_base_coeffs for {} is ' + 'undefined'.format(self.__cls__.__name__)) + + def translate(self, vector, clone=True): + """Translate surface in given direction + + Parameters + ---------- + vector : iterable of float + Direction in which surface should be translated + clone : bool + Whether to return a clone of the Surface or the Surface itself. + Defaults to True + + Returns + ------- + openmc.Surface + Translated surface + + """ + + if clone: + surf = self.clone() + else: + surf = self + + surf.x0 += vector[0] + surf.y0 += vector[1] + surf.z0 += vector[2] + + return surf class XCone(QuadricMixin, Surface): From f1b14c95e2e1968b803765150d59836c0df70c28 Mon Sep 17 00:00:00 2001 From: Ethan Peterson Date: Thu, 13 Feb 2020 13:43:42 -0500 Subject: [PATCH 17/33] implemented individual translate and evaulate methods and removed extraneous docstrings --- openmc/surface.py | 711 +++++------------- .../cmfd_feed_2g/results_true.dat | 2 +- tests/regression_tests/triso/results_true.dat | 2 +- 3 files changed, 202 insertions(+), 513 deletions(-) diff --git a/openmc/surface.py b/openmc/surface.py index ba6ec15f79..064ac7242f 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -7,7 +7,7 @@ from warnings import warn import numpy as np -from openmc.checkvalue import check_type, check_value, check_length +from openmc.checkvalue import check_type, check_value from openmc.region import Region, Intersection, Union from openmc.mixin import IDManagerMixin @@ -330,29 +330,6 @@ class PlaneMixin(metaclass=ABCMeta): """ pass - @abstractmethod - def _update_from_base_coeffs(self, coeffs): - """Update the current plane from coefficients representing a general - plane of the form :math:`ax + by + cz = d`. - - Parameters - ---------- - coeffs : tuple - Tuple of general coefficients (a, b, c, d) representing the plane - - """ - pass - - def _neg_bounds(self): - """Return the lower and upper bounds of the negative half space""" - return (np.array([-np.inf, -np.inf, -np.inf]), - np.array([np.inf, np.inf, np.inf])) - - def _pos_bounds(self): - """Return the lower and upper bounds of the positive half space""" - return (np.array([-np.inf, -np.inf, -np.inf]), - np.array([np.inf, np.inf, np.inf])) - def evaluate(self, point): """Evaluate the surface equation at a given point. @@ -373,6 +350,7 @@ class PlaneMixin(metaclass=ABCMeta): a, b, c, d = self._get_base_coeffs() return a*x + b*y + c*z - d + @abstractmethod def translate(self, vector, clone=True): """Translate surface in given direction @@ -390,46 +368,7 @@ class PlaneMixin(metaclass=ABCMeta): Translated surface """ - vx, vy, vz = vector - a, b, c, d = self._get_base_coeffs() - d = d + a*vx + b*vy + c*vz - - if clone: - surf = self.clone() - else: - surf = self - - surf._update_from_base_coeffs((a, b, c, d)) - return surf - - 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.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 - - """ - - if side == '-': - return self._neg_bounds() - elif side == '+': - return self._pos_bounds() + pass def to_xml_element(self): """Return XML representation of the surface @@ -566,25 +505,18 @@ class Plane(PlaneMixin, Surface): self._coefficients['d'] = d def _get_base_coeffs(self): - """Return coefficients a, b, c, d representing a general plane of the - form :math:`ax + by + cz = d`. - - """ return (self.a, self.b, self.c, self.d) - def _update_from_base_coeffs(self, coeffs): - """Update the current plane from coefficients representing a general - plane of the form :math:`ax + by + cz = d`. - - Parameters - ---------- - coeffs : tuple - Tuple of general coefficients (a, b, c, d) representing the plane - - """ - - for key, val in zip(self._coeff_keys, coeffs): - setattr(self, key, val) + def translate(self, vector, clone=True): + vx, vy, vz = vector + a, b, c, d = self._get_base_coeffs() + d = d + a*vx + b*vy + c*vz + if clone: + surf = self.clone() + else: + surf = self + surf.d = d + return surf @classmethod def from_points(cls, p1, p2, p3, **kwargs): @@ -683,34 +615,26 @@ class XPlane(PlaneMixin, Surface): self._coefficients['x0'] = x0 def _get_base_coeffs(self): - """Return coefficients a, b, c, d representing a general plane of the - form :math:`ax + by + cz = d`. - - """ return (1., 0., 0., self.x0) - def _update_from_base_coeffs(self, coeffs): - """Update the current plane from coefficients representing a general - plane of the form :math:`ax + by + cz = d`. + 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])) - Parameters - ---------- - coeffs : tuple - Tuple of general coefficients (a, b, c, d) representing the plane + def evaluate(self, point): + return point[0] - self.x0 - """ - a, b, c, d = coeffs - self.x0 = d / a - - def _neg_bounds(self): - """Return the lower and upper bounds of the negative half space""" - return (np.array([-np.inf, -np.inf, -np.inf]), - np.array([self.x0, np.inf, np.inf])) - - def _pos_bounds(self): - """Return the lower and upper bounds of the positive half space""" - return (np.array([self.x0, -np.inf, -np.inf]), - np.array([np.inf, np.inf, np.inf])) + def translate(self, vector, clone=True): + if clone: + surf = self.clone() + else: + surf = self + surf.x0 += vector[0] + return surf Plane.register(XPlane) @@ -780,34 +704,26 @@ class YPlane(PlaneMixin, Surface): self._coefficients['y0'] = y0 def _get_base_coeffs(self): - """Return coefficients a, b, c, d representing a general plane of the - form :math:`ax + by + cz = d`. - - """ return (0., 1., 0., self.y0) - def _update_from_base_coeffs(self, coeffs): - """Update the current plane from coefficients representing a general - plane of the form :math:`ax + by + cz = d`. + 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])) - Parameters - ---------- - coeffs : tuple - Tuple of general coefficients (a, b, c, d) representing the plane + def evaluate(self, point): + return point[1] - self.y0 - """ - a, b, c, d = coeffs - self.y0 = d / b - - def _neg_bounds(self): - """Return the lower and upper bounds of the negative half space""" - return (np.array([-np.inf, -np.inf, -np.inf]), - np.array([np.inf, self.y0, np.inf])) - - def _pos_bounds(self): - """Return the lower and upper bounds of the positive half space""" - return (np.array([-np.inf, self.y0, -np.inf]), - np.array([np.inf, np.inf, np.inf])) + def translate(self, vector, clone=True): + if clone: + surf = self.clone() + else: + surf = self + surf.y0 += vector[1] + return surf Plane.register(YPlane) @@ -877,34 +793,26 @@ class ZPlane(PlaneMixin, Surface): self._coefficients['z0'] = z0 def _get_base_coeffs(self): - """Return coefficients a, b, c, d, e, f, g, h, j, k representing a - general plane form :math:`ax + by + cz = d`. - - """ return (0., 0., 1., self.z0) - def _update_from_base_coeffs(self, coeffs): - """Update the current plane from coefficients representing a general - plane of the form :math:`ax + by + cz = d`. + 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])) - Parameters - ---------- - coeffs : tuple - Tuple of general coefficients (a, b, c, d) representing the plane + def evaluate(self, point): + return point[2] - self.z0 - """ - a, b, c, d = coeffs - self.z0 = d / c - - def _neg_bounds(self): - """Return the lower and upper bounds of the negative half space""" - return (np.array([-np.inf, -np.inf, -np.inf]), - np.array([np.inf, np.inf, self.z0])) - - def _pos_bounds(self): - """Return the lower and upper bounds of the positive half space""" - return (np.array([-np.inf, -np.inf, self.z0]), - np.array([np.inf, np.inf, np.inf])) + def translate(self, vector, clone=True): + if clone: + surf = self.clone() + else: + surf = self + surf.z0 += vector[2] + return surf Plane.register(ZPlane) @@ -924,31 +832,6 @@ class QuadricMixin(metaclass=ABCMeta): """ pass - @abstractmethod - def _update_from_base_coeffs(self, coeffs): - """Update the current surface from coefficients representing a general - quadric surface of the form: - :math:`ax^2 + by^2 + cz^2 + dxy + eyz + fxz+ gx + hy + jz +k = 0`. - - Parameters - ---------- - coeffs : tuple - Tuple of general coefficients (a, b, c, d, e, f, g, h, j, k) - representing the quadric surface - - """ - pass - - def _neg_bounds(self): - """Return the lower and upper bounds of the negative half space""" - return (np.array([-np.inf, -np.inf, -np.inf]), - np.array([np.inf, np.inf, np.inf])) - - def _pos_bounds(self): - """Return the lower and upper bounds of the positive half space""" - return (np.array([-np.inf, -np.inf, -np.inf]), - np.array([np.inf, np.inf, np.inf])) - def get_Abc(self, coeffs=None): """Compute matrix, vector, and scalar coefficients for this surface or for a specified set of coefficients. @@ -991,34 +874,6 @@ class QuadricMixin(metaclass=ABCMeta): """ return np.linalg.eigh(self.get_Abc(coeffs=coeffs)[0]) - 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. To represent infinity, - numpy.inf is used. - - 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 - - """ - - if side == '-': - return self._neg_bounds() - elif side == '+': - return self._pos_bounds() - def evaluate(self, point): """Evaluate the surface equation at a given point. @@ -1039,6 +894,7 @@ class QuadricMixin(metaclass=ABCMeta): A, b, c = self.get_Abc() return np.matmul(x.T, np.matmul(A, x)) + np.matmul(b.T, x) + c + @abstractmethod def translate(self, vector, clone=True): """Translate surface in given direction @@ -1056,22 +912,7 @@ class QuadricMixin(metaclass=ABCMeta): Translated surface """ - vector = np.asarray(vector) - A, bvec, cnst = self.get_Abc() - - a, b, c, d, e, f = self._get_base_coeffs()[0:6] - 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) - - if clone: - surf = self.clone() - else: - surf = self - - surf._update_from_base_coeffs((a, b, c, d, e, f, g, h, j, k)) - - return surf + pass class Cylinder(QuadricMixin, Surface): @@ -1212,11 +1053,6 @@ class Cylinder(QuadricMixin, Surface): self._coefficients['dz'] = dz def _get_base_coeffs(self): - """Return coefficients a, b, c, d, e, f, g, h, j, k representing a - general quadric surface of the form: - :math:`ax^2 + by^2 + cz^2 + dxy + eyz + fxz+ gx + hy + jz + k = 0`. - - """ x0, y0, z0, r = self.x0, self.y0, self.z0, self.r dx, dy, dz = self.dx, self.dy, self.dz dx2, dy2, dz2 = dx**2, dy**2, dz**2 @@ -1235,39 +1071,7 @@ class Cylinder(QuadricMixin, Surface): return (a, b, c, d, e, f, g, h, j, k) - def _update_from_base_coeffs(self, coeffs): - """Update the current surface from coefficients representing a general - quadric surface of the form: - :math:`ax^2 + by^2 + cz^2 + dxy + eyz + fxz+ gx + hy + jz + k = 0`. - - Parameters - ---------- - coeffs : tuple - Tuple of general coefficients (a, b, c, d, e, f, g, h, j, k) - representing the quadric surface - - """ - raise NotImplementedError('_update_from_base_coeffs for {} is ' - 'undefined'.format(self.__cls__.__name__)) - def translate(self, vector, clone=True): - """Translate surface in given direction - - Parameters - ---------- - vector : iterable of float - Direction in which surface should be translated - clone : bool - Whether to return a clone of the Surface or the Surface itself. - Defaults to True - - Returns - ------- - openmc.Surface - Translated surface - - """ - if clone: surf = self.clone() else: @@ -1279,7 +1083,6 @@ class Cylinder(QuadricMixin, Surface): return surf - @classmethod def from_points(cls, p1, p2, r=1., **kwargs): """Return a cylinder given points that define the axis and a radius. @@ -1396,11 +1199,6 @@ class XCylinder(QuadricMixin, Surface): self._coefficients['r'] = r def _get_base_coeffs(self): - """Return coefficients a, b, c, d, e, f, g, h, j, k representing a - general quadric surface of the form: - :math:`ax^2 + by^2 + cz^2 + dxy + eyz + fxz+ gx + hy + jz +k = 0`. - - """ y0, z0, r = self.y0, self.z0, self.r a = d = e = f = g = 0. @@ -1409,34 +1207,27 @@ class XCylinder(QuadricMixin, Surface): return (a, b, c, d, e, f, g, h, j, k) - def _update_from_base_coeffs(self, coeffs): - """Update the current surface from coefficients representing a general - quadric surface of the form: - :math:`ax^2 + by^2 + cz^2 + dxy + eyz + fxz+ gx + hy + jz +k = 0`. + 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])) + elif side == '+': + return (np.array([-np.inf, -np.inf, -np.inf]), + np.array([np.inf, np.inf, np.inf])) - Parameters - ---------- - coeffs : tuple - Tuple of general coefficients (a, b, c, d, e, f, g, h, j, k) - representing the quadric surface + def evaluate(self, point): + y = point[1] - self.y0 + z = point[2] - self.z0 + return y**2 + z**2 - self.r**2 - """ - a, b, c, d, e, f, g, h, j, k = coeffs - y0, z0 = -h / 2, -j / 2 - r = np.sqrt(y0**2 + z0**2 - k) - - for key, val in zip(self._coeff_keys, (y0, z0, r)): - setattr(self, key, val) - - def _neg_bounds(self): - """Return the lower and upper bounds of the negative half space""" - 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])) - - def _pos_bounds(self): - """Return the lower and upper bounds of the positive half space""" - return (np.array([-np.inf, -np.inf, -np.inf]), - np.array([np.inf, np.inf, np.inf])) + def translate(self, vector, clone=True): + if clone: + surf = self.clone() + else: + surf = self + surf.y0 += vector[1] + surf.z0 += vector[2] + return surf Cylinder.register(XCylinder) @@ -1529,11 +1320,6 @@ class YCylinder(QuadricMixin, Surface): self._coefficients['r'] = r def _get_base_coeffs(self): - """Return coefficients a, b, c, d, e, f, g, h, j, k representing a - general quadric surface of the form: - :math:`ax^2 + by^2 + cz^2 + dxy + eyz + fxz+ gx + hy + jz +k = 0`. - - """ x0, z0, r = self.x0, self.z0, self.r b = d = e = f = h = 0. @@ -1542,34 +1328,27 @@ class YCylinder(QuadricMixin, Surface): return (a, b, c, d, e, f, g, h, j, k) - def _update_from_base_coeffs(self, coeffs): - """Update the current surface from coefficients representing a general - quadric surface of the form: - :math:`ax^2 + by^2 + cz^2 + dxy + eyz + fxz+ gx + hy + jz +k = 0`. + 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])) + elif side == '+': + return (np.array([-np.inf, -np.inf, -np.inf]), + np.array([np.inf, np.inf, np.inf])) - Parameters - ---------- - coeffs : tuple - Tuple of general coefficients (a, b, c, d, e, f, g, h, j, k) - representing the quadric surface + def evaluate(self, point): + x = point[0] - self.x0 + z = point[2] - self.z0 + return x**2 + z**2 - self.r**2 - """ - a, b, c, d, e, f, g, h, j, k = coeffs - x0, z0 = -g / 2, -j / 2 - r = np.sqrt(x0**2 + z0**2 - k) - - for key, val in zip(self._coeff_keys, (x0, z0, r)): - setattr(self, key, val) - - def _neg_bounds(self): - """Return the lower and upper bounds of the negative half space""" - 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])) - - def _pos_bounds(self): - """Return the lower and upper bounds of the positive half space""" - return (np.array([-np.inf, -np.inf, -np.inf]), - np.array([np.inf, np.inf, np.inf])) + def translate(self, vector, clone=True): + if clone: + surf = self.clone() + else: + surf = self + surf.x0 += vector[0] + surf.z0 += vector[2] + return surf Cylinder.register(YCylinder) @@ -1662,11 +1441,6 @@ class ZCylinder(QuadricMixin, Surface): self._coefficients['r'] = r def _get_base_coeffs(self): - """Return coefficients a, b, c, d, e, f, g, h, j, k representing a - general quadric surface of the form: - :math:`ax^2 + by^2 + cz^2 + dxy + eyz + fxz+ gx + hy + jz +k = 0`. - - """ x0, y0, r = self.x0, self.y0, self.r c = d = e = f = j = 0. @@ -1675,34 +1449,27 @@ class ZCylinder(QuadricMixin, Surface): return (a, b, c, d, e, f, g, h, j, k) - def _update_from_base_coeffs(self, coeffs): - """Update the current surface from coefficients representing a general - quadric surface of the form: - :math:`ax^2 + by^2 + cz^2 + dxy + eyz + fxz+ gx + hy + jz +k = 0`. + 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])) + elif side == '+': + return (np.array([-np.inf, -np.inf, -np.inf]), + np.array([np.inf, np.inf, np.inf])) - Parameters - ---------- - coeffs : tuple - Tuple of general coefficients (a, b, c, d, e, f, g, h, j, k) - representing the quadric surface + def evaluate(self, point): + x = point[0] - self.x0 + y = point[1] - self.y0 + return x**2 + y**2 - self.r**2 - """ - a, b, c, d, e, f, g, h, j, k = coeffs - x0, y0 = -g / 2, -h / 2 - r = np.sqrt(x0**2 + y0**2 - k) - - for key, val in zip(self._coeff_keys, (x0, y0, r)): - setattr(self, key, val) - - def _neg_bounds(self): - """Return the lower and upper bounds of the negative half space""" - 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])) - - def _pos_bounds(self): - """Return the lower and upper bounds of the positive half space""" - return (np.array([-np.inf, -np.inf, -np.inf]), - np.array([np.inf, np.inf, np.inf])) + def translate(self, vector, clone=True): + if clone: + surf = self.clone() + else: + surf = self + surf.x0 += vector[0] + surf.y0 += vector[1] + return surf Cylinder.register(ZCylinder) @@ -1806,11 +1573,6 @@ class Sphere(QuadricMixin, Surface): self._coefficients['r'] = r def _get_base_coeffs(self): - """Return coefficients a, b, c, d, e, f, g, h, j, k representing a - general quadric surface of the form: - :math:`ax^2 + by^2 + cz^2 + dxy + eyz + fxz+ gx + hy + jz +k = 0`. - - """ x0, y0, z0, r = self.x0, self.y0, self.z0, self.r a = b = c = 1. d = e = f = 0. @@ -1819,36 +1581,31 @@ class Sphere(QuadricMixin, Surface): return (a, b, c, d, e, f, g, h, j, k) - def _update_from_base_coeffs(self, coeffs): - """Update the current surface from coefficients representing a general - quadric surface of the form: - :math:`ax^2 + by^2 + cz^2 + dxy + eyz + fxz+ gx + hy + jz +k = 0`. + 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])) + elif side == '+': + return (np.array([-np.inf, -np.inf, -np.inf]), + np.array([np.inf, np.inf, np.inf])) - Parameters - ---------- - coeffs : tuple - Tuple of general coefficients (a, b, c, d, e, f, g, h, j, k) - representing the quadric surface + def evaluate(self, point): + x = point[0] - self.x0 + y = point[1] - self.y0 + z = point[2] - self.z0 + return x**2 + y**2 + z**2 - self.r**2 - """ - a, b, c, d, e, f, g, h, j, k = coeffs - x0, y0, z0 = -g / 2, -h / 2, -j / 2 - r = np.sqrt(x0**2 + y0**2 + z0**2 - k) - - for key, val in zip(self._coeff_keys, (x0, y0, z0, r)): - setattr(self, key, val) - - def _neg_bounds(self): - """Return the lower and upper bounds of the negative half space""" - 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])) - - def _pos_bounds(self): - """Return the lower and upper bounds of the positive half space""" - return (np.array([-np.inf, -np.inf, -np.inf]), - np.array([np.inf, np.inf, np.inf])) + def translate(self, vector, clone=True): + if clone: + surf = self.clone() + else: + surf = self + surf.x0 += vector[0] + surf.y0 += vector[1] + surf.z0 += vector[2] + return surf class Cone(QuadricMixin, Surface): @@ -1991,11 +1748,6 @@ class Cone(QuadricMixin, Surface): self._coefficients['dz'] = dz def _get_base_coeffs(self): - """Return coefficients a, b, c, d, e, f, g, h, j, k representing a - general quadric surface of the form: - :math:`ax^2 + by^2 + cz^2 + dxy + eyz + fxz+ gx + hy + jz +k = 0`. - - """ x0, y0, z0, r2 = self.x0, self.y0, self.z0, self.r2 dx, dy, dz = self.dx, self.dy, self.dz cos2 = 1 / (1 + r2) @@ -2014,48 +1766,14 @@ class Cone(QuadricMixin, Surface): return (a, b, c, d, e, f, g, h, j, k) - def _update_from_base_coeffs(self, coeffs): - """Update the current surface from coefficients representing a general - quadric surface of the form: - :math:`ax^2 + by^2 + cz^2 + dxy + eyz + fxz+ gx + hy + jz +k = 0`. - - Parameters - ---------- - coeffs : tuple - Tuple of general coefficients (a, b, c, d, e, f, g, h, j, k) - representing the quadric surface - - """ - raise NotImplementedError('_update_from_base_coeffs for {} is ' - 'undefined'.format(self.__cls__.__name__)) - def translate(self, vector, clone=True): - """Translate surface in given direction - - Parameters - ---------- - vector : iterable of float - Direction in which surface should be translated - clone : bool - Whether to return a clone of the Surface or the Surface itself. - Defaults to True - - Returns - ------- - openmc.Surface - Translated surface - - """ - if clone: surf = self.clone() else: surf = self - surf.x0 += vector[0] surf.y0 += vector[1] surf.z0 += vector[2] - return surf @@ -2158,11 +1876,6 @@ class XCone(QuadricMixin, Surface): self._coefficients['r2'] = r2 def _get_base_coeffs(self): - """Return coefficients a, b, c, d, e, f, g, h, j, k representing a - general quadric surface of the form: - :math:`ax^2 + by^2 + cz^2 + dxy + eyz + fxz+ gx + hy + jz +k = 0`. - - """ x0, y0, z0, r2 = self.x0, self.y0, self.z0, self.r2 a = -r2 @@ -2173,25 +1886,21 @@ class XCone(QuadricMixin, Surface): return (a, b, c, d, e, f, g, h, j, k) - def _update_from_base_coeffs(self, coeffs): - """Update the current surface from coefficients representing a general - quadric surface of the form: - :math:`ax^2 + by^2 + cz^2 + dxy + eyz + fxz+ gx + hy + jz +k = 0`. + def evaluate(self, point): + x = point[0] - self.x0 + y = point[1] - self.y0 + z = point[2] - self.z0 + return y**2 + z**2 - self.r2*x**2 - Parameters - ---------- - coeffs : tuple - Tuple of general coefficients (a, b, c, d, e, f, g, h, j, k) - representing the quadric surface - - """ - a, b, c, d, e, f, g, h, j, k = coeffs - - r2 = -a - x0, y0, z0 = g / (2*r2), -h / 2, -j / 2 - - for key, val in zip(self._coeff_keys, (x0, y0, z0, r2)): - setattr(self, key, val) + def translate(self, vector, clone=True): + if clone: + surf = self.clone() + else: + surf = self + surf.x0 += vector[0] + surf.y0 += vector[1] + surf.z0 += vector[2] + return surf Cone.register(XCone) @@ -2296,11 +2005,6 @@ class YCone(QuadricMixin, Surface): self._coefficients['r2'] = r2 def _get_base_coeffs(self): - """Return coefficients a, b, c, d, e, f, g, h, j, k representing a - general quadric surface of the form: - :math:`ax^2 + by^2 + cz^2 + dxy + eyz + fxz+ gx + hy + jz +k = 0`. - - """ x0, y0, z0, r2 = self.x0, self.y0, self.z0, self.r2 b = -r2 @@ -2311,25 +2015,21 @@ class YCone(QuadricMixin, Surface): return (a, b, c, d, e, f, g, h, j, k) - def _update_from_base_coeffs(self, coeffs): - """Update the current surface from coefficients representing a general - quadric surface of the form: - :math:`ax^2 + by^2 + cz^2 + dxy + eyz + fxz+ gx + hy + jz +k = 0`. + def evaluate(self, point): + x = point[0] - self.x0 + y = point[1] - self.y0 + z = point[2] - self.z0 + return x**2 + z**2 - self.r2*y**2 - Parameters - ---------- - coeffs : tuple - Tuple of general coefficients (a, b, c, d, e, f, g, h, j, k) - representing the quadric surface - - """ - a, b, c, d, e, f, g, h, j, k = coeffs - - r2 = -b - x0, y0, z0 = -g / 2, h / (2*r2), -j / 2 - - for key, val in zip(self._coeff_keys, (x0, y0, z0, r2)): - setattr(self, key, val) + def translate(self, vector, clone=True): + if clone: + surf = self.clone() + else: + surf = self + surf.x0 += vector[0] + surf.y0 += vector[1] + surf.z0 += vector[2] + return surf Cone.register(YCone) @@ -2434,11 +2134,6 @@ class ZCone(QuadricMixin, Surface): self._coefficients['r2'] = r2 def _get_base_coeffs(self): - """Return coefficients a, b, c, d, e, f, g, h, j, k representing a - general quadric surface of the form: - :math:`ax^2 + by^2 + cz^2 + dxy + eyz + fxz+ gx + hy + jz +k = 0`. - - """ x0, y0, z0, r2 = self.x0, self.y0, self.z0, self.r2 c = -r2 @@ -2449,25 +2144,21 @@ class ZCone(QuadricMixin, Surface): return (a, b, c, d, e, f, g, h, j, k) - def _update_from_base_coeffs(self, coeffs): - """Update the current surface from coefficients representing a general - quadric surface of the form: - :math:`ax^2 + by^2 + cz^2 + dxy + eyz + fxz+ gx + hy + jz +k = 0`. + def evaluate(self, point): + x = point[0] - self.x0 + y = point[1] - self.y0 + z = point[2] - self.z0 + return x**2 + y**2 - self.r2*z**2 - Parameters - ---------- - coeffs : tuple - Tuple of general coefficients (a, b, c, d, e, f, g, h, j, k) - representing the quadric surface - - """ - a, b, c, d, e, f, g, h, j, k = coeffs - - r2 = -c - x0, y0, z0 = -g / 2, -h / 2, j / (2*r2) - - for key, val in zip(self._coeff_keys, (x0, y0, z0, r2)): - setattr(self, key, val) + def translate(self, vector, clone=True): + if clone: + surf = self.clone() + else: + surf = self + surf.x0 += vector[0] + surf.y0 += vector[1] + surf.z0 += vector[2] + return surf Cone.register(ZCone) @@ -2611,27 +2302,25 @@ class Quadric(QuadricMixin, Surface): self._coefficients['k'] = k def _get_base_coeffs(self): - """Return coefficients a, b, c, d, e, f, g, h, j, k representing a - general quadric surface of the form: - :math:`ax^2 + by^2 + cz^2 + dxy + eyz + fxz+ gx + hy + jz +k = 0`. - - """ return tuple(getattr(self, c) for c in self._coeff_keys) - def _update_from_base_coeffs(self, coeffs): - """Update the current surface from coefficients representing a general - quadric surface of the form: - :math:`ax^2 + by^2 + cz^2 + dxy + eyz + fxz+ gx + hy + jz +k = 0`. + def translate(self, vector, clone=True): + vector = np.asarray(vector) + A, bvec, cnst = self.get_Abc() - Parameters - ---------- - coeffs : tuple - Tuple of general coefficients (a, b, c, d, e, f, g, h, j, k) - representing the quadric surface + 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) - """ - for key, val in zip(self._coeff_keys, coeffs): - setattr(self, key, val) + if clone: + surf = self.clone() + else: + surf = self + + for key, val in zip(('g', 'h', 'j', 'k'), (g, h, j, k)): + setattr(surf, key, val) + + return surf class Halfspace(Region): diff --git a/tests/regression_tests/cmfd_feed_2g/results_true.dat b/tests/regression_tests/cmfd_feed_2g/results_true.dat index d928b33e8a..dc1751574a 100644 --- a/tests/regression_tests/cmfd_feed_2g/results_true.dat +++ b/tests/regression_tests/cmfd_feed_2g/results_true.dat @@ -381,7 +381,7 @@ cmfd balance 2.68486E-04 4.84991E-04 1.08402E-03 -1.09177E-03 +1.09178E-03 5.45977E-04 4.45554E-04 4.01147E-04 diff --git a/tests/regression_tests/triso/results_true.dat b/tests/regression_tests/triso/results_true.dat index c8e6832345..0c5186e45b 100644 --- a/tests/regression_tests/triso/results_true.dat +++ b/tests/regression_tests/triso/results_true.dat @@ -1,2 +1,2 @@ k-combined: -1.701412E+00 3.180881E-02 +1.701412E+00 3.180877E-02 From bed26993ab1bb0e1225370e4bbcc8bfac82755d7 Mon Sep 17 00:00:00 2001 From: Ethan Peterson Date: Thu, 13 Feb 2020 16:17:30 -0500 Subject: [PATCH 18/33] fixed pytest results after ci failure --- tests/regression_tests/cmfd_feed_2g/results_true.dat | 2 +- tests/regression_tests/triso/results_true.dat | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/regression_tests/cmfd_feed_2g/results_true.dat b/tests/regression_tests/cmfd_feed_2g/results_true.dat index dc1751574a..d928b33e8a 100644 --- a/tests/regression_tests/cmfd_feed_2g/results_true.dat +++ b/tests/regression_tests/cmfd_feed_2g/results_true.dat @@ -381,7 +381,7 @@ cmfd balance 2.68486E-04 4.84991E-04 1.08402E-03 -1.09178E-03 +1.09177E-03 5.45977E-04 4.45554E-04 4.01147E-04 diff --git a/tests/regression_tests/triso/results_true.dat b/tests/regression_tests/triso/results_true.dat index 0c5186e45b..c8e6832345 100644 --- a/tests/regression_tests/triso/results_true.dat +++ b/tests/regression_tests/triso/results_true.dat @@ -1,2 +1,2 @@ k-combined: -1.701412E+00 3.180877E-02 +1.701412E+00 3.180881E-02 From 38640840c1107bfd9fcdf94b6541cd6e4d533643 Mon Sep 17 00:00:00 2001 From: Ethan Peterson Date: Fri, 14 Feb 2020 09:16:14 -0500 Subject: [PATCH 19/33] applied some suggestions from code review --- openmc/surface.py | 201 ++++++++++++++++++++++------------------------ 1 file changed, 96 insertions(+), 105 deletions(-) diff --git a/openmc/surface.py b/openmc/surface.py index 064ac7242f..0337c7297b 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -190,14 +190,17 @@ class Surface(IDManagerMixin, metaclass=ABCMeta): Parameters ---------- coeffs : tuple, optional - Tuple of surface coefficients to normalize. If none are supplied, - the coefficients from the current surface will be used. + Tuple of surface coefficients to normalize. Defaults to None. If no + coefficients are supplied then the coefficients will be taken from + the current Surface. Returns ------- tuple of normalized coefficients """ + if coeffs is None: + coeffs = self._get_base_coeffs() coeffs = np.asarray(coeffs) nonzeros = ~np.isclose(coeffs, 0., rtol=0., atol=self._atol) norm_factor = coeffs[nonzeros][0] @@ -213,8 +216,8 @@ class Surface(IDManagerMixin, metaclass=ABCMeta): surface """ - coeffs1 = normalize(self._get_base_coeffs()) - coeffs2 = normalize(other._get_base_coeffs()) + coeffs1 = self.normalize(self._get_base_coeffs()) + coeffs2 = self.normalize(other._get_base_coeffs()) return np.all(np.isclose(coeffs1, coeffs2, rtol=0., atol=self._atol)) @@ -351,16 +354,16 @@ class PlaneMixin(metaclass=ABCMeta): return a*x + b*y + c*z - d @abstractmethod - def translate(self, vector, clone=True): + def translate(self, vector, inplace=False): """Translate surface in given direction Parameters ---------- vector : iterable of float Direction in which surface should be translated - clone : boolean + inplace : boolean Whether or not to return a new instance of a Plane or to modify the - coefficients of this plane. Defaults to True + coefficients of this plane. Defaults to False Returns ------- @@ -447,8 +450,7 @@ class Plane(PlaneMixin, Surface): oldkwargs = deepcopy(kwargs) # work around for accepting Surface kwargs as positional parameters # until they are deprecated - argsdict = {k: v for k, v in zip(('boundary_type', 'name', - 'surface_id'), args)} + argsdict = dict(zip(('boundary_type', 'name', 'surface_id'), args)) for k, v in argsdict.items(): warn(_WARNING_KWARGS.format(type(self).__name__, k), FutureWarning) @@ -507,14 +509,14 @@ class Plane(PlaneMixin, Surface): def _get_base_coeffs(self): return (self.a, self.b, self.c, self.d) - def translate(self, vector, clone=True): + def translate(self, vector, inplace=False): vx, vy, vz = vector a, b, c, d = self._get_base_coeffs() d = d + a*vx + b*vy + c*vz - if clone: - surf = self.clone() - else: + if inplace: surf = self + else: + surf = self.clone() surf.d = d return surf @@ -628,11 +630,11 @@ class XPlane(PlaneMixin, Surface): def evaluate(self, point): return point[0] - self.x0 - def translate(self, vector, clone=True): - if clone: - surf = self.clone() - else: + def translate(self, vector, inplace=False): + if inplace: surf = self + else: + surf = self.clone() surf.x0 += vector[0] return surf @@ -717,11 +719,11 @@ class YPlane(PlaneMixin, Surface): def evaluate(self, point): return point[1] - self.y0 - def translate(self, vector, clone=True): - if clone: - surf = self.clone() - else: + def translate(self, vector, inplace=False): + if inplace: surf = self + else: + surf = self.clone() surf.y0 += vector[1] return surf @@ -806,11 +808,11 @@ class ZPlane(PlaneMixin, Surface): def evaluate(self, point): return point[2] - self.z0 - def translate(self, vector, clone=True): - if clone: - surf = self.clone() - else: + def translate(self, vector, inplace=False): + if inplace: surf = self + else: + surf = self.clone() surf.z0 += vector[2] return surf @@ -895,16 +897,16 @@ class QuadricMixin(metaclass=ABCMeta): return np.matmul(x.T, np.matmul(A, x)) + np.matmul(b.T, x) + c @abstractmethod - def translate(self, vector, clone=True): + def translate(self, vector, inplace=False): """Translate surface in given direction Parameters ---------- vector : iterable of float Direction in which surface should be translated - clone : bool + inplace : boolean Whether to return a clone of the Surface or the Surface itself. - Defaults to True + Defaults to False Returns ------- @@ -1055,7 +1057,7 @@ class Cylinder(QuadricMixin, Surface): def _get_base_coeffs(self): x0, y0, z0, r = self.x0, self.y0, self.z0, self.r dx, dy, dz = self.dx, self.dy, self.dz - dx2, dy2, dz2 = dx**2, dy**2, dz**2 + dx2, dy2, dz2 = dx*dx, dy*dy, dz*dz a = dy2 + dz2 b = dx2 + dz2 @@ -1066,16 +1068,16 @@ class Cylinder(QuadricMixin, Surface): g = 2*(dx*(z0*dz + y0*dy) - x0*(dy2 + dz2)) h = 2*(dy*(z0*dz + x0*dx) - y0*(dx2 + dz2)) j = 2*(dz*(y0*dy + x0*dx) - z0*(dx2 + dy2)) - k = x0**2*(dy2 + dz2) + y0**2*(dx2 + dz2) + z0**2*(dx2 + dy2) \ - + e*y0*z0 + f*x0*z0 + d*x0*y0 - r**2*(dx2 + dy2 + dz2) + k = x0*x0*(dy2 + dz2) + y0*y0*(dx2 + dz2) + z0*z0*(dx2 + dy2) \ + + e*y0*z0 + f*x0*z0 + d*x0*y0 - r*r*(dx2 + dy2 + dz2) return (a, b, c, d, e, f, g, h, j, k) - def translate(self, vector, clone=True): - if clone: - surf = self.clone() - else: + def translate(self, vector, inplace=False): + if inplace: surf = self + else: + surf = self.clone() surf.x0 += vector[0] surf.y0 += vector[1] @@ -1203,7 +1205,7 @@ class XCylinder(QuadricMixin, Surface): a = d = e = f = g = 0. b = c = 1. - h, j, k = -2*y0, -2*z0, y0**2 + z0**2 - r**2 + h, j, k = -2*y0, -2*z0, y0*y0 + z0*z0 - r*r return (a, b, c, d, e, f, g, h, j, k) @@ -1218,13 +1220,13 @@ class XCylinder(QuadricMixin, Surface): def evaluate(self, point): y = point[1] - self.y0 z = point[2] - self.z0 - return y**2 + z**2 - self.r**2 + return y*y + z*z - self.r**2 - def translate(self, vector, clone=True): - if clone: - surf = self.clone() - else: + def translate(self, vector, inplace=False): + if inplace: surf = self + else: + surf = self.clone() surf.y0 += vector[1] surf.z0 += vector[2] return surf @@ -1324,7 +1326,7 @@ class YCylinder(QuadricMixin, Surface): b = d = e = f = h = 0. a = c = 1. - g, j, k = -2*x0, -2*z0, x0**2 + z0**2 - r**2 + g, j, k = -2*x0, -2*z0, x0*x0 + z0*z0 - r*r return (a, b, c, d, e, f, g, h, j, k) @@ -1339,13 +1341,13 @@ class YCylinder(QuadricMixin, Surface): def evaluate(self, point): x = point[0] - self.x0 z = point[2] - self.z0 - return x**2 + z**2 - self.r**2 + return x*x + z*z - self.r**2 - def translate(self, vector, clone=True): - if clone: - surf = self.clone() - else: + def translate(self, vector, inplace=False): + if inplace: surf = self + else: + surf = self.clone() surf.x0 += vector[0] surf.z0 += vector[2] return surf @@ -1445,7 +1447,7 @@ class ZCylinder(QuadricMixin, Surface): c = d = e = f = j = 0. a = b = 1. - g, h, k = -2*x0, -2*y0, x0**2 + y0**2 - r**2 + g, h, k = -2*x0, -2*y0, x0*x0 + y0*y0 - r*r return (a, b, c, d, e, f, g, h, j, k) @@ -1460,13 +1462,13 @@ class ZCylinder(QuadricMixin, Surface): def evaluate(self, point): x = point[0] - self.x0 y = point[1] - self.y0 - return x**2 + y**2 - self.r**2 + return x*x + y*y - self.r**2 - def translate(self, vector, clone=True): - if clone: - surf = self.clone() - else: + def translate(self, vector, inplace=False): + if inplace: surf = self + else: + surf = self.clone() surf.x0 += vector[0] surf.y0 += vector[1] return surf @@ -1577,7 +1579,7 @@ class Sphere(QuadricMixin, Surface): a = b = c = 1. d = e = f = 0. g, h, j = -2*x0, -2*y0, -2*z0 - k = x0**2 + y0**2 + z0**2 - r**2 + k = x0*x0 + y0*y0 + z0*z0 - r*r return (a, b, c, d, e, f, g, h, j, k) @@ -1595,13 +1597,13 @@ class Sphere(QuadricMixin, Surface): x = point[0] - self.x0 y = point[1] - self.y0 z = point[2] - self.z0 - return x**2 + y**2 + z**2 - self.r**2 + return x*x + y*y + z*z - self.r**2 - def translate(self, vector, clone=True): - if clone: - surf = self.clone() - else: + def translate(self, vector, inplace=False): + if inplace: surf = self + else: + surf = self.clone() surf.x0 += vector[0] surf.y0 += vector[1] surf.z0 += vector[2] @@ -1751,26 +1753,26 @@ class Cone(QuadricMixin, Surface): x0, y0, z0, r2 = self.x0, self.y0, self.z0, self.r2 dx, dy, dz = self.dx, self.dy, self.dz cos2 = 1 / (1 + r2) - c1 = (dx**2 + dy**2 + dz**2)*cos2 + c1 = (dx*dx + dy*dy + dz*dz)*cos2 - a = dx**2 - c1 - b = dy**2 - c1 - c = dz**2 - c1 + a = dx*dx - c1 + b = dy*dy - c1 + c = dz*dz - c1 d = 2*dx*dy e = 2*dy*dz f = 2*dx*dz - g = -2*(dx**2*x0 + dx*dy*y0 + dx*dz*z0 - c1) - h = -2*(dy**2*y0 + dx*dy*x0 + dy*dz*z0 - c1) - j = -2*(dz**2*y0 + dx*dz*x0 + dy*dz*y0 - c1) - k = (dx*x0 + dy*y0 + dz*z0)**2 - c1*(x0**2 + y0**2 + z0**2) + g = -2*(dx*dx*x0 + dx*dy*y0 + dx*dz*z0 - c1) + h = -2*(dy*dy*y0 + dx*dy*x0 + dy*dz*z0 - c1) + j = -2*(dz*dz*y0 + dx*dz*x0 + dy*dz*y0 - c1) + k = (dx*x0 + dy*y0 + dz*z0)**2 - c1*(x0*x0 + y0*y0 + z0*z0) return (a, b, c, d, e, f, g, h, j, k) - def translate(self, vector, clone=True): - if clone: - surf = self.clone() - else: + def translate(self, vector, inplace=False): + if inplace: surf = self + else: + surf = self.clone() surf.x0 += vector[0] surf.y0 += vector[1] surf.z0 += vector[2] @@ -1882,7 +1884,7 @@ class XCone(QuadricMixin, Surface): b = c = 1. d = e = f = 0. g, h, j = 2*x0*r2, -2*y0, -2*z0 - k = y0**2 + z0**2 - r2*x0**2 + k = y0*y0 + z0*z0 - r2*x0*x0 return (a, b, c, d, e, f, g, h, j, k) @@ -1890,13 +1892,13 @@ class XCone(QuadricMixin, Surface): x = point[0] - self.x0 y = point[1] - self.y0 z = point[2] - self.z0 - return y**2 + z**2 - self.r2*x**2 + return y*y + z*z - self.r2*x*x - def translate(self, vector, clone=True): - if clone: - surf = self.clone() - else: + def translate(self, vector, inplace=False): + if inplace: surf = self + else: + surf = self.clone() surf.x0 += vector[0] surf.y0 += vector[1] surf.z0 += vector[2] @@ -2011,7 +2013,7 @@ class YCone(QuadricMixin, Surface): a = c = 1. d = e = f = 0. g, h, j = -2*x0, 2*y0*r2, -2*z0 - k = x0**2 + z0**2 - r2*y0**2 + k = x0*x0 + z0*z0 - r2*y0*y0 return (a, b, c, d, e, f, g, h, j, k) @@ -2019,13 +2021,13 @@ class YCone(QuadricMixin, Surface): x = point[0] - self.x0 y = point[1] - self.y0 z = point[2] - self.z0 - return x**2 + z**2 - self.r2*y**2 + return x*x + z*z - self.r2*y*y - def translate(self, vector, clone=True): - if clone: - surf = self.clone() - else: + def translate(self, vector, inplace=False): + if inplace: surf = self + else: + surf = self.clone() surf.x0 += vector[0] surf.y0 += vector[1] surf.z0 += vector[2] @@ -2140,7 +2142,7 @@ class ZCone(QuadricMixin, Surface): a = b = 1. d = e = f = 0. g, h, j = -2*x0, -2*y0, 2*z0*r2 - k = x0**2 + y0**2 - r2*z0**2 + k = x0*x0 + y0*y0 - r2*z0*z0 return (a, b, c, d, e, f, g, h, j, k) @@ -2148,13 +2150,13 @@ class ZCone(QuadricMixin, Surface): x = point[0] - self.x0 y = point[1] - self.y0 z = point[2] - self.z0 - return x**2 + y**2 - self.r2*z**2 + return x*x + y*y - self.r2*z*z - def translate(self, vector, clone=True): - if clone: - surf = self.clone() - else: + def translate(self, vector, inplace=False): + if inplace: surf = self + else: + surf = self.clone() surf.x0 += vector[0] surf.y0 += vector[1] surf.z0 += vector[2] @@ -2304,7 +2306,7 @@ class Quadric(QuadricMixin, Surface): def _get_base_coeffs(self): return tuple(getattr(self, c) for c in self._coeff_keys) - def translate(self, vector, clone=True): + def translate(self, vector, inplace=False): vector = np.asarray(vector) A, bvec, cnst = self.get_Abc() @@ -2312,10 +2314,10 @@ class Quadric(QuadricMixin, Surface): k = cnst + np.matmul(vector.T, np.matmul(A, vector)) \ - np.matmul(bvec.T, vector) - if clone: - surf = self.clone() - else: + if inplace: surf = self + else: + surf = self.clone() for key, val in zip(('g', 'h', 'j', 'k'), (g, h, j, k)): setattr(surf, key, val) @@ -2511,17 +2513,6 @@ class Halfspace(Region): return type(self)(memo[key], self.side) -def _get_subclasses(cls): - """Recursively find all subclasses of this class""" - return set(cls.__subclasses__()).union([s for c in cls.__subclasses__() - for s in _get_subclasses(c)]) - - -def _get_subclass_map(cls): - """Generate mapping of class _type attributes to classes""" - return {c._type: c for c in _get_subclasses(cls)} - - -_SURFACE_CLASSES = _get_subclass_map(Surface) +_SURFACE_CLASSES = {cls._type: cls for cls in Surface.__subclasses__()} From 33c46c78c87c6b1b8400d5bab5a44b1bf6164048 Mon Sep 17 00:00:00 2001 From: Ethan Peterson Date: Fri, 14 Feb 2020 09:55:57 -0500 Subject: [PATCH 20/33] refactored translate method back out to mixins --- openmc/surface.py | 181 +++++++++------------------------------------- 1 file changed, 35 insertions(+), 146 deletions(-) diff --git a/openmc/surface.py b/openmc/surface.py index 0337c7297b..d45c8387f3 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -353,7 +353,6 @@ class PlaneMixin(metaclass=ABCMeta): a, b, c, d = self._get_base_coeffs() return a*x + b*y + c*z - d - @abstractmethod def translate(self, vector, inplace=False): """Translate surface in given direction @@ -371,7 +370,18 @@ class PlaneMixin(metaclass=ABCMeta): Translated surface """ - pass + vx, vy, vz = vector + a, b, c, d = self._get_base_coeffs() + d = d + a*vx + b*vy + c*vz + + if inplace: + surf = self + else: + surf = self.clone() + + setattr(surf, self._coeff_keys[-1], d) + + return surf def to_xml_element(self): """Return XML representation of the surface @@ -509,17 +519,6 @@ class Plane(PlaneMixin, Surface): def _get_base_coeffs(self): return (self.a, self.b, self.c, self.d) - def translate(self, vector, inplace=False): - vx, vy, vz = vector - a, b, c, d = self._get_base_coeffs() - d = d + a*vx + b*vy + c*vz - if inplace: - surf = self - else: - surf = self.clone() - surf.d = d - return surf - @classmethod def from_points(cls, p1, p2, p3, **kwargs): """Return a plane given three points that pass through it. @@ -630,14 +629,6 @@ class XPlane(PlaneMixin, Surface): def evaluate(self, point): return point[0] - self.x0 - def translate(self, vector, inplace=False): - if inplace: - surf = self - else: - surf = self.clone() - surf.x0 += vector[0] - return surf - Plane.register(XPlane) @@ -719,14 +710,6 @@ class YPlane(PlaneMixin, Surface): def evaluate(self, point): return point[1] - self.y0 - def translate(self, vector, inplace=False): - if inplace: - surf = self - else: - surf = self.clone() - surf.y0 += vector[1] - return surf - Plane.register(YPlane) @@ -808,14 +791,6 @@ class ZPlane(PlaneMixin, Surface): def evaluate(self, point): return point[2] - self.z0 - def translate(self, vector, inplace=False): - if inplace: - surf = self - else: - surf = self.clone() - surf.z0 += vector[2] - return surf - Plane.register(ZPlane) @@ -896,7 +871,6 @@ class QuadricMixin(metaclass=ABCMeta): A, b, c = self.get_Abc() return np.matmul(x.T, np.matmul(A, x)) + np.matmul(b.T, x) + c - @abstractmethod def translate(self, vector, inplace=False): """Translate surface in given direction @@ -914,7 +888,29 @@ class QuadricMixin(metaclass=ABCMeta): Translated surface """ - pass + vector = np.asarray(vector) + + if inplace: + surf = self + else: + surf = self.clone() + + if set(('x0', 'y0', 'z0')).intersection(set(surf._coeff_keys)): + for vi, xi in zip(vector, ('x0', 'y0', 'z0')): + val = getattr(surf, xi, None) + if val is not None: + setattr(surf, xi, val + vi) + else: + 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) + + for key, val in zip(('g', 'h', 'j', 'k'), (g, h, j, k)): + setattr(surf, key, val) + + return surf class Cylinder(QuadricMixin, Surface): @@ -1073,18 +1069,6 @@ class Cylinder(QuadricMixin, Surface): return (a, b, c, d, e, f, g, h, j, k) - def translate(self, vector, inplace=False): - if inplace: - surf = self - else: - surf = self.clone() - - surf.x0 += vector[0] - surf.y0 += vector[1] - surf.z0 += vector[2] - - return surf - @classmethod def from_points(cls, p1, p2, r=1., **kwargs): """Return a cylinder given points that define the axis and a radius. @@ -1222,15 +1206,6 @@ class XCylinder(QuadricMixin, Surface): z = point[2] - self.z0 return y*y + z*z - self.r**2 - def translate(self, vector, inplace=False): - if inplace: - surf = self - else: - surf = self.clone() - surf.y0 += vector[1] - surf.z0 += vector[2] - return surf - Cylinder.register(XCylinder) @@ -1343,15 +1318,6 @@ class YCylinder(QuadricMixin, Surface): z = point[2] - self.z0 return x*x + z*z - self.r**2 - def translate(self, vector, inplace=False): - if inplace: - surf = self - else: - surf = self.clone() - surf.x0 += vector[0] - surf.z0 += vector[2] - return surf - Cylinder.register(YCylinder) @@ -1464,15 +1430,6 @@ class ZCylinder(QuadricMixin, Surface): y = point[1] - self.y0 return x*x + y*y - self.r**2 - def translate(self, vector, inplace=False): - if inplace: - surf = self - else: - surf = self.clone() - surf.x0 += vector[0] - surf.y0 += vector[1] - return surf - Cylinder.register(ZCylinder) @@ -1599,16 +1556,6 @@ class Sphere(QuadricMixin, Surface): z = point[2] - self.z0 return x*x + y*y + z*z - self.r**2 - def translate(self, vector, inplace=False): - if inplace: - surf = self - else: - surf = self.clone() - surf.x0 += vector[0] - surf.y0 += vector[1] - surf.z0 += vector[2] - return surf - class Cone(QuadricMixin, Surface): """A conical surface parallel to the x-, y-, or z-axis. @@ -1768,16 +1715,6 @@ class Cone(QuadricMixin, Surface): return (a, b, c, d, e, f, g, h, j, k) - def translate(self, vector, inplace=False): - if inplace: - surf = self - else: - surf = self.clone() - surf.x0 += vector[0] - surf.y0 += vector[1] - surf.z0 += vector[2] - return surf - class XCone(QuadricMixin, Surface): """A cone parallel to the x-axis of the form :math:`(y - y_0)^2 + (z - z_0)^2 = @@ -1894,16 +1831,6 @@ class XCone(QuadricMixin, Surface): z = point[2] - self.z0 return y*y + z*z - self.r2*x*x - def translate(self, vector, inplace=False): - if inplace: - surf = self - else: - surf = self.clone() - surf.x0 += vector[0] - surf.y0 += vector[1] - surf.z0 += vector[2] - return surf - Cone.register(XCone) @@ -2023,16 +1950,6 @@ class YCone(QuadricMixin, Surface): z = point[2] - self.z0 return x*x + z*z - self.r2*y*y - def translate(self, vector, inplace=False): - if inplace: - surf = self - else: - surf = self.clone() - surf.x0 += vector[0] - surf.y0 += vector[1] - surf.z0 += vector[2] - return surf - Cone.register(YCone) @@ -2152,16 +2069,6 @@ class ZCone(QuadricMixin, Surface): z = point[2] - self.z0 return x*x + y*y - self.r2*z*z - def translate(self, vector, inplace=False): - if inplace: - surf = self - else: - surf = self.clone() - surf.x0 += vector[0] - surf.y0 += vector[1] - surf.z0 += vector[2] - return surf - Cone.register(ZCone) @@ -2306,24 +2213,6 @@ class Quadric(QuadricMixin, Surface): def _get_base_coeffs(self): return tuple(getattr(self, c) for c in self._coeff_keys) - def translate(self, vector, inplace=False): - vector = np.asarray(vector) - 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) - - if inplace: - surf = self - else: - surf = self.clone() - - for key, val in zip(('g', 'h', 'j', 'k'), (g, h, j, k)): - setattr(surf, key, val) - - return surf - class Halfspace(Region): """A positive or negative half-space region. From 35fb12b3619b8d2614ac007c59c1f6a1584152db Mon Sep 17 00:00:00 2001 From: Ethan Peterson Date: Fri, 14 Feb 2020 11:37:26 -0500 Subject: [PATCH 21/33] implemented Plane.__new__ for returning simplified versions --- openmc/surface.py | 95 ++++++++++++++++++++++++++++++++++++----------- 1 file changed, 73 insertions(+), 22 deletions(-) diff --git a/openmc/surface.py b/openmc/surface.py index d45c8387f3..2e75945d37 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -379,7 +379,7 @@ class PlaneMixin(metaclass=ABCMeta): else: surf = self.clone() - setattr(surf, self._coeff_keys[-1], d) + setattr(surf, surf._coeff_keys[-1], d) return surf @@ -410,9 +410,9 @@ class Plane(PlaneMixin, Surface): a : float, optional The 'A' parameter for the plane. Defaults to 1. b : float, optional - The 'B' parameter for the plane. Defaults to 0. + The 'B' parameter for the plane. Defaults to 1. c : float, optional - The 'C' parameter for the plane. Defaults to 0. + The 'C' parameter for the plane. Defaults to 1. d : float, optional The 'D' parameter for the plane. Defaults to 0. boundary_type : {'transmission, 'vacuum', 'reflective', 'white'}, optional @@ -455,30 +455,82 @@ class Plane(PlaneMixin, Surface): _type = 'plane' _coeff_keys = ('a', 'b', 'c', 'd') - def __init__(self, a=1., b=0., c=0., d=0., *args, **kwargs): - # work around until capital letter kwargs are deprecated - oldkwargs = deepcopy(kwargs) - # work around for accepting Surface kwargs as positional parameters - # until they are deprecated - argsdict = dict(zip(('boundary_type', 'name', 'surface_id'), args)) - for k, v in argsdict.items(): - warn(_WARNING_KWARGS.format(type(self).__name__, k), FutureWarning) - - kwargs.update(argsdict) - - for k in 'ABCD': - kwargs.pop(k, None) + def __new__(cls, *args, **kwargs): + # Intercept Plane creation to return a simpler form if warranted + new_cls, kwargs = cls._process_args(*args, **kwargs) + obj = super().__new__(new_cls) + if new_cls is not cls: + obj.__init__(**kwargs) + return obj + def __init__(self, a=1., b=1., c=1., d=0., *args, **kwargs): super().__init__(**kwargs) for key, val in zip(self._coeff_keys, (a, b, c, d)): setattr(self, key, val) - for k, v in oldkwargs.items(): - if k in 'ABCD': - warn(_WARNING_UPPER.format(type(self).__name__, k.lower(), k), + @classmethod + def _process_args(cls, *args, **kwargs): + """ Process arguments to Plane class to determine which type of Plane + should be instantiated. + + Returns + ------- + tuple : tuple of type and dict + Returns a tuple where the first element is the class (XPlane, + YPlane, ZPlane, or Plane) that should be instantiated. The second + element is the dictionary of keyword arguments that should be + passed to that object's constructor. + + """ + # *args should ultimately be limited to a, b, c, d as specified in + # __init__, but to preserve the API it is allowed to accept Surface + # parameters for now, but will raise warnings if this is done. + argtup = ('a', 'b', 'c', 'd', 'boundary_type', 'name', 'surface_id') + kwargs.update(dict(zip(argtup, args))) + + # Warn if capital letter arguments are passed + for k in 'ABCD': + val = kwargs.pop(k, None) + if val is not None: + warn(_WARNING_UPPER.format(cls.__name__, k.lower(), k), FutureWarning) - setattr(self, k.lower(), v) + kwargs[k.lower()] = val + + # Warn if Surface parameters are passed by position, not by keyword + for k in ('boundary_type', 'name', 'surface_id'): + val = kwargs.get(k, None) + if val is not None: + warn(_WARNING_KWARGS.format(cls.__name__, k), + FutureWarning) + + # Determine if a simpler version of a Plane should be returned and + # modify kwargs accordingly + a = kwargs.get('a', 1.) + b = kwargs.get('b', 1.) + c = kwargs.get('c', 1.) + + # If two of a, b, or c are zero, a simpler version should be returned + if np.sum(np.isclose((a, b, c), 0., rtol=0., atol=cls._atol)) >= 2: + a = kwargs.pop('a', 1.) + b = kwargs.pop('b', 1.) + c = kwargs.pop('c', 1.) + d = kwargs.pop('d', 0.) + # Return XPlane class and kwargs + if ~np.isclose(a, 0.0, rtol=0., atol=cls._atol): + kwargs['x0'] = d/a + return XPlane, kwargs + # Return YPlane class and kwargs + elif ~np.isclose(b, 0.0, rtol=0., atol=cls._atol): + kwargs['y0'] = d/b + return YPlane, kwargs + # Return ZPlane class and kwargs + elif ~np.isclose(c, 0.0, rtol=0., atol=cls._atol): + kwargs['z0'] = d/c + return ZPlane, kwargs + else: + # Return Plane class and kwargs + return cls, kwargs @property def a(self): @@ -597,8 +649,7 @@ class XPlane(PlaneMixin, Surface): def __init__(self, x0=0., *args, **kwargs): # work around for accepting Surface kwargs as positional parameters # until they are deprecated - argsdict = {k: v for k, v in zip(('boundary_type', 'name', - 'surface_id'), args)} + argsdict = dict(zip(('boundary_type', 'name', 'surface_id'), args)) for k, v in argsdict.items(): warn(_WARNING_KWARGS.format(type(self).__name__, k), FutureWarning) kwargs.update(argsdict) From ae696503b6af3d09505556965e38e86881fa5d9b Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 17 Feb 2020 07:52:52 -0600 Subject: [PATCH 22/33] Avoid copysign in instances where an int return value is needed --- src/cell.cpp | 6 +++--- src/output.cpp | 2 +- src/particle.cpp | 4 ++-- 3 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/cell.cpp b/src/cell.cpp index cd6199b755..1d175f05d2 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -377,7 +377,7 @@ CSGCell::CSGCell(pugi::xml_node cell_node) throw std::runtime_error{"Invalid surface ID " + std::to_string(abs(r)) + " specified in region for cell " + std::to_string(id_) + "."}; } - r = copysign(it->second + 1, r); + r = (r > 0) ? it->second + 1 : -(it->second + 1); } } @@ -536,8 +536,8 @@ CSGCell::to_hdf5(hid_t cell_group) const region_spec << " |"; } else { // Note the off-by-one indexing - region_spec << " " - << copysign(model::surfaces[abs(token)-1]->id_, token); + auto surf_id = model::surfaces[abs(token)-1]->id_; + region_spec << " " << ((token > 0) ? surf_id : -surf_id); } } write_string(group, "region", region_spec.str(), false); diff --git a/src/output.cpp b/src/output.cpp index d98cfd4a1d..18b94be9be 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -192,7 +192,7 @@ extern "C" void print_particle(Particle* p) // Display miscellaneous info. if (p->surface_ != 0) { const Surface& surf {*model::surfaces[std::abs(p->surface_)-1]}; - fmt::print(" Surface = {}\n", std::copysign(surf.id_, p->surface_)); + fmt::print(" Surface = {}\n", (p->surface_ > 0) ? surf.id_ : -surf.id_); } fmt::print(" Weight = {}\n", p->wgt_); if (settings::run_CE) { diff --git a/src/particle.cpp b/src/particle.cpp index d2d8b2e1b3..ea0d03cd01 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -1,7 +1,7 @@ #include "openmc/particle.h" #include // copy, min -#include // log, abs, copysign +#include // log, abs #include @@ -539,7 +539,7 @@ Particle::cross_surface() // TODO: off-by-one surface_ = rotational ? surf_p->i_periodic_ + 1 : - std::copysign(surf_p->i_periodic_ + 1, surface_); + ((surface_ > 0) ? surf_p->i_periodic_ + 1 : -(surf_p->i_periodic_ + 1)); // Figure out what cell particle is in now n_coord_ = 1; From 2332a679fadbf2142ae21a1a0fbe46fd90b62f48 Mon Sep 17 00:00:00 2001 From: Ethan Peterson Date: Tue, 18 Feb 2020 10:08:22 -0500 Subject: [PATCH 23/33] reverted simplifying Plane constructor and default args --- openmc/surface.py | 98 +++++-------------- .../deplete/example_geometry.py | 1 + 2 files changed, 25 insertions(+), 74 deletions(-) diff --git a/openmc/surface.py b/openmc/surface.py index 2e75945d37..c73ce6fc87 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -203,7 +203,7 @@ class Surface(IDManagerMixin, metaclass=ABCMeta): coeffs = self._get_base_coeffs() coeffs = np.asarray(coeffs) nonzeros = ~np.isclose(coeffs, 0., rtol=0., atol=self._atol) - norm_factor = coeffs[nonzeros][0] + norm_factor = np.abs(coeffs[nonzeros][0]) return tuple([c/norm_factor for c in coeffs]) def is_equal(self, other): @@ -410,9 +410,9 @@ class Plane(PlaneMixin, Surface): a : float, optional The 'A' parameter for the plane. Defaults to 1. b : float, optional - The 'B' parameter for the plane. Defaults to 1. + The 'B' parameter for the plane. Defaults to 0. c : float, optional - The 'C' parameter for the plane. Defaults to 1. + The 'C' parameter for the plane. Defaults to 0. d : float, optional The 'D' parameter for the plane. Defaults to 0. boundary_type : {'transmission, 'vacuum', 'reflective', 'white'}, optional @@ -455,82 +455,34 @@ class Plane(PlaneMixin, Surface): _type = 'plane' _coeff_keys = ('a', 'b', 'c', 'd') - def __new__(cls, *args, **kwargs): - # Intercept Plane creation to return a simpler form if warranted - new_cls, kwargs = cls._process_args(*args, **kwargs) - obj = super().__new__(new_cls) - if new_cls is not cls: - obj.__init__(**kwargs) - return obj - - def __init__(self, a=1., b=1., c=1., d=0., *args, **kwargs): - super().__init__(**kwargs) - - for key, val in zip(self._coeff_keys, (a, b, c, d)): - setattr(self, key, val) - - @classmethod - def _process_args(cls, *args, **kwargs): - """ Process arguments to Plane class to determine which type of Plane - should be instantiated. - - Returns - ------- - tuple : tuple of type and dict - Returns a tuple where the first element is the class (XPlane, - YPlane, ZPlane, or Plane) that should be instantiated. The second - element is the dictionary of keyword arguments that should be - passed to that object's constructor. - - """ + def __init__(self, a=1., b=0., c=0., d=0., *args, **kwargs): # *args should ultimately be limited to a, b, c, d as specified in # __init__, but to preserve the API it is allowed to accept Surface # parameters for now, but will raise warnings if this is done. argtup = ('a', 'b', 'c', 'd', 'boundary_type', 'name', 'surface_id') kwargs.update(dict(zip(argtup, args))) + # Warn if Surface parameters are passed by position, not by keyword + superkwargs = {} + for k in ('boundary_type', 'name', 'surface_id'): + val = kwargs.get(k, None) + if val is not None: + superkwargs[k] = val + warn(_WARNING_KWARGS.format(type(self), k), + FutureWarning) + + super().__init__(**superkwargs) + + for key, val in zip(self._coeff_keys, (a, b, c, d)): + setattr(self, key, val) + # Warn if capital letter arguments are passed for k in 'ABCD': val = kwargs.pop(k, None) if val is not None: - warn(_WARNING_UPPER.format(cls.__name__, k.lower(), k), + warn(_WARNING_UPPER.format(type(self), k.lower(), k), FutureWarning) - kwargs[k.lower()] = val - - # Warn if Surface parameters are passed by position, not by keyword - for k in ('boundary_type', 'name', 'surface_id'): - val = kwargs.get(k, None) - if val is not None: - warn(_WARNING_KWARGS.format(cls.__name__, k), - FutureWarning) - - # Determine if a simpler version of a Plane should be returned and - # modify kwargs accordingly - a = kwargs.get('a', 1.) - b = kwargs.get('b', 1.) - c = kwargs.get('c', 1.) - - # If two of a, b, or c are zero, a simpler version should be returned - if np.sum(np.isclose((a, b, c), 0., rtol=0., atol=cls._atol)) >= 2: - a = kwargs.pop('a', 1.) - b = kwargs.pop('b', 1.) - c = kwargs.pop('c', 1.) - d = kwargs.pop('d', 0.) - # Return XPlane class and kwargs - if ~np.isclose(a, 0.0, rtol=0., atol=cls._atol): - kwargs['x0'] = d/a - return XPlane, kwargs - # Return YPlane class and kwargs - elif ~np.isclose(b, 0.0, rtol=0., atol=cls._atol): - kwargs['y0'] = d/b - return YPlane, kwargs - # Return ZPlane class and kwargs - elif ~np.isclose(c, 0.0, rtol=0., atol=cls._atol): - kwargs['z0'] = d/c - return ZPlane, kwargs - else: - # Return Plane class and kwargs - return cls, kwargs + setattr(self, k.lower(), val) @property def a(self): @@ -729,9 +681,8 @@ class YPlane(PlaneMixin, Surface): def __init__(self, y0=0., *args, **kwargs): # work around for accepting Surface kwargs as positional parameters # until they are deprecated - argsdict = {k: v for k, v in zip(('boundary_type', 'name', - 'surface_id'), args)} - for k, v in argsdict.items(): + argsdict = dict(zip(('boundary_type', 'name', 'surface_id'), args)) + for k in argsdict.keys(): warn(_WARNING_KWARGS.format(type(self).__name__, k), FutureWarning) kwargs.update(argsdict) @@ -810,9 +761,8 @@ class ZPlane(PlaneMixin, Surface): def __init__(self, z0=0., *args, **kwargs): # work around for accepting Surface kwargs as positional parameters # until they are deprecated - argsdict = {k: v for k, v in zip(('boundary_type', 'name', - 'surface_id'), args)} - for k, v in argsdict.items(): + argsdict = dict(zip(('boundary_type', 'name', 'surface_id'), args)) + for k in argsdict.keys(): warn(_WARNING_KWARGS.format(type(self).__name__, k), FutureWarning) kwargs.update(argsdict) diff --git a/tests/regression_tests/deplete/example_geometry.py b/tests/regression_tests/deplete/example_geometry.py index f79044558c..1275e23b26 100644 --- a/tests/regression_tests/deplete/example_geometry.py +++ b/tests/regression_tests/deplete/example_geometry.py @@ -188,6 +188,7 @@ def segment_pin(n_rings, n_wedges, r_fuel, r_gap, r_clad): fuel_wedges = [openmc.Plane(a=math.cos(theta[i]), b=math.sin(theta[i])) for i in range(n_wedges)] + print(fuel_wedges) gap_ring = openmc.ZCylinder(x0=0, y0=0, r=r_gap) clad_ring = openmc.ZCylinder(x0=0, y0=0, r=r_clad) From 03f33943ad46d1a7bfdc17b707f40bb66fd24455 Mon Sep 17 00:00:00 2001 From: Ethan Peterson Date: Tue, 18 Feb 2020 11:54:33 -0500 Subject: [PATCH 24/33] removed debug statement --- tests/regression_tests/deplete/example_geometry.py | 1 - 1 file changed, 1 deletion(-) diff --git a/tests/regression_tests/deplete/example_geometry.py b/tests/regression_tests/deplete/example_geometry.py index 1275e23b26..f79044558c 100644 --- a/tests/regression_tests/deplete/example_geometry.py +++ b/tests/regression_tests/deplete/example_geometry.py @@ -188,7 +188,6 @@ def segment_pin(n_rings, n_wedges, r_fuel, r_gap, r_clad): fuel_wedges = [openmc.Plane(a=math.cos(theta[i]), b=math.sin(theta[i])) for i in range(n_wedges)] - print(fuel_wedges) gap_ring = openmc.ZCylinder(x0=0, y0=0, r=r_gap) clad_ring = openmc.ZCylinder(x0=0, y0=0, r=r_clad) From dd74b2f43f1d2060486de2823ee30058b8971dbd Mon Sep 17 00:00:00 2001 From: Sam Powell-Gill Date: Wed, 19 Feb 2020 09:26:13 +0000 Subject: [PATCH 25/33] Fixed bug which prevents creating photon data with no atomic relax --- openmc/data/photon.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/openmc/data/photon.py b/openmc/data/photon.py index 32ed3010a4..9cc697b782 100644 --- a/openmc/data/photon.py +++ b/openmc/data/photon.py @@ -784,8 +784,9 @@ class IncidentPhoton(EqualityMixin): sub_group = shell_group.create_group(key) # Write atomic relaxation - if key in self.atomic_relaxation.subshells: - self.atomic_relaxation.to_hdf5(sub_group, key) + if self.atomic_relaxation is not None: + if key in self.atomic_relaxation.subshells: + self.atomic_relaxation.to_hdf5(sub_group, key) else: continue From 6fd08f5aecd835e695b65657598c8a7e7c64e26b Mon Sep 17 00:00:00 2001 From: Sam Powell-Gill Date: Wed, 19 Feb 2020 11:59:35 +0000 Subject: [PATCH 26/33] Unit test for exporting photon data with no atomic relax --- tests/unit_tests/test_data_photon.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/tests/unit_tests/test_data_photon.py b/tests/unit_tests/test_data_photon.py index c767d19e6c..02635bcc9a 100644 --- a/tests/unit_tests/test_data_photon.py +++ b/tests/unit_tests/test_data_photon.py @@ -144,3 +144,11 @@ def test_export_to_hdf5(tmpdir, element): element2.bremsstrahlung['electron_energy']).all() # Export to hdf5 again element2.export_to_hdf5(filename, 'w') + +def test_photodat_only(tmpdir): + endf_data = os.environ['OPENMC_ENDF_DATA'] + filename = str(tmpdir.join('tmp.h5')) + p_file = 'photoat-{:03}_{}_000.endf'.format(1, 'H') + p_path = os.path.join(endf_data, 'photoat', p_file) + data=openmc.data.IncidentPhoton.from_endf(p_path) + data.export_to_hdf5(filename, 'w') \ No newline at end of file From 9a21eabd55d1353ca8241abbd91578b4ede84502 Mon Sep 17 00:00:00 2001 From: Ethan Peterson Date: Wed, 19 Feb 2020 08:46:50 -0500 Subject: [PATCH 27/33] Apply suggestions from code review Co-Authored-By: Paul Romano --- openmc/surface.py | 23 +++++++++++------------ 1 file changed, 11 insertions(+), 12 deletions(-) diff --git a/openmc/surface.py b/openmc/surface.py index c73ce6fc87..6d08507d69 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -19,7 +19,7 @@ _WARNING_UPPER = """\ will not accept the capitalized version.\ """ -_WARNING_KWARGS = """ +_WARNING_KWARGS = """\ "{}(...) accepts keyword arguments only for '{}'. Future versions of OpenMC \ will not accept positional parameters for superclass arguments.\ """ @@ -219,7 +219,7 @@ class Surface(IDManagerMixin, metaclass=ABCMeta): coeffs1 = self.normalize(self._get_base_coeffs()) coeffs2 = self.normalize(other._get_base_coeffs()) - return np.all(np.isclose(coeffs1, coeffs2, rtol=0., atol=self._atol)) + return np.allclose(coeffs1, coeffs2, rtol=0., atol=self._atol) @abstractmethod def evaluate(self, point): @@ -599,10 +599,10 @@ class XPlane(PlaneMixin, Surface): _coeff_keys = ('x0',) def __init__(self, x0=0., *args, **kwargs): - # work around for accepting Surface kwargs as positional parameters + # work around for accepting Surface kwargs as positional parameters # until they are deprecated argsdict = dict(zip(('boundary_type', 'name', 'surface_id'), args)) - for k, v in argsdict.items(): + for k in argsdict: warn(_WARNING_KWARGS.format(type(self).__name__, k), FutureWarning) kwargs.update(argsdict) @@ -679,10 +679,10 @@ class YPlane(PlaneMixin, Surface): _coeff_keys = ('y0',) def __init__(self, y0=0., *args, **kwargs): - # work around for accepting Surface kwargs as positional parameters + # work around for accepting Surface kwargs as positional parameters # until they are deprecated argsdict = dict(zip(('boundary_type', 'name', 'surface_id'), args)) - for k in argsdict.keys(): + for k in argsdict: warn(_WARNING_KWARGS.format(type(self).__name__, k), FutureWarning) kwargs.update(argsdict) @@ -759,10 +759,10 @@ class ZPlane(PlaneMixin, Surface): _coeff_keys = ('z0',) def __init__(self, z0=0., *args, **kwargs): - # work around for accepting Surface kwargs as positional parameters + # work around for accepting Surface kwargs as positional parameters # until they are deprecated argsdict = dict(zip(('boundary_type', 'name', 'surface_id'), args)) - for k in argsdict.keys(): + for k in argsdict: warn(_WARNING_KWARGS.format(type(self).__name__, k), FutureWarning) kwargs.update(argsdict) @@ -825,8 +825,8 @@ class QuadricMixin(metaclass=ABCMeta): else: a, b, c, d, e, f, g, h, j, k = coeffs - A = np.asarray([[a, d/2, f/2], [d/2, b, e/2], [f/2, e/2, c]]) - bvec = np.asarray([g, h, j]) + A = np.array([[a, d/2, f/2], [d/2, b, e/2], [f/2, e/2, c]]) + bvec = np.array([g, h, j]) return A, bvec, k @@ -979,7 +979,7 @@ class Cylinder(QuadricMixin, Surface): """ _type = 'cylinder' - _coeff_keys = ('x0', 'y0', 'z0', 'r', 'dx', 'dy','dz') + _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): @@ -2405,4 +2405,3 @@ class Halfspace(Region): _SURFACE_CLASSES = {cls._type: cls for cls in Surface.__subclasses__()} - From 44d2527faf77df9f5dea25b7ae26fd4484e70133 Mon Sep 17 00:00:00 2001 From: Ethan Peterson Date: Wed, 19 Feb 2020 09:56:46 -0500 Subject: [PATCH 28/33] applied more suggestions from code review and added coefficient cache for general cones and cylinders --- openmc/surface.py | 209 ++++++++++++++++++++++++++++++++++------------ 1 file changed, 157 insertions(+), 52 deletions(-) diff --git a/openmc/surface.py b/openmc/surface.py index c73ce6fc87..918be232b7 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -187,6 +187,7 @@ class Surface(IDManagerMixin, metaclass=ABCMeta): def normalize(self, coeffs=None): """Normalize coefficients by first nonzero value + Parameters ---------- coeffs : tuple, optional @@ -222,12 +223,47 @@ class Surface(IDManagerMixin, metaclass=ABCMeta): return np.all(np.isclose(coeffs1, coeffs2, rtol=0., atol=self._atol)) @abstractmethod - def evaluate(self, point): - pass + def _get_base_coeffs(self): + """Return polynomial coefficients representing the implicit surface + equation. + + """ @abstractmethod - def translate(self, vector): - pass + def evaluate(self, point): + """Evaluate the surface equation at a given point. + + Parameters + ---------- + point : 3-tuple of float + The Cartesian coordinates, :math:`(x',y',z')`, at which the surface + equation should be evaluated. + + Returns + ------- + float + Evaluation of the surface polynomial at point :math:`(x',y',z')` + + """ + + @abstractmethod + def translate(self, vector, inplace=False): + """Translate surface in given direction + + Parameters + ---------- + vector : iterable of float + Direction in which surface should be translated + inplace : boolean + Whether or not to return a new instance of this Surface or to + modify the coefficients of this Surface. Defaults to False + + Returns + ------- + instance of openmc.Surface + Translated surface + + """ def to_xml_element(self): """Return XML representation of the surface @@ -325,13 +361,11 @@ class PlaneMixin(metaclass=ABCMeta): self._periodic_surface = periodic_surface periodic_surface._periodic_surface = self - @abstractmethod def _get_base_coeffs(self): """Return coefficients a, b, c, d representing a general plane of the form :math:`ax + by + cz = d`. """ - pass def evaluate(self, point): """Evaluate the surface equation at a given point. @@ -374,10 +408,7 @@ class PlaneMixin(metaclass=ABCMeta): a, b, c, d = self._get_base_coeffs() d = d + a*vx + b*vy + c*vz - if inplace: - surf = self - else: - surf = self.clone() + surf = self if inplace else self.clone() setattr(surf, surf._coeff_keys[-1], d) @@ -798,17 +829,13 @@ Plane.register(ZPlane) class QuadricMixin(metaclass=ABCMeta): """A Mixin class implementing common functionality for quadric surfaces""" - def __init__(self, **kwargs): - super().__init__(**kwargs) - @abstractmethod def _get_base_coeffs(self): """Return coefficients a, b, c, d, e, f, g, h, j, k representing a general quadric surface of the form: :math:`ax^2 + by^2 + cz^2 + dxy + eyz + fxz+ gx + hy + jz +k = 0`. """ - pass def get_Abc(self, coeffs=None): """Compute matrix, vector, and scalar coefficients for this surface or @@ -891,10 +918,7 @@ class QuadricMixin(metaclass=ABCMeta): """ vector = np.asarray(vector) - if inplace: - surf = self - else: - surf = self.clone() + surf = self if inplace else self.clone() if set(('x0', 'y0', 'z0')).intersection(set(surf._coeff_keys)): for vi, xi in zip(vector, ('x0', 'y0', 'z0')): @@ -916,8 +940,8 @@ class QuadricMixin(metaclass=ABCMeta): class Cylinder(QuadricMixin, Surface): """A cylinder with radius r, centered on the point (x0, y0, z0) with an - axis specified by the line through points (x0, y0, z0) and (x0+u, y0+v, - z0+w) + axis specified by the line through points (x0, y0, z0) and (x0+dx, y0+dy, + z0+dz) Parameters ---------- @@ -982,6 +1006,10 @@ 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) @@ -1020,55 +1048,90 @@ class Cylinder(QuadricMixin, Surface): def x0(self, x0): check_type('x0 coefficient', x0, Real) self._coefficients['x0'] = x0 + self._out_of_date = True @y0.setter def y0(self, y0): check_type('y0 coefficient', y0, Real) self._coefficients['y0'] = y0 + self._out_of_date = True @z0.setter def z0(self, z0): check_type('z0 coefficient', z0, Real) self._coefficients['z0'] = z0 + self._out_of_date = True @r.setter def r(self, r): check_type('r coefficient', r, Real) self._coefficients['r'] = r + self._out_of_date = True @dx.setter def dx(self, dx): check_type('dx coefficient', dx, Real) self._coefficients['dx'] = dx + self._out_of_date = True @dy.setter def dy(self, dy): check_type('dy coefficient', dy, Real) self._coefficients['dy'] = dy + self._out_of_date = True @dz.setter def dz(self, dz): check_type('dz coefficient', dz, Real) self._coefficients['dz'] = dz + self._out_of_date = True + + def _set_base_coeffs(self, coeffs): + """Set quadric coefficients for quick access + + Parameters + ---------- + coeffs : tuple + Tuple of Quadric coefficients (a, b, c, d, e, f, g, h, j, k) + """ + + for k, c in zip(Quadric._coeff_keys, coeffs): + setattr(self, '_' + k, c) + self._out_of_date = False def _get_base_coeffs(self): - x0, y0, z0, r = self.x0, self.y0, self.z0, self.r - dx, dy, dz = self.dx, self.dy, self.dz - dx2, dy2, dz2 = dx*dx, dy*dy, dz*dz + if self._out_of_date: + # 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 + r = self.r - a = dy2 + dz2 - b = dx2 + dz2 - c = dx2 + dy2 - d = -2*dx*dy - e = -2*dy*dz - f = -2*dx*dz - g = 2*(dx*(z0*dz + y0*dy) - x0*(dy2 + dz2)) - h = 2*(dy*(z0*dz + x0*dx) - y0*(dx2 + dz2)) - j = 2*(dz*(y0*dy + x0*dx) - z0*(dx2 + dy2)) - k = x0*x0*(dy2 + dz2) + y0*y0*(dx2 + dz2) + z0*z0*(dx2 + dy2) \ - + e*y0*z0 + f*x0*z0 + d*x0*y0 - r*r*(dx2 + dy2 + dz2) + # 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 - return (a, b, c, d, e, f, g, h, j, k) + # 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. + a = dy*dy + dz*dz + b = dx*dx + dz*dz + c = dx*dx + dy*dy + d = -2*dx*dy + e = -2*dy*dz + f = -2*dx*dz + g = 2*(cy*dz - cz*dy) + h = 2*(cz*dx - cx*dz) + j = 2*(cx*dy - cy*dx) + k = cx*cx + cy*cy + cz*cz - (dx*dx + dy*dy + dz*dz)*r*r + self._set_base_coeffs((a, b, c, d, e, f, g, h, j, k)) + + return tuple(getattr(self, '_' + k) for k in Quadric._coeff_keys) @classmethod def from_points(cls, p1, p2, r=1., **kwargs): @@ -1090,6 +1153,10 @@ 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) @@ -1624,6 +1691,9 @@ 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'), @@ -1666,55 +1736,92 @@ class Cone(QuadricMixin, Surface): def x0(self, x0): check_type('x0 coefficient', x0, Real) self._coefficients['x0'] = x0 + self._out_of_date = True @y0.setter def y0(self, y0): check_type('y0 coefficient', y0, Real) self._coefficients['y0'] = y0 + self._out_of_date = True @z0.setter def z0(self, z0): check_type('z0 coefficient', z0, Real) self._coefficients['z0'] = z0 + self._out_of_date = True @r2.setter def r2(self, r2): check_type('r^2 coefficient', r2, Real) self._coefficients['r2'] = r2 + self._out_of_date = True @dx.setter def dx(self, dx): check_type('dx coefficient', dx, Real) self._coefficients['dx'] = dx + self._out_of_date = True @dy.setter def dy(self, dy): check_type('dy coefficient', dy, Real) self._coefficients['dy'] = dy + self._out_of_date = True @dz.setter def dz(self, dz): check_type('dz coefficient', dz, Real) self._coefficients['dz'] = dz + self._out_of_date = True + + def _set_base_coeffs(self, coeffs): + """Set quadric coefficients for quick access + + Parameters + ---------- + coeffs : tuple + Tuple of Quadric coefficients (a, b, c, d, e, f, g, h, j, k) + """ + + for k, c in zip(Quadric._coeff_keys, coeffs): + setattr(self, '_' + k, c) + self._out_of_date = False def _get_base_coeffs(self): - x0, y0, z0, r2 = self.x0, self.y0, self.z0, self.r2 - dx, dy, dz = self.dx, self.dy, self.dz - cos2 = 1 / (1 + r2) - c1 = (dx*dx + dy*dy + dz*dz)*cos2 + # The equation for a general cone with vertex at point p = (x0, y0, z0) + # and axis specified by the unit vector d = (dx, dy, dz) and opening + # half angle theta can be described by the equation + # + # (d*(r - p))^2 - (r - p)*(r - p)cos^2(theta) = 0 + # + # where * is the dot product and the vector r is the evaulation point + # r = (x, y, z) + # + # The argument r2 for cones is actually tan^2(theta) so that + # cos^2(theta) = 1 / (1 + r2) - a = dx*dx - c1 - b = dy*dy - c1 - c = dz*dz - c1 - d = 2*dx*dy - e = 2*dy*dz - f = 2*dx*dz - g = -2*(dx*dx*x0 + dx*dy*y0 + dx*dz*z0 - c1) - h = -2*(dy*dy*y0 + dx*dy*x0 + dy*dz*z0 - c1) - j = -2*(dz*dz*y0 + dx*dz*x0 + dy*dz*y0 - c1) - k = (dx*x0 + dy*y0 + dz*z0)**2 - c1*(x0*x0 + y0*y0 + z0*z0) + if self._out_of_date: + 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) - return (a, b, c, d, e, f, g, h, j, k) + 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*y0 + dx*dz*x0 + dy*dz*y0 - cos2) + k = (dx*x0 + dy*y0 + dz*z0)**2 - cos2*(x0*x0 + y0*y0 + z0*z0) + self._set_base_coeffs((a, b, c, d, e, f, g, h, j, k)) + + return tuple(getattr(self, '_' + k) for k in Quadric._coeff_keys) class XCone(QuadricMixin, Surface): @@ -2404,5 +2511,3 @@ class Halfspace(Region): _SURFACE_CLASSES = {cls._type: cls for cls in Surface.__subclasses__()} - - From 6c2ef8b7f0590da2fbd5a20c6747d69ac541bb73 Mon Sep 17 00:00:00 2001 From: Ethan Peterson Date: Wed, 19 Feb 2020 14:23:12 -0500 Subject: [PATCH 29/33] removed coefficient caching --- openmc/surface.py | 147 ++++++++++++++-------------------------------- 1 file changed, 45 insertions(+), 102 deletions(-) diff --git a/openmc/surface.py b/openmc/surface.py index 578e902ed4..92c828d65f 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -361,12 +361,6 @@ class PlaneMixin(metaclass=ABCMeta): self._periodic_surface = periodic_surface periodic_surface._periodic_surface = self - def _get_base_coeffs(self): - """Return coefficients a, b, c, d representing a general plane of the - form :math:`ax + by + cz = d`. - - """ - def evaluate(self, point): """Evaluate the surface equation at a given point. @@ -830,13 +824,6 @@ Plane.register(ZPlane) class QuadricMixin(metaclass=ABCMeta): """A Mixin class implementing common functionality for quadric surfaces""" - def _get_base_coeffs(self): - """Return coefficients a, b, c, d, e, f, g, h, j, k representing a - general quadric surface of the form: - :math:`ax^2 + by^2 + cz^2 + dxy + eyz + fxz+ gx + hy + jz +k = 0`. - - """ - def get_Abc(self, coeffs=None): """Compute matrix, vector, and scalar coefficients for this surface or for a specified set of coefficients. @@ -1048,90 +1035,68 @@ class Cylinder(QuadricMixin, Surface): def x0(self, x0): check_type('x0 coefficient', x0, Real) self._coefficients['x0'] = x0 - self._out_of_date = True @y0.setter def y0(self, y0): check_type('y0 coefficient', y0, Real) self._coefficients['y0'] = y0 - self._out_of_date = True @z0.setter def z0(self, z0): check_type('z0 coefficient', z0, Real) self._coefficients['z0'] = z0 - self._out_of_date = True @r.setter def r(self, r): check_type('r coefficient', r, Real) self._coefficients['r'] = r - self._out_of_date = True @dx.setter def dx(self, dx): check_type('dx coefficient', dx, Real) self._coefficients['dx'] = dx - self._out_of_date = True @dy.setter def dy(self, dy): check_type('dy coefficient', dy, Real) self._coefficients['dy'] = dy - self._out_of_date = True @dz.setter def dz(self, dz): check_type('dz coefficient', dz, Real) self._coefficients['dz'] = dz - self._out_of_date = True - - def _set_base_coeffs(self, coeffs): - """Set quadric coefficients for quick access - - Parameters - ---------- - coeffs : tuple - Tuple of Quadric coefficients (a, b, c, d, e, f, g, h, j, k) - """ - - for k, c in zip(Quadric._coeff_keys, coeffs): - setattr(self, '_' + k, c) - self._out_of_date = False def _get_base_coeffs(self): - if self._out_of_date: - # 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 - r = self.r + # 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 + r = self.r - # 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 + # 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. - a = dy*dy + dz*dz - b = dx*dx + dz*dz - c = dx*dx + dy*dy - d = -2*dx*dy - e = -2*dy*dz - f = -2*dx*dz - g = 2*(cy*dz - cz*dy) - h = 2*(cz*dx - cx*dz) - j = 2*(cx*dy - cy*dx) - k = cx*cx + cy*cy + cz*cz - (dx*dx + dy*dy + dz*dz)*r*r - self._set_base_coeffs((a, b, c, d, e, f, g, h, j, k)) + # 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. + a = dy*dy + dz*dz + b = dx*dx + dz*dz + c = dx*dx + dy*dy + d = -2*dx*dy + e = -2*dy*dz + f = -2*dx*dz + g = 2*(cy*dz - cz*dy) + h = 2*(cz*dx - cx*dz) + j = 2*(cx*dy - cy*dx) + k = cx*cx + cy*cy + cz*cz - (dx*dx + dy*dy + dz*dz)*r*r - return tuple(getattr(self, '_' + k) for k in Quadric._coeff_keys) + return (a, b, c, d, e, f, g, h, j, k) @classmethod def from_points(cls, p1, p2, r=1., **kwargs): @@ -1736,56 +1701,36 @@ class Cone(QuadricMixin, Surface): def x0(self, x0): check_type('x0 coefficient', x0, Real) self._coefficients['x0'] = x0 - self._out_of_date = True @y0.setter def y0(self, y0): check_type('y0 coefficient', y0, Real) self._coefficients['y0'] = y0 - self._out_of_date = True @z0.setter def z0(self, z0): check_type('z0 coefficient', z0, Real) self._coefficients['z0'] = z0 - self._out_of_date = True @r2.setter def r2(self, r2): check_type('r^2 coefficient', r2, Real) self._coefficients['r2'] = r2 - self._out_of_date = True @dx.setter def dx(self, dx): check_type('dx coefficient', dx, Real) self._coefficients['dx'] = dx - self._out_of_date = True @dy.setter def dy(self, dy): check_type('dy coefficient', dy, Real) self._coefficients['dy'] = dy - self._out_of_date = True @dz.setter def dz(self, dz): check_type('dz coefficient', dz, Real) self._coefficients['dz'] = dz - self._out_of_date = True - - def _set_base_coeffs(self, coeffs): - """Set quadric coefficients for quick access - - Parameters - ---------- - coeffs : tuple - Tuple of Quadric coefficients (a, b, c, d, e, f, g, h, j, k) - """ - - for k, c in zip(Quadric._coeff_keys, coeffs): - setattr(self, '_' + k, c) - self._out_of_date = False def _get_base_coeffs(self): # The equation for a general cone with vertex at point p = (x0, y0, z0) @@ -1800,28 +1745,26 @@ class Cone(QuadricMixin, Surface): # The argument r2 for cones is actually tan^2(theta) so that # cos^2(theta) = 1 / (1 + r2) - if self._out_of_date: - 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, 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) - 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*y0 + dx*dz*x0 + dy*dz*y0 - cos2) - k = (dx*x0 + dy*y0 + dz*z0)**2 - cos2*(x0*x0 + y0*y0 + z0*z0) - self._set_base_coeffs((a, b, c, d, e, f, g, h, j, k)) + 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*y0 + dx*dz*x0 + dy*dz*y0 - cos2) + k = (dx*x0 + dy*y0 + dz*z0)**2 - cos2*(x0*x0 + y0*y0 + z0*z0) - return tuple(getattr(self, '_' + k) for k in Quadric._coeff_keys) + return (a, b, c, d, e, f, g, h, j, k) class XCone(QuadricMixin, Surface): From 6bd4c5cd900ab38088faa68024d29f46e5c72552 Mon Sep 17 00:00:00 2001 From: Sam Powell-Gill Date: Thu, 20 Feb 2020 08:48:26 +0000 Subject: [PATCH 30/33] Added suggested changes for #1489 --- tests/unit_tests/test_data_photon.py | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/tests/unit_tests/test_data_photon.py b/tests/unit_tests/test_data_photon.py index 02635bcc9a..2fa0f85a7b 100644 --- a/tests/unit_tests/test_data_photon.py +++ b/tests/unit_tests/test_data_photon.py @@ -2,6 +2,7 @@ from collections.abc import Mapping, Callable import os +from pathlib import Path import numpy as np import pandas as pd @@ -145,10 +146,8 @@ def test_export_to_hdf5(tmpdir, element): # Export to hdf5 again element2.export_to_hdf5(filename, 'w') -def test_photodat_only(tmpdir): - endf_data = os.environ['OPENMC_ENDF_DATA'] - filename = str(tmpdir.join('tmp.h5')) - p_file = 'photoat-{:03}_{}_000.endf'.format(1, 'H') - p_path = os.path.join(endf_data, 'photoat', p_file) - data=openmc.data.IncidentPhoton.from_endf(p_path) - data.export_to_hdf5(filename, 'w') \ No newline at end of file +def test_photodat_only(run_in_tmpdir): + endf_dir = Path(os.environ['OPENMC_ENDF_DATA']) + photoatomic_file = endf_dir / 'photoat/photoat-001_H_000.endf' + data=openmc.data.IncidentPhoton.from_endf(photoatomic_file) + data.export_to_hdf5('tmp.h5', 'w') \ No newline at end of file From 0fbfad4d49c543f9203967bb7f9ec5a5225d3734 Mon Sep 17 00:00:00 2001 From: Sam Powell-Gill Date: Thu, 20 Feb 2020 16:29:08 +0000 Subject: [PATCH 31/33] Made requested formatting change --- tests/unit_tests/test_data_photon.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/unit_tests/test_data_photon.py b/tests/unit_tests/test_data_photon.py index 2fa0f85a7b..f7274e9c19 100644 --- a/tests/unit_tests/test_data_photon.py +++ b/tests/unit_tests/test_data_photon.py @@ -148,6 +148,6 @@ def test_export_to_hdf5(tmpdir, element): def test_photodat_only(run_in_tmpdir): endf_dir = Path(os.environ['OPENMC_ENDF_DATA']) - photoatomic_file = endf_dir / 'photoat/photoat-001_H_000.endf' - data=openmc.data.IncidentPhoton.from_endf(photoatomic_file) + photoatomic_file = endf_dir / 'photoat' / 'photoat-001_H_000.endf' + data = openmc.data.IncidentPhoton.from_endf(photoatomic_file) data.export_to_hdf5('tmp.h5', 'w') \ No newline at end of file From 48ee3e78ec9e82919d2e3cf11b7ba2d8dbfbe20c Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 20 Feb 2020 13:09:15 -0600 Subject: [PATCH 32/33] Make sure DAGMC is discoverable if linking against OpenMC --- CMakeLists.txt | 24 ++++++++++++------- cmake/Modules/FindDAGMC.cmake | 10 ++++---- ...enMCConfig.cmake => OpenMCConfig.cmake.in} | 7 ++++++ 3 files changed, 28 insertions(+), 13 deletions(-) rename cmake/{OpenMCConfig.cmake => OpenMCConfig.cmake.in} (58%) diff --git a/CMakeLists.txt b/CMakeLists.txt index ec6b500b45..c776a7747c 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -9,6 +9,11 @@ set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/bin) # Set module path set(CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/cmake/Modules) +# Allow user to specify _ROOT variables +if (NOT (CMAKE_VERSION VERSION_LESS 3.12)) + cmake_policy(SET CMP0074 NEW) +endif() + #=============================================================================== # Command line options #=============================================================================== @@ -35,18 +40,18 @@ endif() #=============================================================================== if(dagmc) find_package(DAGMC REQUIRED) - link_directories(${DAGMC_LIBRARY_DIRS}) + add_library(dagmc-imported INTERFACE IMPORTED) + target_link_directories(dagmc-imported INTERFACE ${DAGMC_LIBRARY_DIRS}) + target_link_libraries(dagmc-imported INTERFACE ${DAGMC_LIBRARIES}) + target_include_directories(dagmc-imported INTERFACE ${DAGMC_INCLUDE_DIRS}) +else() + set(DAGMC_FOUND false) endif() #=============================================================================== # HDF5 for binary output #=============================================================================== -# Allow user to specify HDF5_ROOT -if (NOT (CMAKE_VERSION VERSION_LESS 3.12)) - cmake_policy(SET CMP0074 NEW) -endif() - # Unfortunately FindHDF5.cmake will always prefer a serial HDF5 installation # over a parallel installation if both appear on the user's PATH. To get around # this, we check for the environment variable HDF5_ROOT and if it exists, use it @@ -357,8 +362,7 @@ target_link_libraries(libopenmc ${ldflags} ${HDF5_LIBRARIES} ${HDF5_HL_LIBRARIES if(dagmc) target_compile_definitions(libopenmc PRIVATE DAGMC) - target_link_libraries(libopenmc ${DAGMC_LIBRARIES}) - target_include_directories(libopenmc PRIVATE ${DAGMC_INCLUDE_DIRS}) + target_link_libraries(libopenmc dagmc-imported) endif() #=============================================================================== @@ -388,6 +392,8 @@ add_custom_command(TARGET libopenmc POST_BUILD # Install executable, scripts, manpage, license #=============================================================================== +configure_file(cmake/OpenMCConfig.cmake.in "${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/OpenMCConfig.cmake" @ONLY) + set(INSTALL_CONFIGDIR ${CMAKE_INSTALL_LIBDIR}/cmake/OpenMC) install(TARGETS openmc libopenmc faddeeva EXPORT openmc-targets @@ -401,7 +407,7 @@ install(EXPORT openmc-targets DESTINATION ${INSTALL_CONFIGDIR}) install(DIRECTORY src/relaxng DESTINATION ${CMAKE_INSTALL_DATADIR}/openmc) -install(FILES cmake/OpenMCConfig.cmake DESTINATION ${INSTALL_CONFIGDIR}) +install(FILES "${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/OpenMCConfig.cmake" DESTINATION ${INSTALL_CONFIGDIR}) install(FILES man/man1/openmc.1 DESTINATION ${CMAKE_INSTALL_MANDIR}/man1) install(FILES LICENSE DESTINATION "${CMAKE_INSTALL_DOCDIR}" RENAME copyright) install(DIRECTORY include/ DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}) diff --git a/cmake/Modules/FindDAGMC.cmake b/cmake/Modules/FindDAGMC.cmake index c644c8886e..bb9bc56ddb 100644 --- a/cmake/Modules/FindDAGMC.cmake +++ b/cmake/Modules/FindDAGMC.cmake @@ -12,7 +12,9 @@ find_path(DAGMC_CMAKE_CONFIG NAMES DAGMCConfig.cmake PATHS ENV LD_LIBRARY_PATH PATH_SUFFIXES lib Lib cmake lib/cmake NO_DEFAULT_PATH) - -message(STATUS "Found DAGMC in ${DAGMC_CMAKE_CONFIG}") - -include(${DAGMC_CMAKE_CONFIG}/DAGMCConfig.cmake) +if(DAGMC_CMAKE_CONFIG) + message(STATUS "Found DAGMC in ${DAGMC_CMAKE_CONFIG}") + include(${DAGMC_CMAKE_CONFIG}/DAGMCConfig.cmake) +else() + message(WARNING "Cound not find DAGMC") +endif() diff --git a/cmake/OpenMCConfig.cmake b/cmake/OpenMCConfig.cmake.in similarity index 58% rename from cmake/OpenMCConfig.cmake rename to cmake/OpenMCConfig.cmake.in index 29a0e4542f..ad00a28616 100644 --- a/cmake/OpenMCConfig.cmake +++ b/cmake/OpenMCConfig.cmake.in @@ -5,6 +5,13 @@ find_package(gsl-lite REQUIRED HINTS ${OpenMC_CMAKE_DIR}/../gsl-lite) find_package(pugixml REQUIRED HINTS ${OpenMC_CMAKE_DIR}/../pugixml) find_package(xtl REQUIRED HINTS ${OpenMC_CMAKE_DIR}/../xtl) find_package(xtensor REQUIRED HINTS ${OpenMC_CMAKE_DIR}/../xtensor) +if(@DAGMC_FOUND@) + find_package(DAGMC REQUIRED HINTS @DAGMC_LIBRARY_DIRS@) + add_library(dagmc-imported INTERFACE IMPORTED) + target_link_directories(dagmc-imported INTERFACE ${DAGMC_LIBRARY_DIRS}) + target_link_libraries(dagmc-imported INTERFACE ${DAGMC_LIBRARIES}) + target_include_directories(dagmc-imported INTERFACE ${DAGMC_INCLUDE_DIRS}) +endif() if(NOT TARGET OpenMC::libopenmc) include("${OpenMC_CMAKE_DIR}/OpenMCTargets.cmake") From c7c3db82467869c85a5249e5567ee12d0cd2c620 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sun, 23 Feb 2020 22:15:28 -0600 Subject: [PATCH 33/33] Don't use target_link_directorires (only supported in cmake 3.13+) --- CMakeLists.txt | 2 +- cmake/OpenMCConfig.cmake.in | 4 +++- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index c776a7747c..29eed4fb04 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -41,7 +41,7 @@ endif() if(dagmc) find_package(DAGMC REQUIRED) add_library(dagmc-imported INTERFACE IMPORTED) - target_link_directories(dagmc-imported INTERFACE ${DAGMC_LIBRARY_DIRS}) + link_directories(${DAGMC_LIBRARY_DIRS}) target_link_libraries(dagmc-imported INTERFACE ${DAGMC_LIBRARIES}) target_include_directories(dagmc-imported INTERFACE ${DAGMC_INCLUDE_DIRS}) else() diff --git a/cmake/OpenMCConfig.cmake.in b/cmake/OpenMCConfig.cmake.in index ad00a28616..1c6a488f54 100644 --- a/cmake/OpenMCConfig.cmake.in +++ b/cmake/OpenMCConfig.cmake.in @@ -8,8 +8,10 @@ find_package(xtensor REQUIRED HINTS ${OpenMC_CMAKE_DIR}/../xtensor) if(@DAGMC_FOUND@) find_package(DAGMC REQUIRED HINTS @DAGMC_LIBRARY_DIRS@) add_library(dagmc-imported INTERFACE IMPORTED) - target_link_directories(dagmc-imported INTERFACE ${DAGMC_LIBRARY_DIRS}) target_link_libraries(dagmc-imported INTERFACE ${DAGMC_LIBRARIES}) + foreach(dir ${DAGMC_LIBRARY_DIRS}) + target_link_libraries(dagmc-imported INTERFACE "-L${dir}") + endforeach() target_include_directories(dagmc-imported INTERFACE ${DAGMC_INCLUDE_DIRS}) endif()