From 86fc40a6456e789676eec3410afff38e340ea8f8 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 19 Aug 2024 10:07:28 -0500 Subject: [PATCH] Add OrthogonalBox composite surface (#3118) --- docs/source/pythonapi/model.rst | 1 + openmc/model/surface_composite.py | 80 ++++++++++++++++++++++ openmc/surface.py | 13 +++- tests/unit_tests/test_surface_composite.py | 54 +++++++++++++++ 4 files changed, 145 insertions(+), 3 deletions(-) diff --git a/docs/source/pythonapi/model.rst b/docs/source/pythonapi/model.rst index 99459f64d09..21944018e7d 100644 --- a/docs/source/pythonapi/model.rst +++ b/docs/source/pythonapi/model.rst @@ -26,6 +26,7 @@ Composite Surfaces openmc.model.CylinderSector openmc.model.HexagonalPrism openmc.model.IsogonalOctagon + openmc.model.OrthogonalBox openmc.model.Polygon openmc.model.RectangularParallelepiped openmc.model.RectangularPrism diff --git a/openmc/model/surface_composite.py b/openmc/model/surface_composite.py index 29411fe4608..a2cb0243849 100644 --- a/openmc/model/surface_composite.py +++ b/openmc/model/surface_composite.py @@ -645,6 +645,86 @@ def __pos__(self): return +self.xmax | -self.xmin | +self.ymax | -self.ymin | +self.zmax | -self.zmin +class OrthogonalBox(CompositeSurface): + """Arbitrarily oriented orthogonal box + + This composite surface is composed of four or six planar surfaces that form + an arbitrarily oriented orthogonal box when combined. + + Parameters + ---------- + v : iterable of float + (x,y,z) coordinates of a corner of the box + a1 : iterable of float + Vector of first side starting from ``v`` + a2 : iterable of float + Vector of second side starting from ``v`` + a3 : iterable of float, optional + Vector of third side starting from ``v``. When not specified, it is + assumed that the box will be infinite along the vector normal to the + plane specified by ``a1`` and ``a2``. + **kwargs + Keyword arguments passed to underlying plane classes + + Attributes + ---------- + ax1_min, ax1_max : openmc.Plane + Planes representing minimum and maximum along first axis + ax2_min, ax2_max : openmc.Plane + Planes representing minimum and maximum along second axis + ax3_min, ax3_max : openmc.Plane + Planes representing minimum and maximum along third axis + + """ + _surface_names = ('ax1_min', 'ax1_max', 'ax2_min', 'ax2_max', 'ax3_min', 'ax3_max') + + def __init__(self, v, a1, a2, a3=None, **kwargs): + v = np.array(v) + a1 = np.array(a1) + a2 = np.array(a2) + if has_a3 := a3 is not None: + a3 = np.array(a3) + else: + a3 = np.cross(a1, a2) # normal to plane specified by a1 and a2 + + # Generate corners of box + p1 = v + p2 = v + a1 + p3 = v + a2 + p4 = v + a3 + p5 = v + a1 + a2 + p6 = v + a2 + a3 + p7 = v + a1 + a3 + + # Generate 6 planes of box + self.ax1_min = openmc.Plane.from_points(p1, p3, p4, **kwargs) + self.ax1_max = openmc.Plane.from_points(p2, p5, p7, **kwargs) + self.ax2_min = openmc.Plane.from_points(p1, p4, p2, **kwargs) + self.ax2_max = openmc.Plane.from_points(p3, p6, p5, **kwargs) + if has_a3: + self.ax3_min = openmc.Plane.from_points(p1, p2, p3, **kwargs) + self.ax3_max = openmc.Plane.from_points(p4, p7, p6, **kwargs) + + # Make sure a point inside the box produces the correct senses. If not, + # flip the plane coefficients so it does. + mid_point = v + (a1 + a2 + a3)/2 + nums = (1, 2, 3) if has_a3 else (1, 2) + for num in nums: + min_surf = getattr(self, f'ax{num}_min') + max_surf = getattr(self, f'ax{num}_max') + if mid_point in -min_surf: + min_surf.flip_normal() + if mid_point in +max_surf: + max_surf.flip_normal() + + def __neg__(self): + region = (+self.ax1_min & -self.ax1_max & + +self.ax2_min & -self.ax2_max) + if hasattr(self, 'ax3_min'): + region &= (+self.ax3_min & -self.ax3_max) + return region + + class XConeOneSided(CompositeSurface): """One-sided cone parallel the x-axis diff --git a/openmc/surface.py b/openmc/surface.py index 0fa2ae43949..2f10750a89b 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -802,6 +802,13 @@ def from_points(cls, p1, p2, p3, **kwargs): d = np.dot(n, p1) return cls(a=a, b=b, c=c, d=d, **kwargs) + def flip_normal(self): + """Modify plane coefficients to reverse the normal vector.""" + self.a = -self.a + self.b = -self.b + self.c = -self.c + self.d = -self.d + class XPlane(PlaneMixin, Surface): """A plane perpendicular to the x axis of the form :math:`x - x_0 = 0` @@ -1762,7 +1769,7 @@ class Cone(QuadricMixin, Surface): z-coordinate of the apex in [cm]. Defaults to 0. r2 : float, optional Parameter related to the aperture [:math:`\\rm cm^2`]. - It can be interpreted as the increase in the radius squared per cm along + It can be interpreted as the increase in the radius squared per cm along the cone's axis of revolution. dx : float, optional x-component of the vector representing the axis of the cone. @@ -1920,7 +1927,7 @@ class XCone(QuadricMixin, Surface): z-coordinate of the apex in [cm]. Defaults to 0. r2 : float, optional Parameter related to the aperture [:math:`\\rm cm^2`]. - It can be interpreted as the increase in the radius squared per cm along + It can be interpreted as the increase in the radius squared per cm along the cone's axis of revolution. boundary_type : {'transmission, 'vacuum', 'reflective', 'white'}, optional Boundary condition that defines the behavior for particles hitting the @@ -2021,7 +2028,7 @@ class YCone(QuadricMixin, Surface): z-coordinate of the apex in [cm]. Defaults to 0. r2 : float, optional Parameter related to the aperture [:math:`\\rm cm^2`]. - It can be interpreted as the increase in the radius squared per cm along + It can be interpreted as the increase in the radius squared per cm along the cone's axis of revolution. boundary_type : {'transmission, 'vacuum', 'reflective', 'white'}, optional Boundary condition that defines the behavior for particles hitting the diff --git a/tests/unit_tests/test_surface_composite.py b/tests/unit_tests/test_surface_composite.py index 01c827223fd..d862ae6b0de 100644 --- a/tests/unit_tests/test_surface_composite.py +++ b/tests/unit_tests/test_surface_composite.py @@ -498,3 +498,57 @@ def test_cruciform_prism(axis): openmc.model.CruciformPrism([1.0, 0.5, 2.0, 3.0]) with pytest.raises(ValueError): openmc.model.CruciformPrism([3.0, 2.0, 0.5, 1.0]) + + +def test_box(): + v = (-1.0, -1.0, -2.5) + a1 = (2.0, -1.0, 0.0) + a2 = (1.0, 2.0, 0.0) + a3 = (0.0, 0.0, 5.0) + s = openmc.model.OrthogonalBox(v, a1, a2, a3) + for num in (1, 2, 3): + assert isinstance(getattr(s, f'ax{num}_min'), openmc.Plane) + assert isinstance(getattr(s, f'ax{num}_max'), openmc.Plane) + + # Make sure boundary condition propagates + s.boundary_type = 'reflective' + assert s.boundary_type == 'reflective' + for num in (1, 2, 3): + assert getattr(s, f'ax{num}_min').boundary_type == 'reflective' + assert getattr(s, f'ax{num}_max').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 ll[2] == pytest.approx(-2.5) + assert ur[2] == pytest.approx(2.5) + + # __contains__ on associated half-spaces + assert (0., 0., 0.) in -s + assert (-2., 0., 0.) not in -s + assert (0., 0.9, 0.) in -s + assert (0., 0., -3.) in +s + assert (0., 0., 3.) in +s + + # translate method + s_t = s.translate((1., 1., 0.)) + assert (-0.01, 0., 0.) in +s_t + assert (0.01, 0., 0.) in -s_t + + # Make sure repr works + repr(s) + + # Version with infinite 3rd dimension + s = openmc.model.OrthogonalBox(v, a1, a2) + assert not hasattr(s, 'ax3_min') + assert not hasattr(s, 'ax3_max') + ll, ur = (-s).bounding_box + assert np.all(np.isinf(ll)) + assert np.all(np.isinf(ur)) + assert (0., 0., 0.) in -s + assert (-2., 0., 0.) not in -s + assert (0., 0.9, 0.) in -s + assert (0., 0., -3.) not in +s + assert (0., 0., 3.) not in +s