From 132f85c8de301cc2853437dc9d23a901a68c5613 Mon Sep 17 00:00:00 2001 From: Peter Willendrup Date: Tue, 20 Oct 2026 13:46:08 +0200 Subject: [PATCH] display: match McCode for hollow box/cylinder, box normal and annulus width The MCDISPLAY primitives drew some shapes differently from the McCode runtime (mccode-r.c) and from McCode's own mcdisplay tools, so geometry shown or exported by mccode-antlr did not match the component physics. - box: the thickness cavity was w-t (half-sizes d-t/2) instead of the closed w-2t x h-2t x l-2t cavity mcdis_box draws, and the normal was ignored outside brep. The local y axis is now rotated onto n by the smallest rotation (half turn about z for n = -y), as mcdis_rotated_box, in to_polylines, to_mesh and brep. brep also had the axes of a default-normal box scrambled (0.1x0.2x0.3 came out 0.2x0.3x0.1). - box to_mesh and brep now include the cavity. - cylinder: thickness adds the closed r-t, h-2t cavity, as mcdis_cylinder, in to_polylines, to_mesh and brep. - annulus: the 5th argument is the ring width, as McStas and McXtrace Annulus.comp pass it, so the inner radius is max(0, r_outer - width). The field is renamed r_inner -> thickness to say so. Cavities are drawn only when all their sizes are positive, as in mccode-r.c. Fixes #398 Co-Authored-By: Claude Opus 5.5 --- src/mccode_antlr/display/primitives.py | 114 +++++++++++++++++------- src/mccode_antlr/display/render/brep.py | 22 +++-- tests/test_display.py | 64 +++++++++++++ tests/test_display_brep.py | 35 +++++++- 4 files changed, 196 insertions(+), 39 deletions(-) diff --git a/src/mccode_antlr/display/primitives.py b/src/mccode_antlr/display/primitives.py index 18c667de..77ef94e4 100644 --- a/src/mccode_antlr/display/primitives.py +++ b/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,23 +323,33 @@ 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 @@ -327,7 +357,8 @@ 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,15 +462,24 @@ 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 @@ -441,13 +487,15 @@ def to_polylines(self, params=None): 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,16 +661,18 @@ 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)) @@ -630,7 +680,7 @@ class Annulus(Primitive): 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] diff --git a/src/mccode_antlr/display/render/brep.py b/src/mccode_antlr/display/render/brep.py index 20c8dc1f..ccc6f7f6 100644 --- a/src/mccode_antlr/display/render/brep.py +++ b/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) diff --git a/tests/test_display.py b/tests/test_display.py index d3d5a0a5..bd5a521f 100644 --- a/tests/test_display.py +++ b/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) diff --git a/tests/test_display_brep.py b/tests/test_display_brep.py index 6b1adcee..12e2c4a7 100644 --- a/tests/test_display_brep.py +++ b/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)