From 3aceb2dc94651f670e97f28d6f9e261aa8e7d310 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 26 Feb 2020 12:07:02 -0600 Subject: [PATCH 1/9] Make sure gets installed. Closes #1496 --- CMakeLists.txt | 2 ++ 1 file changed, 2 insertions(+) diff --git a/CMakeLists.txt b/CMakeLists.txt index 29eed4fb04..f251d75751 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -170,6 +170,8 @@ target_link_libraries(xtensor INTERFACE xtl) # GSL header-only library #=============================================================================== +set(GSL_LITE_OPT_INSTALL_COMPAT_HEADER ON CACHE BOOL + "Install MS-GSL compatibility header ") add_subdirectory(vendor/gsl-lite) # Make sure contract violations throw exceptions From a1f7cb7387e9920edf3e101d4b0e6ef36b934caf Mon Sep 17 00:00:00 2001 From: Ethan Peterson Date: Fri, 21 Feb 2020 20:49:06 -0500 Subject: [PATCH 2/9] Implementing Surface Rotations --- openmc/surface.py | 226 ++++++++++++++++++++++++++++++- tests/unit_tests/test_surface.py | 31 +++++ 2 files changed, 254 insertions(+), 3 deletions(-) diff --git a/openmc/surface.py b/openmc/surface.py index 567968b2c8..dd2db91f27 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -1,5 +1,6 @@ from abc import ABCMeta, abstractmethod from collections import OrderedDict +from collections.abc import Iterable from copy import deepcopy from numbers import Real from xml.etree import ElementTree as ET @@ -8,7 +9,7 @@ import math 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, IDWarning @@ -26,6 +27,62 @@ will not accept positional parameters for superclass arguments.\ """ +def get_rotation_matrix(rotation, order='xyz'): + """Generate a 3x3 rotation matrix from input angles + + Parameters + ---------- + rotation : 3-tuple of float + A 3-tuple of angles (psi, theta, phi) in degrees where the first + element is the rotation about the x-axis in the fixed laboratory frame, + the second element is the rotation about the y-axis in the fixed + laboratory frame, and the third element is the rotation about the + z-axis in the fixed laboratory frame. The rotations are active + rotations. + order : str, optinoal + A string of 'x', 'y', and 'z' in some order specifying which rotation + to perform first, second, and third. Defaults to 'zyx' which means, the + rotation by angle phi about z will be applied first, followed by theta + about y and then psi about x. + + """ + check_type('surface rotation', rotation, Iterable, Real) + check_length('surface rotation', rotation, 3) + + # Calculate rotation matrix from angles psi, theta, phi, and the order + # of the desired rotations - order= 'xyz' applies the rotation matrix + # R = Rz(phi) @ Ry(theta) @ Rx(psi) to a vector from the left + # order is the extrinsic order (actual multiplication of matrices). This + # means that the associated intrinsic rotation swaps the order of + # multiplication i.e. an explicit Rz(psi) @ Ry(theta) @ Rx(phi) is known + # as an x-y-z rotation and its associated intrinsic rotation is referred + # to as a z-y'x'' rotation. + # + # Default order produces the same active rotation as the cell.rotation + # method except that method introduces a minus sign to all the angles after + # they are read in which is bad + phi, theta, psi = np.array(rotation)*(np.pi/180.) + cx, sx = np.cos(phi), np.sin(phi) + cy, sy = np.cos(theta), np.sin(theta) + cz, sz = np.cos(psi), np.sin(psi) + Rx = np.array([[1., 0., 0.], [0., cx, -sx], [0., sx, cx]]) + Ry = np.array([[cy, 0., sy], [0., 1., 0.], [-sy, 0., cy]]) + Rz = np.array([[cz, -sz, 0.], [sz, cz, 0.], [0., 0., 1.]]) + angledct = {'x': Rx, 'y': Ry, 'z': Rz} + R1, R2, R3 = tuple(angledct[i] for i in order) + return R3 @ R2 @ R1 + + +def get_rotation_matrix_2(rotation): + phi, theta, psi = np.array(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) + return 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]]) + + class Surface(IDManagerMixin, metaclass=ABCMeta): """An implicit surface with an associated boundary condition. @@ -266,6 +323,28 @@ class Surface(IDManagerMixin, metaclass=ABCMeta): """ + @abstractmethod + def rotate(self, rotation, pivot=(0., 0., 0.), order='zyx', inplace=False): + """Rotate surface by given Tait-Bryan angles + + Parameters + ---------- + rotation : iterable of float + Intrinsic Tait-Bryan angles in degrees used to rotate the surface + pivot : iterable of float, optional + (x, y, z) coordinates for the point to rotate about. Defaults to + (0., 0., 0.) + inplace : boolean + Whether or not to return a new instance of Surface or to modify the + coefficients of this Surface in place. Defaults to False. + + Returns + ------- + openmc.Surface + Rotated surface + + """ + def to_xml_element(self): """Return XML representation of the surface @@ -414,6 +493,33 @@ class PlaneMixin(metaclass=ABCMeta): ur = np.array([v if not np.isnan(v) else np.inf for v in vals]) return (ll, ur) +======= + a, b, c = self._get_base_coeffs()[0:3] + return np.array((a, b, c)) / np.sqrt(a*a + b*b + c*c) + + def bounding_box(self, side): + # Compute the bounding box based on the normal vector to the plane + nhat = self._get_normal() + lb = np.array([-np.inf, -np.inf, -np.inf]) + ub = np.array([np.inf, np.inf, np.inf]) + # If the plane is axis aligned, find the proper bounding box + if np.any(np.isclose(np.abs(nhat), 1., rtol=0., atol=self._atol)): + sign = np.dot(np.array([1., 1., 1.]), nhat) + a, b, c, d = self._get_base_coeffs() + vals = tuple(d/i if round(i) != 0 else np.nan for i in (a, b, c)) + if side == '-': + if sign > 0: + ub = np.array([v if ~np.isnan(v) else np.inf for v in vals]) + else: + lb = np.array([v if ~np.isnan(v) else -np.inf for v in vals]) + elif side == '+': + if sign > 0: + lb = np.array([v if ~np.isnan(v) else -np.inf for v in vals]) + else: + ub = np.array([v if ~np.isnan(v) else np.inf for v in vals]) + + return (lb, ub) +>>>>>>> Implementing Surface Rotations def evaluate(self, point): """Evaluate the surface equation at a given point. @@ -452,9 +558,11 @@ class PlaneMixin(metaclass=ABCMeta): Translated surface """ - vx, vy, vz = vector + if np.allclose(vector, 0., rtol=0., atol=self._atol): + return self + a, b, c, d = self._get_base_coeffs() - d = d + a*vx + b*vy + c*vz + d = d + np.array([a, b, c]) @ np.array(vector) surf = self if inplace else self.clone() @@ -462,6 +570,43 @@ class PlaneMixin(metaclass=ABCMeta): return surf + def rotate(self, rotation, pivot=(0., 0., 0.), order='zyx', inplace=False): + """Rotate surface by given Tait-Bryan angles + + Parameters + ---------- + rotation : iterable of float + Intrinsic Tait-Bryan angles in degrees used to rotate the surface + pivot : iterable of float, optional + (x, y, z) coordinates for the point to rotate about. Defaults to + (0., 0., 0.) + inplace : boolean + Whether or not to return a new instance of a Plane or to modify the + coefficients of this plane in place. Defaults to False. + + Returns + ------- + openmc.Plane + Rotated surface + + """ + + pivot = np.asarray(pivot) + surf = self.translate(-pivot, inplace=inplace) + + Rmat = get_rotation_matrix(rotation, order=order) + a, b, c, d = surf._get_base_coeffs() + # Compute new rotated coefficients a, b, c + a, b, c = Rmat @ np.array([a, b, c]) + + kwargs = {'boundary_type': surf.boundary_type, 'name': surf.name} + if inplace: + kwargs['sufrace_id'] = surf.id + + surf = Plane(a=a, b=b, c=c, d=d, **kwargs) + + return surf.translate(pivot, inplace=inplace) + def to_xml_element(self): """Return XML representation of the surface @@ -976,6 +1121,8 @@ class QuadricMixin(metaclass=ABCMeta): """ vector = np.asarray(vector) + if np.allclose(vector, 0., rtol=0., atol=self._atol): + return self surf = self if inplace else self.clone() @@ -999,6 +1146,72 @@ class QuadricMixin(metaclass=ABCMeta): return surf + def rotate(self, rotation, pivot=(0., 0., 0.), order='zyx', inplace=False): + """Rotate surface by given Tait-Bryan angles + + Parameters + ---------- + rotation : iterable of float + Intrinsic Tait-Bryan angles in degrees used to rotate the surface + pivot : iterable of float, optional + (x, y, z) coordinates for the point to rotate about. Defaults to + (0., 0., 0.) + inplace : boolean + Whether or not to return a new instance of this Surface or to modify the + coefficients of this Surface in place. Defaults to False. + + Returns + ------- + openmc.Surface + Rotated surface + + """ + + pivot = np.asarray(pivot) + surf = self.translate(-pivot, inplace=inplace) + Rmat = get_rotation_matrix(rotation, order=order) + + kwargs = {'boundary_type': self.boundary_type, 'name': self.name} + if inplace: + kwargs['sufrace_id'] = self.id + + # Handle rotations for classes with well defined origins and axes + attrdct = {} + if any((isinstance(surf, cls) for cls in (Cylinder, Sphere, Cone))): + if isinstance(surf, Sphere): + surf = surf if inplace else surf.clone() + else: + # Get generic form of Surface + surf = _VIRTUAL_BASES[surf._type](**kwargs) + # get the axis and origin of this Surface + axis = np.array([surf.dx, surf.dy, surf.dz]) + # Effectively clone all attributes that arent in axis or origin + attrdct.update({k: getattr(surf, k) for k in surf._coeff_keys}) + # update axis attributes from rotation + attrdct.update(dict(zip(('dx', 'dy', 'dz'), Rmat @ axis))) + # update origin attributes from rotation + origin = np.array([surf.x0, surf.y0, surf.z0]) + attrdct.update(dict(zip(('x0', 'y0', 'z0'), Rmat @ origin))) + + # Handle rotations for general Quadric types + elif isinstance(surf, Quadric): + surf = surf if inplace else surf.clone() + + A, bvec, k = surf.get_Abc() + Arot = Rmat @ A @ Rmat.T + + a, b, c = np.diagonal(Arot) + d, e, f = 2*Arot[0, 1], 2*Arot[1, 2], 2*Arot[0, 2] + g, h, j = Rmat @ bvec + attrdct.update(dict(zip(surf._coeff_keys, + (a, b, c, d, e, f, g, h, j, k)))) + + # Set all rotated attributes + for key, val in attrdct.items(): + setattr(surf, key, val) + + return surf.translate(pivot, inplace=inplace) + class Cylinder(QuadricMixin, Surface): """A cylinder with radius r, centered on the point (x0, y0, z0) with an @@ -2627,3 +2840,10 @@ class Halfspace(Region): return type(self)(memo[key], self.side) _SURFACE_CLASSES = {cls._type: cls for cls in Surface.__subclasses__()} +_VIRTUAL_BASES = {'{}plane'.format(s): Plane for s in ('x-', 'y-', 'z-', '')} +_VIRTUAL_BASES.update({'{}cylinder'.format(s): Cylinder + for s in ('x-', 'y-', 'z-','')}) +_VIRTUAL_BASES.update({'{}cone'.format(s): Cone + for s in ('x-', 'y-', 'z-', '')}) +_VIRTUAL_BASES['sphere'] = Sphere +_VIRTUAL_BASES['quadric'] = Quadric diff --git a/tests/unit_tests/test_surface.py b/tests/unit_tests/test_surface.py index 749ae96212..8142debd0b 100644 --- a/tests/unit_tests/test_surface.py +++ b/tests/unit_tests/test_surface.py @@ -37,6 +37,18 @@ def test_plane(): 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) + # rotate method + yp = openmc.YPlane(np.abs(s.d)/np.sqrt(s.a**2 + s.b**2 + s.c**2)) + phi = np.rad2deg(np.arctan2(1, 2)) + theta = 0. + psi = np.rad2deg(np.arctan2(1, np.sqrt(5))) + sr = s.rotate((psi, theta, phi)) + assert yp.normalize() == pytest.approx(sr.normalize()) + # test rotation ordering + phi = np.rad2deg(np.arctan2(1, np.sqrt(2))) + sr = s.rotate((0., -45., phi), order='xyz') + assert yp.normalize() == pytest.approx(sr.normalize()) + # Make sure repr works repr(s) @@ -81,6 +93,14 @@ def test_xplane(): st = s.translate((1.0, 0.0, 0.0)) assert st.x0 == s.x0 + 1 + # rotate method + # make sure rotating around x axis does nothing to coefficients + sr = s.rotate((37.4, 0., 0.)) + assert s._get_base_coeffs() == pytest.approx(sr._get_base_coeffs()) + # rotating around z by 90 deg then x by -90 deg should give negative z-plane + sr = s.rotate((-90., 0., 90)) + assert (0., 0., -1., 3.) == pytest.approx(sr._get_base_coeffs()) + # Make sure repr works repr(s) @@ -110,6 +130,17 @@ def test_yplane(): st = s.translate((0.0, 1.0, 0.0)) assert st.y0 == s.y0 + 1 + # rotate method + # make sure rotating around y axis does nothing to coefficients + sr = s.rotate((0., -12.4, 0.), order='yxz') + assert s._get_base_coeffs() == pytest.approx(sr._get_base_coeffs()) + # rotating around z by 90 deg then x by -90 deg should give negative z-plane + #sr = s.rotate((-90., 0., 90)) + #assert (0., 0., -1., 3.) == pytest.approx(sr._get_base_coeffs()) + + # Make sure repr works + repr(s) + def test_zplane(): s = openmc.ZPlane(z0=3.) From cb92d025376e5c637b05fbe6a46f5ab87a78802e Mon Sep 17 00:00:00 2001 From: Ethan Peterson Date: Tue, 25 Feb 2020 17:56:12 -0500 Subject: [PATCH 3/9] initial unit tests working --- openmc/surface.py | 212 ++++++++++--------------------- tests/unit_tests/test_surface.py | 29 +++-- 2 files changed, 90 insertions(+), 151 deletions(-) diff --git a/openmc/surface.py b/openmc/surface.py index dd2db91f27..03be4dbb00 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -28,61 +28,41 @@ will not accept positional parameters for superclass arguments.\ def get_rotation_matrix(rotation, order='xyz'): - """Generate a 3x3 rotation matrix from input angles + r"""Generate a 3x3 rotation matrix from input angles Parameters ---------- rotation : 3-tuple of float - A 3-tuple of angles (psi, theta, phi) in degrees where the first - element is the rotation about the x-axis in the fixed laboratory frame, - the second element is the rotation about the y-axis in the fixed + A 3-tuple of angles :math:`(\phi, \theta, \psi)` in degrees where the + first element is the rotation about the x-axis in the fixed laboratory + frame, the second element is the rotation about the y-axis in the fixed laboratory frame, and the third element is the rotation about the z-axis in the fixed laboratory frame. The rotations are active rotations. order : str, optinoal A string of 'x', 'y', and 'z' in some order specifying which rotation - to perform first, second, and third. Defaults to 'zyx' which means, the - rotation by angle phi about z will be applied first, followed by theta - about y and then psi about x. + to perform first, second, and third. Defaults to 'xyz' which means, the + rotation by angle phi about x will be applied first, followed by theta + about y and then phi about z. This corresponds to an x-y-z extrinsic + rotation as well as a z-y'-x'' intrinsic rotation using Tait-Bryan + angles :math:`(\phi, \theta, \psi)`. """ check_type('surface rotation', rotation, Iterable, Real) check_length('surface rotation', rotation, 3) - # Calculate rotation matrix from angles psi, theta, phi, and the order - # of the desired rotations - order= 'xyz' applies the rotation matrix - # R = Rz(phi) @ Ry(theta) @ Rx(psi) to a vector from the left - # order is the extrinsic order (actual multiplication of matrices). This - # means that the associated intrinsic rotation swaps the order of - # multiplication i.e. an explicit Rz(psi) @ Ry(theta) @ Rx(phi) is known - # as an x-y-z rotation and its associated intrinsic rotation is referred - # to as a z-y'x'' rotation. - # - # Default order produces the same active rotation as the cell.rotation - # method except that method introduces a minus sign to all the angles after - # they are read in which is bad phi, theta, psi = np.array(rotation)*(np.pi/180.) cx, sx = np.cos(phi), np.sin(phi) cy, sy = np.cos(theta), np.sin(theta) cz, sz = np.cos(psi), np.sin(psi) - Rx = np.array([[1., 0., 0.], [0., cx, -sx], [0., sx, cx]]) - Ry = np.array([[cy, 0., sy], [0., 1., 0.], [-sy, 0., cy]]) - Rz = np.array([[cz, -sz, 0.], [sz, cz, 0.], [0., 0., 1.]]) - angledct = {'x': Rx, 'y': Ry, 'z': Rz} - R1, R2, R3 = tuple(angledct[i] for i in order) + R = {} + R['x'] = np.array([[1., 0., 0.], [0., cx, -sx], [0., sx, cx]]) + R['y'] = np.array([[cy, 0., sy], [0., 1., 0.], [-sy, 0., cy]]) + R['z'] = np.array([[cz, -sz, 0.], [sz, cz, 0.], [0., 0., 1.]]) + R1, R2, R3 = tuple(R[xi] for xi in order) return R3 @ R2 @ R1 -def get_rotation_matrix_2(rotation): - phi, theta, psi = np.array(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) - return 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]]) - - class Surface(IDManagerMixin, metaclass=ABCMeta): """An implicit surface with an associated boundary condition. @@ -324,8 +304,8 @@ class Surface(IDManagerMixin, metaclass=ABCMeta): """ @abstractmethod - def rotate(self, rotation, pivot=(0., 0., 0.), order='zyx', inplace=False): - """Rotate surface by given Tait-Bryan angles + def rotate(self, rotation, pivot=(0., 0., 0.), order='xyz', inplace=False): + r"""Rotate surface by given Tait-Bryan angles Parameters ---------- @@ -334,6 +314,13 @@ class Surface(IDManagerMixin, metaclass=ABCMeta): pivot : iterable of float, optional (x, y, z) coordinates for the point to rotate about. Defaults to (0., 0., 0.) + order : str, optional + A string of 'x', 'y', and 'z' in some order specifying which + rotation to perform first, second, and third. Defaults to 'xyz' + which means, the rotation by angle phi about x will be applied + first, followed by theta about y and then phi about z. This + corresponds to an x-y-z extrinsic rotation as well as a z-y'-x'' + intrinsic rotation using Tait-Bryan angles :math:`(\phi, \theta, \psi)`. inplace : boolean Whether or not to return a new instance of Surface or to modify the coefficients of this Surface in place. Defaults to False. @@ -493,33 +480,6 @@ class PlaneMixin(metaclass=ABCMeta): ur = np.array([v if not np.isnan(v) else np.inf for v in vals]) return (ll, ur) -======= - a, b, c = self._get_base_coeffs()[0:3] - return np.array((a, b, c)) / np.sqrt(a*a + b*b + c*c) - - def bounding_box(self, side): - # Compute the bounding box based on the normal vector to the plane - nhat = self._get_normal() - lb = np.array([-np.inf, -np.inf, -np.inf]) - ub = np.array([np.inf, np.inf, np.inf]) - # If the plane is axis aligned, find the proper bounding box - if np.any(np.isclose(np.abs(nhat), 1., rtol=0., atol=self._atol)): - sign = np.dot(np.array([1., 1., 1.]), nhat) - a, b, c, d = self._get_base_coeffs() - vals = tuple(d/i if round(i) != 0 else np.nan for i in (a, b, c)) - if side == '-': - if sign > 0: - ub = np.array([v if ~np.isnan(v) else np.inf for v in vals]) - else: - lb = np.array([v if ~np.isnan(v) else -np.inf for v in vals]) - elif side == '+': - if sign > 0: - lb = np.array([v if ~np.isnan(v) else -np.inf for v in vals]) - else: - ub = np.array([v if ~np.isnan(v) else np.inf for v in vals]) - - return (lb, ub) ->>>>>>> Implementing Surface Rotations def evaluate(self, point): """Evaluate the surface equation at a given point. @@ -570,27 +530,7 @@ class PlaneMixin(metaclass=ABCMeta): return surf - def rotate(self, rotation, pivot=(0., 0., 0.), order='zyx', inplace=False): - """Rotate surface by given Tait-Bryan angles - - Parameters - ---------- - rotation : iterable of float - Intrinsic Tait-Bryan angles in degrees used to rotate the surface - pivot : iterable of float, optional - (x, y, z) coordinates for the point to rotate about. Defaults to - (0., 0., 0.) - inplace : boolean - Whether or not to return a new instance of a Plane or to modify the - coefficients of this plane in place. Defaults to False. - - Returns - ------- - openmc.Plane - Rotated surface - - """ - + def rotate(self, rotation, pivot=(0., 0., 0.), order='xyz', inplace=False): pivot = np.asarray(pivot) surf = self.translate(-pivot, inplace=inplace) @@ -601,7 +541,7 @@ class PlaneMixin(metaclass=ABCMeta): kwargs = {'boundary_type': surf.boundary_type, 'name': surf.name} if inplace: - kwargs['sufrace_id'] = surf.id + kwargs['surface_id'] = surf.id surf = Plane(a=a, b=b, c=c, d=d, **kwargs) @@ -1146,70 +1086,46 @@ class QuadricMixin(metaclass=ABCMeta): return surf - def rotate(self, rotation, pivot=(0., 0., 0.), order='zyx', inplace=False): - """Rotate surface by given Tait-Bryan angles - - Parameters - ---------- - rotation : iterable of float - Intrinsic Tait-Bryan angles in degrees used to rotate the surface - pivot : iterable of float, optional - (x, y, z) coordinates for the point to rotate about. Defaults to - (0., 0., 0.) - inplace : boolean - Whether or not to return a new instance of this Surface or to modify the - coefficients of this Surface in place. Defaults to False. - - Returns - ------- - openmc.Surface - Rotated surface - - """ - + def rotate(self, rotation, pivot=(0., 0., 0.), order='xyz', inplace=False): + # Get pivot and rotation matrix pivot = np.asarray(pivot) - surf = self.translate(-pivot, inplace=inplace) Rmat = get_rotation_matrix(rotation, order=order) - kwargs = {'boundary_type': self.boundary_type, 'name': self.name} - if inplace: - kwargs['sufrace_id'] = self.id + # Translate surface to the pivot point + tsurf = self.translate(-pivot, inplace=inplace) - # Handle rotations for classes with well defined origins and axes - attrdct = {} - if any((isinstance(surf, cls) for cls in (Cylinder, Sphere, Cone))): - if isinstance(surf, Sphere): - surf = surf if inplace else surf.clone() - else: - # Get generic form of Surface - surf = _VIRTUAL_BASES[surf._type](**kwargs) - # get the axis and origin of this Surface - axis = np.array([surf.dx, surf.dy, surf.dz]) - # Effectively clone all attributes that arent in axis or origin - attrdct.update({k: getattr(surf, k) for k in surf._coeff_keys}) - # update axis attributes from rotation - attrdct.update(dict(zip(('dx', 'dy', 'dz'), Rmat @ axis))) - # update origin attributes from rotation - origin = np.array([surf.x0, surf.y0, surf.z0]) - attrdct.update(dict(zip(('x0', 'y0', 'z0'), Rmat @ origin))) - - # Handle rotations for general Quadric types - elif isinstance(surf, Quadric): - surf = surf if inplace else surf.clone() + # If the surface is already generalized just clone it + if type(tsurf) is tsurf._virtual_base: + surf = tsurf if inplace else tsurf.clone() + else: + base_cls = type(tsurf)._virtual_base + # Copy necessary surface attributes to new kwargs dictionary + kwargs = {'boundary_type': tsurf.boundary_type, 'name': tsurf.name} + if inplace: + kwargs['surface_id'] = tsurf.id + kwargs.update({k: getattr(tsurf, k) for k in base_cls._coeff_keys}) + # Create new instance of the virtual base class + surf = base_cls(**kwargs) + # Perform rotations on axis, origin, or quadric coefficients + if hasattr(surf, 'dx'): + for key, val in zip(('dx', 'dy', 'dz'), Rmat @ tsurf._axis): + setattr(surf, key, val) + if hasattr(surf, 'x0'): + for key, val in zip(('x0', 'y0', 'z0'), Rmat @ tsurf._origin): + setattr(surf, key, val) + else: A, bvec, k = surf.get_Abc() Arot = Rmat @ A @ Rmat.T a, b, c = np.diagonal(Arot) d, e, f = 2*Arot[0, 1], 2*Arot[1, 2], 2*Arot[0, 2] g, h, j = Rmat @ bvec - attrdct.update(dict(zip(surf._coeff_keys, - (a, b, c, d, e, f, g, h, j, k)))) - # Set all rotated attributes - for key, val in attrdct.items(): - setattr(surf, key, val) + for key, val in zip(surf._coeff_keys, (a, b, c, d, e, f, g, h, j, k)): + setattr(surf, key, val) + # translate back to the original frame and return the surface return surf.translate(pivot, inplace=inplace) @@ -2840,10 +2756,20 @@ class Halfspace(Region): return type(self)(memo[key], self.side) _SURFACE_CLASSES = {cls._type: cls for cls in Surface.__subclasses__()} -_VIRTUAL_BASES = {'{}plane'.format(s): Plane for s in ('x-', 'y-', 'z-', '')} -_VIRTUAL_BASES.update({'{}cylinder'.format(s): Cylinder - for s in ('x-', 'y-', 'z-','')}) -_VIRTUAL_BASES.update({'{}cone'.format(s): Cone - for s in ('x-', 'y-', 'z-', '')}) -_VIRTUAL_BASES['sphere'] = Sphere -_VIRTUAL_BASES['quadric'] = Quadric + + +# Set virtual base classes for "casting" up the heirarchy +Plane._virtual_base = Plane +XPlane._virtual_base = Plane +YPlane._virtual_base = Plane +ZPlane._virtual_base = Plane +Cylinder._virtual_base = Cylinder +XCylinder._virtual_base = Cylinder +YCylinder._virtual_base = Cylinder +ZCylinder._virtual_base = Cylinder +Cone._virtual_base = Cone +XCone._virtual_base = Cone +YCone._virtual_base = Cone +ZCone._virtual_base = Cone +Sphere._virtual_base = Sphere +Quadric._virtual_base = Quadric diff --git a/tests/unit_tests/test_surface.py b/tests/unit_tests/test_surface.py index 8142debd0b..4b0da6f4b0 100644 --- a/tests/unit_tests/test_surface.py +++ b/tests/unit_tests/test_surface.py @@ -39,10 +39,9 @@ def test_plane(): # rotate method yp = openmc.YPlane(np.abs(s.d)/np.sqrt(s.a**2 + s.b**2 + s.c**2)) - phi = np.rad2deg(np.arctan2(1, 2)) - theta = 0. - psi = np.rad2deg(np.arctan2(1, np.sqrt(5))) - sr = s.rotate((psi, theta, phi)) + psi = np.rad2deg(np.arctan2(1, 2)) + phi = np.rad2deg(np.arctan2(1, np.sqrt(5))) + sr = s.rotate((phi, 0., psi), order='zyx') assert yp.normalize() == pytest.approx(sr.normalize()) # test rotation ordering phi = np.rad2deg(np.arctan2(1, np.sqrt(2))) @@ -98,7 +97,7 @@ def test_xplane(): sr = s.rotate((37.4, 0., 0.)) assert s._get_base_coeffs() == pytest.approx(sr._get_base_coeffs()) # rotating around z by 90 deg then x by -90 deg should give negative z-plane - sr = s.rotate((-90., 0., 90)) + sr = s.rotate((-90., 0., 90), order='zyx') assert (0., 0., -1., 3.) == pytest.approx(sr._get_base_coeffs()) # Make sure repr works @@ -134,9 +133,9 @@ def test_yplane(): # make sure rotating around y axis does nothing to coefficients sr = s.rotate((0., -12.4, 0.), order='yxz') assert s._get_base_coeffs() == pytest.approx(sr._get_base_coeffs()) - # rotating around z by 90 deg then x by -90 deg should give negative z-plane - #sr = s.rotate((-90., 0., 90)) - #assert (0., 0., -1., 3.) == pytest.approx(sr._get_base_coeffs()) + # rotate around x by -90 deg and y by 90 deg should give negative x-plane + sr = s.rotate((-90, 90, 0.)) + assert (-1, 0., 0., 3.) == pytest.approx(sr._get_base_coeffs()) # Make sure repr works repr(s) @@ -167,6 +166,14 @@ def test_zplane(): st = s.translate((0.0, 0.0, 1.0)) assert st.z0 == s.z0 + 1 + # rotate method + # make sure rotating around z axis does nothing to coefficients + sr = s.rotate((0., 0., 123), order='zxy') + assert s._get_base_coeffs() == pytest.approx(sr._get_base_coeffs()) + # rotate around x by -90 deg and y by 90 deg should give negative x-plane + sr = s.rotate((-90, 0., 90.)) + assert (-1., 0., 0., 3.) == pytest.approx(sr._get_base_coeffs()) + # Make sure repr works repr(s) @@ -212,6 +219,12 @@ def test_cylinder(): assert st.dz == s.dz assert st.r == s.r + # rotate method + sr = s.rotate((90, 90, 90)) + R = np.array([[0, 0, 1], [0, 1, 0], [-1, 0, 0]]) + assert sr._origin == pytest.approx(R @ s._origin) + assert sr._axis == pytest.approx(R @ s._axis) + # Make sure repr works repr(s) From 7d69ff6de18476ef7194abc6264c967f4d22f389 Mon Sep 17 00:00:00 2001 From: Ethan Peterson Date: Wed, 26 Feb 2020 12:33:51 -0500 Subject: [PATCH 4/9] applied suggestions from code review --- openmc/surface.py | 48 ++++++++++++++++++++++---------- tests/unit_tests/test_surface.py | 8 +++--- 2 files changed, 38 insertions(+), 18 deletions(-) diff --git a/openmc/surface.py b/openmc/surface.py index 03be4dbb00..562b4c859d 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -51,15 +51,17 @@ def get_rotation_matrix(rotation, order='xyz'): check_type('surface rotation', rotation, Iterable, Real) check_length('surface rotation', rotation, 3) - phi, theta, psi = np.array(rotation)*(np.pi/180.) - cx, sx = np.cos(phi), np.sin(phi) - cy, sy = np.cos(theta), np.sin(theta) - cz, sz = np.cos(psi), np.sin(psi) - R = {} - R['x'] = np.array([[1., 0., 0.], [0., cx, -sx], [0., sx, cx]]) - R['y'] = np.array([[cy, 0., sy], [0., 1., 0.], [-sy, 0., cy]]) - R['z'] = np.array([[cz, -sz, 0.], [sz, cz, 0.], [0., 0., 1.]]) - R1, R2, R3 = tuple(R[xi] for xi in order) + phi, theta, psi = np.array(rotation)*(math.pi/180.) + cx, sx = math.cos(phi), math.sin(phi) + cy, sy = math.cos(theta), math.sin(theta) + cz, sz = math.cos(psi), math.sin(psi) + R = { + 'x': np.array([[1., 0., 0.], [0., cx, -sx], [0., sx, cx]]), + 'y': np.array([[cy, 0., sy], [0., 1., 0.], [-sy, 0., cy]]), + 'z': np.array([[cz, -sz, 0.], [sz, cz, 0.], [0., 0., 1.]]), + } + + R1, R2, R3 = (R[xi] for xi in order) return R3 @ R2 @ R1 @@ -522,7 +524,7 @@ class PlaneMixin(metaclass=ABCMeta): return self a, b, c, d = self._get_base_coeffs() - d = d + np.array([a, b, c]) @ np.array(vector) + d = d + np.dot([a, b, c], vector) surf = self if inplace else self.clone() @@ -532,12 +534,22 @@ class PlaneMixin(metaclass=ABCMeta): def rotate(self, rotation, pivot=(0., 0., 0.), order='xyz', inplace=False): pivot = np.asarray(pivot) + rotation = np.asarray(rotation) + + # Allow rotaiton matrix to be passed in directly, otherwise build it + if np.rank(rotation) == 2: + check_type('surface rotation', rotation, Iterable, Real) + check_length('surface rotation', rotation.ravel(), 9) + Rmat = rotation + else: + Rmat = get_rotation_matrix(rotation, order=order) + + # Translate surface to pivot surf = self.translate(-pivot, inplace=inplace) - Rmat = get_rotation_matrix(rotation, order=order) a, b, c, d = surf._get_base_coeffs() # Compute new rotated coefficients a, b, c - a, b, c = Rmat @ np.array([a, b, c]) + a, b, c = Rmat @ [a, b, c] kwargs = {'boundary_type': surf.boundary_type, 'name': surf.name} if inplace: @@ -1089,14 +1101,22 @@ class QuadricMixin(metaclass=ABCMeta): def rotate(self, rotation, pivot=(0., 0., 0.), order='xyz', inplace=False): # Get pivot and rotation matrix pivot = np.asarray(pivot) - Rmat = get_rotation_matrix(rotation, order=order) + rotation = np.asarray(rotation) + + # Allow rotaiton matrix to be passed in directly, otherwise build it + if np.rank(rotation) == 2: + check_type('surface rotation', rotation, Iterable, Real) + check_length('surface rotation', rotation.ravel(), 9) + Rmat = rotation + else: + Rmat = get_rotation_matrix(rotation, order=order) # Translate surface to the pivot point tsurf = self.translate(-pivot, inplace=inplace) # If the surface is already generalized just clone it if type(tsurf) is tsurf._virtual_base: - surf = tsurf if inplace else tsurf.clone() + surf = tsurf if inplace else tsurf.clone() else: base_cls = type(tsurf)._virtual_base # Copy necessary surface attributes to new kwargs dictionary diff --git a/tests/unit_tests/test_surface.py b/tests/unit_tests/test_surface.py index 4b0da6f4b0..da764e3bbc 100644 --- a/tests/unit_tests/test_surface.py +++ b/tests/unit_tests/test_surface.py @@ -38,13 +38,13 @@ def test_plane(): assert (st.a, st.b, st.c, st.d) == (s.a, s.b, s.c, 4) # rotate method - yp = openmc.YPlane(np.abs(s.d)/np.sqrt(s.a**2 + s.b**2 + s.c**2)) - psi = np.rad2deg(np.arctan2(1, 2)) - phi = np.rad2deg(np.arctan2(1, np.sqrt(5))) + yp = openmc.YPlane(abs(s.d)/math.sqrt(s.a**2 + s.b**2 + s.c**2)) + psi = math.degrees(np.arctan2(1, 2)) + phi = math.degrees(np.arctan2(1, np.sqrt(5))) sr = s.rotate((phi, 0., psi), order='zyx') assert yp.normalize() == pytest.approx(sr.normalize()) # test rotation ordering - phi = np.rad2deg(np.arctan2(1, np.sqrt(2))) + phi = math.degrees(np.arctan2(1, np.sqrt(2))) sr = s.rotate((0., -45., phi), order='xyz') assert yp.normalize() == pytest.approx(sr.normalize()) From a25882e76f8ae39b577c29767e1d454412815e55 Mon Sep 17 00:00:00 2001 From: Ethan Peterson Date: Wed, 26 Feb 2020 17:42:01 -0500 Subject: [PATCH 5/9] more unit tests --- openmc/surface.py | 4 +-- tests/unit_tests/test_surface.py | 61 ++++++++++++++++++++++++++++++++ 2 files changed, 63 insertions(+), 2 deletions(-) diff --git a/openmc/surface.py b/openmc/surface.py index 562b4c859d..eedbd13170 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -537,7 +537,7 @@ class PlaneMixin(metaclass=ABCMeta): rotation = np.asarray(rotation) # Allow rotaiton matrix to be passed in directly, otherwise build it - if np.rank(rotation) == 2: + if rotation.ndim == 2: check_type('surface rotation', rotation, Iterable, Real) check_length('surface rotation', rotation.ravel(), 9) Rmat = rotation @@ -1104,7 +1104,7 @@ class QuadricMixin(metaclass=ABCMeta): rotation = np.asarray(rotation) # Allow rotaiton matrix to be passed in directly, otherwise build it - if np.rank(rotation) == 2: + if rotation.ndim == 2: check_type('surface rotation', rotation, Iterable, Real) check_length('surface rotation', rotation.ravel(), 9) Rmat = rotation diff --git a/tests/unit_tests/test_surface.py b/tests/unit_tests/test_surface.py index da764e3bbc..77c8cbe5a5 100644 --- a/tests/unit_tests/test_surface.py +++ b/tests/unit_tests/test_surface.py @@ -224,6 +224,9 @@ def test_cylinder(): R = np.array([[0, 0, 1], [0, 1, 0], [-1, 0, 0]]) assert sr._origin == pytest.approx(R @ s._origin) assert sr._axis == pytest.approx(R @ s._axis) + # test passing in rotation matrix + sr2 = s.rotate(R) + assert sr2.is_equal(sr) # Make sure repr works repr(s) @@ -253,6 +256,15 @@ def test_xcylinder(): assert st.z0 == s.z0 + 1 assert st.r == s.r + # rotate method + sr = s.rotate((90, 90, 90)) + R = np.array([[0, 0, 1], [0, 1, 0], [-1, 0, 0]]) + assert sr._origin == pytest.approx(R @ s._origin) + assert sr._axis == pytest.approx(R @ s._axis) + # test passing in rotation matrix + sr2 = s.rotate(R) + assert sr2._get_base_coeffs() == pytest.approx(sr._get_base_coeffs()) + # Make sure repr works repr(s) @@ -290,6 +302,18 @@ def test_ycylinder(): assert st.z0 == s.z0 + 1 assert st.r == s.r + # rotate method + sr = s.rotate((90, 90, 90)) + R = np.array([[0, 0, 1], [0, 1, 0], [-1, 0, 0]]) + assert sr._origin == pytest.approx(R @ s._origin) + assert sr._axis == pytest.approx(R @ s._axis) + # test passing in rotation matrix + sr2 = s.rotate(R) + assert sr2._get_base_coeffs() == pytest.approx(sr._get_base_coeffs()) + + # Make sure repr works + repr(s) + def test_zcylinder(): x, y, r = 3, 5, 2 @@ -315,6 +339,15 @@ def test_zcylinder(): assert st.y0 == s.y0 + 1 assert st.r == s.r + # rotate method + sr = s.rotate((90, 90, 90)) + R = np.array([[0, 0, 1], [0, 1, 0], [-1, 0, 0]]) + assert sr._origin == pytest.approx(R @ s._origin) + assert sr._axis == pytest.approx(R @ s._axis) + # test passing in rotation matrix + sr2 = s.rotate(R) + assert sr2._get_base_coeffs() == pytest.approx(sr._get_base_coeffs()) + # Make sure repr works repr(s) @@ -345,6 +378,15 @@ def test_sphere(): assert st.z0 == s.z0 + 1 assert st.r == s.r + # rotate method + pivot = np.array([1, -2, 3]) + sr = s.rotate((90, 90, 90), pivot=pivot) + R = np.array([[0, 0, 1], [0, 1, 0], [-1, 0, 0]]) + assert sr._origin == pytest.approx((R @ (s._origin - pivot)) + pivot) + # test passing in rotation matrix + sr2 = s.rotate(R, pivot=pivot) + assert sr2._get_base_coeffs() == pytest.approx(sr._get_base_coeffs()) + # Make sure repr works repr(s) @@ -370,6 +412,15 @@ def cone_common(apex, r2, cls): assert st.z0 == s.z0 + 1 assert st.r2 == s.r2 + # rotate method + sr = s.rotate((90, 90, 90)) + R = np.array([[0, 0, 1], [0, 1, 0], [-1, 0, 0]]) + assert sr._origin == pytest.approx(R @ s._origin) + assert sr._axis == pytest.approx(R @ s._axis) + # test passing in rotation matrix + sr2 = s.rotate(R) + assert sr2._get_base_coeffs() == pytest.approx(sr._get_base_coeffs()) + # Make sure repr works repr(s) @@ -422,6 +473,15 @@ def test_cone(): assert st.dz == s.dz assert st.r2 == s.r2 + # rotate method + sr = s.rotate((90, 90, 90)) + R = np.array([[0, 0, 1], [0, 1, 0], [-1, 0, 0]]) + assert sr._origin == pytest.approx(R @ s._origin) + assert sr._axis == pytest.approx(R @ s._axis) + # test passing in rotation matrix + sr2 = s.rotate(R) + assert sr2._get_base_coeffs() == pytest.approx(sr._get_base_coeffs()) + # Make sure repr works repr(s) @@ -453,6 +513,7 @@ def test_quadric(): assert s.b == coeffs['b'] assert s.c == coeffs['c'] assert s.k == coeffs['k'] + assert openmc.Sphere(r=10).is_equal(s) # All other coeffs should be zero for coeff in ('d', 'e', 'f', 'g', 'h', 'j'): From 98f914b4c26405bc6f9204917f051d08921ebb32 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 27 Feb 2020 09:38:15 -0600 Subject: [PATCH 6/9] Avoid target_link_libraries on imported target for CMake <3.12 compatibility --- CMakeLists.txt | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index f251d75751..b761929523 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -42,8 +42,9 @@ if(dagmc) find_package(DAGMC REQUIRED) add_library(dagmc-imported INTERFACE IMPORTED) link_directories(${DAGMC_LIBRARY_DIRS}) - target_link_libraries(dagmc-imported INTERFACE ${DAGMC_LIBRARIES}) - target_include_directories(dagmc-imported INTERFACE ${DAGMC_INCLUDE_DIRS}) + set_target_properties(dagmc-imported PROPERTIES + INTERFACE_LINK_LIBRARIES "${DAGMC_LIBRARIES}" + INTERFACE_INCLUDE_DIRECTORIES "${DAGMC_INCLUDE_DIRS}") else() set(DAGMC_FOUND false) endif() From c88bd495976aefed998680046221df4645224e36 Mon Sep 17 00:00:00 2001 From: Ethan Peterson Date: Thu, 27 Feb 2020 12:31:15 -0500 Subject: [PATCH 7/9] finished unit tests --- openmc/region.py | 186 ++++++++++--------------------- openmc/surface.py | 96 +++++++++++++--- tests/unit_tests/test_region.py | 22 ++++ tests/unit_tests/test_surface.py | 9 ++ 4 files changed, 174 insertions(+), 139 deletions(-) diff --git a/openmc/region.py b/openmc/region.py index 2b2fce14d0..22718ac7b2 100644 --- a/openmc/region.py +++ b/openmc/region.py @@ -225,7 +225,6 @@ class Region(metaclass=ABCMeta): # at the end return output[0] - @abstractmethod def clone(self, memo=None): """Create a copy of this region - each of the surfaces in the region's nodes will be cloned and will have new unique IDs. @@ -242,9 +241,14 @@ class Region(metaclass=ABCMeta): The clone of this region """ - pass - @abstractmethod + if memo is None: + memo = {} + + clone = deepcopy(self) + clone[:] = [n.clone(memo) for n in self] + return clone + def translate(self, vector, memo=None): """Translate region in given direction @@ -262,7 +266,52 @@ class Region(metaclass=ABCMeta): Translated region """ - pass + + if memo is None: + memo = {} + return type(self)(n.translate(vector, memo) for n in self) + + def rotate(self, rotation, pivot=(0., 0., 0.), order='xyz', inplace=False, + memo=None): + r"""Rotate surface by angles provided or by applying matrix directly. + + Parameters + ---------- + rotation : 3-tuple of float, or 3x3 iterable + A 3-tuple of angles :math:`(\phi, \theta, \psi)` in degrees where + the first element is the rotation about the x-axis in the fixed + laboratory frame, the second element is the rotation about the + y-axis in the fixed laboratory frame, and the third element is the + rotation about the z-axis in the fixed laboratory frame. The + rotations are active rotations. Additionally a 3x3 rotation matrix + can be specified directly either as a nested iterable or array. + pivot : iterable of float, optional + (x, y, z) coordinates for the point to rotate about. Defaults to + (0., 0., 0.) + order : str, optional + A string of 'x', 'y', and 'z' in some order specifying which + rotation to perform first, second, and third. Defaults to 'xyz' + which means, the rotation by angle :math:`\phi` about x will be + applied first, followed by :math:`\theta` about y and then + :math:`\psi` about z. This corresponds to an x-y-z extrinsic + rotation as well as a z-y'-x'' intrinsic rotation using Tait-Bryan + angles :math:`(\phi, \theta, \psi)`. + inplace : boolean + Whether or not to return a new instance of Surface or to modify the + coefficients of this Surface in place. Defaults to False. + memo : dict or None + Dictionary used for memoization + + Returns + ------- + openmc.Region + Translated region + + """ + if memo is None: + memo = {} + return type(self)(n.rotate(rotation, pivot=pivot, order=order, + inplace=inplace, memo=memo) for n in self) class Intersection(Region, MutableSequence): @@ -354,51 +403,6 @@ class Intersection(Region, MutableSequence): upper_right[:] = np.minimum(upper_right, upper_right_n) return lower_left, upper_right - def clone(self, memo=None): - """Create a copy of this region - each of the surfaces in the - intersection's nodes will be cloned and will have new unique IDs. - - Parameters - ---------- - memo : dict or None - A nested dictionary of previously cloned objects. This parameter - is used internally and should not be specified by the user. - - Returns - ------- - clone : openmc.Intersection - The clone of this intersection - - """ - - if memo is None: - memo = {} - - clone = deepcopy(self) - 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. This parameter is used internally - and should not be specified by the user. - - 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. @@ -487,51 +491,6 @@ class Union(Region, MutableSequence): upper_right[:] = np.maximum(upper_right, upper_right_n) return lower_left, upper_right - def clone(self, memo=None): - """Create a copy of this region - each of the surfaces in the - union's nodes will be cloned and will have new unique IDs. - - Parameters - ---------- - memo : dict or None - A nested dictionary of previously cloned objects. This parameter - is used internally and should not be specified by the user. - - Returns - ------- - clone : openmc.Union - The clone of this union - - """ - - if memo is None: - memo = {} - - clone = deepcopy(self) - 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. This parameter is used internally - and should not be specified by the user. - - 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. @@ -642,22 +601,6 @@ class Complement(Region): region.remove_redundant_surfaces(redundant_surfaces) def clone(self, memo=None): - """Create a copy of this region - each of the surfaces in the - complement's node will be cloned and will have new unique IDs. - - Parameters - ---------- - memo : dict or None - A nested dictionary of previously cloned objects. This parameter - is used internally and should not be specified by the user. - - Returns - ------- - clone : openmc.Complement - The clone of this complement - - """ - if memo is None: memo = {} @@ -666,22 +609,13 @@ class Complement(Region): 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. This parameter is used internally - and should not be specified by the user. - - Returns - ------- - openmc.Complement - Translated region - - """ if memo is None: memo = {} return type(self)(self.node.translate(vector, memo)) + + def rotate(self, rotation, pivot=(0., 0., 0.), order='xyz', inplace=False, + memo=None): + if memo is None: + memo = {} + return type(self)(self.node.rotate(rotation, pivot=pivot, order=order, + inplace=inplace, memo=memo)) diff --git a/openmc/surface.py b/openmc/surface.py index eedbd13170..f0198e02bc 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -39,13 +39,13 @@ def get_rotation_matrix(rotation, order='xyz'): laboratory frame, and the third element is the rotation about the z-axis in the fixed laboratory frame. The rotations are active rotations. - order : str, optinoal + order : str, optional A string of 'x', 'y', and 'z' in some order specifying which rotation to perform first, second, and third. Defaults to 'xyz' which means, the - rotation by angle phi about x will be applied first, followed by theta - about y and then phi about z. This corresponds to an x-y-z extrinsic - rotation as well as a z-y'-x'' intrinsic rotation using Tait-Bryan - angles :math:`(\phi, \theta, \psi)`. + rotation by angle :math:`\phi` about x will be applied first, followed + by :math:`\theta` about y and then :math:`\psi` about z. This + corresponds to an x-y-z extrinsic rotation as well as a z-y'-x'' + intrinsic rotation using Tait-Bryan angles :math:`(\phi, \theta, \psi)`. """ check_type('surface rotation', rotation, Iterable, Real) @@ -307,22 +307,29 @@ class Surface(IDManagerMixin, metaclass=ABCMeta): @abstractmethod def rotate(self, rotation, pivot=(0., 0., 0.), order='xyz', inplace=False): - r"""Rotate surface by given Tait-Bryan angles + r"""Rotate surface by angles provided or by applying matrix directly. Parameters ---------- - rotation : iterable of float - Intrinsic Tait-Bryan angles in degrees used to rotate the surface + rotation : 3-tuple of float, or 3x3 iterable + A 3-tuple of angles :math:`(\phi, \theta, \psi)` in degrees where + the first element is the rotation about the x-axis in the fixed + laboratory frame, the second element is the rotation about the + y-axis in the fixed laboratory frame, and the third element is the + rotation about the z-axis in the fixed laboratory frame. The + rotations are active rotations. Additionally a 3x3 rotation matrix + can be specified directly either as a nested iterable or array. pivot : iterable of float, optional (x, y, z) coordinates for the point to rotate about. Defaults to (0., 0., 0.) order : str, optional A string of 'x', 'y', and 'z' in some order specifying which rotation to perform first, second, and third. Defaults to 'xyz' - which means, the rotation by angle phi about x will be applied - first, followed by theta about y and then phi about z. This - corresponds to an x-y-z extrinsic rotation as well as a z-y'-x'' - intrinsic rotation using Tait-Bryan angles :math:`(\phi, \theta, \psi)`. + which means, the rotation by angle :math:`\phi` about x will be + applied first, followed by :math:`\theta` about y and then + :math:`\psi` about z. This corresponds to an x-y-z extrinsic + rotation as well as a z-y'-x'' intrinsic rotation using Tait-Bryan + angles :math:`(\phi, \theta, \psi)`. inplace : boolean Whether or not to return a new instance of Surface or to modify the coefficients of this Surface in place. Defaults to False. @@ -1291,6 +1298,19 @@ class Cylinder(QuadricMixin, Surface): check_type('dz coefficient', dz, Real) self._coefficients['dz'] = dz + def bounding_box(self, side): + if side == '-': + r = self.r + ll = [xi - r if np.isclose(dxi, 0., rtol=0., atol=self._atol) + else -np.inf for xi, dxi in zip(self._origin, self._axis)] + ur = [xi + r if np.isclose(dxi, 0., rtol=0., atol=self._atol) + else np.inf for xi, dxi in zip(self._origin, self._axis)] + return (np.array(ll), np.array(ur)) + + elif side == '+': + return (np.array([-np.inf, -np.inf, -np.inf]), + np.array([np.inf, np.inf, np.inf])) + def _get_base_coeffs(self): # Get x, y, z coordinates of two points x1, y1, z1 = self._origin @@ -2772,9 +2792,59 @@ class Halfspace(Region): if key not in memo: memo[key] = self.surface.translate(vector) - # Return translated surface + # Return translated half-space return type(self)(memo[key], self.side) + def rotate(self, rotation, pivot=(0., 0., 0.), order='xyz', inplace=False, + memo=None): + r"""Rotate surface by angles provided or by applying matrix directly. + + Parameters + ---------- + rotation : 3-tuple of float, or 3x3 iterable + A 3-tuple of angles :math:`(\phi, \theta, \psi)` in degrees where + the first element is the rotation about the x-axis in the fixed + laboratory frame, the second element is the rotation about the + y-axis in the fixed laboratory frame, and the third element is the + rotation about the z-axis in the fixed laboratory frame. The + rotations are active rotations. Additionally a 3x3 rotation matrix + can be specified directly either as a nested iterable or array. + pivot : iterable of float, optional + (x, y, z) coordinates for the point to rotate about. Defaults to + (0., 0., 0.) + order : str, optional + A string of 'x', 'y', and 'z' in some order specifying which + rotation to perform first, second, and third. Defaults to 'xyz' + which means, the rotation by angle :math:`\phi` about x will be + applied first, followed by :math:`\theta` about y and then + :math:`\psi` about z. This corresponds to an x-y-z extrinsic + rotation as well as a z-y'-x'' intrinsic rotation using Tait-Bryan + angles :math:`(\phi, \theta, \psi)`. + inplace : boolean + Whether or not to return a new instance of Surface or to modify the + coefficients of this Surface in place. Defaults to False. + memo : dict or None + Dictionary used for memoization + + Returns + ------- + openmc.Halfspace + Translated half-space + + """ + if memo is None: + memo = {} + + # If rotated surface not in memo, add it + key = (self.surface, tuple(rotation), tuple(pivot), order, inplace) + if key not in memo: + memo[key] = self.surface.rotate(rotation, pivot=pivot, order=order, + inplace=inplace) + + # Return rotated half-space + return type(self)(memo[key], self.side) + + _SURFACE_CLASSES = {cls._type: cls for cls in Surface.__subclasses__()} diff --git a/tests/unit_tests/test_region.py b/tests/unit_tests/test_region.py index f9656ebb4b..8b2e06ad1b 100644 --- a/tests/unit_tests/test_region.py +++ b/tests/unit_tests/test_region.py @@ -40,6 +40,13 @@ def test_union(reset): assert (6, 0, 0) not in regt assert (8, 0, 0) in regt + # rotate method + regr = region.rotate((0., 90., 0.)) + assert (-4, 0, 0) not in regr + assert (0, 0, 6) in regr + assert (0, 0, -6) in regr + assert (0, 0, 3) not in regr + def test_intersection(reset): s1 = openmc.XPlane(x0=5, surface_id=1) @@ -73,6 +80,13 @@ def test_intersection(reset): assert (6, 0, 0) in regt assert (8, 0, 0) not in regt + # rotate method + regr = region.rotate((0., 90., 0.)) + assert (-4, 0, 0) in regr + assert (0, 0, 6) not in regr + assert (0, 0, -6) not in regr + assert (0, 0, 3) in regr + def test_complement(reset): zcyl = openmc.ZCylinder(r=1., surface_id=1) @@ -106,6 +120,14 @@ def test_complement(reset): assert ll == pytest.approx((0., 0., -4.)) assert ur == pytest.approx((2., 2., 6.)) + # rotate method + inside_r = inside.rotate((90., 0., 0.)) + ll, ur = inside_r.bounding_box + assert (.5, 2, 0) in inside_r + assert (0, 0, 6) not in inside_r + assert ll == pytest.approx((-1., -5., -1.)) + assert ur == pytest.approx((1., 5., 1.)) + def test_get_surfaces(): s1 = openmc.XPlane() diff --git a/tests/unit_tests/test_surface.py b/tests/unit_tests/test_surface.py index 77c8cbe5a5..b5e4015494 100644 --- a/tests/unit_tests/test_surface.py +++ b/tests/unit_tests/test_surface.py @@ -533,6 +533,15 @@ def test_quadric(): assert (st.g, st.h, st.j) == (-2, -2, -2) assert st.k == s.k + 3 + # rotate method + x0, y0, z0, r2 = 2, 3, 4, 4 + dx, dy, dz = 1, -1, 1 + s = openmc.Cone(x0=x0, y0=y0, z0=z0, dx=dx, dy=dy, dz=dz, r2=r2) + q = openmc.Quadric(*s._get_base_coeffs()) + qr = q.rotate((45, 60, 30)) + sr = s.rotate((45, 60, 30)) + assert qr.is_equal(sr) + def test_cylinder_from_points(): seed(1) # Make random numbers reproducible From c29ffa947dd08ab4410c1927729b5d3d1f5195f5 Mon Sep 17 00:00:00 2001 From: Ethan Peterson Date: Fri, 28 Feb 2020 08:03:36 -0500 Subject: [PATCH 8/9] Apply suggestions from code review Co-Authored-By: Paul Romano --- openmc/surface.py | 6 +++--- tests/unit_tests/test_surface.py | 6 +++--- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/openmc/surface.py b/openmc/surface.py index f0198e02bc..0909b02a87 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -541,9 +541,9 @@ class PlaneMixin(metaclass=ABCMeta): def rotate(self, rotation, pivot=(0., 0., 0.), order='xyz', inplace=False): pivot = np.asarray(pivot) - rotation = np.asarray(rotation) + rotation = np.asarray(rotation, dtype=float) - # Allow rotaiton matrix to be passed in directly, otherwise build it + # Allow rotation matrix to be passed in directly, otherwise build it if rotation.ndim == 2: check_type('surface rotation', rotation, Iterable, Real) check_length('surface rotation', rotation.ravel(), 9) @@ -1108,7 +1108,7 @@ class QuadricMixin(metaclass=ABCMeta): def rotate(self, rotation, pivot=(0., 0., 0.), order='xyz', inplace=False): # Get pivot and rotation matrix pivot = np.asarray(pivot) - rotation = np.asarray(rotation) + rotation = np.asarray(rotation, dtype=float) # Allow rotaiton matrix to be passed in directly, otherwise build it if rotation.ndim == 2: diff --git a/tests/unit_tests/test_surface.py b/tests/unit_tests/test_surface.py index b5e4015494..b37dc8f310 100644 --- a/tests/unit_tests/test_surface.py +++ b/tests/unit_tests/test_surface.py @@ -39,12 +39,12 @@ def test_plane(): # rotate method yp = openmc.YPlane(abs(s.d)/math.sqrt(s.a**2 + s.b**2 + s.c**2)) - psi = math.degrees(np.arctan2(1, 2)) - phi = math.degrees(np.arctan2(1, np.sqrt(5))) + psi = math.degrees(math.atan2(1, 2)) + phi = math.degrees(math.atan2(1, math.sqrt(5))) sr = s.rotate((phi, 0., psi), order='zyx') assert yp.normalize() == pytest.approx(sr.normalize()) # test rotation ordering - phi = math.degrees(np.arctan2(1, np.sqrt(2))) + phi = math.degrees(math.atan2(1, math.sqrt(2))) sr = s.rotate((0., -45., phi), order='xyz') assert yp.normalize() == pytest.approx(sr.normalize()) From 248ac357ac89da44f397ede5136ea83f43bfe5cf Mon Sep 17 00:00:00 2001 From: Ethan Peterson Date: Fri, 28 Feb 2020 09:21:32 -0500 Subject: [PATCH 9/9] fixed future kwarg warning issue --- openmc/surface.py | 91 +++++++++++++++++++++++------------------------ 1 file changed, 45 insertions(+), 46 deletions(-) diff --git a/openmc/surface.py b/openmc/surface.py index 0909b02a87..021f134b20 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -27,6 +27,15 @@ will not accept positional parameters for superclass arguments.\ """ +def _future_kwargs_warning_helper(cls, *args, **kwargs): + # Warn if Surface parameters are passed by position, not by keyword + argsdict = dict(zip(('boundary_type', 'name', 'surface_id'), args)) + for k in argsdict: + warn(_WARNING_KWARGS.format(cls.__name__, k), FutureWarning) + kwargs.update(argsdict) + return kwargs + + def get_rotation_matrix(rotation, order='xyz'): r"""Generate a 3x3 rotation matrix from input angles @@ -545,7 +554,6 @@ class PlaneMixin(metaclass=ABCMeta): # Allow rotation matrix to be passed in directly, otherwise build it if rotation.ndim == 2: - check_type('surface rotation', rotation, Iterable, Real) check_length('surface rotation', rotation.ravel(), 9) Rmat = rotation else: @@ -642,30 +650,23 @@ class Plane(PlaneMixin, Surface): # *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) - + kwargs = _future_kwargs_warning_helper(type(self), *args, **kwargs) # Warn if capital letter arguments are passed + capdict = {} for k in 'ABCD': val = kwargs.pop(k, None) if val is not None: warn(_WARNING_UPPER.format(type(self), k.lower(), k), FutureWarning) - setattr(self, k.lower(), val) + capdict[k.lower()] = val + + super().__init__(**kwargs) + + for key, val in zip(self._coeff_keys, (a, b, c, d)): + setattr(self, key, val) + + for key, val in capdict.items(): + setattr(self, key, val) @classmethod def __subclasshook__(cls, c): @@ -787,11 +788,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 = dict(zip(('boundary_type', 'name', 'surface_id'), args)) - for k in argsdict: - warn(_WARNING_KWARGS.format(type(self).__name__, k), FutureWarning) - kwargs.update(argsdict) - + kwargs = _future_kwargs_warning_helper(type(self), *args, **kwargs) super().__init__(**kwargs) self.x0 = x0 @@ -869,11 +866,7 @@ 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 = dict(zip(('boundary_type', 'name', 'surface_id'), args)) - for k in argsdict: - warn(_WARNING_KWARGS.format(type(self).__name__, k), FutureWarning) - kwargs.update(argsdict) - + kwargs = _future_kwargs_warning_helper(type(self), *args, **kwargs) super().__init__(**kwargs) self.y0 = y0 @@ -951,11 +944,7 @@ 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 = dict(zip(('boundary_type', 'name', 'surface_id'), args)) - for k in argsdict: - warn(_WARNING_KWARGS.format(type(self).__name__, k), FutureWarning) - kwargs.update(argsdict) - + kwargs = _future_kwargs_warning_helper(type(self), *args, **kwargs) super().__init__(**kwargs) self.z0 = z0 @@ -1112,7 +1101,6 @@ class QuadricMixin(metaclass=ABCMeta): # Allow rotaiton matrix to be passed in directly, otherwise build it if rotation.ndim == 2: - check_type('surface rotation', rotation, Iterable, Real) check_length('surface rotation', rotation.ravel(), 9) Rmat = rotation else: @@ -1223,7 +1211,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): + def __init__(self, x0=0., y0=0., z0=0., r=1., dx=0., dy=0., dz=1., *args, + **kwargs): + kwargs = _future_kwargs_warning_helper(type(self), *args, **kwargs) super().__init__(**kwargs) for key, val in zip(self._coeff_keys, (x0, y0, z0, r, dx, dy, dz)): @@ -1438,12 +1428,13 @@ class XCylinder(QuadricMixin, Surface): _type = 'x-cylinder' _coeff_keys = ('y0', 'z0', 'r') - def __init__(self, y0=0., z0=0., r=1., **kwargs): + def __init__(self, y0=0., z0=0., r=1., *args, **kwargs): R = kwargs.pop('R', None) if R is not None: warn(_WARNING_UPPER.format(type(self).__name__, 'r', 'R'), FutureWarning) r = R + kwargs = _future_kwargs_warning_helper(type(self), *args, **kwargs) super().__init__(**kwargs) for key, val in zip(self._coeff_keys, (y0, z0, r)): @@ -1563,12 +1554,13 @@ class YCylinder(QuadricMixin, Surface): _type = 'y-cylinder' _coeff_keys = ('x0', 'z0', 'r') - def __init__(self, x0=0., z0=0., r=1., **kwargs): + def __init__(self, x0=0., z0=0., r=1., *args, **kwargs): R = kwargs.pop('R', None) if R is not None: warn(_WARNING_UPPER.format(type(self).__name__, 'r', 'R'), FutureWarning) r = R + kwargs = _future_kwargs_warning_helper(type(self), *args, **kwargs) super().__init__(**kwargs) for key, val in zip(self._coeff_keys, (x0, z0, r)): @@ -1688,12 +1680,13 @@ class ZCylinder(QuadricMixin, Surface): _type = 'z-cylinder' _coeff_keys = ('x0', 'y0', 'r') - def __init__(self, x0=0., y0=0., r=1., **kwargs): + def __init__(self, x0=0., y0=0., r=1., *args, **kwargs): R = kwargs.pop('R', None) if R is not None: warn(_WARNING_UPPER.format(type(self).__name__, 'r', 'R'), FutureWarning) r = R + kwargs = _future_kwargs_warning_helper(type(self), *args, **kwargs) super().__init__(**kwargs) for key, val in zip(self._coeff_keys, (x0, y0, r)): @@ -1815,12 +1808,13 @@ class Sphere(QuadricMixin, Surface): _type = 'sphere' _coeff_keys = ('x0', 'y0', 'z0', 'r') - def __init__(self, x0=0., y0=0., z0=0., r=1., **kwargs): + def __init__(self, x0=0., y0=0., z0=0., r=1., *args, **kwargs): R = kwargs.pop('R', None) if R is not None: warn(_WARNING_UPPER.format(type(self).__name__, 'r', 'R'), FutureWarning) r = R + kwargs = _future_kwargs_warning_helper(type(self), *args, **kwargs) super().__init__(**kwargs) for key, val in zip(self._coeff_keys, (x0, y0, z0, r)): @@ -1953,12 +1947,14 @@ class Cone(QuadricMixin, Surface): _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): + def __init__(self, x0=0., y0=0., z0=0., r2=1., dx=0., dy=0., dz=1., *args, + **kwargs): R2 = kwargs.pop('R2', None) if R2 is not None: warn(_WARNING_UPPER.format(type(self).__name__, 'r2', 'R2'), FutureWarning) r2 = R2 + kwargs = _future_kwargs_warning_helper(type(self), *args, **kwargs) super().__init__(**kwargs) for key, val in zip(self._coeff_keys, (x0, y0, z0, r2, dx, dy, dz)): @@ -2134,12 +2130,13 @@ class XCone(QuadricMixin, Surface): _type = 'x-cone' _coeff_keys = ('x0', 'y0', 'z0', 'r2') - def __init__(self, x0=0., y0=0., z0=0., r2=1., **kwargs): + def __init__(self, x0=0., y0=0., z0=0., r2=1., *args, **kwargs): R2 = kwargs.pop('R2', None) if R2 is not None: warn(_WARNING_UPPER.format(type(self).__name__, 'r2', 'R2'), FutureWarning) r2 = R2 + kwargs = _future_kwargs_warning_helper(type(self), *args, **kwargs) super().__init__(**kwargs) for key, val in zip(self._coeff_keys, (x0, y0, z0, r2)): @@ -2262,12 +2259,13 @@ class YCone(QuadricMixin, Surface): _type = 'y-cone' _coeff_keys = ('x0', 'y0', 'z0', 'r2') - def __init__(self, x0=0., y0=0., z0=0., r2=1., **kwargs): + def __init__(self, x0=0., y0=0., z0=0., r2=1., *args, **kwargs): R2 = kwargs.pop('R2', None) if R2 is not None: warn(_WARNING_UPPER.format(type(self).__name__, 'r2', 'R2'), FutureWarning) r2 = R2 + kwargs = _future_kwargs_warning_helper(type(self), *args, **kwargs) super().__init__(**kwargs) for key, val in zip(self._coeff_keys, (x0, y0, z0, r2)): @@ -2390,12 +2388,13 @@ class ZCone(QuadricMixin, Surface): _type = 'z-cone' _coeff_keys = ('x0', 'y0', 'z0', 'r2') - def __init__(self, x0=0., y0=0., z0=0., r2=1., **kwargs): + def __init__(self, x0=0., y0=0., z0=0., r2=1., *args, **kwargs): R2 = kwargs.pop('R2', None) if R2 is not None: warn(_WARNING_UPPER.format(type(self).__name__, 'r2', 'R2'), FutureWarning) r2 = R2 + kwargs = _future_kwargs_warning_helper(type(self), *args, **kwargs) super().__init__(**kwargs) for key, val in zip(self._coeff_keys, (x0, y0, z0, r2)): @@ -2507,8 +2506,8 @@ class Quadric(QuadricMixin, 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., **kwargs): - + k=0., *args, **kwargs): + kwargs = _future_kwargs_warning_helper(type(self), *args, **kwargs) super().__init__(**kwargs) for key, val in zip(self._coeff_keys, (a, b, c, d, e, f, g, h, j, k)):