From 95e3f94630f017decf7d5ebc9cd4bcdfd6fb5875 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 4 Apr 2023 22:49:11 +0200 Subject: [PATCH] Implement CruciformPrism composite surface --- docs/source/pythonapi/model.rst | 3 +- openmc/model/surface_composite.py | 74 +++++++++++++++++++++- tests/unit_tests/test_surface_composite.py | 62 ++++++++++++++++++ 3 files changed, 137 insertions(+), 2 deletions(-) diff --git a/docs/source/pythonapi/model.rst b/docs/source/pythonapi/model.rst index cdabee396..c9c6e117a 100644 --- a/docs/source/pythonapi/model.rst +++ b/docs/source/pythonapi/model.rst @@ -24,14 +24,15 @@ Composite Surfaces :nosignatures: :template: myclass.rst + openmc.model.CruciformPrism openmc.model.CylinderSector openmc.model.IsogonalOctagon + openmc.model.Polygon openmc.model.RectangularParallelepiped openmc.model.RightCircularCylinder openmc.model.XConeOneSided openmc.model.YConeOneSided openmc.model.ZConeOneSided - openmc.model.Polygon TRISO Fuel Modeling ------------------- diff --git a/openmc/model/surface_composite.py b/openmc/model/surface_composite.py index 25826723c..6acfdde4a 100644 --- a/openmc/model/surface_composite.py +++ b/openmc/model/surface_composite.py @@ -630,7 +630,7 @@ class ZConeOneSided(CompositeSurface): class Polygon(CompositeSurface): - """Create a polygon composite surface from a path of closed points. + """Polygon formed from a path of closed points. .. versionadded:: 0.13.3 @@ -1079,3 +1079,75 @@ class Polygon(CompositeSurface): disp_vec = distance / costheta * unit_nvec return type(self)(self.points + disp_vec, basis=self.basis) + + +class CruciformPrism(CompositeSurface): + """Generalized cruciform prism + + This surface represents a prism parallel to an axis formed by planes at + multiple distances from the center. Equivalent to the 'gcross' derived + surface in Serpent. + + .. versionadded:: 0.13.4 + + Parameters + ---------- + distances : iterable of float + A monotonically increasing (or decreasing) iterable of distances in [cm] + that form the planes of the generalized cruciform. + center : iterable of float + The center of the prism in the two non-parallel axes (e.g., (x, y) when + axis is 'z') in [cm] + axis : {'x', 'y', 'z'} + Axis to which the prism is parallel + **kwargs + Keyword arguments passed to underlying plane classes + + """ + + def __init__(self, distances, center=(0., 0.), axis='z', **kwargs): + x0, y0 = center + self.distances = distances + + if axis == 'x': + cls_horizontal = openmc.YPlane + cls_vertical = openmc.ZPlane + elif axis == 'y': + cls_horizontal = openmc.XPlane + cls_vertical = openmc.ZPlane + elif axis == 'z': + cls_horizontal = openmc.XPlane + cls_vertical = openmc.YPlane + else: + raise ValueError("axis must be 'x', 'y', or 'z'") + + # Create each planar surface + surfnames = [] + for i, d in enumerate(distances): + setattr(self, f'hmin{i}', cls_horizontal(x0 - d, **kwargs)) + setattr(self, f'hmax{i}', cls_horizontal(x0 + d, **kwargs)) + setattr(self, f'vmin{i}', cls_vertical(y0 - d, **kwargs)) + setattr(self, f'vmax{i}', cls_vertical(y0 + d, **kwargs)) + surfnames.extend([f'hmin{i}', f'hmax{i}', f'vmin{i}', f'vmax{i}']) + + # Set _surfnames to satisfy CompositeSurface protocol + self._surfnames = tuple(surfnames) + + @property + def _surface_names(self): + return self._surfnames + + def __neg__(self): + n = len(self.distances) + regions = [] + for i in range(n): + regions.append( + +getattr(self, f'hmin{i}') & + -getattr(self, f'hmax{i}') & + +getattr(self, f'vmin{n-1-i}') & + -getattr(self, f'vmax{n-1-i}') + ) + return openmc.Union(regions) + + def __pos__(self): + return ~(-self) diff --git a/tests/unit_tests/test_surface_composite.py b/tests/unit_tests/test_surface_composite.py index 3f94893ed..066dd3307 100644 --- a/tests/unit_tests/test_surface_composite.py +++ b/tests/unit_tests/test_surface_composite.py @@ -315,6 +315,7 @@ def test_isogonal_octagon(axis, plane_tb, plane_lr, axis_idx): # Make sure repr works repr(s) + def test_polygon(): # define a 5 pointed star centered on 1, 1 star = np.array([[1. , 2. ], @@ -393,3 +394,64 @@ def test_polygon(): [6.88, 3.02]]) with pytest.raises(ValueError): openmc.model.Polygon(rz_points) + + +@pytest.mark.parametrize("axis", ["x", "y", "z"]) +def test_cruciform_prism(axis): + center = x0, y0 = (3., 4.) + distances = [2., 3., 5.] + s = openmc.model.CruciformPrism(distances, center, axis=axis) + + if axis == 'x': + i1, i2 = 1, 2 + elif axis == 'y': + i1, i2 = 0, 2 + elif axis == 'z': + i1, i2 = 0, 1 + plane_cls = (openmc.XPlane, openmc.YPlane, openmc.ZPlane) + + # Check type of surfaces + for i in range(3): + assert isinstance(getattr(s, f'hmin{i}'), plane_cls[i1]) + assert isinstance(getattr(s, f'hmax{i}'), plane_cls[i1]) + assert isinstance(getattr(s, f'vmin{i}'), plane_cls[i2]) + assert isinstance(getattr(s, f'vmax{i}'), plane_cls[i2]) + + # Make sure boundary condition propagates + s.boundary_type = 'reflective' + for i in range(3): + assert getattr(s, f'hmin{i}').boundary_type == 'reflective' + assert getattr(s, f'hmax{i}').boundary_type == 'reflective' + assert getattr(s, f'vmin{i}').boundary_type == 'reflective' + assert getattr(s, f'vmax{i}').boundary_type == 'reflective' + + # Check bounding box + ll, ur = (+s).bounding_box + assert np.all(np.isinf(ll)) + assert np.all(np.isinf(ur)) + ll, ur = (-s).bounding_box + assert ur[i1] == pytest.approx(x0 + distances[-1]) + assert ur[i2] == pytest.approx(y0 + distances[-1]) + assert ll[i1] == pytest.approx(x0 - distances[-1]) + assert ll[i2] == pytest.approx(y0 - distances[-1]) + + # __contains__ on associated half-spaces + point_pos, point_neg = np.zeros(3), np.zeros(3) + point_pos[i1] = x0 + 3.1 + point_pos[i2] = y0 + 2.05 + point_neg[i1] = x0 + 3.5 + point_neg[i2] = y0 + 1.99 + assert point_pos in +s + assert point_pos not in -s + assert point_neg in -s + assert point_neg not in +s + + # translate method + t = uniform(-5.0, 5.0) + s_t = s.translate((t, t, t)) + ll_t, ur_t = (-s_t).bounding_box + assert ur_t == pytest.approx(ur + t) + assert ll_t == pytest.approx(ll + t) + + # Make sure repr works + repr(s)