diff --git a/openmc/model/surface_composite.py b/openmc/model/surface_composite.py index 6f9123311fb..1fa4234fbf7 100644 --- a/openmc/model/surface_composite.py +++ b/openmc/model/surface_composite.py @@ -90,7 +90,7 @@ class CylinderSector(CompositeSurface): counterclockwise direction with respect to the first basis axis (+y, +z, or +x). Must be greater than :attr:`theta1`. center : iterable of float - Coordinate for central axes of cylinders in the (y, z), (z, x), or (x, y) + Coordinate for central axes of cylinders in the (y, z), (x, z), or (x, y) basis. Defaults to (0,0). axis : {'x', 'y', 'z'} Central axis of the cylinders defining the inner and outer surfaces of @@ -133,13 +133,13 @@ def __init__(self, phi2 = pi / 180 * theta2 # Coords for axis-perpendicular planes - p1 = np.array([0., 0., 1.]) + p1 = np.array([center[0], center[1], 1.]) - p2_plane1 = np.array([r1 * cos(phi1), r1 * sin(phi1), 0.]) - p3_plane1 = np.array([r2 * cos(phi1), r2 * sin(phi1), 0.]) + p2_plane1 = np.array([r1 * cos(phi1) + center[0], r1 * sin(phi1) + center[1], 0.]) + p3_plane1 = np.array([r2 * cos(phi1) + center[0], r2 * sin(phi1) + center[1], 0.]) - p2_plane2 = np.array([r1 * cos(phi2), r1 * sin(phi2), 0.]) - p3_plane2 = np.array([r2 * cos(phi2), r2 * sin(phi2), 0.]) + p2_plane2 = np.array([r1 * cos(phi2) + center[0], r1 * sin(phi2)+ center[1], 0.]) + p3_plane2 = np.array([r2 * cos(phi2) + center[0], r2 * sin(phi2)+ center[1], 0.]) points = [p1, p2_plane1, p3_plane1, p2_plane2, p3_plane2] if axis == 'z': @@ -147,7 +147,7 @@ def __init__(self, self.inner_cyl = openmc.ZCylinder(*center, r1, **kwargs) self.outer_cyl = openmc.ZCylinder(*center, r2, **kwargs) elif axis == 'y': - coord_map = [1, 2, 0] + coord_map = [0, 2, 1] self.inner_cyl = openmc.YCylinder(*center, r1, **kwargs) self.outer_cyl = openmc.YCylinder(*center, r2, **kwargs) elif axis == 'x': @@ -155,6 +155,7 @@ def __init__(self, self.inner_cyl = openmc.XCylinder(*center, r1, **kwargs) self.outer_cyl = openmc.XCylinder(*center, r2, **kwargs) + # Reorder the points to correspond to the correct central axis for p in points: p[:] = p[coord_map] @@ -192,8 +193,8 @@ def from_theta_alpha(cls, with respect to the first basis axis (+y, +z, or +x). Note that negative values translate to an offset in the clockwise direction. center : iterable of float - Coordinate for central axes of cylinders in the (y, z), (z, x), or (x, y) - basis. Defaults to (0,0). + Coordinate for central axes of cylinders in the (y, z), (x, z), or + (x, y) basis. Defaults to (0,0). axis : {'x', 'y', 'z'} Central axis of the cylinders defining the inner and outer surfaces of the sector. Defaults to 'z'. @@ -216,10 +217,16 @@ def from_theta_alpha(cls, return cls(r1, r2, theta1, theta2, center=center, axis=axis, **kwargs) def __neg__(self): - return -self.outer_cyl & +self.inner_cyl & -self.plane1 & +self.plane2 + if isinstance(self.inner_cyl, openmc.YCylinder): + return -self.outer_cyl & +self.inner_cyl & +self.plane1 & -self.plane2 + else: + return -self.outer_cyl & +self.inner_cyl & -self.plane1 & +self.plane2 def __pos__(self): - return +self.outer_cyl | -self.inner_cyl | +self.plane1 | -self.plane2 + if isinstance(self.inner_cyl, openmc.YCylinder): + return +self.outer_cyl | -self.inner_cyl | -self.plane1 | +self.plane2 + else: + return +self.outer_cyl | -self.inner_cyl | +self.plane1 | -self.plane2 class IsogonalOctagon(CompositeSurface): @@ -284,11 +291,11 @@ def __init__(self, center, r1, r2, axis='z', **kwargs): c1, c2 = center # Coords for axis-perpendicular planes - ctop = c1 + r1 - cbottom = c1 - r1 + cright = c1 + r1 + cleft = c1 - r1 - cright = c2 + r1 - cleft = c2 - r1 + ctop = c2 + r1 + cbottom = c2 - r1 # Side lengths if r2 > r1 * sqrt(2): @@ -309,7 +316,16 @@ def __init__(self, center, r1, r2, axis='z', **kwargs): p2_lr = np.array([L_basis_ax, -r1, 0.]) p3_lr = np.array([L_basis_ax, -r1, 1.]) - points = [p1_ur, p2_ur, p3_ur, p1_lr, p2_lr, p3_lr] + p1_ll = -p1_ur + p2_ll = -p2_ur + p3_ll = -p3_ur + + p1_ul = -p1_lr + p2_ul = -p2_lr + p3_ul = -p3_lr + + points = [p1_ur, p2_ur, p3_ur, p1_lr, p2_lr, p3_lr, + p1_ll, p2_ll, p3_ll, p1_ul, p2_ul, p3_ul] # Orientation specific variables if axis == 'z': @@ -331,17 +347,19 @@ def __init__(self, center, r1, r2, axis='z', **kwargs): self.right = openmc.YPlane(cright, **kwargs) self.left = openmc.YPlane(cleft, **kwargs) - # Put our coordinates in (x,y,z) order + # Put our coordinates in (x,y,z) order and add the offset for p in points: + p[0] += c1 + p[1] += c2 p[:] = p[coord_map] self.upper_right = openmc.Plane.from_points(p1_ur, p2_ur, p3_ur, **kwargs) self.lower_right = openmc.Plane.from_points(p1_lr, p2_lr, p3_lr, **kwargs) - self.lower_left = openmc.Plane.from_points(-p1_ur, -p2_ur, -p3_ur, + self.lower_left = openmc.Plane.from_points(p1_ll, p2_ll, p3_ll, **kwargs) - self.upper_left = openmc.Plane.from_points(-p1_lr, -p2_lr, -p3_lr, + self.upper_left = openmc.Plane.from_points(p1_ul, p2_ul, p3_ul, **kwargs) def __neg__(self): diff --git a/tests/unit_tests/test_surface_composite.py b/tests/unit_tests/test_surface_composite.py index 19221212e06..62ec18b3261 100644 --- a/tests/unit_tests/test_surface_composite.py +++ b/tests/unit_tests/test_surface_composite.py @@ -158,18 +158,23 @@ def test_cone_one_sided(axis, point_pos, point_neg, ll_true): @pytest.mark.parametrize( - "axis, indices", [ - ("X", [0, 1, 2]), - ("Y", [1, 2, 0]), - ("Z", [2, 0, 1]), + "axis, indices, center", [ + ("X", [2, 0, 1], (0., 0.)), + ("Y", [0, 2, 1], (0., 0.)), + ("Z", [0, 1, 2], (0., 0.)), + ("X", [2, 0, 1], (10., 5.)), + ("Y", [0, 2, 1], (10., 5.)), + ("Z", [0, 1, 2], (10., 5.)), + ] ) -def test_cylinder_sector(axis, indices): +def test_cylinder_sector(axis, indices, center): + c1, c2 = center r1, r2 = 0.5, 1.5 d = (r2 - r1) / 2 phi1 = -60. phi2 = 60 - s = openmc.model.CylinderSector(r1, r2, phi1, phi2, + s = openmc.model.CylinderSector(r1, r2, phi1, phi2, center=center, axis=axis.lower()) assert isinstance(s.outer_cyl, getattr(openmc, axis + "Cylinder")) assert isinstance(s.inner_cyl, getattr(openmc, axis + "Cylinder")) @@ -189,16 +194,18 @@ def test_cylinder_sector(axis, indices): assert np.all(np.isinf(ll)) assert np.all(np.isinf(ur)) ll, ur = (-s).bounding_box - assert ll == pytest.approx(np.roll([-np.inf, -r2, -r2], indices[0])) - assert ur == pytest.approx(np.roll([np.inf, r2, r2], indices[0])) + test_point_ll = np.array([-r2 + c1, -r2 + c2, -np.inf]) + assert ll == pytest.approx(test_point_ll[indices]) + test_point_ur = np.array([r2 + c1, r2 + c2, np.inf]) + assert ur == pytest.approx(test_point_ur[indices]) # __contains__ on associated half-spaces - point_pos = np.roll([0, r2 + 1, 0], indices[0]) - assert point_pos in +s - assert point_pos not in -s - point_neg = np.roll([0, r1 + d, r1 + d], indices[0]) - assert point_neg in -s - assert point_neg not in +s + point_pos = np.array([0 + c1, r2 + 1 + c2, 0]) + assert point_pos[indices] in +s + assert point_pos[indices] not in -s + point_neg = np.array([r1 + d + c1, r1 + d + c2, 0]) + assert point_neg[indices] in -s + assert point_neg[indices] not in +s # translate method t = uniform(-5.0, 5.0) @@ -399,7 +406,7 @@ def test_polygon(): with pytest.raises(ValueError): openmc.model.Polygon(rz_points) - # Test "M" shaped polygon + # Test "M" shaped polygon points = np.array([[8.5151581, -17.988337], [10.381711000000001, -17.988337], [12.744357, -24.288728000000003],