Add Plane.from_points and cylinder_from_points convenience functions

docs/source/pythonapi/model.rst

@@ -11,6 +11,7 @@ Convenience Functions
    :template: myfunction.rst
 
    openmc.model.borated_water
+   openmc.model.cylinder_from_points
    openmc.model.get_hexagonal_prism
    openmc.model.get_rectangular_prism
    openmc.model.subdivide

openmc/model/funcs.py

@@ -3,7 +3,7 @@ from collections.abc import Iterable
 from math import sqrt
 from numbers import Real
 
-from openmc import XPlane, YPlane, Plane, ZCylinder
+from openmc import XPlane, YPlane, Plane, ZCylinder, Quadric
 from openmc.checkvalue import check_type, check_value
 import openmc.data
 
@@ -345,6 +345,54 @@ def get_hexagonal_prism(edge_length=1., orientation='y', origin=(0., 0.),
     return prism
 
 
+def cylinder_from_points(p1, p2, r, **kwargs):
+    """Return cylinder defined by two points passing through its center.
+
+    Parameters
+    ----------
+    p1, p2 : 3-tuples
+        Coordinates of two points that pass through the center of the cylinder
+    r : float
+        Radius of the cylinder
+    kwargs : dict
+        Keyword arguments passed to the :class:`openmc.Quadric` constructor
+
+    Returns
+    -------
+    openmc.Quadric
+        Quadric surface representing the cylinder.
+
+    """
+    # Get x, y, z coordinates of two points
+    x1, y1, z1 = p1
+    x2, y2, z2 = p2
+
+    # Define intermediate terms
+    dx = x2 - x1
+    dy = y2 - y1
+    dz = z2 - z1
+    cx = y1*z2 + y2*z1
+    cy = -(x1*z2 + x2*z1)
+    cz = x1*y2 + x2*y1
+
+    # Given p=(x,y,z), p1=(x1, y1, z1), p2=(x2, y2, z2), the equation for the
+    # cylinder can be derived as r = |(p - p1) ⨯ (p - p2)| / |p2 - p1|.
+    # Expanding out all terms and grouping according to what Quadric expects
+    # gives the following coefficients.
+    kwargs['a'] = dy*dy + dz*dz
+    kwargs['b'] = dx*dx + dz*dz
+    kwargs['c'] = dx*dx + dy*dy
+    kwargs['d'] = -2*dx*dy
+    kwargs['e'] = -2*dy*dz
+    kwargs['f'] = -2*dx*dz
+    kwargs['g'] = cy*dz - cz*dy
+    kwargs['h'] = cz*dx - cx*dz
+    kwargs['j'] = cx*dy - cy*dx
+    kwargs['k'] = -(dx*dx + dy*dy + dz*dz)*r*r
+
+    return openmc.Quadric(**kwargs)
+
+
 def subdivide(surfaces):
     """Create regions separated by a series of surfaces.
 

openmc/surface.py

@@ -456,6 +456,38 @@ class Plane(Surface):
                 element.set("periodic_surface_id", str(self.periodic_surface.id))
         return element
 
+    @classmethod
+    def from_points(cls, p1, p2, p3, **kwargs):
+        """Return a plane given three points that pass through it.
+
+        Parameters
+        ----------
+        p1, p2, p3 : 3-tuples
+            Points that pass through the plane
+        kwargs : dict
+            Keyword arguments passed to the :class:`Plane` constructor
+
+        Returns
+        -------
+        Plane
+            Plane that passes through the three points
+
+        """
+        # Convert to numpy arrays
+        p1 = np.asarray(p1)
+        p2 = np.asarray(p2)
+        p3 = np.asarray(p3)
+
+        # Find normal vector to plane by taking cross product of two vectors
+        # connecting p1->p2 and p1->p3
+        n = np.cross(p2 - p1, p3 - p1)
+
+        # The equation of the plane will by n·(<x,y,z> - p1) = 0. Determine
+        # coefficients A, B, C, and D based on that
+        A, B, C = n
+        D = np.dot(n, p1)
+        return cls(A=A, B=B, C=C, D=D, **kwargs)
+
 
 class XPlane(Plane):
     """A plane perpendicular to the x axis of the form :math:`x - x_0 = 0`

tests/unit_tests/test_surface.py

@@ -33,6 +33,20 @@ def test_plane():
     repr(s)
 
 
+def test_plane_from_points():
+    # Generate the plane x - y = 1 given three points
+    p1 = (0, -1, 0)
+    p2 = (1, 0, 0)
+    p3 = (1, 0, 1)
+    s = openmc.Plane.from_points(p1, p2, p3)
+
+    # Confirm correct coefficients
+    assert s.a == 1.0
+    assert s.b == -1.0
+    assert s.c == 0.0
+    assert s.d == 1.0
+
+
 def test_xplane():
     s = openmc.XPlane(x0=3., boundary_type='reflective')
     assert s.x0 == 3.
@@ -257,3 +271,17 @@ def test_quadric():
     # evaluate method
     assert s.evaluate((0., 0., 0.)) == pytest.approx(coeffs['k'])
     assert s.evaluate((1., 1., 1.)) == pytest.approx(3 + coeffs['k'])
+
+
+def test_cylinder_from_points():
+    # Generate 45-degree rotated cylinder in x-y plane with radius 1
+    p1 = (0, 0, 0)
+    p2 = (1, 1, 0)
+    s = openmc.model.cylinder_from_points(p1, p2, 1)
+
+    # Points p1 and p2 need to be inside cylinder
+    assert p1 in -s
+    assert p2 in -s
+    assert (-1, 1, 0) in +s
+    assert (1, -1, 0) in +s
+    assert (0, 0, 1.5) in +s