display: match McCode for hollow box/cylinder, box normal and annulus width

src/mccode_antlr/display/primitives.py

@@ -97,6 +97,26 @@ def _circle_with_normal(cx, cy, cz, r, nx, ny, nz, n: int = 24) -> np.ndarray:
     return pts
 
 
+def _y_onto(nx, ny, nz) -> np.ndarray:
+    """Rotation matrix taking local y onto (nx,ny,nz) by the smallest rotation.
+
+    As ``mcdis_rotated_box`` in McCode's ``mccode-r.c``: a half turn about z when
+    the normal is -y, and no rotation for a zero normal.
+    """
+    n = np.array([nx, ny, nz], dtype=float)
+    norm = np.linalg.norm(n)
+    if norm < 1e-10:
+        return np.eye(3)
+    n /= norm
+    a = np.cross([0.0, 1.0, 0.0], n)
+    s = np.linalg.norm(a)
+    if s < 1e-12:
+        return np.eye(3) if n[1] > 0 else np.diag([-1.0, -1.0, 1.0])
+    a /= s
+    k = np.array([[0, -a[2], a[1]], [a[2], 0, -a[0]], [-a[1], a[0], 0]])
+    return np.eye(3) + s * k + (1 - n[1]) * (k @ k)  # Rodrigues, sin = s, cos = n[1]
+
+
 # ---------------------------------------------------------------------------
 # Base primitive
 # ---------------------------------------------------------------------------
@@ -278,10 +298,10 @@ class Box(Primitive):
     """``box(cx,cy,cz, xw,yh,zd[, thickness[, nx,ny,nz]])`` — rectangular box.
 
     The new McCode overhaul adds an optional *thickness* (hollow-wall depth) and an
-    orientation normal (*nx*, *ny*, *nz*).  When *thickness* is zero or absent the
-    box is rendered as 12 wire edges, otherwise a second inner box is also drawn.
-    The normal is currently used only by the 3-D renderer; the wire representation
-    always uses the legacy axis-aligned form.
+    orientation normal (*nx*, *ny*, *nz*), onto which the box's local y axis is
+    rotated by the smallest rotation.  A *thickness* > 0 adds a closed inner cavity
+    ``xw-2t × yh-2t × zd-2t`` with the same centre and orientation, drawn only when
+    all three are positive, as ``mcdis_box`` in McCode's ``mccode-r.c``.
     """
     cx: Expr = field(default_factory=Expr._null)
     cy: Expr = field(default_factory=Expr._null)
@@ -303,31 +323,42 @@ def evaluate(self, params):
             condition=self.condition,
         )
 
