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,4 +1,5 @@
 import numpy as np
+from random import random
 import openmc
 import pytest
 
@@ -329,14 +330,28 @@ 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)
+    for _ in range(10):
+        # Generate cylinder in random direction
+        p1 = np.array([random(), random(), random()])
+        p2 = np.array([random(), random(), random()])
+        r = random()
+        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 = 100*random() - 200
+        p = p1 + t*(p2 - p1)
+        assert p in -s
+
+        # Check that a point outside cylinder is in positive 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 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