Fix RCC with negative height (#45)

Author: paulromano

src/openmc_mcnp_adapter/openmc_conversion.py

@@ -2,7 +2,7 @@
 # SPDX-License-Identifier: MIT
 
 import argparse
-from math import pi
+from math import pi, isclose
 import os
 import re
 import tempfile
@@ -72,26 +72,39 @@ def rotation_matrix(v1, v2):
     3x3 rotation matrix
 
     """
-    # Normalize vectors
+    # Normalize vectors and compute cosine
     u1 = v1 / np.linalg.norm(v1)
     u2 = v2 / np.linalg.norm(v2)
-
-    # Calculate axis of rotation
-    axis = np.cross(u1, u2)
-    axis /= np.linalg.norm(axis)
+    cos_angle = float(np.clip(np.dot(u1, u2), -1.0, 1.0))
 
     I = np.identity(3)
 
     # Handle special case where vectors are parallel or anti-parallel
-    if abs(np.dot(u2, axis) - 1.0) < 1e-8:
-        return I
-    elif abs(np.dot(u2, axis) + 1.0) < 1e-8:
-        return -I
+    if isclose(abs(cos_angle), 1.0, rel_tol=1e-8):
+        if cos_angle > 0.0:
+            return I
+        else:
+            # Proper 180° rotation: rotate about any axis that is orthogonal
+            # with u1. Because |k| = 1 and cos(180°) = -1, the rotation matrix
+            # is simply K = I + 2 * (k k^T - I) = 2 k k^T - I
+
+            # Choose reference vector not parallel to u1
+            ref = np.array([1.0, 0.0, 0.0]) if abs(u1[0]) < 0.9 else np.array([0.0, 1.0, 0.0])
+
+            # Create orthogonal unit vector
+            k = np.cross(u1, ref)
+            k /= np.linalg.norm(k)
+
+            # Create rotation matrix
+            return 2.0 * np.outer(k, k) - I
     else:
         # Calculate rotation angle
-        cos_angle = np.dot(u1, u2)
         sin_angle = np.sqrt(1 - cos_angle*cos_angle)
 
+        # Calculate axis of rotation
+        axis = np.cross(u1, u2)
+        axis /= np.linalg.norm(axis)
+
         # Create cross-product matrix K
         kx, ky, kz = axis
         K = np.array([
@@ -326,20 +339,11 @@ def flip_sense(surf):
                 raise NotImplementedError(f"{s['mnemonic']} surface with {len(coeffs)} parameters")
         elif s['mnemonic'] == 'rcc':
             vx, vy, vz, hx, hy, hz, r = coeffs
-            if hx == 0.0 and hy == 0.0:
-                if hz < 0.0:
-                    vz += hz
-                    hz = -hz
+            if hx == 0.0 and hy == 0.0 and hz > 0.0:
                 surf = RCC((vx, vy, vz), hz, r, axis='z')
-            elif hy == 0.0 and hz == 0.0:
-                if hx < 0.0:
-                    vx += hx
-                    hx = -hx
+            elif hy == 0.0 and hz == 0.0 and hx > 0.0:
                 surf = RCC((vx, vy, vz), hx, r, axis='x')
-            elif hx == 0.0 and hz == 0.0:
-                if hy < 0.0:
-                    vy += hy
-                    hy = -hy
+            elif hx == 0.0 and hz == 0.0 and hy > 0.0:
                 surf = RCC((vx, vy, vz), hy, r, axis='y')
             else:
                 # Create vectors for Z-axis and cylinder orientation

tests/test_rotation_matrix.py

@@ -0,0 +1,87 @@
+import numpy as np
+import pytest
+
+from openmc_mcnp_adapter import rotation_matrix
+
+
+def is_orthogonal(R: np.ndarray, atol: float = 1e-12) -> bool:
+    """Check if the matrix R is orthogonal"""
+    I = np.identity(3)
+    return np.allclose(R.T @ R, I, atol=atol) and np.allclose(R @ R.T, I, atol=atol)
+
+
+@pytest.mark.parametrize(
+    "v1, v2",
+    [
+        # Same direction, different magnitudes
+        (np.array([1.0, 2.0, 3.0]), np.array([2.0, 4.0, 6.0])),
+        (np.array([0.0, 0.0, 1.0]), np.array([0.0, 0.0, 10.0])),
+    ],
+)
+def test_rotation_parallel(v1, v2):
+    """Test rotation_matrix for parallel vectors"""
+    R = rotation_matrix(v1, v2)
+    # Should be exactly identity from the parallel branch
+    assert np.allclose(R, np.identity(3))
+    assert is_orthogonal(R)
+    assert np.isclose(np.linalg.det(R), 1.0)
+
+
+@pytest.mark.parametrize(
+    "v1, v2",
+    [
+        (np.array([0.0, 0.0, 1.0]), np.array([0.0, 0.0, -5.0])),  # z -> -z
+        (np.array([1.0, 0.0, 0.0]), np.array([-2.0, 0.0, 0.0])),  # x -> -x
+    ],
+)
+def test_rotation_antiparallel(v1, v2):
+    """Test rotation_matrix for anti-parallel vectors"""
+    R = rotation_matrix(v1, v2)
+
+    # Maps v1 direction to v2 direction
+    v2_hat = v2 / np.linalg.norm(v2)
+    mapped = R @ (v1 / np.linalg.norm(v1))
+    assert np.allclose(mapped, v2_hat, atol=1e-12)
+
+    # No NaNs, orthogonal and det +1
+    assert not np.isnan(R).any()
+    assert is_orthogonal(R)
+    assert np.isclose(np.linalg.det(R), 1.0, atol=1e-12)
+    # 180-degree rotation has trace -1
+    assert np.isclose(np.trace(R), -1.0, atol=1e-12)
+
+
+@pytest.mark.parametrize(
+    "v1, v2",
+    [
+        (np.array([0.0, 0.0, 1.0]), np.array([1.0, 2.0, 3.0])),
+        (np.array([0.0, 1.0, 0.0]), np.array([3.0, -1.0, 2.0])),
+    ],
+)
+def test_rotation_general(v1, v2):
+    """Test rotation_matrix for general vectors"""
+    R = rotation_matrix(v1, v2)
+
+    # Maps v1 to v2 direction
+    v2_hat = v2 / np.linalg.norm(v2)
+    mapped = R @ (v1 / np.linalg.norm(v1))
+    assert np.allclose(mapped, v2_hat, atol=1e-12)
+
+    # Orthogonal and proper rotation
+    assert is_orthogonal(R)
+    assert np.isclose(np.linalg.det(R), 1.0, atol=1e-12)
+
+    # A vector perpendicular to v1 remains perpendicular to R v1 (i.e., to v2_hat)
+    # Use a simple perpendicular vector: pick the least-aligned basis axis and project out
+    basis = [
+        np.array([1.0, 0.0, 0.0]),
+        np.array([0.0, 1.0, 0.0]),
+        np.array([0.0, 0.0, 1.0])
+    ]
+    u1_hat = v1 / np.linalg.norm(v1)
+    e = min(basis, key=lambda b: abs(np.dot(b, u1_hat)))
+    x = e - np.dot(e, u1_hat) * u1_hat
+    x /= np.linalg.norm(x)
+
+    # Should be perpendicular to v2_hat
+    assert np.isclose(np.dot(R @ x, v2_hat), 0.0, atol=1e-12)

tests/test_surfaces.py

@@ -5,7 +5,6 @@
 from openmc.model.surface_composite import OrthogonalBox, \
     RectangularParallelepiped, RightCircularCylinder, ConicalFrustum
 from openmc_mcnp_adapter import mcnp_str_to_model, get_openmc_surfaces
-import pytest
 from pytest import approx, mark
 
 
@@ -239,21 +238,25 @@ def test_rpp_macrobody():
 
 
 @mark.parametrize(
-    "coeffs, expected_bottom_z, expected_top_z, r",
+    "coeffs, expected_bottom, expected_top, r, coeff",
     [
-        # Base at (0,0,0), height +5 along z
-        ((0.0, 0.0, 0.0, 0.0, 0.0, 5.0, 1.5), 0.0, 5.0, 1.5),
-        # Negative height vector should be flipped internally (known failing behavior)
-        pytest.param((0.0, 0.0, 0.0, 0.0, 0.0, -5.0, 1.0), 0.0, -5.0, 1.0,
-                     marks=pytest.mark.xfail(reason="Negative height handling currently broken", strict=False)),
+        # Base at (0,0,0), positive/negative height along x
+        ((0.0, 0.0, 0.0, 5.0, 0.0, 0.0, 1.5), 0.0, 5.0, 1.5, 'a'),
+        ((0.0, 0.0, 0.0, -5.0, 0.0, 0.0, 1.0), 0.0, -5.0, 1.0, 'a'),
+        # Base at (0,0,0), positive/negative height along y
+        ((0.0, 0.0, 0.0, 0.0, 5.0, 0.0, 1.5), 0.0, 5.0, 1.5, 'b'),
+        ((0.0, 0.0, 0.0, 0.0, -5.0, 0.0, 1.0), 0.0, -5.0, 1.0, 'b'),
+        # Base at (0,0,0), positive/negative height along z
+        ((0.0, 0.0, 0.0, 0.0, 0.0, 5.0, 1.5), 0.0, 5.0, 1.5, 'c'),
+        ((0.0, 0.0, 0.0, 0.0, 0.0, -5.0, 1.0), 0.0, -5.0, 1.0, 'c'),
     ],
 )
-def test_rcc_macrobody(coeffs, expected_bottom_z, expected_top_z, r):
+def test_rcc_macrobody(coeffs, expected_bottom, expected_top, r, coeff):
     surf = convert_surface("rcc", coeffs)
     assert isinstance(surf, RightCircularCylinder)
     assert surf.cyl.r == approx(r)
-    assert surf.bottom.d / surf.bottom.c == approx(expected_bottom_z)
-    assert surf.top.d / surf.top.c == approx(expected_top_z)
+    assert surf.bottom.d / getattr(surf.bottom, coeff) == approx(expected_bottom)
+    assert surf.top.d / getattr(surf.top, coeff) == approx(expected_top)
 
 
 def test_trc_macrobody():