From cb92d025376e5c637b05fbe6a46f5ab87a78802e Mon Sep 17 00:00:00 2001 From: Ethan Peterson Date: Tue, 25 Feb 2020 17:56:12 -0500 Subject: [PATCH] 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)