crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-27 21:55:41 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

initial unit tests working

Browse source
This commit is contained in:
Ethan Peterson 2020-02-25 17:56:12 -05:00
parent a1f7cb7387
commit cb92d02537
2 changed files with 90 additions and 151 deletions
Download patch file Download diff file Expand all files Collapse all files

212
openmc/surface.py
View file

@ -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

29
tests/unit_tests/test_surface.py
View file

@ -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)