From 394ac75560646ac8214583a3fb581f8f66e6516f Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 6 Mar 2019 14:27:40 -0600 Subject: [PATCH] Add translate methods to surfaces and regions --- openmc/region.py | 85 +++++++++- openmc/surface.py | 263 ++++++++++++++++++++++++++++++- tests/unit_tests/test_region.py | 18 +++ tests/unit_tests/test_surface.py | 55 +++++++ 4 files changed, 409 insertions(+), 12 deletions(-) diff --git a/openmc/region.py b/openmc/region.py index ceffadbe36..ffcbeea600 100644 --- a/openmc/region.py +++ b/openmc/region.py @@ -229,14 +229,27 @@ class Region(metaclass=ABCMeta): clone : openmc.Region The clone of this region - Raises - ------ - NotImplementedError - This method is not implemented for the abstract region class. + """ + pass + + @abstractmethod + def translate(self, vector, memo=None): + """Translate region in given direction + + Parameters + ---------- + vector : iterable of float + Direction in which region should be translated + memo : dict or None + Dictionary used for memoization + + Returns + ------- + openmc.Region + Translated region """ - raise NotImplementedError('The clone method is not implemented for ' - 'the abstract region class.') + pass class Intersection(Region, MutableSequence): @@ -352,6 +365,26 @@ class Intersection(Region, MutableSequence): clone[:] = [n.clone(memo) for n in self] return clone + def translate(self, vector, memo=None): + """Translate region in given direction + + Parameters + ---------- + vector : iterable of float + Direction in which region should be translated + memo : dict or None + Dictionary used for memoization + + Returns + ------- + openmc.Intersection + Translated region + + """ + if memo is None: + memo = {} + return type(self)(n.translate(vector, memo) for n in self) + class Union(Region, MutableSequence): r"""Union of two or more regions. @@ -464,6 +497,26 @@ class Union(Region, MutableSequence): clone[:] = [n.clone(memo) for n in self] return clone + def translate(self, vector, memo=None): + """Translate region in given direction + + Parameters + ---------- + vector : iterable of float + Direction in which region should be translated + memo : dict or None + Dictionary used for memoization + + Returns + ------- + openmc.Union + Translated region + + """ + if memo is None: + memo = {} + return type(self)(n.translate(vector, memo) for n in self) + class Complement(Region): """Complement of a region. @@ -584,3 +637,23 @@ class Complement(Region): clone = deepcopy(self) clone.node = self.node.clone(memo) return clone + + def translate(self, vector, memo=None): + """Translate region in given direction + + Parameters + ---------- + vector : iterable of float + Direction in which region should be translated + memo : dict or None + Dictionary used for memoization + + Returns + ------- + openmc.Complement + Translated region + + """ + if memo is None: + memo = {} + return type(self)(self.node.translate(vector, memo)) diff --git a/openmc/surface.py b/openmc/surface.py index 7efa04e384..306fb6d735 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -1,4 +1,4 @@ -from abc import ABCMeta +from abc import ABCMeta, abstractmethod from collections import OrderedDict from copy import deepcopy from functools import partial @@ -15,7 +15,7 @@ from openmc.mixin import IDManagerMixin _BOUNDARY_TYPES = ['transmission', 'vacuum', 'reflective', 'periodic'] -class Surface(IDManagerMixin): +class Surface(IDManagerMixin, metaclass=ABCMeta): """An implicit surface with an associated boundary condition. An implicit surface is defined as the set of zeros of a function of the @@ -142,7 +142,6 @@ class Surface(IDManagerMixin): desired half-space """ - return (np.array([-np.inf, -np.inf, -np.inf]), np.array([np.inf, np.inf, np.inf])) @@ -175,6 +174,14 @@ class Surface(IDManagerMixin): return memo[self] + @abstractmethod + def evaluate(self, point): + pass + + @abstractmethod + def translate(self, vector): + pass + def to_xml_element(self): """Return XML representation of the surface @@ -432,13 +439,34 @@ class Plane(Surface): Returns ------- float - :math:`Ax' + By' + Cz' - d` + :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 @@ -593,6 +621,26 @@ class XPlane(Plane): """ 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) + class YPlane(Plane): """A plane perpendicular to the y axis of the form :math:`y - y_0 = 0` @@ -699,6 +747,26 @@ class YPlane(Plane): """ 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) + class ZPlane(Plane): """A plane perpendicular to the z axis of the form :math:`z - z_0 = 0` @@ -805,8 +873,28 @@ class ZPlane(Plane): """ return point[2] - self.z0 + def translate(self, vector): + """Translate surface in given direction -class Cylinder(Surface, metaclass=ABCMeta): + Parameters + ---------- + vector : iterable of float + Direction in which surface should be translated + + Returns + ------- + openmc.ZPlane + Translated surface + + """ + vz = vector[2] + if vz == 0.0: + return self + else: + return type(self)(z0=self.z0 + vz) + + +class Cylinder(Surface): """A cylinder whose length is parallel to the x-, y-, or z-axis. Parameters @@ -977,6 +1065,28 @@ class XCylinder(Cylinder): 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): """An infinite cylinder whose length is parallel to the y-axis of the form @@ -1099,6 +1209,28 @@ class YCylinder(Cylinder): 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): """An infinite cylinder whose length is parallel to the z-axis of the form @@ -1221,6 +1353,28 @@ class ZCylinder(Cylinder): 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): """A sphere of the form :math:`(x - x_0)^2 + (y - y_0)^2 + (z - z_0)^2 = R^2`. @@ -1369,8 +1523,31 @@ class Sphere(Surface): z = point[2] - self.z0 return x**2 + y**2 + z**2 - self.r**2 + def translate(self, vector): + """Translate surface in given direction -class Cone(Surface, metaclass=ABCMeta): + 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): """A conical surface parallel to the x-, y-, or z-axis. Parameters @@ -1463,6 +1640,29 @@ class Cone(Surface, metaclass=ABCMeta): check_type('R^2 coefficient', R2, Real) self._coefficients['R2'] = R2 + def translate(self, vector): + """Translate surface in given direction + + Parameters + ---------- + vector : iterable of float + Direction in which surface should be translated + + Returns + ------- + openmc.Cone + 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, R2=self.r2) + class XCone(Cone): """A cone parallel to the x-axis of the form :math:`(y - y_0)^2 + (z - z_0)^2 = @@ -1842,6 +2042,30 @@ class Quadric(Surface): 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. @@ -1987,3 +2211,30 @@ class Halfspace(Region): clone = deepcopy(self) clone.surface = self.surface.clone(memo) return clone + + def translate(self, vector, memo=None): + """Translate half-space in given direction + + Parameters + ---------- + vector : iterable of float + Direction in which region should be translated + memo : dict or None + Dictionary used for memoization + + Returns + ------- + openmc.Halfspace + Translated half-space + + """ + if memo is None: + memo = {} + + # If translated surface not in memo, add it + key = (self.surface, tuple(vector)) + if key not in memo: + memo[key] = self.surface.translate(vector) + + # Return translated surface + return type(self)(memo[key], self.side) diff --git a/tests/unit_tests/test_region.py b/tests/unit_tests/test_region.py index 377c8fbd58..b623cd6061 100644 --- a/tests/unit_tests/test_region.py +++ b/tests/unit_tests/test_region.py @@ -34,6 +34,12 @@ def test_union(reset): assert (6, -1, 0) not in reg2 assert str(reg2) == '((1 | -2) 3)' + # translate method + regt = region.translate((2.0, 0.0, 0.0)) + assert (-4, 0, 0) in regt + assert (6, 0, 0) not in regt + assert (8, 0, 0) in regt + def test_intersection(reset): s1 = openmc.XPlane(surface_id=1, x0=5) @@ -61,6 +67,12 @@ def test_intersection(reset): assert (-6, -2, 0) not in reg2 assert str(reg2) == '((-1 2) | 3)' + # translate method + regt = region.translate((2.0, 0.0, 0.0)) + assert (-4, 0, 0) not in regt + assert (6, 0, 0) in regt + assert (8, 0, 0) not in regt + def test_complement(reset): zcyl = openmc.ZCylinder(surface_id=1, R=1.) @@ -88,6 +100,12 @@ def test_complement(reset): assert (0, 0, 6) not in inside assert (0, 0, 6) in outside + # translate method + inside_t = inside.translate((1.0, 1.0, 1.0)) + ll, ur = inside_t.bounding_box + assert ll == pytest.approx((0., 0., -4.)) + assert ur == pytest.approx((2., 2., 6.)) + def test_get_surfaces(): s1 = openmc.XPlane() diff --git a/tests/unit_tests/test_surface.py b/tests/unit_tests/test_surface.py index 01a7cf7108..45fecb4936 100644 --- a/tests/unit_tests/test_surface.py +++ b/tests/unit_tests/test_surface.py @@ -29,6 +29,10 @@ def test_plane(): x, y, z = (4, 3, 6) assert s.evaluate((x, y, z)) == pytest.approx(s.a*x + s.b*y + s.c*z - s.d) + # translate method + st = s.translate((1.0, 0.0, 0.0)) + assert (st.a, st.b, st.c, st.d) == (s.a, s.b, s.c, 4) + # Make sure repr works repr(s) @@ -69,6 +73,10 @@ def test_xplane(): # evaluate method assert s.evaluate((5., 0., 0.)) == pytest.approx(2.) + # translate method + st = s.translate((1.0, 0.0, 0.0)) + assert st.x0 == s.x0 + 1 + # Make sure repr works repr(s) @@ -94,6 +102,10 @@ def test_yplane(): # evaluate method assert s.evaluate((0., 0., 0.)) == pytest.approx(-3.) + # translate method + st = s.translate((0.0, 1.0, 0.0)) + assert st.y0 == s.y0 + 1 + def test_zplane(): s = openmc.ZPlane(z0=3.) @@ -116,6 +128,10 @@ def test_zplane(): # evaluate method assert s.evaluate((0., 0., 10.)) == pytest.approx(7.) + # translate method + st = s.translate((0.0, 0.0, 1.0)) + assert st.z0 == s.z0 + 1 + # Make sure repr works repr(s) @@ -138,6 +154,12 @@ def test_xcylinder(): # evaluate method assert s.evaluate((0, y, z)) == pytest.approx(-r**2) + # translate method + st = s.translate((1.0, 1.0, 1.0)) + assert st.y0 == s.y0 + 1 + assert st.z0 == s.z0 + 1 + assert st.r == s.r + # Make sure repr works repr(s) @@ -169,6 +191,12 @@ def test_ycylinder(): # evaluate method assert s.evaluate((x, 0, z)) == pytest.approx(-r**2) + # translate method + st = s.translate((1.0, 1.0, 1.0)) + assert st.x0 == s.x0 + 1 + assert st.z0 == s.z0 + 1 + assert st.r == s.r + def test_zcylinder(): x, y, r = 3, 5, 2 @@ -188,6 +216,12 @@ def test_zcylinder(): # evaluate method assert s.evaluate((x, y, 0)) == pytest.approx(-r**2) + # translate method + st = s.translate((1.0, 1.0, 1.0)) + assert st.x0 == s.x0 + 1 + assert st.y0 == s.y0 + 1 + assert st.r == s.r + # Make sure repr works repr(s) @@ -211,6 +245,13 @@ def test_sphere(): # evaluate method assert s.evaluate((x, y, z)) == pytest.approx(-r**2) + # translate method + st = s.translate((1.0, 1.0, 1.0)) + assert st.x0 == s.x0 + 1 + assert st.y0 == s.y0 + 1 + assert st.z0 == s.z0 + 1 + assert st.r == s.r + # Make sure repr works repr(s) @@ -229,6 +270,13 @@ def cone_common(apex, r2, cls): # evaluate method -- should be zero at apex assert s.evaluate((x, y, z)) == pytest.approx(0.0) + # translate method + st = s.translate((1.0, 1.0, 1.0)) + assert st.x0 == s.x0 + 1 + assert st.y0 == s.y0 + 1 + assert st.z0 == s.z0 + 1 + assert st.r2 == s.r2 + # Make sure repr works repr(s) @@ -272,6 +320,13 @@ def test_quadric(): assert s.evaluate((0., 0., 0.)) == pytest.approx(coeffs['k']) assert s.evaluate((1., 1., 1.)) == pytest.approx(3 + coeffs['k']) + # translate method + st = s.translate((1.0, 1.0, 1.0)) + for coeff in 'abcdef': + assert getattr(s, coeff) == getattr(st, coeff) + assert (st.g, st.h, st.j) == (-2, -2, -2) + assert st.k == s.k + 3 + def test_cylinder_from_points(): # Generate 45-degree rotated cylinder in x-y plane with radius 1