Merge pull request #1285 from paulromano/cylinder-points-fix

Fix cylinder_from_points function

openmc/model/funcs.py

@@ -395,9 +395,9 @@ def cylinder_from_points(p1, p2, r, **kwargs):
     dx = x2 - x1
     dy = y2 - y1
     dz = z2 - z1
-    cx = y1*z2 + y2*z1
-    cy = -(x1*z2 + x2*z1)
-    cz = x1*y2 + x2*y1
+    cx = y1*z2 - y2*z1
+    cy = x2*z1 - x1*z2
+    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|.

tests/unit_tests/test_surface.py

@@ -1,3 +1,6 @@
+from functools import partial
+from random import uniform, seed
+
 import numpy as np
 import openmc
 import pytest
@@ -329,14 +332,33 @@ def test_quadric():
 
 
 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)
+    seed(1)  # Make random numbers reproducible
+    for _ in range(100):
+        # Generate cylinder in random direction
+        xi = partial(uniform, -10.0, 10.0)
+        p1 = np.array([xi(), xi(), xi()])
+        p2 = np.array([xi(), xi(), xi()])
+        r = uniform(1.0, 100.0)
+        s = openmc.model.cylinder_from_points(p1, p2, r)
 
-    # 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
+        # Points p1 and p2 need to be inside cylinder
+        assert p1 in -s
+        assert p2 in -s
+
+        # Points further along the line should be inside cylinder as well
+        t = uniform(-100.0, 100.0)
+        p = p1 + t*(p2 - p1)
+        assert p in -s
+
+        # Check that points outside cylinder are in positive half-space and
+        # inside are in negative half-space. We do this by constructing a plane
+        # that includes the cylinder's axis, finding the normal to the plane,
+        # and using it to find a point slightly more/less than one radius away
+        # from the axis.
+        plane = openmc.Plane.from_points(p1, p2, (0., 0., 0.))
+        n = np.array([plane.a, plane.b, plane.c])
+        n /= np.linalg.norm(n)
+        assert p1 + 1.1*r*n in +s
+        assert p2 + 1.1*r*n in +s
+        assert p1 + 0.9*r*n in -s
+        assert p2 + 0.9*r*n in -s