diff --git a/CHANGELOG.md b/CHANGELOG.md index 056d19f..67896b6 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,15 @@ format follows [Keep a Changelog](https://keepachangelog.com/en/1.1.0/). ### Added +- `Helix(start, radius, pitch, turns, wire_radius, axis=..., left_handed=...)`: a helical + wire of round cross-section with hemispherical ends, e.g. the winding of an RF coil. Its + level set is the distance from the centreline (Newton steps along it from the + linearised minima of the nearest turns), exact within a wire radius of the wire and a + bound farther away; `centreline(t)`, `wire_length`, volume and surface area. +- `PolarRepeat(source, count, axis=..., origin=...)`: ``count`` copies of a 3D shape about + an axis, evaluated by folding the point's angle into the first sector, so its cost does + not grow with ``count`` (the slits of a Faraday shield: ``tube - PolarRepeat(slit, 24)``; + the cut cells of 24 slits as fast as of one). - A Numerics section in the docs (header and sidebar): an overview of the geometry (`docs/src/content/docs/numerics/index.mdx`, with the main cards) and one page per module, from the shapes to the grids (starlight-theme's `methodcard: { pages: true }`). diff --git a/docs/src/content/docs/getting-started/geometry-api.md b/docs/src/content/docs/getting-started/geometry-api.md index 2bccc0c..afb6795 100644 --- a/docs/src/content/docs/getting-started/geometry-api.md +++ b/docs/src/content/docs/getting-started/geometry-api.md @@ -146,6 +146,8 @@ and where `level_set_exact` is true it is the Euclidean signed distance. | `Aperture2D` | exact while the plates do not touch | | `Cylinder`, `Tube`, `Cone`, `Torus` (any axis, with or without fillets) | exact | | `SlantedTube` | the larger of the uncut tube's and the cut plane's: exact on the faces, a bound near the cut's rim | +| `Helix` | the distance from the centreline minus the wire's radius: exact within a wire radius of the wire, a bound farther away | +| `PolarRepeat` (copies about an axis) | the source's, the least over the point's sector and its two neighbours: a bound | | `a \| b`, `a & b`, `a - b` | `min(a, b)`, `max(a, b)`, `max(a, -b)`: a bound, not exact | | rotations, translations, reflections, uniform scaling | exact if the source is | | nonuniform scaling, shears | the source's times the smallest stretch: a bound | diff --git a/docs/src/content/docs/numerics/index.mdx b/docs/src/content/docs/numerics/index.mdx index f1429d0..df326b6 100644 --- a/docs/src/content/docs/numerics/index.mdx +++ b/docs/src/content/docs/numerics/index.mdx @@ -42,6 +42,8 @@ distances of boxes, segments, capped cones and tori: [Tube](/simgeometry/numerics/objects-tube/), [Cone](/simgeometry/numerics/objects-cone/), [Torus](/simgeometry/numerics/objects-torus/), and the [SlantedTube](/simgeometry/numerics/objects-slanted-tube/), a bound near its cut; +- the [Helix](/simgeometry/numerics/objects-helix/), a helical wire (an RF coil), its + distance found by Newton steps along the centreline; - [Booleans, affine transforms and sections](/simgeometry/numerics/operations/), whose level sets are bounds, and [collections](/simgeometry/numerics/collection/) of labelled objects. diff --git a/src/simgeometry/__init__.py b/src/simgeometry/__init__.py index 8a19279..6f547a6 100644 --- a/src/simgeometry/__init__.py +++ b/src/simgeometry/__init__.py @@ -16,6 +16,7 @@ Cylinder, FunctionObject2D, FunctionObject3D, + Helix, Line, Polygon, Polyline, @@ -26,6 +27,7 @@ Triangle, Tube, ) +from simgeometry.operations import PolarRepeat from simgeometry.simulation_domain import SimulationDomain # simgeometry's regions stay off unless the application enables the "simgeometry" tag. @@ -38,6 +40,7 @@ __all__ = [ "GeometryCollection", "Grid", + "PolarRepeat", "SimulationDomain", "Aperture2D", "Aperture3D", @@ -47,6 +50,7 @@ "Cylinder", "FunctionObject2D", "FunctionObject3D", + "Helix", "Line", "Polygon", "Polyline", diff --git a/src/simgeometry/objects/__init__.py b/src/simgeometry/objects/__init__.py index 563a8fa..29e96dc 100644 --- a/src/simgeometry/objects/__init__.py +++ b/src/simgeometry/objects/__init__.py @@ -9,6 +9,7 @@ from simgeometry.objects.cuboid import Cuboid from simgeometry.objects.cylinder import Cylinder from simgeometry.objects.function_object import FunctionObject2D, FunctionObject3D +from simgeometry.objects.helix import Helix from simgeometry.objects.line import Line from simgeometry.objects.polygon import Polygon from simgeometry.objects.polyline import Polyline @@ -31,6 +32,7 @@ "Cylinder", "FunctionObject2D", "FunctionObject3D", + "Helix", # 1D objects "Line", "Polygon", diff --git a/src/simgeometry/objects/helix.py b/src/simgeometry/objects/helix.py new file mode 100644 index 0000000..9558aaf --- /dev/null +++ b/src/simgeometry/objects/helix.py @@ -0,0 +1,335 @@ +"""Helical wires: the points within a radius of a helix, e.g. the winding of an RF coil.""" + +from __future__ import annotations + +import cunumpy as xp +import numpy as np +from methodcard import method + +from simgeometry._types import Array, BoundsLike, ConfigValue +from simgeometry.objects.axisymmetric import _vector +from simgeometry.shapes import Shape3D, _component_major + + +def _reference(axis: np.ndarray) -> np.ndarray: + """The unit vector normal to ``axis`` closest to x (to y for an axis along x).""" + helper = np.eye(3)[0] if abs(axis[0]) < 0.9 else np.eye(3)[1] + first = helper - (helper @ axis) * axis + return first / np.linalg.norm(first) + + +@method( + equations={ + "centreline": r"\mathbf c(t) = \mathbf p_0 + R\,(\cos t\,\mathbf u + \sin t\,\mathbf v) + \frac{p\,t}{2\pi}\,\mathbf a, \qquad 0 \le t \le 2\pi N", + "level_set": r"\phi(\mathbf x) = \min_{0 \le t \le 2\pi N} |\mathbf x - \mathbf c(t)| - r_w", + "descent": r"t \leftarrow t - \frac{f'(t)}{f''(t) \text{ or } 2\rho R + 2b^2}, \qquad f(t) = \rho^2 + R^2 - 2\rho R\cos(t - \varphi) + (s - b t)^2, \quad b = \frac{p}{2\pi}", + "bound": r"|\mathbf x - \mathbf c(t)| \ge \sqrt{(\rho - R)^2 + \max(-s, s - pN, 0)^2}", + }, + assumptions=( + r"$\mathbf a$ the unit axis, $\mathbf u$ the direction of the start point from the axis, $\mathbf v = \pm\mathbf a\times\mathbf u$ (right- or left-handed); $(s, \rho, \varphi)$ the point's cylindrical coordinates in that frame", + r"the wire's ends are hemispheres (the points within $r_w$ of the end points of the centreline)", + r"within $2 r_w$ of the cylinder the centreline lies on: the minimum over $t$ from four starts, the two ends and the linearised minima of the two turns around $s / b$, by six Newton steps (descent steps on the majorant $f'' \le 2\rho R + 2b^2$ where $f''$ is not clearly positive), clipped to $[0, 2\pi N]$; exact there (tested against dense sampling)", + r"farther away $\phi$ is the bound (the distance from the cylinder, minus $r_w$): a level set that never exceeds the distance", + ), +) +class Helix(Shape3D): + """A helical wire of round cross-section: the points within ``wire_radius`` of a helix. + + The centreline winds ``turns`` times around ``axis`` at ``radius``, advancing + ``pitch`` along the axis per turn, from its start point at + ``start + radius * u`` (``u`` the unit vector normal to the axis in the + direction of ``start_direction``, by default the one closest to x) to its end + ``pitch * turns`` further along the axis. Right-handed by default (counter- + clockwise seen from the end the axis points to); ``left_handed`` reverses it. The + wire's ends are hemispherical caps. The `level_set` is the signed distance to + the centreline minus the wire's radius, found by Newton steps on the squared + distance along the centreline, within a wire radius of the wire; farther away it is + a lower bound (``level_set_exact`` is False). + + Example + ------- + The six windings of an RF coil, radius 0.14 m, 2 cm pitch, a copper tube of 6 mm + outer diameter, along z from z = 0.02:: + + coil = Helix((0, 0, 0.02), 0.14, 0.02, 6, 0.003) + """ + + has_level_set = True + level_set_exact = False + + def _bounds(self) -> tuple[tuple[float, ...], tuple[float, ...]]: + """The box of the cylinder of radius ``radius + wire_radius`` around the wire.""" + extent = (self.radius + self.wire_radius) * np.sqrt( + np.clip(1.0 - self._axis**2, 0.0, None) + ) + first = self._start - self.wire_radius * self._axis + last = self._start + (self.length + self.wire_radius) * self._axis + lower = np.minimum(first, last) - extent + upper = np.maximum(first, last) + extent + return tuple(float(c) for c in lower), tuple(float(c) for c in upper) + + def __init__( + self, + start: BoundsLike, + radius: float, + pitch: float, + turns: float, + wire_radius: float, + axis: BoundsLike = (0.0, 0.0, 1.0), + start_direction: BoundsLike | None = None, + left_handed: bool = False, + bc: ConfigValue = None, + inverted: bool = False, + ) -> None: + """ + Initialize a helical wire in 3D. + + Parameters: + start (list or tuple): The point on the axis level with the wire's start. + radius (float): The radius of the centreline. + pitch (float): The advance along the axis per turn, positive. + turns (float): The number of turns, positive (need not be whole). + wire_radius (float): The wire's radius, below ``radius``. + axis (list or tuple, optional): The direction the helix advances along. + start_direction (list or tuple, optional): The direction from the axis to the + start point (its part normal to the axis); default the one closest to x. + left_handed (bool, optional): Wind clockwise seen from the end ``axis`` + points to. + bc (optional): Boundary condition information. + inverted (bool, optional): Cover everything except the wire. + """ + axis = _vector(axis, "axis") + if np.linalg.norm(axis) == 0: + raise ValueError("The axis must have a nonzero length") + axis = axis / np.linalg.norm(axis) + if not 0 < wire_radius < radius: + raise ValueError("A helix needs 0 < wire_radius < radius") + if pitch <= 0 or turns <= 0: + raise ValueError("A helix needs pitch > 0 and turns > 0") + if start_direction is None: + u = _reference(axis) + else: + direction = _vector(start_direction, "start_direction") + u = direction - (direction @ axis) * axis + if np.linalg.norm(u) < 1e-12: + raise ValueError("start_direction must not be along the axis") + u = u / np.linalg.norm(u) + v = np.cross(axis, u) + self._start = _vector(start, "start") + self._axis = axis + self._u = u + self._v = -v if left_handed else v + self.radius = float(radius) + self.pitch = float(pitch) + self.turns = float(turns) + self.wire_radius = float(wire_radius) + self._bc = bc + self._inverted = bool(inverted) + super().__init__(*self._bounds()) + + @property + def start(self) -> tuple[float, float, float]: + """The point on the axis level with the wire's start.""" + return tuple(float(c) for c in self._start) + + @property + def axis(self) -> tuple[float, float, float]: + """The unit vector the helix advances along.""" + return tuple(float(c) for c in self._axis) + + @property + def start_direction(self) -> tuple[float, float, float]: + """The unit vector from the axis to the start point.""" + return tuple(float(c) for c in self._u) + + @property + def left_handed(self) -> bool: + """Whether the helix winds clockwise seen from the end the axis points to.""" + return bool(np.cross(self._axis, self._u) @ self._v < 0) + + @property + def length(self) -> float: + """The axial length of the centreline, ``pitch * turns``.""" + return self.pitch * self.turns + + @property + def wire_length(self) -> float: + """The length of the centreline, ``turns * sqrt((2 pi radius)^2 + pitch^2)``.""" + return float(self.turns * np.hypot(2 * np.pi * self.radius, self.pitch)) + + def centreline(self, t: Array) -> Array: + """The points ``(3, n)`` of the centreline at the angles ``t`` (``0`` to ``2 pi N``).""" + t = xp.asarray(t, dtype=float) + b = self.pitch / (2 * np.pi) + return ( + xp.asarray(self._start)[:, None] + + self.radius * xp.cos(t)[None, :] * xp.asarray(self._u)[:, None] + + self.radius * xp.sin(t)[None, :] * xp.asarray(self._v)[:, None] + + b * t[None, :] * xp.asarray(self._axis)[:, None] + ) + + def copy(self) -> Helix: + """Create a copy of the helix.""" + return Helix( + self.start, + self.radius, + self.pitch, + self.turns, + self.wire_radius, + axis=self.axis, + start_direction=self.start_direction, + left_handed=self.left_handed, + bc=self._bc, + inverted=self._inverted, + ) + + def _frame_coordinates(self, points: Array): + """Axial position ``s``, distance ``rho`` from the axis and angle ``phi``.""" + relative = points - xp.asarray(self._start)[:, None] + s = xp.sum(relative * xp.asarray(self._axis)[:, None], axis=0) + x = xp.sum(relative * xp.asarray(self._u)[:, None], axis=0) + y = xp.sum(relative * xp.asarray(self._v)[:, None], axis=0) + return s, xp.sqrt(x * x + y * y), xp.arctan2(y, x) + + def _bound(self, s: Array, rho: Array) -> Array: + """A lower bound of the distance from the centreline: from its cylinder.""" + beyond = xp.maximum(xp.maximum(-s, s - self.length), 0.0) + return xp.sqrt((rho - self.radius) ** 2 + beyond * beyond) + + def _nearest(self, s: Array, rho: Array, phi: Array) -> Array: + """The distance from the centreline, by descent along it from four starts.""" + radius, b = self.radius, self.pitch / (2 * np.pi) + top, turn = 2 * np.pi * self.turns, 2 * np.pi + coupling = rho * radius + majorant = 2 * coupling + 2 * b * b + k = xp.floor((s / b - phi) / turn) + starts = [] + for j in (0, 1): + # the minimum of the cosine's basin around phi + 2 pi (k + j), linearised: + # rho R sin(tau) = b (s - b t) with sin(tau) = tau + base = phi + turn * (k + j) + # moved by whole turns into the range (points beyond the ends) + base = xp.where( + base > top, base - turn * xp.ceil((base - top) / turn), base + ) + base = xp.where(base < 0, base + turn * xp.ceil(-base / turn), base) + starts.append(base + b * (s - b * base) / (coupling + b * b)) + + def descend(t, s, phi, coupling, majorant): + t = xp.clip(t, 0.0, top) + for _ in range(5): + angle = t - phi + slope = 2 * coupling * xp.sin(angle) - 2 * b * (s - b * t) + curvature = 2 * coupling * xp.cos(angle) + 2 * b * b + # Newton where the curvature is clearly positive, else a descent + # step on the majorant f'' <= 2 rho R + 2 b^2 + divisor = xp.where(curvature > 0.05 * majorant, curvature, majorant) + t = xp.clip(t - slope / divisor, 0.0, top) + return t + + def squared(t, s, rho, phi): + return ( + rho * rho + + radius**2 + - 2 * rho * radius * xp.cos(t - phi) + + (s - b * t) ** 2 + ) + + best = None + for t in starts: + t = descend(t, s, phi, coupling, majorant) + value = squared(t, s, rho, phi) + best = value if best is None else xp.minimum(best, value) + # the two ends, for the points within a pitch of them (axially) + margin = self.pitch + self.wire_radius + for end in (0.0, top): + near = xp.flatnonzero(xp.abs(s - b * end) < margin) + if near.shape[0]: + sub = (s[near], phi[near], coupling[near], majorant[near]) + t = descend(xp.full(near.shape, end), *sub) + best[near] = xp.minimum( + best[near], squared(t, s[near], rho[near], phi[near]) + ) + return xp.sqrt(xp.maximum(best, 0.0)) + + def _distance(self, points: Array, reach: float) -> Array: + """The distance from the centreline: exact within ``reach`` of its cylinder. + + Farther away the lower bound of `_bound` (the distance from the cylinder the + centreline lies on). + """ + s, rho, phi = self._frame_coordinates(points) + distance = self._bound(s, rho) + near = xp.flatnonzero(distance < reach) + if near.shape[0]: + distance[near] = self._nearest(s[near], rho[near], phi[near]) + return distance + + def _level_set(self, points: Array) -> Array: + """The distance from the centreline minus the wire's radius (a bound beyond it).""" + return self._distance(points, 2 * self.wire_radius) - self.wire_radius + + def contains_points(self, points: Array) -> Array: + """Vectorised `contains`: a mask of the ``(3, n)`` points inside.""" + points = _component_major(points, 3) + inside = self._distance(points, self.wire_radius) <= self.wire_radius + return ~inside if self.inverted else inside + + def contains(self, coords: Array) -> bool: + """ + Check if the given coordinates are inside the wire. + + Parameters: + coords (list or tuple): Coordinates (x, y, z). + + Returns: + bool: True if the point is inside the object, False otherwise. + """ + point = xp.asarray(coords, dtype=float).reshape(3, 1) + return bool(xp.to_numpy(self.contains_points(point))[0]) + + def translate(self, offset: tuple[float, float, float]) -> None: + """ + Translate the helix by the given offset. + + Parameters: + offset (tuple): (dx, dy, dz) translation vector. + """ + self._start = self._start + _vector(offset, "offset") + self.set_bounding_box(*self._bounds()) + + def mirror(self, axis: int, position: float) -> None: + """ + Mirror the helix across the plane where coordinate ``axis`` is ``position``. + + A mirror image of a right-handed helix is left-handed. + + Parameters: + axis (int): Axis index (0 for x-axis, 1 for y-axis, 2 for z-axis). + position (float): Position along the axis to mirror. + """ + if axis < 0 or axis > 2: + raise ValueError("Axis must be 0 (x-axis), 1 (y-axis), or 2 (z-axis).") + flip = np.ones(3) + flip[axis] = -1.0 + self._start = self._start * flip + self._start[axis] += 2 * position + self._axis = self._axis * flip + self._u = self._u * flip + self._v = self._v * flip + self.set_bounding_box(*self._bounds()) + + def volume(self) -> float: + """The volume: the swept disc along the centreline and the two hemispheres. + + Exact (Pappus for a curve) as long as the wire does not overlap itself + (``pitch >= 2 * wire_radius`` and ``wire_radius`` below the radius of + curvature). + """ + disc = np.pi * self.wire_radius**2 + return float(disc * self.wire_length + 4 / 3 * np.pi * self.wire_radius**3) + + def surface_area(self) -> float: + """The surface area: the swept circle and the two hemispheres (no overlap).""" + r = self.wire_radius + return float(2 * np.pi * r * self.wire_length + 4 * np.pi * r**2) diff --git a/src/simgeometry/operations.py b/src/simgeometry/operations.py index 37ec1df..6cb27db 100644 --- a/src/simgeometry/operations.py +++ b/src/simgeometry/operations.py @@ -397,6 +397,106 @@ def evaluate(self, *args, quantity=None, **kwargs): return super().evaluate(*args, quantity=quantity, **kwargs) +@method( + equations={ + "fold": r"\varphi' = \varphi - k\,\Delta, \quad k = \operatorname{round}(\varphi / \Delta), \quad \Delta = \frac{2\pi}{N}", + "level_set": r"\phi_N(p) = \min_{j \in \{k - 1, k, k + 1\}} \phi\big(R_{-j\Delta}\, p\big)", + }, + assumptions=( + r"$N$ copies of the source rotated about the axis by multiples of $\Delta$; $\varphi$ the point's angle about the axis from the reference direction (the unit vector normal to the axis closest to x)", + r"the source lies within its sector $|\varphi| < \Delta / 2$ (copies do not overlap), so membership is the source's at the point rotated into the first sector", + r"the level set is the smallest over the point's sector and its two neighbours: a bound on the distance where the nearest copy is among them (copies away from the axis)", + ), +) +class PolarRepeat(Shape): + """``count`` copies of a 3D shape rotated about an axis, evaluated as one. + + Copy ``j`` is the source rotated by ``j * 360 / count`` degrees about ``axis`` + through ``origin``. A point is mapped into the first copy's sector (its angle + folded to ``|phi| < 180 / count`` degrees) and tested against the source alone, + so the cost does not grow with ``count`` as a union of rotated copies would: the + 24 slits of a Faraday shield are ``tube - PolarRepeat(slit, 24)``. The source must + lie within its sector, centred on the reference direction (the unit vector + normal to the axis closest to x, as `Helix` uses). + """ + + level_set_exact = False + + def __init__(self, source, count, *, axis=(0.0, 0.0, 1.0), origin=None): + if source.ndim != 3: + raise ValueError("PolarRepeat needs a 3D shape") + if int(count) != count or count < 1: + raise ValueError("count must be a positive integer") + axis = np.asarray(axis, float) + if axis.shape != (3,) or not np.all(np.isfinite(axis)) or not np.any(axis): + raise ValueError("The axis must be a finite nonzero 3-vector") + self.source = copy.deepcopy(source) + self.count = int(count) + self.axis = axis / np.linalg.norm(axis) + self.origin = np.zeros(3) if origin is None else np.asarray(origin, float) + helper = np.eye(3)[0] if abs(self.axis[0]) < 0.9 else np.eye(3)[1] + u = helper - (helper @ self.axis) * self.axis + self._u = u / np.linalg.norm(u) + self._v = np.cross(self.axis, self._u) + self._inverted = source.inverted + self._bc = getattr(source, "_bc", None) + relative = corners(source) - self.origin[:, None] + along = self.axis @ relative + radius = float( + np.max(np.linalg.norm(relative - np.outer(self.axis, along), axis=0)) + ) + extent = radius * np.sqrt(np.clip(1.0 - self.axis**2, 0.0, None)) + ends = self.origin[:, None] + np.outer(self.axis, [along.min(), along.max()]) + super().__init__(3, ends.min(axis=1) - extent, ends.max(axis=1) + extent) + + def _folded(self, points, shift=0): + """The points rotated by ``-(k + shift)`` sector angles, ``k`` their sector.""" + relative = points - xp.asarray(self.origin)[:, None] + axis, u, v = (xp.asarray(w)[:, None] for w in (self.axis, self._u, self._v)) + s = xp.sum(relative * axis, axis=0) + x = xp.sum(relative * u, axis=0) + y = xp.sum(relative * v, axis=0) + sector = 2 * np.pi / self.count + phi = xp.arctan2(y, x) + phi = phi - sector * (xp.round(phi / sector) + shift) + rho = xp.sqrt(x * x + y * y) + return ( + xp.asarray(self.origin)[:, None] + + s[None, :] * axis + + (rho * xp.cos(phi))[None, :] * u + + (rho * xp.sin(phi))[None, :] * v + ) + + def contains_points(self, points): + """Membership of the ``(3, n)`` points: the source's in the first sector.""" + points = _component_major(points, 3) + return self.source.contains_points(self._folded(points)) + + def contains(self, coords): + """Whether the point is inside one of the copies.""" + point = xp.asarray(coords, dtype=float).reshape(3, 1) + return bool(xp.to_numpy(self.contains_points(point))[0]) + + @property + def has_level_set(self): + """Whether the source has a level set.""" + return bool(getattr(self.source, "has_level_set", False)) + + def level_set(self, points): + """The smallest of the source's level set in the point's and the next sectors.""" + if not self.has_level_set: + raise NotImplementedError(f"{type(self.source).__name__} has no level set") + points = _component_major(points, 3) + values = [ + self.source.level_set(self._folded(points, shift)) for shift in (-1, 0, 1) + ] + return xp.minimum(xp.minimum(values[0], values[1]), values[2]) + + def volume(self): + """``count`` times the source's volume (the copies do not overlap).""" + return self.count * self.source.volume() + + def rotation(shape, angle, *, axis=(0, 0, 1), origin=None): angle = np.deg2rad(float(angle)) c, s = np.cos(angle), np.sin(angle) diff --git a/src/simgeometry/tests/test_helix.py b/src/simgeometry/tests/test_helix.py new file mode 100644 index 0000000..4e11771 --- /dev/null +++ b/src/simgeometry/tests/test_helix.py @@ -0,0 +1,154 @@ +"""Helix: the level set against dense sampling, volume, bounds, frames, crossings, Yee fractions.""" + +from __future__ import annotations + +import numpy as np +import pytest + +from simgeometry import Grid, Helix +from simgeometry.crossings import crossings + +RNG = np.random.default_rng(5) + +SHAPES = { + "coil": Helix((0, 0, 0.02), 0.14, 0.02, 6, 0.003), + "oblique left": Helix( + (0.1, -0.2, 0.3), 0.05, 0.03, 2.5, 0.008, axis=(1, 1, 0.5), left_handed=True + ), + "steep": Helix((0, 0, 0), 0.02, 0.2, 1.3, 0.005, start_direction=(0, 1, 0)), +} + + +def _box_points(shape, n, pad=0.0): + lower = np.asarray(shape.lower_bound) - pad + upper = np.asarray(shape.upper_bound) + pad + return RNG.uniform(lower[:, None], upper[:, None], (3, n)), np.prod(upper - lower) + + +def _sampled_distance(shape, points, samples=100_000): + t = np.linspace(0, 2 * np.pi * shape.turns, samples) + curve = np.asarray(shape.centreline(t)) + return np.array( + [np.min(np.linalg.norm(curve - p[:, None], axis=0)) for p in points.T] + ) + + +@pytest.mark.parametrize("name", SHAPES) +def test_level_set_is_the_distance_to_the_centreline(name): + shape = SHAPES[name] + box, _ = _box_points(shape, 300, pad=0.03) + along = np.asarray(shape.centreline(RNG.uniform(0, 2 * np.pi * shape.turns, 300))) + points = np.concatenate( + [box, along + RNG.normal(0, shape.wire_radius, along.shape)], axis=1 + ) + expected = _sampled_distance(shape, points) - shape.wire_radius + values = np.asarray(shape.level_set(points)) + # exact within a wire radius of the wire (the dense sampling's own error is + # below 1e-6 here), a lower bound farther away + near = expected < shape.wire_radius + assert near.sum() > 100 + np.testing.assert_allclose(values[near], expected[near], atol=2e-6) + assert np.all(values <= expected + 2e-6) + inside = np.asarray(shape.contains_points(points)) + np.testing.assert_array_equal(inside, values <= 0) + assert shape.has_level_set and not shape.level_set_exact + for point in points[:, :5].T: + assert shape.contains(point) == bool(inside[np.argmax(points[0] == point[0])]) + + +@pytest.mark.parametrize("name", SHAPES) +def test_volume_and_bounds(name): + shape = SHAPES[name] + points, box = _box_points(shape, 300_000) + inside = np.asarray(shape.contains_points(points)) + assert np.mean(inside) * box == pytest.approx(shape.volume(), rel=0.03) + outside, _ = _box_points(shape, 20_000, pad=0.02) + beyond = np.any( + (outside < np.asarray(shape.lower_bound)[:, None]) + | (outside > np.asarray(shape.upper_bound)[:, None]), + axis=0, + ) + assert not np.asarray(shape.contains_points(outside[:, beyond])).any() + + +def test_frame_handedness_and_ends(): + coil = SHAPES["coil"] + first, last = np.asarray(coil.centreline([0.0, 2 * np.pi * coil.turns])).T + np.testing.assert_allclose(first, (0.14, 0, 0.02), atol=1e-15) + np.testing.assert_allclose(last, (0.14, 0, 0.14), atol=1e-12) + # right-handed: a quarter turn goes from +x to +y + np.testing.assert_allclose( + np.asarray(coil.centreline([np.pi / 2]))[:, 0], (0, 0.14, 0.025), atol=1e-12 + ) + assert not coil.left_handed and SHAPES["oblique left"].left_handed + assert coil.length == pytest.approx(0.12) + assert coil.wire_length == pytest.approx(6 * np.hypot(2 * np.pi * 0.14, 0.02)) + assert coil.surface_area() == pytest.approx( + 2 * np.pi * 0.003 * coil.wire_length + 4 * np.pi * 0.003**2 + ) + # hemispherical caps: just beyond the start, along the wire's tangent + assert coil.contains((0.14, -0.0029, 0.02)) + assert not coil.contains((0.14, -0.0031, 0.02)) + + +def test_copy_translate_and_mirror(): + shape = SHAPES["oblique left"] + points, _ = _box_points(shape, 2000, pad=0.02) + copied = shape.copy() + np.testing.assert_allclose(copied.level_set(points), shape.level_set(points)) + assert copied.left_handed and copied.start_direction == shape.start_direction + moved = shape.copy() + moved.translate((0.1, 0.0, -0.2)) + np.testing.assert_allclose( + moved.level_set(points + np.array([[0.1], [0.0], [-0.2]])), + shape.level_set(points), + atol=1e-12, + ) + mirrored = shape.copy() + mirrored.mirror(1, 0.05) + image = points.copy() + image[1] = 0.1 - image[1] + np.testing.assert_allclose( + mirrored.level_set(image), shape.level_set(points), atol=1e-12 + ) + assert not mirrored.left_handed + with pytest.raises(ValueError, match="Axis"): + mirrored.mirror(3, 0.0) + + +def test_crossings_of_the_wire(): + coil = SHAPES["coil"] + # radial segments through the first turn's wire at z = 0.02 + 0.02 * angle / 2 pi + angles = np.linspace(0.1, 2 * np.pi - 0.1, 50) + z = 0.02 + 0.02 * angles / (2 * np.pi) + start = np.stack([0.1 * np.cos(angles), 0.1 * np.sin(angles), z]) + end = np.stack([0.2 * np.cos(angles), 0.2 * np.sin(angles), z]) + result = crossings(coil, start, end) + assert np.all(np.asarray(result.hit)) + radius = 0.1 + 0.1 * np.asarray(result.fraction) + # the wire's inner side, where the plane z = const cuts it (to 1e-3 of the pitch) + np.testing.assert_allclose(radius, 0.14 - 0.003, atol=2e-5) + + +def test_yee_fractions_of_the_wire(): + coil = Helix((0, 0, -0.01), 0.03, 0.02, 1, 0.005) + grid = Grid(3, (-0.04, -0.04, -0.02), (0.04, 0.04, 0.02), (33, 33, 17)) + fractions = coil.yee_fractions(grid, samples=16) + h = 0.08 / 32 + # each family of edges, its inside lengths times h^2, is the wire's volume + for family in fractions.edges: + assert np.sum(family) * h**3 == pytest.approx(coil.volume(), rel=0.05) + # the projected traces times h over all three families: twice the surface's area + area = sum(np.sum(t) * h for t in fractions.projected_traces) + assert area == pytest.approx(2 * coil.surface_area(), rel=0.05) + + +def test_invalid_input(): + with pytest.raises(ValueError, match="wire_radius"): + Helix((0, 0, 0), 0.1, 0.02, 3, 0.2) + with pytest.raises(ValueError, match="pitch"): + Helix((0, 0, 0), 0.1, 0.0, 3, 0.01) + with pytest.raises(ValueError, match="nonzero"): + Helix((0, 0, 0), 0.1, 0.02, 3, 0.01, axis=(0, 0, 0)) + with pytest.raises(ValueError, match="start_direction"): + Helix((0, 0, 0), 0.1, 0.02, 3, 0.01, start_direction=(0, 0, 2)) diff --git a/src/simgeometry/tests/test_polar_repeat.py b/src/simgeometry/tests/test_polar_repeat.py new file mode 100644 index 0000000..7c0b5bc --- /dev/null +++ b/src/simgeometry/tests/test_polar_repeat.py @@ -0,0 +1,67 @@ +"""PolarRepeat: the copies of a shape about an axis, against their explicit union.""" + +from __future__ import annotations + +import numpy as np +import pytest + +from simgeometry import Circle, Cuboid, Grid, PolarRepeat, Sphere, Tube + +RNG = np.random.default_rng(3) + + +def _union(shape, count, axis=(0, 0, 1), origin=None): + union = shape + for k in range(1, count): + union = union | shape.rotated(360.0 * k / count, axis=axis, origin=origin) + return union + + +@pytest.mark.parametrize( + ("source", "count", "axis", "origin"), + [ + (Cuboid((0.12, 0, 0.08), (0.02, 0.002, 0.15)), 24, (0, 0, 1), None), + (Sphere((0.2, 0.0, 0.1), 0.05), 5, (0, 0, 1), (0.0, 0.0, 0.1)), + (Sphere((0.0, 0.3, 0.0), 0.1), 3, (1, 0, 0), (0.1, 0.0, 0.0)), + ], +) +def test_membership_and_level_set_are_those_of_the_union(source, count, axis, origin): + repeat = PolarRepeat(source, count, axis=axis, origin=origin) + union = _union(source, count, axis, origin) + lower = np.asarray(repeat.lower_bound) - 0.02 + upper = np.asarray(repeat.upper_bound) + 0.02 + points = RNG.uniform(lower[:, None], upper[:, None], (3, 50_000)) + np.testing.assert_array_equal( + np.asarray(repeat.contains_points(points)), + np.asarray(union.contains_points(points)), + ) + np.testing.assert_allclose( + repeat.level_set(points), union.level_set(points), atol=1e-12 + ) + assert repeat.has_level_set and not repeat.level_set_exact + inside = np.asarray(union.contains_points(points)) + assert np.all((points[:, inside] >= lower[:, None] + 0.02 - 1e-12).all(axis=0)) + assert np.all((points[:, inside] <= upper[:, None] - 0.02 + 1e-12).all(axis=0)) + assert repeat.contains(points[:, np.argmax(inside)]) + assert repeat.volume() == pytest.approx(count * source.volume()) + + +def test_slotted_tube_fractions_match_the_union(): + tube = Tube((0, 0, 0), (0, 0, 0.1), 0.03, 0.036) + slit = Cuboid((0.033, 0, 0.05), (0.02, 0.003, 0.08)) + grid = Grid(3, (-0.04, -0.04, 0.0), (0.04, 0.04, 0.1), (17, 17, 11)) + fast = (tube - PolarRepeat(slit, 8)).yee_fractions(grid, samples=16) + slow = (tube - _union(slit, 8)).yee_fractions(grid, samples=16) + for a in range(3): + np.testing.assert_allclose(fast.edges[a], slow.edges[a], atol=1e-9) + np.testing.assert_allclose(fast.faces[a], slow.faces[a], atol=1e-9) + np.testing.assert_allclose(fast.traces[a], slow.traces[a], atol=1e-9) + + +def test_invalid_input(): + with pytest.raises(ValueError, match="3D"): + PolarRepeat(Circle((1, 0), 0.1), 3) + with pytest.raises(ValueError, match="count"): + PolarRepeat(Sphere((1, 0, 0), 0.1), 0) + with pytest.raises(ValueError, match="axis"): + PolarRepeat(Sphere((1, 0, 0), 0.1), 3, axis=(0, 0, 0)) diff --git a/tutorials/shapes_and_node_masks.ipynb b/tutorials/shapes_and_node_masks.ipynb index 5de280a..5b39d56 100644 --- a/tutorials/shapes_and_node_masks.ipynb +++ b/tutorials/shapes_and_node_masks.ipynb @@ -244,7 +244,13 @@ "id": "15", "metadata": {}, "source": [ - "`node_mask()` works the same way on a 3D grid. `show_3d(points=...)` draws the nodes inside a tilted torus as a labelled point set over its translucent surface, in simgeometry's `(3, n)` layout as `node_points()` returns them; rotate the view to see that they fill the tube of the torus and leave its hole empty." + "## A coil and a slotted shield\n", + "\n", + "`Helix(start, radius, pitch, turns, wire_radius)` is a helical wire, such as the copper\n", + "tube of an RF coil: the points within `wire_radius` of a helix around `axis` (z by\n", + "default) that starts at `start + radius * x`. `PolarRepeat(shape, count)` makes `count`\n", + "copies of a shape about an axis and tests a point against one of them, so the 24 slits of\n", + "a Faraday shield cost no more than one." ] }, { @@ -253,6 +259,38 @@ "id": "16", "metadata": {}, "outputs": [], + "source": [ + "from simgeometry import Cuboid, Helix, PolarRepeat\n", + "\n", + "coil = Helix((0, 0, 0.025), 0.1385, 0.02, 6, 0.004) # six turns, an 8 mm tube\n", + "tube = Tube((0, 0, 0), (0, 0, 0.17), 0.1185, 0.1265)\n", + "shield = tube - PolarRepeat(Cuboid((0.1225, 0, 0.085), (0.024, 0.002, 0.15)), 24)\n", + "print(\"wire length:\", round(coil.wire_length, 3), \"m\")\n", + "points = (\n", + " np.random.default_rng(0).uniform(coil.lower_bound, coil.upper_bound, (200_000, 3)).T\n", + ")\n", + "box = np.prod(np.subtract(coil.upper_bound, coil.lower_bound))\n", + "print(\"volume:\", coil.volume(), \"sampled:\", coil.contains_points(points).mean() * box)\n", + "print(\n", + " \"in a slit, in the wall:\",\n", + " shield.contains_points(np.array([[0.1225, 0.1225], [0, 0.003], [0.085, 0.085]])),\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "17", + "metadata": {}, + "source": [ + "`node_mask()` works the same way on a 3D grid. `show_3d(points=...)` draws the nodes inside a tilted torus as a labelled point set over its translucent surface, in simgeometry's `(3, n)` layout as `node_points()` returns them; rotate the view to see that they fill the tube of the torus and leave its hole empty." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "18", + "metadata": {}, + "outputs": [], "source": [ "grid3 = Grid(3, (-1.0, -1.0, -1.0), (1.0, 1.0, 1.0), (21, 21, 21))\n", "torus = Torus((0, 0, 0), 0.6, 0.25, axis=(1, 0, 1))\n",