-    def _corners(self, cx, cy, cz, dx, dy, dz):
-        return np.array([
+    def _corners(self, cx, cy, cz, dx, dy, dz, rot=None):
+        local = np.array([
             (cx-dx, cy-dy, cz-dz), (cx+dx, cy-dy, cz-dz),
             (cx+dx, cy+dy, cz-dz), (cx-dx, cy+dy, cz-dz),
             (cx-dx, cy-dy, cz+dz), (cx+dx, cy-dy, cz+dz),
             (cx+dx, cy+dy, cz+dz), (cx-dx, cy+dy, cz+dz),
         ])
+        if rot is None:
+            return local
+        c = np.array([cx, cy, cz])
+        return c + (local - c) @ rot.T
+
+    def _cavity(self, dx, dy, dz, t):
+        """Half-sizes of the inner cavity, or None for a solid box."""
+        inner = (dx - t, dy - t, dz - t)
+        return inner if t > 0 and min(inner) > 0 else None
 
     def to_polylines(self, params=None):
         p = params or {}
         cx, cy, cz = _eval(self.cx, p), _eval(self.cy, p), _eval(self.cz, p)
         dx, dy, dz = _eval(self.xw, p)/2, _eval(self.yh, p)/2, _eval(self.zd, p)/2
+        rot = _y_onto(_eval(self.nx, p), _eval(self.ny, p), _eval(self.nz, p))
         edges = [(0,1),(1,2),(2,3),(3,0),(4,5),(5,6),(6,7),(7,4),(0,4),(1,5),(2,6),(3,7)]
-        result = [self._corners(cx, cy, cz, dx, dy, dz)[list(e)] for e in edges]
-        t = _eval(self.thickness, p)
-        if t:
-            inner = self._corners(cx, cy, cz, max(0, dx-t/2), max(0, dy-t/2), max(0, dz-t/2))
+        result = [self._corners(cx, cy, cz, dx, dy, dz, rot)[list(e)] for e in edges]
+        cavity = self._cavity(dx, dy, dz, _eval(self.thickness, p))
+        if cavity:
+            inner = self._corners(cx, cy, cz, *cavity, rot)
             result += [inner[list(e)] for e in edges]
         return result
 
     def to_mesh(self, params=None):
         p = params or {}
         cx, cy, cz = _eval(self.cx, p), _eval(self.cy, p), _eval(self.cz, p)
         dx, dy, dz = _eval(self.xw, p)/2, _eval(self.yh, p)/2, _eval(self.zd, p)/2
-        verts = self._corners(cx, cy, cz, dx, dy, dz).astype(np.float32)
+        rot = _y_onto(_eval(self.nx, p), _eval(self.ny, p), _eval(self.nz, p))
+        verts = self._corners(cx, cy, cz, dx, dy, dz, rot)
         # 6 faces, each split into 2 triangles
         faces = np.array([
             [0,1,2],[0,2,3],  # -z face
@@ -337,7 +368,11 @@ def to_mesh(self, params=None):
             [0,3,7],[0,7,4],  # -x face
             [1,5,6],[1,6,2],  # +x face
         ], dtype=np.int32)
-        return verts, faces
+        cavity = self._cavity(dx, dy, dz, _eval(self.thickness, p))
+        if cavity:  # inner shell, wound the other way round
+            verts = np.vstack([verts, self._corners(cx, cy, cz, *cavity, rot)])
+            faces = np.vstack([faces, faces[:, ::-1] + 8])
+        return verts.astype(np.float32), faces
 
 
 @dataclass
@@ -394,7 +429,9 @@ def to_mesh(self, params=None):
 class Cylinder(Primitive):
     """``cylinder(cx,cy,cz, r,h[, thickness[, nx,ny,nz]])`` — cylinder.
 
-    The new McCode overhaul adds an optional *thickness* for hollow cylinders.
+    The new McCode overhaul adds an optional *thickness* for hollow cylinders: a
+    closed inner cavity of radius ``r-t`` and height ``h-2t``, drawn only when both
+    are positive, as ``mcdis_cylinder`` in McCode's ``mccode-r.c``.
     The axis orientation (*nx*, *ny*, *nz*) was already supported.
     """
     cx: Expr = field(default_factory=Expr._null)
@@ -425,29 +462,40 @@ def _axis_frame(self, nx, ny, nz):
             normal /= norm
         return normal
 
+    def _resolve(self, p):
+        """Centre, axis and the (r, h) of the outer surface plus the cavity, if any."""
+        centre = np.array([_eval(self.cx, p), _eval(self.cy, p), _eval(self.cz, p)])
+        normal = self._axis_frame(_eval(self.nx, p), _eval(self.ny, p), _eval(self.nz, p))
+        r, h, t = _eval(self.r, p), _eval(self.h, p), _eval(self.thickness, p)
+        shells = [(r, h)]
+        if t > 0 and r - t > 0 and h - 2 * t > 0:
+            shells.append((r - t, h - 2 * t))
+        return centre, normal, shells
+
     def to_polylines(self, params=None):
-        p = params or {}
-        cx, cy, cz = _eval(self.cx, p), _eval(self.cy, p), _eval(self.cz, p)
-        r, h = _eval(self.r, p), _eval(self.h, p)
-        nx, ny, nz = _eval(self.nx, p), _eval(self.ny, p), _eval(self.nz, p)
-        normal = self._axis_frame(nx, ny, nz)
+        centre, normal, shells = self._resolve(params or {})
+        return [pl for r, h in shells for pl in self._wire(centre, normal, r, h)]
+
+    def _wire(self, centre, normal, r, h):
         half = 0.5 * h * normal
-        c1 = np.array([cx, cy, cz]) - half
-        c2 = np.array([cx, cy, cz]) + half
+        c1 = centre - half
+        c2 = centre + half
         cap1 = _circle_with_normal(*c1, r, *normal)
         cap2 = _circle_with_normal(*c2, r, *normal)
         n_pts = cap1.shape[0] - 1
         lines = [np.array([cap1[i * n_pts // 4], cap2[i * n_pts // 4]]) for i in range(4)]
         return [cap1, cap2] + lines
 
     def to_mesh(self, params=None):
-        p = params or {}
-        cx, cy, cz = _eval(self.cx, p), _eval(self.cy, p), _eval(self.cz, p)
-        r, h = _eval(self.r, p), _eval(self.h, p)
-        nx, ny, nz = _eval(self.nx, p), _eval(self.ny, p), _eval(self.nz, p)
-        normal = self._axis_frame(nx, ny, nz)
+        centre, normal, shells = self._resolve(params or {})
+        verts, faces = self._solid(centre, normal, *shells[0])
+        if len(shells) > 1:  # inner shell, wound the other way round
+            v, f = self._solid(centre, normal, *shells[1])
+            verts, faces = np.vstack([verts, v]), np.vstack([faces, f[:, ::-1] + len(verts)])
+        return verts, faces
+
+    def _solid(self, centre, normal, r, h):
         half = 0.5 * h * normal
-        centre = np.array([cx, cy, cz])
         c1, c2 = centre - half, centre + half
         slices = 32
         # Build perimeter points for both caps (without the closing duplicate)
@@ -613,24 +661,26 @@ def to_mesh(self, params=None):
 
 @dataclass
 class Annulus(Primitive):
-    """``annulus(x,y,z,r_outer,r_inner,nx,ny,nz)`` — annular disc with arbitrary normal.
+    """``annulus(x,y,z,r_outer,thickness,nx,ny,nz)`` — annular disc with arbitrary normal.
 
-    Corresponds to ``mcdis_annulus`` in C.  Wire representation is two concentric
+    Corresponds to ``mcdis_annulus`` in C.  *thickness* is the ring width, as
+    McStas and McXtrace ``Annulus.comp`` pass it, so the inner radius is
+    ``max(0, r_outer - thickness)``.  Wire representation is two concentric
     circles; mesh representation is a ring of quads.
     """
     cx: Expr = field(default_factory=Expr._null)
     cy: Expr = field(default_factory=Expr._null)
     cz: Expr = field(default_factory=Expr._null)
     r_outer: Expr = field(default_factory=Expr._null)
-    r_inner: Expr = field(default_factory=Expr._null)
+    thickness: Expr = field(default_factory=Expr._null)
     nx: Expr = field(default_factory=lambda: Expr.float(0))
     ny: Expr = field(default_factory=lambda: Expr.float(1))
     nz: Expr = field(default_factory=lambda: Expr.float(0))
 
     def evaluate(self, params):
         return Annulus(
             _eval(self.cx, params), _eval(self.cy, params), _eval(self.cz, params),
-            _eval(self.r_outer, params), _eval(self.r_inner, params),
+            _eval(self.r_outer, params), _eval(self.thickness, params),
             _eval(self.nx, params), _eval(self.ny, params), _eval(self.nz, params),
             condition=self.condition,
         )
@@ -639,7 +689,7 @@ def to_polylines(self, params=None):
         p = params or {}
         cx, cy, cz = _eval(self.cx, p), _eval(self.cy, p), _eval(self.cz, p)
         ro = _eval(self.r_outer, p)
-        ri = _eval(self.r_inner, p)
+        ri = max(0.0, ro - _eval(self.thickness, p))
         nx, ny, nz = _eval(self.nx, p), _eval(self.ny, p), _eval(self.nz, p)
         return [
             _circle_with_normal(cx, cy, cz, ro, nx, ny, nz),
@@ -650,7 +700,7 @@ def to_mesh(self, params=None):
         p = params or {}
         cx, cy, cz = _eval(self.cx, p), _eval(self.cy, p), _eval(self.cz, p)
         ro = _eval(self.r_outer, p)
-        ri = _eval(self.r_inner, p)
+        ri = max(0.0, ro - _eval(self.thickness, p))
         nx, ny, nz = _eval(self.nx, p), _eval(self.ny, p), _eval(self.nz, p)
         slices = 32
         outer = _circle_with_normal(cx, cy, cz, ro, nx, ny, nz, n=slices)[:-1]

src/mccode_antlr/display/render/brep.py

@@ -208,11 +208,18 @@ def _brep_rectangle(prim, p):
 
 def _brep_box(prim, p):
     import build123d as bd
+    from ..primitives import _y_onto
     cx, cy, cz = _eval(prim.cx,p), _eval(prim.cy,p), _eval(prim.cz,p)
     xw, yh, zd = _eval(prim.xw,p), _eval(prim.yh,p), _eval(prim.zd,p)
     nx, ny, nz = _eval(prim.nx,p), _eval(prim.ny,p), _eval(prim.nz,p)
+    t = _eval(prim.thickness,p)
     box = bd.Box(xw, yh, zd)
-    return _orient(box, cx, cy, cz, nx, ny, nz)
+    if t > 0 and min(xw, yh, zd) - 2*t > 0:  # closed inner cavity, as mcdis_box
+        box = box - bd.Box(xw - 2*t, yh - 2*t, zd - 2*t)
+    # local y onto n by the smallest rotation, as mcdis_rotated_box
+    rot = _y_onto(nx, ny, nz)
+    plane = bd.Plane(origin=(cx, cy, cz), x_dir=tuple(rot[:, 0]), z_dir=tuple(rot[:, 2]))
+    return box.locate(bd.Location(plane))
 
 
 def _brep_sphere(prim, p):
@@ -225,7 +232,10 @@ def _brep_cylinder(prim, p):
     import build123d as bd
     cx, cy, cz = _eval(prim.cx,p), _eval(prim.cy,p), _eval(prim.cz,p)
     nx, ny, nz = _eval(prim.nx,p), _eval(prim.ny,p), _eval(prim.nz,p)
-    cyl = bd.Cylinder(_eval(prim.r,p), _eval(prim.h,p))
+    r, h, t = _eval(prim.r,p), _eval(prim.h,p), _eval(prim.thickness,p)
+    cyl = bd.Cylinder(r, h)
+    if t > 0 and r - t > 0 and h - 2*t > 0:  # closed inner cavity, as mcdis_cylinder
+        cyl = cyl - bd.Cylinder(r - t, h - 2*t)
     return _orient(cyl, cx, cy, cz, nx, ny, nz)
 
 
@@ -251,10 +261,10 @@ def _brep_annulus(prim, p):
     cx, cy, cz = _eval(prim.cx,p), _eval(prim.cy,p), _eval(prim.cz,p)
     nx, ny, nz = _eval(prim.nx,p), _eval(prim.ny,p), _eval(prim.nz,p)
     r_outer = _eval(prim.r_outer, p)
-    r_inner = _eval(prim.r_inner, p)
-    outer = bd.Cylinder(r_outer, 1e-3)
-    inner = bd.Cylinder(r_inner, 1e-3)
-    ring = outer - inner
+    r_inner = r_outer - _eval(prim.thickness, p)  # thickness is the ring width
+    ring = bd.Cylinder(r_outer, 1e-3)
+    if r_inner > 0:
+        ring = ring - bd.Cylinder(r_inner, 1e-3)
     return _orient(ring, cx, cy, cz, nx, ny, nz)
 
 

tests/test_display.py

@@ -31,6 +31,7 @@
     Sphere,
     parse_display_source,
 )
+from mccode_antlr.display.primitives import _y_onto
 
 
 # ---------------------------------------------------------------------------
@@ -577,3 +578,66 @@ def test_cone_mesh(self):
         verts, faces = c.to_mesh()
         assert verts.shape[1] == 3
         assert faces.shape[1] == 3
+
+
+# ===========================================================================
+# Hollow box/cylinder, box orientation and annulus width, as McCode's mccode-r.c
+# ===========================================================================
+
+def _extent(polylines):
+    pts = np.vstack(polylines)
+    return pts.max(axis=0) - pts.min(axis=0)
+
+
+class TestHollowShapes:
+    def test_box_normal_rotates_local_y(self):
+        # local y onto +x: a quarter turn about -z, so xwidth lies along y
+        b = Box(0, 0, 0, 0.1, 0.2, 0.3, 0, 1, 0, 0)
+        np.testing.assert_allclose(_extent(b.to_polylines()), [0.2, 0.1, 0.3], atol=1e-12)
+        verts, _ = b.to_mesh()
+        np.testing.assert_allclose(verts.max(axis=0) - verts.min(axis=0), [0.2, 0.1, 0.3], atol=1e-6)
+
+    def test_box_normal_minus_y_is_half_turn_about_z(self):
+        b = Box(0, 0, 0, 0.1, 0.2, 0.3, 0, 0, -1, 0)
+        corners = b._corners(0, 0, 0, 0.05, 0.1, 0.15, _y_onto(0, -1, 0))
+        np.testing.assert_allclose(corners[0], [0.05, 0.1, -0.15])
+
+    def test_box_cavity_is_two_thicknesses_smaller(self):
+        b = Box(0, 0, 0, 0.1, 0.2, 0.3, 0.01)
+        pls = b.to_polylines()
+        assert len(pls) == 24
+        np.testing.assert_allclose(_extent(pls[12:]), [0.08, 0.18, 0.28])
+        verts, faces = b.to_mesh()
+        assert verts.shape == (16, 3) and faces.shape == (24, 3)
+
+    def test_box_too_thick_is_solid(self):
+        b = Box(0, 0, 0, 0.1, 0.2, 0.3, 0.05)  # xwidth - 2t = 0
+        assert len(b.to_polylines()) == 12
+        assert b.to_mesh()[0].shape == (8, 3)
+
+    def test_cylinder_cavity(self):
+        c = Cylinder(0, 0, 0, 0.05, 0.5, 0.01, 0, 1, 0)
+        pls = c.to_polylines()
+        assert len(pls) == 12
+        np.testing.assert_allclose(_extent(pls[6:]), [0.08, 0.48, 0.08], atol=1e-3)
+        solid = Cylinder(0, 0, 0, 0.05, 0.5, 0, 0, 1, 0).to_mesh()
+        verts, faces = c.to_mesh()
+        assert len(verts) == 2 * len(solid[0]) and len(faces) == 2 * len(solid[1])
+
+    def test_cylinder_too_thick_is_solid(self):
+        c = Cylinder(0, 0, 0, 0.05, 0.5, 0.05, 0, 1, 0)  # r - t = 0
+        assert len(c.to_polylines()) == 6
+
+    def test_annulus_thickness_is_ring_width(self):
+        # as Annulus.comp: mcdis_annulus(0,0,0, radius, thickness, 0,1,0)
+        a = Annulus(0, 0, 0, 0.06, 0.02, 0, 1, 0)
+        outer, inner = (np.linalg.norm(pl[:, [0, 2]], axis=1) for pl in a.to_polylines())
+        np.testing.assert_allclose(outer, 0.06)
+        np.testing.assert_allclose(inner, 0.04)
+        verts, _ = a.to_mesh()
+        radii = np.unique(np.round(np.linalg.norm(verts[:, [0, 2]], axis=1), 6))
+        np.testing.assert_allclose(radii, [0.04, 0.06], rtol=1e-6)
+
+    def test_annulus_wider_than_radius_is_a_disc(self):
+        a = Annulus(0, 0, 0, 0.06, 0.1, 0, 1, 0)
+        np.testing.assert_allclose(a.to_polylines()[1], 0, atol=1e-12)

tests/test_display_brep.py

@@ -128,7 +128,7 @@ def test_disc(self):
 
     def test_annulus(self):
         prim = Annulus(cx=E(0), cy=E(0), cz=E(0),
-                       r_outer=E(0.1), r_inner=E(0.05),
+                       r_outer=E(0.1), thickness=E(0.05),
                        nx=E(0), ny=E(1), nz=E(0))
         shape = primitive_to_brep(prim)
         _check_shape(shape)
@@ -315,3 +315,36 @@ def test_cli_missing_instr_returns_nonzero(self, tmp_path):
         from mccode_antlr.cli.step import main
         rc = main([str(tmp_path / "nonexistent.instr")])
         assert rc != 0
+
+
+# ---------------------------------------------------------------------------
+# Hollow box/cylinder, box orientation and annulus width, as McCode's mccode-r.c
+# ---------------------------------------------------------------------------
+
+def _size(shape):
+    bb = shape.bounding_box()
+    return np.array([bb.size.X, bb.size.Y, bb.size.Z])
+
+
+class TestHollowShapes:
+    def test_box_default_normal_keeps_axes(self):
+        shape = primitive_to_brep(Box(E(0), E(0), E(0), E(0.1), E(0.2), E(0.3)))
+        np.testing.assert_allclose(_size(shape), [0.1, 0.2, 0.3], atol=1e-6)
+
+    def test_box_normal_rotates_local_y(self):
+        shape = primitive_to_brep(Box(E(0), E(0), E(0), E(0.1), E(0.2), E(0.3),
+                                      E(0), E(1), E(0), E(0)))
+        np.testing.assert_allclose(_size(shape), [0.2, 0.1, 0.3], atol=1e-6)
+
+    def test_box_cavity(self):
+        shape = primitive_to_brep(Box(E(0), E(0), E(0), E(0.1), E(0.2), E(0.3), E(0.01)))
+        assert shape.volume == pytest.approx(0.1*0.2*0.3 - 0.08*0.18*0.28, rel=1e-6)
+
+    def test_cylinder_cavity(self):
+        shape = primitive_to_brep(Cylinder(E(0), E(0), E(0), E(0.05), E(0.5), E(0.01),
+                                           E(0), E(1), E(0)))
+        assert shape.volume == pytest.approx(np.pi*(0.05**2*0.5 - 0.04**2*0.48), rel=1e-4)
+
+    def test_annulus_thickness_is_ring_width(self):
+        shape = primitive_to_brep(Annulus(E(0), E(0), E(0), E(0.06), E(0.02), E(0), E(1), E(0)))
+        assert shape.volume == pytest.approx(np.pi*(0.06**2 - 0.04**2)*1e-3, rel=1e-4)