diff --git a/docs/architecture/lib-map.md b/docs/architecture/lib-map.md index d6aef129..e6081d02 100644 --- a/docs/architecture/lib-map.md +++ b/docs/architecture/lib-map.md @@ -55,6 +55,7 @@ rendering, slicing support, and verification: - `fit_camera_math.py` - `fit_check.py` - `assembly_cli.py` +- `assembly_pair_cli.py` - `fits.py` - `fits_apply.py` - `fits_registry.py` @@ -67,6 +68,7 @@ rendering, slicing support, and verification: - `gcode_resume.py` - `gcode_labels.py` - `gcode_overhangs.py` +- `gcode_walls.py` - `gcode_mix_gate.py` - `gcode_scan.py` - `stripe_contrast.py` @@ -284,12 +286,14 @@ Every current root Python module is listed here so future additions are intentio - `fit_check.py` — already-printed insert vs its socket: printed-size lock, per-layer clearance (manifold3d slices, shapely distances) and the open path along the insertion direction; backs `3d fit-check`. - `fits.py` — the one fit/tolerance table (`data/fits.yaml`, schema `3d.fits/1`): fit types, aliases, printer x material x nozzle x layer corrections, `resolve()` with provenance and `fit_clearance` calibration readings. - `assembly_cli.py` — the worker behind `3d assembly` (argument parsing, report assembly; geometry in `lib/assembly/`). +- `assembly_pair_cli.py` — the worker behind `3d assembly pair` (two-part fit check; geometry in `lib/assembly/pair.py`, `helical.py`, `printed.py`). - `parts_3mf.py` — reads the `3d:role` / `3d:printable` tags of a project 3MF and guards print boundaries (`arrange`, `retarget`, `prep`) with `lib/parts.py`. - `fits_apply.py` — pure text operations of `3d fits apply`: the OpenSCAD snippet of a resolved clearance and the idempotent rewrite of `// @fit`-marked declared clearances. - `fits_registry.py` — bridge from a free-text printer/material name (fleet key, model, `Elegoo Rapid PLA+`) to the registries `fits.resolve` cites in its provenance. - `gcode_resume.py` — resume a stopped multi-head Orca print above a layer (layer walk, machine state at a byte, resume header from the device's `resume:` profile, per-head retraction tracking, `--no-tower` off-bed purges); stdlib-only; backs `3d gcode-resume`. Its per-line parser (`advance`) and the streamed `state_at` (machine state at a printer's byte cursor or a layer) also give PrintOps incidents their G-code context. - `gcode_labels.py` — the slicer's own labels in a G-code: object labels in both dialects (Klipper `EXCLUDE_OBJECT_*`, Bambu Studio `; start printing object` / `M624`…`M625` blocks named by their `; printing object` comment, `M486`, PrusaSlicer comments), layer-change markers, feature tags and line widths; stdlib-only; the one parser behind head-sequence, overhangs, visible-faces and the print-area plate section. - `gcode_overhangs.py` — extrusion printed on nothing per object and layer of a sliced G-code (numpy/scipy raster of each layer's carrier, sparse reach/bridge/step-out tests on the pixels that stick out, layers split over forked workers); backs `3d gcode-overhangs` and the `overhangs` gate of `print-preflight`. +- `gcode_walls.py` — the outer-wall lines of one object, layer by layer, from a sliced G-code (object / layer labels via `gcode_labels`, width and path height per piece); backs the as-printed numbers of `3d assembly pair`. - `gcode_scan.py` — one streaming pass over a sliced G-code, layer by layer: per-layer Z / height / extrusion runs (type, tool, XY) handed to consumers, plus M600 pauses, tool changes, grams per tool and feature type and footer settings; the whole file is never held; numpy only. - `stripe_contrast.py` — perceptual layer-stripe contrast of a stack of single-colour layers at their real heights: eye-blurred profile vs its period mean, worst CIEDE2000 per layer (the `stripe contrast` check of `3d gcode-mix-gate`). - `gcode_mix_gate.py` — the multi-colour gates of `3d gcode-mix-gate`: colour volumes attributed from the part STLs (pure / mixed per layer on the layer's own outer-wall evidence / accent / thin), pause, heads and prime tower, pocket void, layer heights, footer settings, and the `--mix-teeth` / `--cavity-teeth` controls that must FAIL (a gate that cannot fail is vacuous); numpy/scipy/trimesh. diff --git a/docs/commands/README.md b/docs/commands/README.md index 91c1d17b..ed5b0e7d 100644 --- a/docs/commands/README.md +++ b/docs/commands/README.md @@ -31,7 +31,7 @@ explain why the command exists, what its inputs mean, and how it fits into the l | `convert` | [convert.md](convert.md) | Convert between formats (multicolour scad→Orca project [3MF](GLOSSARY.md#3mf) with per-object extruders, stl/usdz/png, retarget alias). | | `generate` | [generate.md](generate.md) | Text description + named dimensions -> verified parametric `.scad` via an AI backend (generate/validate/render/check/fix loop). | | `fit-niche` | [fit-niche.md](fit-niche.md) | Generate a parametric insert/plug that mates into a described cavity with FDM clearances. | -| `assembly` | [assembly.md](assembly.md) | Whole-assembly checks over declared mates: blocking intersections, per-mate fit clearance, insertion/removal sweeps on any axis, assembly order, finger/tool ergonomics, explode, BOM, 3MF positioning (`check`, `sweep`, `order`, `ergonomics`, `explode`, `bom`, `position`, `show`). | +| `assembly` | [assembly.md](assembly.md) | Whole-assembly checks over declared mates: blocking intersections, per-mate fit clearance, insertion/removal sweeps on any axis, assembly order, finger/tool ergonomics, explode, BOM, 3MF positioning (`check`, `sweep`, `order`, `ergonomics`, `explode`, `bom`, `position`, `show`), plus `pair`: two meshes, interpenetration, minimum gap and per-feature radial clearance of threads and ribs, model and as printed from G-code. | | `fits` | [fits.md](fits.md) | Fit tolerances per fit type (snap/sliding/tight/press/glue) corrected for printer x material: `list`, `show`, `resolve` (with provenance), `apply` (snippet / correct declared clearances), `corrections`. | | `library` | [library.md](library.md) | Local model-library database (searches, sources, designs, makes, photos, publications): `path`, `migrate`, `status`, `export`, offline; read-only `makes`, `publications`, `relations`, `graph DESIGN` (source <- remix <- make). | | `make` | [make.md](make.md) | Record the prints (makes) of a design in the library: `add`, `status` (queued -> printing -> done/failed), `list`; local database only. | diff --git a/docs/commands/assembly.md b/docs/commands/assembly.md index 77d71c12..560a2506 100644 --- a/docs/commands/assembly.md +++ b/docs/commands/assembly.md @@ -7,7 +7,7 @@ hooks and screws, and where each printed part sits on the plates. Umbrella #570, #553 (declared insertion axis) and covers the assembly-motion part of #61. Every subcommand takes `--json` (schema `3d.assembly/1`, `kind` = `check | sweep | sweeps | order | ergonomics | explode | -bom | position | show`). Exit codes: `0` ok, `1` a verdict failed (blocking intersection, blocked path, no access, no +bom | position | show | pair`). Exit codes: `0` ok, `1` a verdict failed (blocking intersection, blocked path, no access, no straight order), `2` bad spec / usage, `127` missing dependency. ## Usage @@ -23,6 +23,13 @@ straight order), `2` bad spec / usage, `127` missing dependency. [--process NAME] [--filament NAME] [--preset-dir DIR]... [--set K=V]... [--include-non-printable] [--reason TEXT] [--yes] [--json] 3d assembly show TARGET [--json] +3d assembly pair A B [--scale-a S] [--scale-b S] [--pose-a POSE] [--pose-b POSE] [--axis AXIS] [--axis-point X,Y,Z] + [--a-side inner|outer --b-side inner|outer [--pitch MM] [--hand auto|right|left] [--a-z LO:HI] [--b-z LO:HI] + [--a-r LO:HI] [--b-r LO:HI] [--feature-name N] | --features FILE] + [--min-radial-gap MM | --fit TYPE [--shape S] [--printer P] [--material M] [--nozzle MM] [--layer MM]] + [--min-clearance MM] [--zone XMIN,YMIN,ZMIN,XMAX,YMAX,ZMAX] [--max-overlap-mm3 V] [--seat-tol MM] + [--screw-to-best] [--gcode-a F --object-a NAME --print-a POSE] [--gcode-b F --object-b NAME --print-b POSE] + [--skip-below MM] [--max-faces N] [--json] ``` `TARGET` is an assembly spec (`assembly.yaml` / `.json`), a project directory holding `assembly.yaml` next to its @@ -38,6 +45,7 @@ straight order), `2` bad spec / usage, `127` missing dependency. 3d assembly position assembly.yaml --3mf out/plates.3mf --base-project print/mole.3mf --yes 3d assembly position assembly.yaml --3mf out/plates.3mf --include-non-printable --reason "fit test" --yes 3d assembly show assembly.yaml +3d assembly pair rim.stl plug.stl --a-side inner --b-side outer --pitch 3.007 --screw-to-best ``` ## The spec (`3d.assembly/1`) @@ -85,6 +93,71 @@ loads: # carried in the report for the S4 solver * **Roles**: `printed | vitamin | reference | stand_in` (`lib/parts.py`); a vitamin or reference model must name its source URL and licence, a stand-in is labelled as invented in the BOM and drawn ghosted. +## `pair` - do two parts mate? (no spec) + +`3d assembly check` judges declared mates of a spec. `pair` takes two meshes (STL / 3MF / OBJ / PLY), each with an optional +scale and pose into the assembled frame, and answers the question for one pair, down to a thread or a rib. A 3MF is read with the +stdlib reader (no lxml); one that holds several objects needs `plate.3mf#ID` (the refusal lists the ids), and its build-item +transform is part of the mesh: + +```bash +# rim (the ball's bore) and plug at the print scale, plug seated under the rim; thread pitch 3.007 mm; model only +3d assembly pair rim.stl plug.stl --scale-a 2.00493695 --scale-b 2.00493695 --pose-b '{translate: [0, 0, -37.782]}' \ + --a-side inner --b-side outer --pitch 3.007 --a-z -36.78:-28.28 --a-r 0:36.4 --b-z -36.78:-30.98 --b-r 0:36 --screw-to-best +# the same, as printed: the plug from its K1 G-code, the new rim from its slice; both placed on the bed automatically +3d assembly pair rim.stl plug.stl ... --gcode-a rim.gcode --object-a "rim.stl" --gcode-b plate.gcode --object-b "Plug x2.stl" +# several features from a file, the clearance the fit table asks for on this printer and plastic +3d assembly pair a.stl b.stl --features features.yaml --fit sliding --printer k1se --material PETG +``` + +* **interpenetration** - exact overlap regions (volume, depth, xyz as in `check`); a region thinner than `--seat-tol` + (0.02 mm) is a seat contact, not a collision. FAILs above `--max-overlap-mm3` (default 0.01): a designed press fit is allowed + only by saying how much. +* **clearance** - the exact minimum gap of the two solids (0 where parts seat on a face); `--zone BOX` measures only the surface + samples inside a box (spacing reported), which leaves a seat face out; `--min-clearance` makes it a verdict. +* **features** - the part of the question a volume cannot answer. A feature is the surface of each part around a COMMON AXIS + (`--axis`, `--axis-point`; z windows are the coordinate along it, r windows the radius from it): `--a-side inner` is the wall + of a bore, `outer` a boss; one of each. The surface radius is read against the PHASE `u` (a thread: `z/pitch -+ theta/2pi`, one + value per `u` for a perfect single-start thread; ribs and keys, `--pitch 0`: the angle), 360 bins, from the model's sections + every 0.12 mm. The bore minus the boss is the radial gap, at the posed phase and at every other relative phase - which is + screwing (a turn = an axial advance of one pitch) or clocking. Reported: the gap at the pose, the best phase and the tightest + radial gap there, the free play (mm of axial travel of a thread, degrees for ribs), male crest -> female root and female crest -> + male root clearances, handedness (detected, a mismatch FAILs: the threads cannot screw together), pitch coherence (a wrong pitch + shows as a large value and a warning), coverage of the phase bins (unmeasured = FAIL under 50 %), and the concentricity of each + feature's own axis with the declared one. FAILs when the best-phase gap is under `--min-radial-gap` (default 0.1) or the + `--fit TYPE` row (`clearance - tol` per face of `3d fits resolve`, printer x material corrected). The posed phase may interfere + while a free phase exists: that is a WARN with the move that fixes it, in turns, degrees and mm of advance. + A bore with no `--a-r` is read from the axis out to the boss's largest radius + 3 mm (else its own outer wall would count). +* **`--screw-to-best`** turns B about the axis to the best phase of the first feature before anything else is measured (a screwed + fit has no unique pose: the phase is free). The turn keeps B's axial position, so a part that seats on a face stays seated. +* **as printed** - `--gcode-a/--gcode-b` (+ `--object-a/b` when the plate holds several objects; the slicer's own object labels - a wrong or + ambiguous name, a name that matches several copies, or none on a multi-object plate lists the file's labels; `"Cup.stl id:0 + copy 1"` picks one copy) read the OUTER WALL paths of the object (`lib/gcode_walls.py`: Orca / Snapmaker Orca / Bambu / + PrusaSlicer, arcs, relative extrusion). The model is put on the bed (`--print-a POSE`: orientation like a pose; without a + `translate` it is dropped on z = 0 and REGISTERED on the walls in XY by a translation-only fit - the slicer's own placement), the + per-layer deviation of the printed wall from the model's section at the layer mid height is reported (the wall centre line is + half a line width inside the surface; layers under `--skip-below` are listed apart: first-layer flare), overall and inside the + feature's window, and the feature clearance is recomputed with the printed surface (`centre line +- width/2` radially). A part + without a G-code is taken from its model, the report says which side is printed. The slicer's compensation + (`xy_hole_compensation`, elephant-foot, line width) therefore shows in the numbers, which a model-only check cannot see. + +`features.yaml` (several features, the single-feature flags are the same keys): + +```yaml +features: + - name: rim_thread + pitch: 3.007 # 0 / absent = ribs or keys + hand: auto # auto | right | left + min_gap_mm: 0.2 # per-feature requirement (else --min-radial-gap / --fit / 0.1) + a: {side: inner, z: [-36.78, -28.28], r: [0, 36.4]} + b: {side: outer, z: [-36.78, -30.98], r: [0, 36.0]} +``` + +A value that starts with `-` may be given after `=` (`--a-z=-36:-28`) or as the next argument (`--a-z -36:-28`). +Exit codes as the other subcommands: 0 ok, 1 a verdict failed, 2 bad input (an unknown object name, a window without surface, two +bores), 127 missing dependency. `--json` document: `kind: pair` with `parts`, `axis`, `interpenetration`, `clearance`, `features` +(each with `model` and, with a G-code, `as_printed`, `printed_vs_model`), `printed`, `fit`, `summary`. + ## What each check does * **check** - per pair of parts, exact overlap regions (manifold3d, coarse to fine: boxes, then a crop to the overlap diff --git a/lib/assembly/helical.py b/lib/assembly/helical.py new file mode 100644 index 00000000..fb1ea8bf --- /dev/null +++ b/lib/assembly/helical.py @@ -0,0 +1,316 @@ +"""helical.py — radial clearance of two coaxial features (threads, ribs, keys, bores) from their surface profiles. + +WHAT: a feature is the surface of one part around a common axis (the wall of a bore, the flank of a boss). Its PROFILE is + the surface radius as a function of the phase `u` (one turn = 1): + thread (pitch p > 0) u = z/p -+ theta/2pi - constant along a helix, so a perfect single-start thread is one + number per u; the sign is the handedness (right: u = z/p - theta/2pi), + ribs / keys (p = 0) u = theta/2pi - the radius against the angle, whatever the height. + The female profile (the bore, material outside the surface) and the male profile (the boss, material inside) give + the radial gap `rf(u) - rm(u)` - at the posed phase, and at every other relative phase `d`, which is exactly the + screwing (rotation = axial advance for a thread; rotation for ribs). The best phase, how much of the turn is free + (axial play = fraction * pitch) and the crest / root clearances are read off that. +WHY: interpenetration volume and the minimum 3D gap say a thread BINDS; they do not say which feature of it, nor + whether it screws at all (wrong pitch or handedness, a root fuller than the crest of the partner). The per-phase + radial clearance does, and it works the same on the model and on the printed wall paths (`printed.py`). +INVARIANTS: the radius of a bin is the MEDIAN of the points in it (robust against the first-layer flare and stray + segments); `coherence_mm` is the 90th percentile of the distance of the points to their bin's median - a wrong pitch + or handedness shows as a large value; empty bins are NaN and ignored by the gap search, and `coverage` says how + many there are. All coordinates are in the AXIS frame (axis = +z, origin on the axis). +""" + +from __future__ import annotations + +import math +from dataclasses import dataclass, field +from typing import Any + +import numpy as np +import numpy.typing as npt + +from gcode_walls import densify_segments + +BINS = 360 +MIN_POINTS_PER_BIN = 3 +TIE_MM = 5e-4 # phases whose gap is within this of the best are equally good +HAND_DECISIVE_RATIO = 0.8 # the wrong handedness must be this much less coherent for the pick to count as decided + +Array = npt.NDArray[np.float64] +Matrix = npt.NDArray[np.float64] + + +def axis_frame(origin: Array, direction: Array) -> Matrix: + """4x4 matrix taking assembly coordinates to the axis frame: `origin` -> 0, `direction` -> +z.""" + z = np.asarray(direction, dtype=float) + z = z / np.linalg.norm(z) + helper = np.array([1.0, 0.0, 0.0]) if abs(z[0]) < 0.9 else np.array([0.0, 1.0, 0.0]) + x = helper - z * float(helper @ z) + x /= np.linalg.norm(x) + y = np.cross(z, x) + rot = np.vstack([x, y, z]) + out = np.eye(4) + out[:3, :3] = rot + out[:3, 3] = -rot @ np.asarray(origin, dtype=float) + return out + + +def apply(matrix: Matrix, pts: Array) -> Array: + return pts @ matrix[:3, :3].T + matrix[:3, 3] + + +def cylindrical(pts: Array) -> tuple[Array, Array, Array]: + """(r, theta, z) of axis-frame points.""" + return np.hypot(pts[:, 0], pts[:, 1]), np.arctan2(pts[:, 1], pts[:, 0]), pts[:, 2] + + +def window( + pts: Array, z_range: tuple[float, float] | None, r_range: tuple[float, float] | None +) -> Array: + """The axis-frame points inside a z window and a radius window (None = unbounded).""" + keep = np.ones(len(pts), dtype=bool) + if z_range is not None: + keep &= (pts[:, 2] >= z_range[0]) & (pts[:, 2] <= z_range[1]) + if r_range is not None: + r = np.hypot(pts[:, 0], pts[:, 1]) + keep &= (r >= r_range[0]) & (r <= r_range[1]) + return pts[keep] + + +def phase(theta: Array, z: Array, pitch: float, hand: int) -> Array: + """The phase u in [0, 1): `hand` +1 right-handed, -1 left-handed; pitch 0 -> the angle alone.""" + turns = theta / (2.0 * math.pi) + if pitch <= 0.0: + return np.asarray(turns % 1.0) + return np.asarray((z / pitch - hand * turns) % 1.0) + + +@dataclass +class Profile: + """Radius against phase for one feature (`radius[k]` for bin k of BINS; NaN = no data).""" + + radius: Array + counts: Array + coherence_mm: float + pitch_mm: float + hand: int # +1 right, -1 left, 0 none (ribs) + + @property + def coverage(self) -> float: + return float(np.mean(np.isfinite(self.radius))) + + @property + def r_min(self) -> float: + return float(np.nanmin(self.radius)) + + @property + def r_max(self) -> float: + return float(np.nanmax(self.radius)) + + +ENVELOPE_PERCENTILE = {"outer": 90.0, "inner": 10.0} # the blocking side of a bin: a boss's largest, a bore's smallest radius + + +def build_profile( + pts: Array, + pitch: float, + hand: int, + bins: int = BINS, + offset: Array | None = None, + side: str | None = None, +) -> Profile: + """Radius per phase bin of axis-frame points; `offset` is added to each point's radius (half line width). + + A bin of a helical surface holds one radius (the points of a perfect thread agree to a few um; the spread left is + mesh faceting and layer noise), so a thread takes the MEDIAN. Ribs and keys (pitch 0) have steep walls, where one + angle bin spans a range of radii: there the envelope that blocks is kept - the 90th percentile of a boss + (`side='outer'`), the 10th of a bore (`'inner'`).""" + r, theta, z = cylindrical(pts) + if offset is not None: + r = r + offset + u = phase(theta, z, pitch, hand) + b = np.minimum((u * bins).astype(int), bins - 1) + order = np.argsort(b, kind="stable") + bs, rs = b[order], r[order] + cuts = np.searchsorted(bs, np.arange(bins + 1)) + radius = np.full(bins, np.nan) + counts = np.zeros(bins) + spread: list[Array] = [] + for k in range(bins): + seg = rs[cuts[k] : cuts[k + 1]] + counts[k] = len(seg) + if len(seg) >= MIN_POINTS_PER_BIN: + med = float(np.median(seg)) + radius[k] = ( + med + if pitch > 0 or side not in ENVELOPE_PERCENTILE + else float(np.percentile(seg, ENVELOPE_PERCENTILE[side])) + ) + spread.append(np.abs(seg - med)) + coherence = float(np.percentile(np.concatenate(spread), 90)) if spread else math.nan + return Profile(radius, counts, coherence, pitch, hand) + + +def pick_handedness( + pts: Array, + pitch: float, + requested: str, + offset: Array | None = None, + side: str | None = None, +) -> tuple[Profile, dict[str, Any]]: + """The profile of the handedness that explains the points best (or the requested one), with how decisive it was.""" + if pitch <= 0.0: + prof = build_profile(pts, 0.0, 0, offset=offset, side=side) + return prof, {"hand": "none", "decided": True, "coherence_mm": {}} + right = build_profile(pts, pitch, +1, offset=offset, side=side) + left = build_profile(pts, pitch, -1, offset=offset, side=side) + spreads = { + "right": round(right.coherence_mm, 4), + "left": round(left.coherence_mm, 4), + } + best = right if right.coherence_mm <= left.coherence_mm else left + other = left if best is right else right + decided = ( + np.isfinite(best.coherence_mm) + and np.isfinite(other.coherence_mm) + and best.coherence_mm < HAND_DECISIVE_RATIO * other.coherence_mm + ) + if requested in ("right", "left"): + chosen = right if requested == "right" else left + else: + chosen = best + info = { + "hand": "right" if chosen.hand > 0 else "left", + "decided": bool(decided), + "coherence_mm": spreads, + "best_fit": "right" if best.hand > 0 else "left", + } + return chosen, info + + +def fit_axis_offset(pts: Array) -> tuple[float, float, float]: + """Least-squares circle through the points' xy: (dx, dy, residual rms) - how far the feature's own axis is from + the declared one (a feature that is not concentric with the axis it is judged on cannot be read radially).""" + x, y = pts[:, 0], pts[:, 1] + a = np.column_stack([2 * x, 2 * y, np.ones(len(x))]) + sol, *_ = np.linalg.lstsq(a, x * x + y * y, rcond=None) + cx, cy = float(sol[0]), float(sol[1]) + r = math.sqrt(max(sol[2] + cx * cx + cy * cy, 0.0)) + rms = float(np.sqrt(np.mean((np.hypot(x - cx, y - cy) - r) ** 2))) + return cx, cy, rms + + +@dataclass +class Gap: + """The radial gap of the female over the male profile at one relative phase.""" + + min_mm: float + mean_mm: float + at_bin: int + + +def gap_at(female: Array, male: Array, shift: int = 0) -> Gap: + """Radial gap female - male with the male turned by `shift` bins; NaN bins on either side are skipped.""" + g = female - np.roll(male, shift) + ok = np.isfinite(g) + if not ok.any(): + return Gap(math.nan, math.nan, -1) + k = int(np.nanargmin(g)) + return Gap(float(g[k]), float(np.nanmean(g)), k) + + +def sliding(female: Array, male: Array) -> Array: + """The minimum radial gap at every relative phase (index = shift in bins).""" + n = len(female) + return np.array([gap_at(female, male, d).min_mm for d in range(n)]) + + +@dataclass +class Clearance: + """The per-feature radial clearance of a female/male profile pair.""" + + pitch_mm: float + at_pose: Gap + best_shift_bins: int + best_gap_mm: float + free_fraction: float + crest_male_to_root_female_mm: float + crest_female_to_root_male_mm: float + by_shift: Array = field(repr=False, default_factory=lambda: np.zeros(0)) + + @property + def free_play(self) -> tuple[float, str]: + """The free travel at the best position: mm along the axis for a thread, degrees for ribs.""" + if self.pitch_mm > 0: + return self.free_fraction * self.pitch_mm, "mm axial" + return self.free_fraction * 360.0, "deg" + + +def clearance(female: Profile, male: Profile) -> Clearance: + by = sliding(female.radius, male.radius) + if np.isfinite(by).any(): + best_gap = float(np.nanmax(by)) + ties = np.flatnonzero(by >= best_gap - TIE_MM) # a symmetric part is best at several phases: take the nearest + best = int(ties[np.argmin(np.minimum(ties, len(by) - ties))]) + free = float(np.nansum(by > 0) / np.isfinite(by).sum()) + else: + best, best_gap, free = 0, math.nan, math.nan + return Clearance( + female.pitch_mm, + gap_at(female.radius, male.radius, 0), + best, + best_gap, + free, + float(np.nanmax(female.radius) - np.nanmax(male.radius)), + float(np.nanmin(female.radius) - np.nanmin(male.radius)), + by, + ) + + +def shift_to_motion(shift_bins: int, pitch: float, hand: int, bins: int = BINS) -> dict[str, float]: + """What moving the MALE part by `shift_bins` bins of phase means, the shorter way round. + + A thread: an axial advance `axial_mm` (u = z/p - hand*theta/2pi: +dz moves a feature by +dz/p turns) or, with the + axial position kept (a part that seats on a face), the turn `rotate_deg` that has the same effect (a turn by phi + moves a feature by -hand*phi/360 turns). Ribs: only the turn (u = theta/2pi). The sign follows `roll(male, shift)`: + bin k of the rolled profile is bin k - shift of the male, i.e. the male profile moves toward higher u.""" + turns = shift_bins / bins + if turns > 0.5: + turns -= 1.0 + if pitch > 0: + return { + "rotate_deg": round(-hand * turns * 360.0, 3) + 0.0, + "axial_mm": round(turns * pitch, 4) + 0.0, + "turns": round(turns, 4) + 0.0, + } + return {"rotate_deg": round(turns * 360.0, 3) + 0.0, "turns": round(turns, 4) + 0.0} + + +def mesh_wall_points( + mesh_axis_frame: Any, + z_range: tuple[float, float], + r_range: tuple[float, float] | None, + step_mm: float, + spacing_mm: float, +) -> Array: + """Surface points of a mesh already in the axis frame: its sections every `step_mm` along z, densified, windowed.""" + from trimesh import intersections + + z0, z1 = z_range + heights = np.arange(z0 + step_mm / 2.0, z1, step_mm) + if len(heights) == 0: + heights = np.array([(z0 + z1) / 2.0]) + lines, to_3d, _ = intersections.mesh_multiplane( + mesh_axis_frame, np.zeros(3), np.array([0.0, 0.0, 1.0]), heights + ) + out: list[Array] = [] + for seg, mat in zip(lines, to_3d, strict=True): + if len(seg) == 0: + continue + seg = np.asarray(seg, dtype=float) + pts2, _ = densify_segments(seg[:, 0, :], seg[:, 1, :], spacing_mm) + pts3 = ( + np.column_stack([pts2, np.zeros(len(pts2)), np.ones(len(pts2))]) + @ np.asarray(mat).T + ) + out.append(pts3[:, :3]) + pts = np.vstack(out) if out else np.zeros((0, 3)) + return window(pts, None, r_range) diff --git a/lib/assembly/pair.py b/lib/assembly/pair.py new file mode 100644 index 00000000..6ec82857 --- /dev/null +++ b/lib/assembly/pair.py @@ -0,0 +1,945 @@ +"""pair.py — `3d assembly pair A B`: do two parts mate? Interpenetration, minimum clearance, per-feature radial clearance. + +WHAT: two meshes (STL / 3MF / OBJ / PLY), each with an optional scale and pose, are put in the assembled frame and judged: + interpenetration exact overlap volume, regions (volume, depth, xyz) - `overlap.pair_regions`; + clearance exact minimum gap between the solids - `overlap.min_gap` (0 when they overlap); + features per coaxial feature (thread, rib, key, bore), the RADIAL clearance read from surface profiles + (`helical`): at the posed phase and at the best phase of screwing / clocking, the free play, crest + and root clearances, handedness and pitch coherence; + as printed with `--gcode-a/--gcode-b` the same feature clearance from the printed outer-wall paths + (`printed`), the model-to-print deviation per layer, and the registration residual. + `--screw-to-best` turns part B about the axis to the best phase of the first feature before everything else is + measured: the pose of a screwed fit is not unique, the phase is the free parameter (the turn keeps the axial position, + so a seat face stays seated; the equivalent axial advance of a thread is reported next to it). +WHY: `3d assembly check` needs an assembly spec of declared mates; the question "do these two parts screw together on + this plastic" is a pair question with a feature in it. A volume and a minimum gap cannot say which flank binds nor + whether the printed walls (compensation, line width) close the gap the model has. +INVARIANTS: nothing is repaired or simplified; a feature that is not measured (empty bins, wrong pitch / handedness) + FAILS rather than passing; the verdicts (`pass | warn | fail`) fold into `ok`; poses and windows are in the assembled + frame, the axis is a line (point + direction) in it, feature z windows are the coordinate along that axis from the + axis point's plane through the world origin (for the default +Z axis: world z). +""" + +from __future__ import annotations + +import math +from dataclasses import dataclass, field, replace +from pathlib import Path +from typing import Any + +import numpy as np +import numpy.typing as npt + +from assembly import helical, printed +from assembly.overlap import min_gap, pair_regions, sampled_gap +from assembly.scene import ( + Box, + box_intersection, + load_mesh_file, + pose_matrix, + solid_box, + solid_from_mesh, + translation_matrix, +) +from assembly.spec import SCHEMA, Pose +from errors import InputNotFound, InvalidArgument + +CMD = "assembly" +DEFAULT_MIN_RADIAL_GAP_MM = 0.1 +SEARCH_MM = 5.0 +EPS_MM3 = 0.01 +MIN_COVERAGE_FAIL = 0.5 +MIN_COVERAGE_WARN = 0.9 +AXIS_OFFSET_WARN_MM = 0.1 +SLICE_STEP_MM = 0.12 + +Array = npt.NDArray[np.float64] +Matrix = npt.NDArray[np.float64] +Range = tuple[float, float] + + +@dataclass(frozen=True) +class FeatureSide: + """Which surface of a part a feature is: `inner` = the wall of a bore (material outside), `outer` = a boss.""" + + side: str + z: Range | None = None + r: Range | None = None + + +@dataclass(frozen=True) +class Feature: + name: str + a: FeatureSide + b: FeatureSide + pitch_mm: float = 0.0 + hand: str = "auto" + min_gap_mm: float | None = None + + +@dataclass +class PartInput: + path: Path + scale: float = 1.0 + pose: Pose = () + print_pose: Pose = () + gcode: Path | None = None + gcode_object: str | None = None + object_id: int | None = None # which object of a 3MF (`plate.3mf#3`) + + +@dataclass +class PairOptions: + a: PartInput + b: PartInput + features: list[Feature] = field(default_factory=list) + axis_point: tuple[float, float, float] | None = None + axis: tuple[float, float, float] = (0.0, 0.0, 1.0) + min_radial_gap_mm: float | None = None + min_clearance_mm: float | None = None + max_overlap_mm3: float = EPS_MM3 + seat_tol_mm: float = 0.02 + zone: Box | None = None + screw_to_best: bool = False + skip_below_mm: float = 0.5 + max_faces: int = 400_000 + fit_resolution: dict[str, Any] | None = None + + +def parse_range(text: str, where: str) -> Range: + """`LO:HI` -> (lo, hi).""" + try: + lo, hi = (float(x) for x in text.split(":")) + except ValueError as exc: + raise InvalidArgument( + where, text, ["LO:HI, e.g. 1.0:9.5"], command=CMD + ) from exc + if hi < lo: + raise InvalidArgument(where, text, ["LO <= HI"], command=CMD) + return (lo, hi) + + +def _unit(v: tuple[float, float, float]) -> Array: + a = np.asarray(v, dtype=float) + n = float(np.linalg.norm(a)) + if n == 0.0: + raise InvalidArgument("--axis", str(v), ["a non-zero direction"], command=CMD) + return a / n + + +def split_part_path(text: str) -> tuple[Path, int | None]: + """`plate.3mf#3` -> (plate.3mf, 3): one object of a 3MF by its id; any other text is a plain path.""" + base, sep, tail = text.rpartition("#") + if sep and tail.isdigit() and base.lower().endswith(".3mf"): + return Path(base), int(tail) + return Path(text), None + + +def _read_mesh(part: PartInput, label: str) -> Any: + """STL / OBJ / PLY through trimesh; a 3MF through the stdlib reader (no lxml): its single object, or `#ID`'s.""" + if part.path.suffix.lower() != ".3mf": + return load_mesh_file(part.path, label) + import zipfile + from xml.etree import ElementTree as ET + + from assembly import threemf + + try: + objects = threemf.object_meshes(part.path) + except (zipfile.BadZipFile, KeyError, ET.ParseError) as exc: + raise InvalidArgument( + f"{label}.mesh", str(part.path), ["a readable 3MF (3D/3dmodel.model)"], command=CMD, extra=str(exc) + ) from exc + if part.object_id is not None: + if part.object_id not in objects: + raise InvalidArgument( + f"{label}.mesh", + f"{part.path}#{part.object_id}", + [f"an object id of the file: {', '.join(str(i) for i in sorted(objects))}"], + command=CMD, + ) + return objects[part.object_id] + if len(objects) != 1: + listing = [ + f"#{oid}: {len(m.faces)} faces, box {[round(float(v), 2) for v in m.extents]} mm" + for oid, m in sorted(objects.items()) + ] + raise InvalidArgument( + f"{label}.mesh", + f"{part.path} holds {len(objects)} objects", + [f"pick one with `{part.path}#ID`", *listing], + command=CMD, + ) + return next(iter(objects.values())) + + +def _load(part: PartInput, label: str, max_faces: int) -> tuple[Any, Matrix]: + """(mesh scaled, matrix mesh-frame -> assembled frame).""" + if not part.path.is_file(): + raise InputNotFound(str(part.path), command=CMD) + mesh = _read_mesh(part, label) + if len(mesh.faces) > max_faces: + raise InvalidArgument( + f"{label}.mesh", + f"{len(mesh.faces)} faces", + [f"at most {max_faces} faces (--max-faces)"], + command=CMD, + ) + if part.scale != 1.0: + mesh.apply_scale(part.scale) + return mesh, pose_matrix(part.pose) + + +def _bbox(mesh: Any) -> dict[str, list[float]]: + return { + "min": [round(float(v), 4) for v in mesh.bounds[0]], + "max": [round(float(v), 4) for v in mesh.bounds[1]], + } + + +def interpenetration(solid_a: Any, solid_b: Any, max_mm3: float, seat_tol_mm: float) -> dict[str, Any]: + """Overlap regions; a region thinner than `seat_tol_mm` is a CONTACT (a seat face, #140), not interpenetration.""" + regions = pair_regions("a", solid_a, "b", solid_b, eps_mm3=EPS_MM3, keep_mesh=False) + real = [r for r in regions if r.depth_mm >= seat_tol_mm] + contact = [r for r in regions if r.depth_mm < seat_tol_mm] + total = sum(r.volume_mm3 for r in real) + return { + "status": "pass" if total <= max_mm3 else "fail", + "volume_mm3": round(total, 4), + "allowed_mm3": max_mm3, + "count": len(real), + "regions": [r.to_dict() for r in real[:10]], + "contacts": { + "count": len(contact), + "volume_mm3": round(sum(r.volume_mm3 for r in contact), 4), + "max_depth_mm": round(max((r.depth_mm for r in contact), default=0.0), 4), + "seat_tol_mm": seat_tol_mm, + }, + } + + +def clearance( + solid_a: Any, solid_b: Any, mesh_a: Any, mesh_b: Any, required: float | None, zone: Box | None +) -> dict[str, Any]: + window = box_intersection(solid_box(solid_a), solid_box(solid_b), SEARCH_MM) + gap = min_gap(solid_a, solid_b, SEARCH_MM, window) + out: dict[str, Any] = { + "min_gap_mm": round(gap, 4), + "capped_at_mm": SEARCH_MM if gap >= SEARCH_MM else None, + "required_mm": required, + } + status = "pass" + if zone is not None: + zgap, spacing = sampled_gap(mesh_a, mesh_b, zone, SEARCH_MM) + out["zone"] = { + "box": [round(float(v), 4) for v in zone], + "min_gap_mm": round(zgap, 4), + "sample_spacing_mm": round(spacing, 4), + } + gap = zgap + if required is not None and gap < required: + status = "fail" + out["status"] = status + return out + + +@dataclass +class Axis: + origin: Array + direction: Array + matrix: Matrix # assembled -> axis frame + + +def make_axis(opts: PairOptions, mesh_a: Any) -> Axis: + d = _unit(opts.axis) + if opts.axis_point is not None: + o = np.asarray(opts.axis_point, dtype=float) + o = o - float(o @ d) * d + else: + c = np.asarray(mesh_a.bounds, dtype=float).mean(axis=0) + o = c - float(c @ d) * d + return Axis(o, d, helical.axis_frame(o, d)) + + +@dataclass +class PrintedPart: + walls: Any + bed_to_assembly: Matrix + info: dict[str, Any] + bed_mesh: Any = field(repr=False, default=None) # the model on the bed, before the registration shift + shift: tuple[float, float] = (0.0, 0.0) # added to the wall points to lay them on `bed_mesh` + skip_below_mm: float = 0.5 + + +def _drop_matrix(mesh: Any, matrix: Matrix) -> Matrix: + """`matrix` followed by the translation that puts the mesh's lowest point on the bed (z = 0).""" + corners = np.array( + [ + [x, y, z] + for x in mesh.bounds[:, 0] + for y in mesh.bounds[:, 1] + for z in mesh.bounds[:, 2] + ] + ) + pts = helical.apply(matrix, corners) + return translation_matrix([0.0, 0.0, -float(pts[:, 2].min())]) @ matrix + + +def load_printed( + part: PartInput, label: str, mesh: Any, pose: Matrix, skip_below: float +) -> PrintedPart: + from gcode_walls import read_walls + + assert part.gcode is not None + if not part.gcode.is_file(): + raise InputNotFound(str(part.gcode), command=CMD) + try: + walls = read_walls(part.gcode, part.gcode_object) + except (KeyError, ValueError) as exc: + raise InvalidArgument( + f"--object-{label}", + str(part.gcode_object), + [str(exc.args[0])], + command=CMD, + ) from exc + explicit = any(kind == "translate" for kind, _ in part.print_pose) + base = pose_matrix(part.print_pose) + p0 = base if explicit else _drop_matrix(mesh, base) + if not explicit and walls.base_z > 0.05: # a raft: the object's z = 0 sits above the bed + p0 = translation_matrix([0.0, 0.0, walls.base_z]) @ p0 + bed_mesh = mesh.copy() + bed_mesh.apply_transform(p0) + shift = (0.0, 0.0) + resid = math.nan + if not explicit: + shift, resid = printed.register_xy(bed_mesh, walls, skip_below) + dev = printed.deviation(bed_mesh, printed.layers_in_window(walls, np.eye(4), None), skip_below, shift) + placed = translation_matrix([-shift[0], -shift[1], 0.0]) @ p0 + info = { + "gcode": str(part.gcode), + "object": part.gcode_object, + "objects_in_file": walls.objects, + "layers": len(walls.layers), + "wall_mm": round(walls.wall_mm, 1), + "overhang_wall_left_out_mm": round(walls.overhang_mm, 1), + "z_lift_mm": round(walls.base_z, 3), + "placement": { + "registered": not explicit, + "shift_applied_mm": [round(-shift[0], 4), round(-shift[1], 4)], + "residual_rms_mm": None if math.isnan(resid) else round(resid, 4), + }, + "deviation": dev, + } + return PrintedPart(walls, pose @ np.linalg.inv(placed), info, bed_mesh, shift, skip_below) + + +@dataclass +class SideProfile: + profile: helical.Profile + hand_info: dict[str, Any] + axis_offset: tuple[float, float, float] + points: int + + +def _model_points(mesh_axis: Any, side: FeatureSide) -> Array: + zr = side.z or (float(mesh_axis.bounds[0][2]), float(mesh_axis.bounds[1][2])) + return helical.mesh_wall_points(mesh_axis, zr, side.r, SLICE_STEP_MM, 0.1) + + +def _printed_points( + pp: PrintedPart, axis: Axis, side: FeatureSide +) -> tuple[Array, Array]: + from assembly.printed import surface_offsets, to_assembly + + pts, width = to_assembly(pp.walls, axis.matrix @ pp.bed_to_assembly) + keep = np.ones(len(pts), dtype=bool) + if side.z is not None: + keep &= (pts[:, 2] >= side.z[0]) & (pts[:, 2] <= side.z[1]) + if side.r is not None: + r = np.hypot(pts[:, 0], pts[:, 1]) + keep &= (r >= side.r[0]) & (r <= side.r[1]) + return pts[keep], surface_offsets(width[keep], side.side) + + +def _side_profile( + pts: Array, offset: Array | None, feature: Feature, hand: str, side: str +) -> SideProfile: + if len(pts) < helical.MIN_POINTS_PER_BIN * 10: + raise InvalidArgument( + f"features.{feature.name}", + f"{len(pts)} surface points in the window", + ["a z / r window that holds the feature's surface"], + command=CMD, + extra="Check the windows (axis frame: z along the axis, r from it) and the axis point.", + ) + prof, info = helical.pick_handedness(pts, feature.pitch_mm, hand, offset, side) + return SideProfile(prof, info, helical.fit_axis_offset(pts), len(pts)) + + +def _common_hand( + pa: SideProfile, pb: SideProfile, requested: str +) -> tuple[str, str | None]: + """The handedness both parts are judged under, and a mismatch message (both decided and different).""" + if requested in ("right", "left"): + return requested, None + ha, hb = pa.hand_info, pb.hand_info + if ha["hand"] == "none": + return "auto", None + if ha["decided"] and hb["decided"] and ha["hand"] != hb["hand"]: + return ha["hand"], ( + f"handedness mismatch: part A is {ha['hand']}-handed, part B {hb['hand']}-handed - they cannot screw together" + ) + if ha["decided"]: + return ha["hand"], None + if hb["decided"]: + return hb["hand"], None + return ha["best_fit"], None + + +def _profile_summary(p: SideProfile, side: FeatureSide) -> dict[str, Any]: + prof = p.profile + return { + "side": side.side, + "points": p.points, + "r_min_mm": round(prof.r_min, 4), + "r_max_mm": round(prof.r_max, 4), + "coverage": round(prof.coverage, 3), + "coherence_mm": round(prof.coherence_mm, 4), + "handedness": p.hand_info, + "axis_offset_mm": [round(p.axis_offset[0], 4), round(p.axis_offset[1], 4)], + "axis_fit_rms_mm": round(p.axis_offset[2], 4), + } + + +def _female_male( + feature: Feature, pa: SideProfile, pb: SideProfile +) -> tuple[SideProfile, SideProfile, bool]: + """(female, male, a_is_female). A feature is a bore against a boss; two bores or two bosses cannot mate.""" + if feature.a.side == feature.b.side: + raise InvalidArgument( + f"features.{feature.name}", + f"both sides are '{feature.a.side}'", + ["one 'inner' (bore) and one 'outer' (boss) side"], + command=CMD, + ) + a_female = feature.a.side == "inner" + return (pa, pb, True) if a_female else (pb, pa, False) + + +def _clearance_dict( + feature: Feature, + cl: helical.Clearance, + hand: int, + a_is_female: bool, +) -> dict[str, Any]: + # `best_shift_bins` rolls the MALE profile: the part to move is whichever of A / B is male + mover = "b" if a_is_female else "a" + motion = helical.shift_to_motion(cl.best_shift_bins, feature.pitch_mm, hand) + if mover == "a": # moving A by +x is moving B by -x for the same relative phase + motion = {k: -v for k, v in motion.items()} + play, unit = cl.free_play + return { + "at_pose_min_gap_mm": _r(cl.at_pose.min_mm), + "at_pose_mean_gap_mm": _r(cl.at_pose.mean_mm), + "best_phase_min_gap_mm": _r(cl.best_gap_mm), + "best_phase_move_b": motion, + "free_play": { + "value": _r(play, 3), + "unit": unit, + "fraction_of_turn": _r(cl.free_fraction, 3), + }, + "male_crest_to_female_root_mm": _r(cl.crest_male_to_root_female_mm), + "female_crest_to_male_root_mm": _r(cl.crest_female_to_root_male_mm), + } + + +def _r(v: float, nd: int = 4) -> float | None: + return None if math.isnan(v) else round(float(v), nd) + + +def _side_findings(key: str, side: dict[str, Any], pitch_mm: float) -> tuple[list[str], list[str]]: + """(fails, warns) of one part's measured feature surface: coverage, handedness decisiveness, concentricity.""" + who = key.upper() + fails: list[str] = [] + warns: list[str] = [] + if side["coverage"] < MIN_COVERAGE_FAIL: + fails.append(f"{who} feature surface only {side['coverage']:.0%} measured (empty phase bins)") + elif side["coverage"] < MIN_COVERAGE_WARN: + warns.append(f"{who} feature surface {side['coverage']:.0%} measured") + if pitch_mm > 0 and not side["handedness"]["decided"]: + warns.append( + f"{who} handedness is not decisive (coherence {side['handedness']['coherence_mm']}): " + "check the pitch and the windows" + ) + off = math.hypot(*side["axis_offset_mm"]) + # a circle fit only means something for a thread; keys and ribs are not circles (the number is still reported) + if pitch_mm > 0 and off > AXIS_OFFSET_WARN_MM: + warns.append(f"{who} feature axis is {off:.3f} mm off the declared axis") + return fails, warns + + +def _clearance_findings(cl: dict[str, Any], required: float) -> tuple[list[str], list[str]]: + best, pose = cl["best_phase_min_gap_mm"], cl["at_pose_min_gap_mm"] + fails: list[str] = [] + warns: list[str] = [] + if best is None: + fails.append("no phase with both profiles measured") + elif best < required: + fails.append(f"best-phase radial gap {best:.3f} mm < required {required:.3f} mm") + if pose is not None and pose < 0 and (best is None or best >= required): + warns.append( + f"the posed phase interferes (radial gap {pose:.3f} mm) though a free phase exists: " + f"move B by {cl['best_phase_move_b']}" + ) + return fails, warns + + +def _status( + feature: Feature, required: float, res: dict[str, Any] +) -> tuple[str, list[str]]: + fails: list[str] = [] + warns: list[str] = [] + for key in ("a", "b"): + f, w = _side_findings(key, res["profiles"][key], feature.pitch_mm) + fails += f + warns += w + if res.get("mismatch"): + fails.append(res["mismatch"]) + if res.get("clearance") is not None: + f, w = _clearance_findings(res["clearance"], required) + fails += f + warns += w + status = "fail" if fails else "warn" if warns else "pass" + return status, fails + warns + + +@dataclass +class _Ctx: + axis: Axis + mesh_a_axis: Any + mesh_b_axis: Any + printed_a: PrintedPart | None + printed_b: PrintedPart | None + default_required: float + + +def _points_for( + ctx: _Ctx, label: str, side: FeatureSide, source: str +) -> tuple[Array, Array | None]: + pp = ctx.printed_a if label == "a" else ctx.printed_b + if source == "as_printed" and pp is not None: + return _printed_points(pp, ctx.axis, side) + mesh = ctx.mesh_a_axis if label == "a" else ctx.mesh_b_axis + return _model_points(mesh, side), None + + +def analyse_feature(feature: Feature, ctx: _Ctx, source: str) -> dict[str, Any]: + pts_a, off_a = _points_for(ctx, "a", feature.a, source) + pts_b, off_b = _points_for(ctx, "b", feature.b, source) + pa = _side_profile(pts_a, off_a, feature, feature.hand, feature.a.side) + pb = _side_profile(pts_b, off_b, feature, feature.hand, feature.b.side) + hand, mismatch = _common_hand(pa, pb, feature.hand) + if feature.pitch_mm > 0 and hand in ("right", "left") and mismatch is None: + pa = _side_profile(pts_a, off_a, feature, hand, feature.a.side) + pb = _side_profile(pts_b, off_b, feature, hand, feature.b.side) + female, male, a_female = _female_male(feature, pa, pb) + res: dict[str, Any] = { + "source": source, + "profiles": { + "a": _profile_summary(pa, feature.a), + "b": _profile_summary(pb, feature.b), + }, + "mismatch": mismatch, + } + if mismatch is None: + cl = helical.clearance(female.profile, male.profile) + res["clearance"] = _clearance_dict( + feature, cl, female.profile.hand or male.profile.hand, a_female + ) + required = ( + feature.min_gap_mm if feature.min_gap_mm is not None else ctx.default_required + ) + res["required_gap_mm"] = required + res["status"], res["reasons"] = _status(feature, required, res) + return res + + +def _printed_vs_model(feature: Feature, ctx: _Ctx) -> dict[str, Any]: + """How far the printed surface of each printed part sits from the model's, on this feature's profile.""" + out: dict[str, Any] = {} + for label, side, pp in ( + ("a", feature.a, ctx.printed_a), + ("b", feature.b, ctx.printed_b), + ): + if pp is None: + continue + model = _side_profile( + _model_points(ctx.mesh_a_axis if label == "a" else ctx.mesh_b_axis, side), + None, + feature, + feature.hand, + side.side, + ) + pts, off = _printed_points(pp, ctx.axis, side) + prin = _side_profile(pts, off, feature, feature.hand, side.side) + diff = ( + prin.profile.radius - model.profile.radius + if model.profile.hand == prin.profile.hand + else None + ) + window_layers = printed.layers_in_window( + pp.walls, ctx.axis.matrix @ pp.bed_to_assembly, side.z + ) + out[label] = { + "wall_vs_model_in_window": printed.deviation( + pp.bed_mesh, window_layers, pp.skip_below_mm, pp.shift + ), + "printed_r_min_mm": round(prin.profile.r_min, 4), + "printed_r_max_mm": round(prin.profile.r_max, 4), + "model_r_min_mm": round(model.profile.r_min, 4), + "model_r_max_mm": round(model.profile.r_max, 4), + "radius_printed_minus_model_mm": None + if diff is None + else { + "mean": _r(float(np.nanmean(diff))), + "min": _r(float(np.nanmin(diff))), + "max": _r(float(np.nanmax(diff))), + }, + } + return out + + +REACH_MM = 3.0 # a bore surface farther than this beyond the boss's largest radius is not part of the fit + + +def resolve_windows(feature: Feature, ctx: _Ctx) -> Feature: + """Fill an omitted radius window of the BORE side: from the axis out to the boss's largest radius + REACH_MM. + + Without it a tube's own outer wall would be read as part of its bore. An omitted window of the boss, or an omitted z + window, means the whole part.""" + bore_label = "a" if feature.a.side == "inner" else "b" + bore = feature.a if bore_label == "a" else feature.b + boss_label = "b" if bore_label == "a" else "a" + boss = feature.b if bore_label == "a" else feature.a + if bore.r is not None or feature.a.side == feature.b.side: + return feature + mesh = ctx.mesh_a_axis if boss_label == "a" else ctx.mesh_b_axis + pts = _model_points(mesh, boss) + if len(pts) == 0: + return feature + limit = float(np.hypot(pts[:, 0], pts[:, 1]).max()) + REACH_MM + bore = replace(bore, r=(0.0, limit)) + return replace(feature, a=bore) if bore_label == "a" else replace(feature, b=bore) + + +def feature_report(feature: Feature, ctx: _Ctx) -> dict[str, Any]: + feature = resolve_windows(feature, ctx) + sources = ["model"] + if ctx.printed_a is not None or ctx.printed_b is not None: + sources.append("as_printed") + doc: dict[str, Any] = { + "name": feature.name, + "kind": "thread" if feature.pitch_mm > 0 else "ribs", + "pitch_mm": feature.pitch_mm, + "windows": { + "a": {"side": feature.a.side, "z": feature.a.z, "r": feature.a.r}, + "b": {"side": feature.b.side, "z": feature.b.z, "r": feature.b.r}, + }, + } + worst = "pass" + for src in sources: + doc[src] = analyse_feature(feature, ctx, src) + worst = _worse(worst, doc[src]["status"]) + if len(sources) > 1: + doc["printed_vs_model"] = _printed_vs_model(feature, ctx) + doc["as_printed_parts"] = [ + k for k, v in (("a", ctx.printed_a), ("b", ctx.printed_b)) if v is not None + ] + doc["status"] = worst + return doc + + +def _worse(a: str, b: str) -> str: + order = {"pass": 0, "warn": 1, "fail": 2} + return a if order[a] >= order[b] else b + + +def _screw_matrix(axis: Axis, motion: dict[str, float]) -> Matrix: + """The pose change of part B that realises `best_phase_move_b`: a turn about the axis (the axial position, e.g. a + seat face, stays).""" + from assembly.scene import axis_rotation + + return axis_rotation(axis.direction, motion["rotate_deg"], axis.origin) + + +def fit_requirement( + fit: str, + shape: str, + printer: str | None, + material: str | None, + nozzle: float | None, + layer: float | None, +) -> dict[str, Any]: + """The radial gap the `3d fits` table asks of this fit on this printer and plastic: clearance per face - tolerance.""" + import fits_sources + + res = fits_sources.resolve(fit, shape, printer, material, nozzle, layer, command=CMD) + return { + "fit_type": res.fit_type, + "shape": shape, + "clearance_mm": round(res.clearance_mm, 4), + "tol_mm": round(res.tol_mm, 4), + "interference_mm": round(res.interference_mm, 4), + "required_gap_mm": round(res.clearance_mm - res.tol_mm, 4), + "fallback": res.fallback, + "source": "fit table (`3d fits resolve`)", + "provenance": [dict(p) for p in res.provenance], + } + + +def _fit_required(opts: PairOptions) -> float: + if opts.min_radial_gap_mm is not None: + return opts.min_radial_gap_mm + if opts.fit_resolution is not None: + return float(opts.fit_resolution["required_gap_mm"]) + return DEFAULT_MIN_RADIAL_GAP_MM + + +@dataclass +class _Setup: + """Everything `run_pair` measures: the posed meshes, the axis, the printed parts, the screw that was applied.""" + + posed_a: Any + posed_b: Any + axis: Axis + printed_a: PrintedPart | None + printed_b: PrintedPart | None + required: float + screw: dict[str, Any] | None + + +def _prepare(opts: PairOptions) -> _Setup: + mesh_a, qa = _load(opts.a, "a", opts.max_faces) + mesh_b, qb = _load(opts.b, "b", opts.max_faces) + posed_a, posed_b = mesh_a.copy(), mesh_b.copy() + posed_a.apply_transform(qa) + posed_b.apply_transform(qb) + axis = make_axis(opts, posed_a) + printed_a = ( + load_printed(opts.a, "a", mesh_a, qa, opts.skip_below_mm) if opts.a.gcode else None + ) + printed_b = ( + load_printed(opts.b, "b", mesh_b, qb, opts.skip_below_mm) if opts.b.gcode else None + ) + required = _fit_required(opts) + screw = None + if opts.screw_to_best and opts.features: + screw, qb, posed_b = _apply_screw(opts, axis, posed_a, posed_b, qb, required) + if printed_b is not None: + printed_b = replace( + printed_b, bed_to_assembly=screw["matrix"] @ printed_b.bed_to_assembly + ) + return _Setup(posed_a, posed_b, axis, printed_a, printed_b, required, screw) + + +def _head(opts: PairOptions, su: _Setup) -> dict[str, Any]: + doc: dict[str, Any] = { + "schema": SCHEMA, + "kind": "pair", + "parts": { + "a": _part_doc(opts.a, su.posed_a), + "b": _part_doc(opts.b, su.posed_b), + }, + "axis": { + "point": [round(float(v), 4) for v in su.axis.origin], + "direction": [round(float(v), 6) for v in su.axis.direction], + }, + } + if su.screw is not None: + doc["screw_to_best"] = {k: v for k, v in su.screw.items() if k != "matrix"} + return doc + + +def _verdict(doc: dict[str, Any]) -> None: + """`ok` and the names of the failed / warned checks, from the statuses already in `doc`.""" + named = _named(doc) + doc["ok"] = all(s != "fail" for _, s in named) + doc["summary"] = { + "fail": [n for n, s in named if s == "fail"], + "warn": [n for n, s in named if s == "warn"], + } + + +def run_pair(opts: PairOptions) -> dict[str, Any]: + """Judge two parts; the returned document is the `--json` payload (`kind: pair`).""" + su = _prepare(opts) + ctx = _Ctx( + su.axis, + _to_axis(su.posed_a, su.axis), + _to_axis(su.posed_b, su.axis), + su.printed_a, + su.printed_b, + su.required, + ) + doc = _head(opts, su) + solid_a = solid_from_mesh(su.posed_a, "a", max_faces=opts.max_faces) + solid_b = solid_from_mesh(su.posed_b, "b", max_faces=opts.max_faces) + doc["interpenetration"] = interpenetration( + solid_a, solid_b, opts.max_overlap_mm3, opts.seat_tol_mm + ) + doc["clearance"] = clearance( + solid_a, solid_b, su.posed_a, su.posed_b, opts.min_clearance_mm, opts.zone + ) + doc["features"] = [feature_report(f, ctx) for f in opts.features] + printed_doc = { + k: v.info for k, v in (("a", su.printed_a), ("b", su.printed_b)) if v is not None + } + if printed_doc: + doc["printed"] = printed_doc + if opts.fit_resolution is not None: + doc["fit"] = opts.fit_resolution + _verdict(doc) + return doc + + +def _named(doc: dict[str, Any]) -> list[tuple[str, str]]: + rows = [ + ("interpenetration", doc["interpenetration"]["status"]), + ("clearance", doc["clearance"]["status"]), + ] + rows.extend((f"feature:{f['name']}", f["status"]) for f in doc["features"]) + return rows + + +def _to_axis(mesh: Any, axis: Axis) -> Any: + out = mesh.copy() + out.apply_transform(axis.matrix) + return out + + +def _part_doc(part: PartInput, posed: Any) -> dict[str, Any]: + return { + "mesh": str(part.path), + "scale": part.scale, + "bbox_assembled": _bbox(posed), + "faces": len(posed.faces), + "gcode": str(part.gcode) if part.gcode else None, + } + + +def _apply_screw( + opts: PairOptions, + axis: Axis, + posed_a: Any, + posed_b: Any, + qb: Matrix, + required: float, +) -> tuple[dict[str, Any], Matrix, Any]: + """Probe the first feature on the MODEL, then move B by the motion that brings it to its best phase.""" + probe = _Ctx( + axis, _to_axis(posed_a, axis), _to_axis(posed_b, axis), None, None, required + ) + feature = resolve_windows(opts.features[0], probe) + res = analyse_feature(feature, probe, "model") + cl = res.get("clearance") + if cl is None: + raise InvalidArgument( + "--screw-to-best", + "no common phase", + [res.get("mismatch") or "a measurable feature"], + command=CMD, + ) + motion = cl["best_phase_move_b"] + m = _screw_matrix(axis, motion) + moved = posed_b.copy() + moved.apply_transform(m) + return ( + {"feature": feature.name, "applied_to_b": motion, "matrix": m}, + m @ qb, + moved, + ) + + +def pose_from_text(text: str | None, where: str) -> Pose: + """A pose written inline (`{translate: [0, 0, 5], rotate: [0, 0, 90]}` or a list of them) -> spec pose ops.""" + if not text: + return () + import yaml + + from assembly.spec import parse_pose + + try: + value = yaml.safe_load(text) + except yaml.YAMLError as exc: + raise InvalidArgument( + where, text, ["a pose mapping, e.g. '{translate: [0, 0, 5]}'"], command=CMD + ) from exc + return parse_pose(value, where) + + +def features_from_file(path: Path) -> list[Feature]: + """`features:` list of a YAML / JSON file; see docs/commands/assembly.md.""" + import yaml + + if not path.is_file(): + raise InputNotFound(str(path), command=CMD) + data = yaml.safe_load(path.read_text()) + rows = data.get("features") if isinstance(data, dict) else data + if not isinstance(rows, list) or not rows: + raise InvalidArgument( + "--features", + str(path), + ["a `features:` list with at least one entry"], + command=CMD, + ) + return [_feature_from_mapping(r, f"{path.name}[{i}]") for i, r in enumerate(rows)] + + +def _side_from_mapping(raw: Any, where: str) -> FeatureSide: + if not isinstance(raw, dict) or raw.get("side") not in ("inner", "outer"): + raise InvalidArgument( + where, + str(raw), + ["{side: inner|outer, z: [lo, hi], r: [lo, hi]}"], + command=CMD, + ) + + def rng(key: str) -> Range | None: + v = raw.get(key) + if v is None: + return None + return parse_range(f"{v[0]}:{v[1]}", f"{where}.{key}") + + return FeatureSide(raw["side"], rng("z"), rng("r")) + + +def _feature_from_mapping(raw: Any, where: str) -> Feature: + if not isinstance(raw, dict): + raise InvalidArgument(where, str(raw), ["a mapping"], command=CMD) + unknown = set(raw) - {"name", "pitch", "hand", "a", "b", "min_gap_mm"} + if unknown: + raise InvalidArgument( + where, + ", ".join(sorted(unknown)), + ["name, pitch, hand, a, b, min_gap_mm"], + command=CMD, + ) + hand = str(raw.get("hand", "auto")) + if hand not in ("auto", "right", "left"): + raise InvalidArgument( + f"{where}.hand", hand, ["auto", "right", "left"], command=CMD + ) + pitch = float(raw.get("pitch", 0.0)) + if pitch < 0: + raise InvalidArgument( + f"{where}.pitch", str(pitch), ["0 (ribs) or a positive pitch"], command=CMD + ) + gap = raw.get("min_gap_mm") + return Feature( + str(raw.get("name", "feature")), + _side_from_mapping(raw.get("a"), f"{where}.a"), + _side_from_mapping(raw.get("b"), f"{where}.b"), + pitch, + hand, + None if gap is None else float(gap), + ) diff --git a/lib/assembly/pair_report.py b/lib/assembly/pair_report.py new file mode 100644 index 00000000..09aaab68 --- /dev/null +++ b/lib/assembly/pair_report.py @@ -0,0 +1,162 @@ +"""pair_report.py — the text of a `3d assembly pair` report (the `--json` document is the source of truth). + +WHAT: a pure function from the `pair.run_pair` document to the lines printed without `--json`. No geometry, no I/O. +""" + +from __future__ import annotations + +from typing import Any + +from assembly.report import MARK + + +def _g(v: float | None) -> str: + return "n/a" if v is None else f"{v:.3f}" + + +def _printed_lines(doc: dict[str, Any]) -> list[str]: + out: list[str] = [] + for key, info in doc.get("printed", {}).items(): + pl, dev = info["placement"], info["deviation"] + how = "registered" if pl["registered"] else "given" + resid = ( + f", residual {_g(pl['residual_rms_mm'])} mm" + if pl["residual_rms_mm"] is not None + else "" + ) + out.append( + f"as printed {key.upper()}: {info['gcode']} (object {info['object'] or 'all'}), placement {how} " + f"(shift {pl['shift_applied_mm']}{resid})" + ) + if info["overhang_wall_left_out_mm"] > 0.05 * max(info["wall_mm"], 1.0): + out.append( + f" note: {info['overhang_wall_left_out_mm']} mm of `Overhang wall` (vs {info['wall_mm']} mm outer wall) " + "was left out; thread flanks that overhang are thinly measured (see coverage)" + ) + if info["z_lift_mm"]: + out.append(f" note: the object starts {info['z_lift_mm']} mm above the bed (raft); the model is lifted by it") + if dev.get("layers"): + first = dev.get("first_layer_max_mm") + out.append( + f" wall vs model: max {dev['max_mm']:.3f} mm at z {dev['max_at_z_mm']:.2f}, mean of layer maxima " + f"{dev['mean_of_layer_max_mm']:.3f} mm over {dev['layers']} layers" + + (f"; first layer {first:.3f} mm" if first is not None else "") + ) + return out + + +def _source_lines(label: str, res: dict[str, Any], feat: dict[str, Any]) -> list[str]: + tag = MARK[res["status"]] + out = [f" [{label}] {tag} required radial gap >= {res['required_gap_mm']:.3f} mm"] + for key in ("a", "b"): + p = res["profiles"][key] + hand = p["handedness"]["hand"] + out.append( + f" {key.upper()} {p['side']:<5} r {p['r_min_mm']:.3f}..{p['r_max_mm']:.3f} mm, coverage {p['coverage']:.0%}, " + f"coherence {p['coherence_mm']:.3f} mm" + + ("" if hand == "none" else f", {hand}-handed") + ) + cl = res.get("clearance") + if cl is not None: + play = cl["free_play"] + move = cl["best_phase_move_b"] + out.append( + f" radial gap at the posed phase: min {_g(cl['at_pose_min_gap_mm'])} mm (mean {_g(cl['at_pose_mean_gap_mm'])})" + ) + out.append( + f" best phase: tightest radial gap {_g(cl['best_phase_min_gap_mm'])} mm (move B {move}); " + f"free play {_g(play['value'])} {play['unit']} ({_g(play['fraction_of_turn'])} of a turn)" + ) + out.append( + f" male crest -> female root {_g(cl['male_crest_to_female_root_mm'])} mm, " + f"female crest -> male root {_g(cl['female_crest_to_male_root_mm'])} mm" + ) + out.extend(f" ! {reason}" for reason in res["reasons"]) + return out + + +def _feature_lines(feat: dict[str, Any]) -> list[str]: + pitch = f", pitch {feat['pitch_mm']:.3f} mm" if feat["pitch_mm"] else "" + out = [f"feature {feat['name']} ({feat['kind']}{pitch}): {MARK[feat['status']]}"] + out += _source_lines("model", feat["model"], feat) + if "as_printed" in feat: + parts = "+".join(p.upper() for p in feat["as_printed_parts"]) + out += _source_lines( + f"as printed ({parts}, rest model)", feat["as_printed"], feat + ) + for key, d in feat.get("printed_vs_model", {}).items(): + diff = d["radius_printed_minus_model_mm"] + tail = ( + "" + if diff is None + else f"; printed - model radius mean {diff['mean']:+.3f}, min {diff['min']:+.3f}, max {diff['max']:+.3f} mm" + ) + wd = d["wall_vs_model_in_window"] + if wd.get("layers"): + out.append( + f" printed {key.upper()} wall vs model inside the feature window: max {wd['max_mm']:.3f} mm " + f"(p95 {wd['p95_of_layer_max_mm']:.3f}) over {wd['layers']} layers" + ) + out.append( + f" printed {key.upper()} r {d['printed_r_min_mm']:.3f}..{d['printed_r_max_mm']:.3f} vs model " + f"{d['model_r_min_mm']:.3f}..{d['model_r_max_mm']:.3f} mm{tail}" + ) + return out + + +def pair_lines(doc: dict[str, Any]) -> list[str]: + a, b = doc["parts"]["a"], doc["parts"]["b"] + inter, clr = doc["interpenetration"], doc["clearance"] + lines = [ + f"assembly pair: A {a['mesh']} (x{a['scale']}) B {b['mesh']} (x{b['scale']}) -> {'OK' if doc['ok'] else 'FAIL'}", + f"axis through {doc['axis']['point']} along {doc['axis']['direction']}", + ] + if "screw_to_best" in doc: + s = doc["screw_to_best"] + lines.append( + f"screw-to-best: B moved by {s['applied_to_b']} to the best phase of feature '{s['feature']}'" + ) + lines.append( + f"interpenetration: {MARK[inter['status']]} {inter['volume_mm3']:.4f} mm3 in {inter['count']} region(s) " + f"(allowed {inter['allowed_mm3']} mm3)" + ) + for i, r in enumerate(inter["regions"][:5], 1): + c = r["center"] + lines.append( + f" [{i}] {r['volume_mm3']:.4f} mm3, depth {r['depth_mm']:.3f} mm at ({c[0]:.2f}, {c[1]:.2f}, {c[2]:.2f})" + ) + ct = inter["contacts"] + if ct["count"]: + lines.append( + f" contact (thinner than {ct['seat_tol_mm']} mm, a seat): {ct['count']} region(s), " + f"{ct['volume_mm3']:.4f} mm3, depth <= {ct['max_depth_mm']:.4f} mm" + ) + req = ( + f" (required >= {clr['required_mm']} mm)" + if clr["required_mm"] is not None + else "" + ) + cap = f" (search capped at {clr['capped_at_mm']} mm)" if clr["capped_at_mm"] else "" + lines.append( + f"minimum clearance: {MARK[clr['status']]} {clr['min_gap_mm']:.4f} mm{cap}{req}" + ) + if "zone" in clr: + z = clr["zone"] + lines.append( + f" inside the zone {z['box']}: {z['min_gap_mm']:.4f} mm " + f"(surface samples every {z['sample_spacing_mm']} mm)" + ) + if "fit" in doc: + f = doc["fit"] + lines.append( + f"fit {f['fit_type']}/{f['shape']}: clearance {f['clearance_mm']:+.3f} mm, tol {f['tol_mm']:.3f} -> " + f"required radial gap {f['required_gap_mm']:.3f} mm ({f['source']})" + ) + lines += _printed_lines(doc) + for feat in doc["features"]: + lines += _feature_lines(feat) + if doc["summary"]["fail"]: + lines.append("FAILED: " + ", ".join(doc["summary"]["fail"])) + if doc["summary"]["warn"]: + lines.append("warnings: " + ", ".join(doc["summary"]["warn"])) + return lines diff --git a/lib/assembly/printed.py b/lib/assembly/printed.py new file mode 100644 index 00000000..06c70bb4 --- /dev/null +++ b/lib/assembly/printed.py @@ -0,0 +1,233 @@ +"""printed.py — the as-printed side of the fit check: G-code wall paths against the model, in the assembly frame. + +WHAT: for one part with a sliced G-code, + 1. places the model on the bed (`PrintPlacement`): its print orientation (rotate / scale, like a pose), dropped onto + z = 0, and the XY position either given or REGISTERED - a translation-only least-squares fit of the model's + layer sections to the object's printed outer walls (the slicer's own placement, usually the project's object + position), so a G-code that was not made from this exact model-to-bed transform still lines up; + 2. measures the per-layer deviation of the printed wall from the model section: the centre line of an outer wall lies + half a line width inside the surface, so `|signed distance to the section - width/2|` is how far the printed + surface is from the model (elephant foot, hole / thread compensation, a wrong model all show here); + 3. returns the printed wall points mapped into the ASSEMBLY frame (bed -> model -> pose), each with the offset to the + printed SURFACE (+-width/2 along the radius), ready for `helical` profiles. +WHY: the slicer moves walls (xy_hole_compensation, elephant-foot compensation, line width): a thread clearance of 0.6 mm + on the model is another number on the plastic; the cactus rim/plug screwing was verified this way only by an ad-hoc script. +INVARIANTS: the deviation uses the layer's mid height (`z - h/2`); layers below `skip_below_mm` are reported separately + (the first layer is flared on purpose); the registration is reported with its residual - a model that is mirrored or + turned against the G-code does not register and the deviation says so. +""" + +from __future__ import annotations + +import math +from dataclasses import dataclass, field +from typing import Any + +import numpy as np +import numpy.typing as npt + +from assembly import helical +from gcode_walls import WallLayer, Walls, densify_segments + +Array = npt.NDArray[np.float64] +Matrix = npt.NDArray[np.float64] + +REGISTER_LAYERS = 12 +REGISTER_ITERATIONS = 12 +SPACING_MM = 0.1 +DEVIATION_SPACING_MM = 0.4 + + +@dataclass +class Section: + """The cross-section of the model at one height: boundary points, the (holed) polygons and the plane frame.""" + + z: float + polygons: list[Any] + boundary: Array # (N, 2) points along the boundary, ~0.05 mm apart + to_3d: Matrix + tree: Any = field(repr=False, default=None) # cKDTree of `boundary` + + +@dataclass +class PrintPlacement: + matrix: Matrix # model (already scaled) -> bed + registered: bool + residual_mm: float + shift_mm: tuple[float, float] + + +def sections(mesh: Any, heights: Array) -> list[Section | None]: + """Sections of a bed-frame mesh at the given z values (None where the plane misses it).""" + import trimesh + from scipy.spatial import cKDTree + from trimesh import intersections + + lines, to_3d, _ = intersections.mesh_multiplane( + mesh, np.zeros(3), np.array([0.0, 0.0, 1.0]), heights + ) + out: list[Section | None] = [] + for h, seg, mat in zip(heights, lines, to_3d, strict=True): + if len(seg) == 0: + out.append(None) + continue + seg = np.asarray(seg, dtype=float) + path = trimesh.load_path(seg) + polys = list(path.polygons_full) + if not polys: + out.append(None) + continue + pts, _ = densify_segments(seg[:, 0, :], seg[:, 1, :], 0.05) + out.append( + Section(float(h), polys, pts, np.asarray(mat, dtype=float), cKDTree(pts)) + ) + return out + + +def signed_distance(sec: Section, plane_xy: Array) -> tuple[Array, Array, Array]: + """(distance to the boundary, inside flag, nearest boundary point) of plane points.""" + import shapely + + dist, idx = sec.tree.query(plane_xy) + inside = np.zeros(len(plane_xy), dtype=bool) + for poly in sec.polygons: + inside |= shapely.contains_xy(poly, plane_xy[:, 0], plane_xy[:, 1]) + return np.asarray(dist), inside, sec.boundary[idx] + + +def _layer_pts(layer: WallLayer, spacing: float) -> tuple[Array, Array]: + pts, width = layer.points(spacing) + return pts, width + + +def _gauss_newton_rows( + data: list[tuple[Array, Array, Section]], shift: Array +) -> tuple[Array, Array]: + """Jacobian rows (the inward boundary normals) and residuals `signed distance - width/2` of the wall points at `shift`.""" + rows_a: list[Array] = [] + rows_b: list[Array] = [] + for pts, width, sec in data: + p = pts + shift + dist, inside, q = signed_distance(sec, p) + sign = np.where(inside, 1.0, -1.0) + n_in = (p - q) * sign[:, None] + n_in = n_in / np.maximum(np.linalg.norm(n_in, axis=1), 1e-9)[:, None] + rows_a.append(n_in) + rows_b.append(sign * dist - width / 2.0) + return np.vstack(rows_a), np.concatenate(rows_b) + + +def register_xy( + mesh_bed: Any, walls: Walls, skip_below_mm: float +) -> tuple[tuple[float, float], float]: + """Translation (dx, dy) that, ADDED TO THE WALL POINTS, lays them best on the model, and the rms residual of + `signed distance - width/2` after it. Gauss-Newton on the signed distance, started at the centroid difference.""" + layers = [ + lyr + for lyr in walls.layers + if lyr.mid_z is not None and lyr.mid_z > skip_below_mm + ] + pick = layers[:: max(1, len(layers) // REGISTER_LAYERS)] + secs = sections(mesh_bed, np.array([lyr.mid_z for lyr in pick], dtype=float)) + data = [ + (*_layer_pts(lyr, 0.3), s) for lyr, s in zip(pick, secs, strict=True) if s is not None + ] + if not data: + return (0.0, 0.0), math.nan + wall_c = np.mean(np.vstack([d[0] for d in data]), axis=0) + bnd_c = np.mean(np.vstack([d[2].boundary for d in data]), axis=0) + shift = np.asarray(bnd_c - wall_c, dtype=float) # added to the wall points it brings the centroids together + for _ in range(REGISTER_ITERATIONS): + a, b = _gauss_newton_rows(data, shift) + step, *_ = np.linalg.lstsq(a, -b, rcond=None) + shift = shift + step + if float(np.linalg.norm(step)) < 1e-4: + break + _, resid = _gauss_newton_rows(data, shift) + return (float(shift[0]), float(shift[1])), float(np.sqrt(np.mean(resid**2))) + + +def layers_in_window( + walls: Walls, bed_to_axis: Matrix, z_range: tuple[float, float] | None +) -> list[WallLayer]: + """The layers whose wall centroid lies inside `z_range` along the axis (all layers for None).""" + layers = [lyr for lyr in walls.layers if lyr.mid_z is not None] + if z_range is None: + return layers + keep: list[WallLayer] = [] + for lyr in layers: + c = lyr.a.mean(axis=0) + z = helical.apply(bed_to_axis, np.array([[c[0], c[1], lyr.mid_z]]))[0, 2] + if z_range[0] <= z <= z_range[1]: + keep.append(lyr) + return keep + + +def deviation( + mesh_bed: Any, + layers: list[WallLayer], + skip_below_mm: float, + shift: tuple[float, float] = (0.0, 0.0), +) -> dict[str, Any]: + """Per-layer deviation of the printed wall from the model section (`shift` is added to the wall points).""" + secs = sections(mesh_bed, np.array([lyr.mid_z for lyr in layers], dtype=float)) + rows: list[tuple[float, float]] = [] + first: float | None = None + for lyr, sec in zip(layers, secs, strict=True): + if sec is None or lyr.mid_z is None: + continue + pts, width = lyr.points(DEVIATION_SPACING_MM) + p = pts + np.asarray(shift) + dist, inside, _ = signed_distance(sec, p) + signed = np.where(inside, dist, -dist) + dev = np.abs(signed - width / 2.0) + if lyr.mid_z <= skip_below_mm: + first = float(dev.max()) if first is None else max(first, float(dev.max())) + continue + rows.append((lyr.mid_z, float(dev.max()))) + if not rows: + return {"layers": 0} + arr = np.array(rows) + worst = int(np.argmax(arr[:, 1])) + return { + "layers": len(rows), + "max_mm": round(float(arr[worst, 1]), 4), + "max_at_z_mm": round(float(arr[worst, 0]), 3), + "mean_of_layer_max_mm": round(float(arr[:, 1].mean()), 4), + "p95_of_layer_max_mm": round(float(np.percentile(arr[:, 1], 95)), 4), + "first_layer_max_mm": None if first is None else round(first, 4), + "skipped_below_mm": skip_below_mm, + } + + +def surface_offsets(width: Array, side: str) -> Array: + """Radial offset from a wall's centre line to its surface: outward surface (boss) = +w/2, inward (bore) = -w/2.""" + return width * 0.5 * (1.0 if side == "outer" else -1.0) + + +@dataclass +class PrintedWalls: + """Printed wall points in the assembly frame plus how they were placed.""" + + points: Array # (N, 3) assembly frame + width: Array + placement: PrintPlacement + deviation: dict[str, Any] = field(default_factory=dict) + + +def to_assembly( + walls: Walls, bed_to_assembly: Matrix, spacing: float = SPACING_MM +) -> tuple[Array, Array]: + """Printed wall points (bed xy at the layer mid height) -> assembly frame, with their line widths.""" + pts_all: list[Array] = [] + w_all: list[Array] = [] + for lyr in walls.layers: + if lyr.mid_z is None: + continue + xy, w = lyr.points(spacing) + p = np.column_stack([xy, np.full(len(xy), lyr.mid_z)]) + pts_all.append(helical.apply(bed_to_assembly, p)) + w_all.append(w) + if not pts_all: + return np.zeros((0, 3)), np.zeros(0) + return np.vstack(pts_all), np.concatenate(w_all) diff --git a/lib/assembly_pair_cli.py b/lib/assembly_pair_cli.py new file mode 100755 index 00000000..d5a43f6e --- /dev/null +++ b/lib/assembly_pair_cli.py @@ -0,0 +1,248 @@ +#!/usr/bin/env python3 +"""assembly_pair_cli.py — the worker behind `3d assembly pair A B` (run through `cli.pyrun` with trimesh / manifold3d / +numpy / scipy / shapely / PyYAML). + +WHAT: argument parsing and printing for the pair check; the geometry lives in `lib/assembly/pair.py` (+ `helical.py`, + `printed.py`, `lib/gcode_walls.py`). `commands/assembly.py` routes `pair` here and everything else to `assembly_cli.py`. + +CONVENTIONS: as `assembly_cli.py` - `--json` prints exactly one JSON document (schema `3d.assembly/1`, `kind: pair`); a failed + verdict exits 1 (`GateFailure`), bad input 2, a missing dependency 127. Errors are rendered here because this runs in a + child process of the dispatcher. A flag value that starts with one `-` (a negative window bound, `--axis -z`) is joined to + its flag with `=` first, which argparse needs. +""" + +from __future__ import annotations + +import argparse +import json +import sys +from pathlib import Path +from typing import Any + +from assembly_cli import CMD, _apply_mem_limit, _exit_fail, _fail, _peak_mb +from errors import ( + InvalidArgument, + MissingDependency, + ThreeDError, + UsageError, +) + +Box = tuple[float, float, float, float, float, float] + + +def _parser() -> argparse.ArgumentParser: + ap = argparse.ArgumentParser(prog="3d assembly pair", add_help=False) + ap.add_argument("sub", choices=("pair",)) + ap.add_argument("target") + ap.add_argument("other", nargs="?") + for part in ("a", "b"): + ap.add_argument(f"--scale-{part}", type=float, default=1.0) + ap.add_argument(f"--pose-{part}") + ap.add_argument(f"--print-{part}") + ap.add_argument(f"--gcode-{part}") + ap.add_argument(f"--object-{part}") + ap.add_argument(f"--{part}-side", choices=("inner", "outer")) + ap.add_argument(f"--{part}-z") + ap.add_argument(f"--{part}-r") + ap.add_argument("--axis", default="z") + ap.add_argument("--axis-point") + ap.add_argument("--pitch", type=float, default=0.0) + ap.add_argument("--hand", choices=("auto", "right", "left"), default="auto") + ap.add_argument("--feature-name", default="feature") + ap.add_argument("--features") + ap.add_argument("--min-radial-gap", type=float) + ap.add_argument("--min-clearance", type=float) + ap.add_argument("--max-overlap-mm3", type=float, default=0.01) + ap.add_argument("--seat-tol", type=float, default=0.02) + ap.add_argument("--zone") + ap.add_argument("--screw-to-best", action="store_true") + ap.add_argument("--skip-below", type=float, default=0.5) + ap.add_argument("--max-faces", type=int, default=400_000) + ap.add_argument("--fit") + ap.add_argument("--shape", default="round") + ap.add_argument("--printer") + ap.add_argument("--material") + ap.add_argument("--nozzle", type=float) + ap.add_argument("--layer", type=float) + ap.add_argument("--mem-limit-mb", type=int) + ap.add_argument("--json", action="store_true") + return ap + + +def _join_dash_values(argv: list[str]) -> list[str]: + """`--a-z -36:-28` / `--axis -z` -> `--a-z=-36:-28`: argparse takes a value that starts with one `-` for a flag.""" + out: list[str] = [] + for i, tok in enumerate(argv): + prev = argv[i - 1] if i else "" + joins = ( + prev.startswith("--") + and "=" not in prev + and tok.startswith("-") + and not tok.startswith("--") + and len(tok) > 1 + ) + if joins: + out[-1] = f"{prev}={tok}" + else: + out.append(tok) + return out + + +def _zone_box(text: str | None) -> Box | None: + if not text: + return None + try: + v = [float(x) for x in text.split(",")] + except ValueError: + v = [] + if len(v) != 6 or v[0] > v[3] or v[1] > v[4] or v[2] > v[5]: + raise InvalidArgument( + "--zone", text, ["XMIN,YMIN,ZMIN,XMAX,YMAX,ZMAX (min <= max)"], command=CMD + ) + return (v[0], v[1], v[2], v[3], v[4], v[5]) + + +def _features(args: argparse.Namespace) -> list[Any]: + """The feature(s) the flags or `--features FILE` describe (none: interpenetration and clearance only).""" + from assembly import pair + + flagged = args.a_side or args.b_side + if args.features: + if flagged or args.pitch: + raise UsageError( + "--features cannot be combined with --a-side/--b-side/--pitch", + command=CMD, + remediation=["Put the feature in the file, or use the flags alone."], + ) + return pair.features_from_file(Path(args.features)) + if not flagged and not args.pitch: + return [] + if not (args.a_side and args.b_side): + raise UsageError( + "a feature needs both --a-side and --b-side (inner = bore wall, outer = boss)", + command=CMD, + remediation=[ + "Example: --a-side inner --b-side outer --pitch 3.007 --a-z 1:9.5 --b-z 11:16.8" + ], + ) + + def side(label: str, kind: str) -> Any: + z, r = getattr(args, f"{label}_z"), getattr(args, f"{label}_r") + return pair.FeatureSide( + kind, + pair.parse_range(z, f"--{label}-z") if z else None, + pair.parse_range(r, f"--{label}-r") if r else None, + ) + + return [ + pair.Feature( + args.feature_name, + side("a", args.a_side), + side("b", args.b_side), + args.pitch, + args.hand, + ) + ] + + +def _part(args: argparse.Namespace, label: str, path: str) -> Any: + from assembly import pair + + gcode = getattr(args, f"gcode_{label}") + mesh_path, object_id = pair.split_part_path(path) + return pair.PartInput( + mesh_path, + getattr(args, f"scale_{label}"), + pair.pose_from_text(getattr(args, f"pose_{label}"), f"--pose-{label}"), + pair.pose_from_text(getattr(args, f"print_{label}"), f"--print-{label}"), + Path(gcode) if gcode else None, + getattr(args, f"object_{label}"), + object_id, + ) + + +def _options(args: argparse.Namespace) -> Any: + from assembly import pair + from assembly.vec import parse_direction, parse_point + + fit = ( + pair.fit_requirement( + args.fit, args.shape, args.printer, args.material, args.nozzle, args.layer + ) + if args.fit + else None + ) + return pair.PairOptions( + _part(args, "a", args.target), + _part(args, "b", args.other), + _features(args), + parse_point(args.axis_point, "--axis-point") if args.axis_point else None, + parse_direction(args.axis, "--axis"), + args.min_radial_gap, + args.min_clearance, + args.max_overlap_mm3, + args.seat_tol, + _zone_box(args.zone), + args.screw_to_best, + args.skip_below, + args.max_faces, + fit, + ) + + +def cmd_pair(args: argparse.Namespace) -> int: + from assembly import pair, pair_report + + if not args.other: + raise UsageError( + "`3d assembly pair` needs two parts: A and B", + command=CMD, + remediation=[ + "Example: 3d assembly pair rim.stl plug.stl --pose-b '{translate: [0, 0, -10]}'" + ], + ) + doc = pair.run_pair(_options(args)) + doc["resources"] = {"peak_rss_mb": _peak_mb()} + if args.json: + print(json.dumps(doc, indent=2)) + else: + print("\n".join(pair_report.pair_lines(doc))) + if doc["ok"]: + return 0 + return _exit_fail("assembly pair failed: " + ", ".join(doc["summary"]["fail"])) + + +def main(argv: list[str]) -> int: + try: + try: + args = _parser().parse_args(_join_dash_values(argv)) + except SystemExit as exc: + raise UsageError( + "bad arguments", command=CMD, remediation=["Run `3d assembly --help`."] + ) from exc + _apply_mem_limit(args.mem_limit_mb) + try: + return cmd_pair(args) + except MemoryError as exc: + raise _fail( + "out of memory while checking the pair", + ["Raise --mem-limit-mb or lower --max-faces."], + ) from exc + except ImportError as exc: + err: ThreeDError = MissingDependency( + str(exc.name or exc), + install="uv sync --extra dev", + degrades="`3d assembly pair` cannot run", + command=CMD, + ) + sys.stderr.write(err.render() + "\n") + return err.exit_code + except ThreeDError as exc: + text = exc.render() + if text: + sys.stderr.write(text + "\n") + return exc.exit_code + + +if __name__ == "__main__": + sys.exit(main(sys.argv[1:])) diff --git a/lib/commands/assembly.py b/lib/commands/assembly.py index 31aac6ed..a430d91f 100644 --- a/lib/commands/assembly.py +++ b/lib/commands/assembly.py @@ -11,7 +11,9 @@ explode exploded view PNG / USDZ, every part in its own colour with a legend; bom parts (role, reference source, size) and fasteners; stand-ins are labelled; position lay the printed parts on plates and write a project 3MF that `3d arrange` and `3d retarget` accept; - show the parsed spec and the colour legend. + show the parsed spec and the colour legend; + pair two parts, no spec: interpenetration volume, minimum clearance, per-feature radial clearance of threads + and ribs (model and, with G-code, as printed). WHY: one insert/socket pair at a time (`fit-check`) cannot tell whether a case, its board and its snap hooks go together, come apart and can be operated by hand. This is the whole-assembly check with the fit table behind every number. @@ -26,10 +28,10 @@ from cli.registry import Command from errors import UsageError -SUBS = ("check", "sweep", "order", "ergonomics", "explode", "bom", "position", "show") -DEPS = "trimesh,manifold3d,numpy,scipy,rtree,pyyaml,pillow,networkx" +SUBS = ("check", "sweep", "order", "ergonomics", "explode", "bom", "position", "show", "pair") +DEPS = "trimesh,manifold3d,numpy,scipy,rtree,pyyaml,pillow,networkx,shapely" -USAGE = """3d assembly TARGET [options] +USAGE = """3d assembly TARGET [options] Whole-assembly checks over a declared spec. TARGET is an assembly spec (assembly.yaml / .json), a project directory holding assembly.yaml, or (check / sweep / order / ergonomics / explode / bom / show) a 3MF written by `position`, which is then verified as printed. Every subcommand takes --json (schema `3d.assembly/1`, `kind` @@ -72,6 +74,17 @@ show TARGET the parsed spec, resolved context and the colour legend. Example: 3d assembly show assembly.yaml + pair A B [--scale-a S] [--scale-b S] [--pose-a POSE] [--pose-b POSE] [feature flags | --features FILE] [--zone BOX] + [--min-radial-gap MM | --fit TYPE] [--min-clearance MM] [--max-overlap-mm3 V] [--screw-to-best] + [--gcode-a F --object-a NAME --print-a POSE] [--gcode-b F --object-b NAME --print-b POSE] + do two parts mate? A and B are meshes (STL / 3MF / OBJ / PLY; a 3MF with several objects: `plate.3mf#ID` picks + one, the build-item transform is applied), each with an optional scale and pose into the assembled frame. Reports the exact interpenetration volume (regions, depth, xyz; a region thinner than + --seat-tol is a seat contact), the exact minimum clearance, and per FEATURE (a thread, ribs, a bore) the radial + clearance read from surface profiles around an axis: the gap at the posed phase and at the best phase of + screwing / clocking, the free play, crest and root clearances, handedness and pitch coherence. With a G-code + per part the same numbers from the printed outer-wall paths (the slicer's compensation included), the + layer-by-layer wall-vs-model deviation, and the registration of the model on the printed object. + Example: 3d assembly pair rim.stl plug.stl --scale-a 2.00493695 --scale-b 2.00493695 --pose-b '{translate: [0, 0, -37.78]}' --a-side inner --b-side outer --pitch 3.007 --a-z=-36.78:-28.28 --b-z=-36.78:-30.98 --screw-to-best Options: --mate NAME restrict to one mate (repeatable). Example: 3d assembly check assembly.yaml --mate lid @@ -109,9 +122,68 @@ Example: 3d assembly position a.yaml --3mf p.3mf --include-non-printable --reason "fit test" --yes --reason TEXT why the override (logged with it). Example: 3d assembly position a.yaml --3mf p.3mf --include-non-printable --reason "fit test" --yes --yes write the 3MF (position is a dry run without it). Example: 3d assembly position a.yaml --3mf p.3mf --yes + --scale-a S / --scale-b S + scale applied to the mesh first (a print scale such as 2.00493695). + Example: 3d assembly pair rim.stl plug.stl --scale-a 2 --scale-b 2 + --pose-a POSE / --pose-b POSE + place the part in the assembled frame: a spec pose, `{scale, rotate (deg), translate}` or a list of them. + Example: 3d assembly pair rim.stl plug.stl --pose-b '{rotate: [0, 0, 90], translate: [0, 0, -10]}' + --axis AXIS the common axis of the features: x | y | z | -z | X,Y,Z (default z). + Example: 3d assembly pair a.stl b.stl --a-side inner --b-side outer --axis x + --axis-point X,Y,Z a point on the axis in the assembled frame (default: centre of A's box). Feature z windows are + the coordinate along the axis. Example: 3d assembly pair a.stl b.stl --a-side inner --b-side outer --axis-point 110,110,0 + --a-side inner|outer / --b-side inner|outer + the feature's surface on each part: inner = the wall of a bore, outer = a boss; one of each. + Example: 3d assembly pair rim.stl plug.stl --a-side inner --b-side outer + --a-z LO:HI / --b-z LO:HI + the feature's z window along the axis (default: the whole part); a negative bound needs `=`. + Example: 3d assembly pair rim.stl plug.stl --a-side inner --b-side outer --a-z 1:9.5 --b-z=-9:-2 + --a-r LO:HI / --b-r LO:HI + the feature's radius window from the axis (default: a bore reads from the axis to the boss's largest + radius + 3 mm). Example: 3d assembly pair rim.stl plug.stl --a-side inner --b-side outer --a-r 0:36.4 --b-r 0:36 + --pitch MM thread pitch of the feature; 0 or omitted = ribs / keys (profile against the angle only). Handedness + is detected and a mismatch FAILs. Example: 3d assembly pair rim.stl plug.stl --a-side inner --b-side outer --pitch 3.007 + --hand auto|right|left + force the handedness (default auto). Example: 3d assembly pair rim.stl plug.stl --a-side inner --b-side outer --pitch 3 --hand right + --feature-name N the feature's name in the report. Example: 3d assembly pair a.stl b.stl --a-side inner --b-side outer --feature-name rim_thread + --features FILE several features from a YAML / JSON list (name, pitch, hand, min_gap_mm, a/b {side, z, r}). + Example: 3d assembly pair rim.stl plug.stl --features features.yaml + --min-radial-gap MM + radial gap every feature must keep at its best phase (default 0.1). + Example: 3d assembly pair rim.stl plug.stl --a-side inner --b-side outer --pitch 3 --min-radial-gap 0.2 + --fit TYPE take the required radial gap from the fit table (clearance - tol of `3d fits resolve`), with --shape + (default round), --printer, --material, --nozzle, --layer. + Example: 3d assembly pair rim.stl plug.stl --a-side inner --b-side outer --fit sliding --printer k1se --material PETG + --shape S / --nozzle MM / --layer MM + the fit lookup's shape, nozzle and layer height. Example: 3d assembly pair a.stl b.stl --fit sliding --shape round --nozzle 0.4 --layer 0.2 + --min-clearance MM the exact minimum gap the parts must keep (0 where parts seat). Example: 3d assembly pair a.stl b.stl --min-clearance 0.3 + --zone BOX measure the minimum gap only among surface samples inside XMIN,YMIN,ZMIN,XMAX,YMAX,ZMAX (a seat + face outside it does not count). Example: 3d assembly pair a.stl b.stl --zone=-40,-40,1,40,40,9.5 + --max-overlap-mm3 V + interpenetration volume that is still allowed (a designed press fit; default 0.01). + Example: 3d assembly pair a.stl b.stl --max-overlap-mm3 4 + --seat-tol MM an overlap thinner than this is a seat contact, not a collision (default 0.02). + Example: 3d assembly pair a.stl b.stl --seat-tol 0.05 + --screw-to-best turn B about the axis to the best phase of the first feature before measuring (a screwed fit has no + unique pose; the turn keeps the axial position). Example: 3d assembly pair rim.stl plug.stl --a-side inner --b-side outer --pitch 3 --screw-to-best + --gcode-a F / --gcode-b F + the sliced G-code that printed the part; adds the as-printed numbers. + Example: 3d assembly pair rim.stl plug.stl --a-side inner --b-side outer --pitch 3 --gcode-b plug.gcode --object-b "Plug x2.stl" + --object-a NAME / --object-b NAME + which object of the G-code (the slicer's object label; a wrong or ambiguous name lists the file's + labels, `"Cup.stl id:0 copy 1"` picks one copy; omitted = the plate's only object). + Example: 3d assembly pair rim.stl plug.stl --gcode-a plate.gcode --object-a "Pot Pepper.stl" + --print-a POSE / --print-b POSE + how the (scaled) model sits on the bed: rotate / scale like a pose. Without a translate it is + dropped on z = 0 and registered on the printed walls in XY; with one it is used as given. + Example: 3d assembly pair rim.stl plug.stl --gcode-a rim.gcode --print-a '{rotate: [180, 0, 0]}' + --skip-below MM layers below this bed height are left out of the wall deviation (first-layer flare; default 0.5). + Example: 3d assembly pair rim.stl plug.stl --gcode-a rim.gcode --skip-below 0.8 + --max-faces N refuse meshes with more triangles (default 400000). Example: 3d assembly pair a.stl b.stl --max-faces 200000 --json one JSON document, schema `3d.assembly/1`. Example: 3d assembly check assembly.yaml --json | jq .summary -Spec: docs/commands/assembly.md (parts, mates, snap hooks, access, zones, loads).""" +Spec: docs/commands/assembly.md (parts, mates, snap hooks, access, zones, loads; `pair` features and G-code). +Feature flags of `pair` take values that start with a minus after `=` or on their own (`--a-z -36:-28` works).""" def run(argv: list[str]) -> int: @@ -136,6 +208,8 @@ def run(argv: list[str]) -> int: command="assembly", remediation=[f"Example: 3d assembly {argv[0]} assembly.yaml"], ) + if argv[0] == "pair": + return run_tool(DEPS, "assembly_pair_cli.py", argv) return run_tool(DEPS, "assembly_cli.py", argv) diff --git a/lib/gcode_walls.py b/lib/gcode_walls.py new file mode 100644 index 00000000..c2d10395 --- /dev/null +++ b/lib/gcode_walls.py @@ -0,0 +1,417 @@ +"""gcode_walls.py — the printed wall lines of one object, layer by layer, straight from a sliced G-code. + +WHAT: streams an Orca / Snapmaker-Orca / Bambu / PrusaSlicer / Marlin G-code once and returns, for the chosen object, + every layer's perimeter lines (default: the outer wall, `WALL_FEATURES`) as straight pieces with the extrusion + width and the layer height the slicer wrote for them. Arcs (G2/G3, centre form) are followed as chords. The object is + named the way the slicer labels it (`gcode_labels`: EXCLUDE_OBJECT, Bambu blocks, M486, PrusaSlicer comments); a name + that matches several copies or none, an omitted name on a plate with several objects, and a file with no such wall at + all are refused with the labels the file has. +WHY: the as-printed fit check (`3d assembly pair --gcode-a ...`) needs where the nozzle really put the wall: the + slicer's compensation (elephant foot, xy_hole_compensation, line width) moves it away from the model, and a thread or + a rib clearance is a few tenths of a millimetre. Object and layer labels come from `gcode_labels` (the one label + reader); the move parsing is kept here on purpose: `gcode_overhangs.parse` merges every feature kind of an object + into one line set, this reader must keep the feature type, the line width and the path height of each piece. +INVARIANTS: + - Coordinates are the bed frame of the G-code. XY is absolute unless G91; E is relative after M83, absolute after M82 + (G92 E resets it); the position is unknown (NaN) until the first absolute move and again after G28. + - A layer starts at the layer marker (Orca / Prusa `;LAYER_CHANGE`, Bambu `; CHANGE_LAYER`); its Z is the `;Z:` / + `; Z_HEIGHT:` value, the TOP of the layer. The first `;HEIGHT:` / `; LAYER_HEIGHT:` after the marker is only the step + since the previous layer of the whole plate; the object's own layer height is the `;HEIGHT:` in effect when it + draws, which is what each piece records. So the mid height `z - h/2`, where the slicer cuts the model, comes from the + object's own first piece of the layer. + - `base_z` is the bottom of the object's lowest layer (`z - h`): 0 on the bed, the raft height plus its gap on a raft. + - A piece with no `;WIDTH:` before it takes the footer's line width for its kind (outer or inner wall). + - Overhanging stretches of the outer loop are tagged `Overhang wall` by Orca / Bambu (inner loops are too, so the tag + cannot just be added to the kinds): their length is counted in `Walls.overhang_mm` so the caller can say how much + of the wall was left out; `kinds` can ask for them. + - A crafted file cannot exhaust memory: arcs with a radius beyond `MAX_ARC_RADIUS_MM`, pieces beyond `MAX_PIECE_MM` + and more than `MAX_PIECES` pieces per file are refused with `ValueError`; R-form arcs are refused too. +Stdlib + numpy only; imported lazily from the geometry worker. +""" + +from __future__ import annotations + +import contextlib +import math +import re +from dataclasses import dataclass, field +from pathlib import Path + +import numpy as np +import numpy.typing as npt + +from gcode_labels import ( + LAYER_MARKERS, + ObjectLabels, + feature, + line_width, + undecorated, +) +from printer_ops.move_bounds import _MOVE, _WORD, Move +from slicing.mix_gcode_faces import DEFAULT_WIDTH_MM, WALL_FEATURES + +MAX_WIDTH_MM = ( + 3.0 # a configured line width beyond this is a percentage or a typo, not a nozzle +) +MAX_ARC_RADIUS_MM = 500.0 +MAX_PIECE_MM = 2000.0 +MAX_PIECES = 3_000_000 +ARC_CHORD_MM = 0.2 +OVERHANG = "Overhang wall" +_BAMBU_BLOCK = "; start printing object, unique label id:" +_CONFIG_WIDTH = re.compile( + r"^;\s*(outer_wall_line_width|inner_wall_line_width|external_perimeter_extrusion_width" + r"|perimeter_extrusion_width)\s*=\s*([\d.]+)\s*$" +) +_OUTER_KEYS = ("outer_wall_line_width", "external_perimeter_extrusion_width") +_KLIPPER_LABEL = re.compile(r"_id_(\d+)_copy_(\d+)$") + +Points = npt.NDArray[np.float64] +Vector = npt.NDArray[np.float64] +Piece = tuple[ + float, float, float, float, float, float +] # x0, y0, x1, y1, width, path height (NaN: unknown) + + +def densify_segments( + a: Points, b: Points, spacing_mm: float +) -> tuple[Points, npt.NDArray[np.intp]]: + """Points every ~`spacing_mm` along the segments a[i] -> b[i] (each start included, each end left to the next + segment) and, per point, the index of the segment it lies on.""" + if len(a) == 0: + return np.zeros((0, a.shape[1] if a.ndim == 2 else 2)), np.zeros( + 0, dtype=np.intp + ) + length = np.linalg.norm(b - a, axis=1) + count = np.maximum(1, np.ceil(length / spacing_mm).astype(int)) + owner = np.repeat(np.arange(len(count)), count) + first = np.cumsum(count) - count + t = ((np.arange(len(owner)) - first[owner]) / count[owner])[:, None] + return a[owner] + (b[owner] - a[owner]) * t, owner + + +@dataclass +class WallLayer: + """One layer of one object: the straight pieces `a -> b` of every chosen feature line and each piece's width.""" + + z: float # the top of the layer on the bed + h: float # the object's own layer height (the path height of its first piece) + a: Points = field(default_factory=lambda: np.zeros((0, 2))) + b: Points = field(default_factory=lambda: np.zeros((0, 2))) + width: Vector = field(default_factory=lambda: np.zeros(0)) + + @property + def mid_z(self) -> float: + """Where the slicer cuts the model for this layer: half a layer below its top.""" + return self.z - self.h / 2.0 + + def points(self, spacing_mm: float) -> tuple[Points, Vector]: + """Centre-line points every ~`spacing_mm` along the pieces, with the width of the piece each lies on.""" + pts, owner = densify_segments(self.a, self.b, spacing_mm) + return pts, self.width[owner] + + +@dataclass +class Walls: + layers: list[WallLayer] + # every object label the file has, as written (id and copy included) + objects: list[str] + # bottom of the object's lowest layer on the bed: 0, or the raft height plus its gap + base_z: float + # length of the chosen wall lines, and of the `Overhang wall` lines that were left out + wall_mm: float + overhang_mm: float + + +def _norm(label: str) -> str: + """Case, `_` / space and the Klipper `_id_N_copy_M` spelling made comparable.""" + label = _KLIPPER_LABEL.sub(r" id:\1 copy \2", label) + return re.sub(r"[\s_]+", " ", label).strip().casefold() + + +def _bare(label: str) -> str: + """The object's own name: the id / copy suffix cut, and Bambu's `NAME — NAME` doubling folded to `NAME`.""" + name = undecorated(label) + parts = [p.strip() for p in name.split("\u2014")] + return parts[0] if len(parts) > 1 and len(set(parts)) == 1 else name + + +def _matches(label: str, wanted: str) -> bool: + """A project object name (or the slicer's full label with its id and copy) against a G-code object label.""" + suffix = label[len(undecorated(label)) :] + candidates = (label, undecorated(label), _bare(label), _bare(label) + suffix) + return _norm(wanted) in {_norm(c) for c in candidates} + + +def _width(text: str) -> float | None: + """A finite line width in (0, MAX_WIDTH_MM], else None (a percentage, zero, `nan`, a typo).""" + try: + v = float(text) + except ValueError: + return None + return v if 0.0 < v <= MAX_WIDTH_MM else None + + +def _number_after_colon(line: str) -> float | None: + """`;Z:0.2` / `; LAYER_HEIGHT: 0.200000` -> the number (linear time, no regex on a hostile line).""" + try: + v = float(line.partition(":")[2].strip()) + except ValueError: + return None + return v if math.isfinite(v) else None + + +@dataclass +class _Reader: + """Mutable state of the pass; one small method per line kind.""" + + wanted: str | None + kinds: tuple[str, ...] + labels: ObjectLabels = field(default_factory=ObjectLabels) + layers: list[WallLayer] = field(default_factory=list) + objects: list[str] = field(default_factory=list) + # distinct (label, Bambu block id) that matched and printed a wall + printed_by: list[tuple[str, int]] = field(default_factory=list) + config_widths: dict[str, float] = field(default_factory=dict) + # NaN until the first `;WIDTH:`; the footer's width fills it in at the end + width: float = math.nan + path_h: float = math.nan + kind: str = "" + block: int = 0 + x: float = math.nan + y: float = math.nan + rel_xy: bool = False + rel_e: bool = False + e: float = 0.0 + seg: list[Piece] = field(default_factory=list) + z: float | None = None + started: bool = False + pieces: int = 0 + wall_mm: float = 0.0 + overhang_mm: float = 0.0 + # the object printing now is the wanted one (cached when the label changes) + owner_ok: bool = False + + def __post_init__(self) -> None: + self.owner_ok = self.wanted is None + + def open_layer(self) -> None: + self.close_layer() + self.started = True + self.z = None + self.path_h = math.nan + self.seg = [] + + def close_layer(self) -> None: + if not self.started or not self.seg or self.z is None: + return + a = np.asarray(self.seg, dtype=float) + known = a[~np.isnan(a[:, 5]), 5] + h = float(known[0]) if len(known) else 0.0 + self.layers.append(WallLayer(self.z, h, a[:, 0:2], a[:, 2:4], a[:, 4])) + + def comment(self, line: str) -> None: + body = line[1:].lstrip() + if line.startswith(LAYER_MARKERS): + self.open_layer() + elif body.startswith(("Z:", "Z_HEIGHT:")): + self.z = _number_after_colon(line) + elif body.startswith(("HEIGHT:", "LAYER_HEIGHT:")): + self.path_h = _number_after_colon(line) or self.path_h + elif (tag := feature(line)) is not None: + self.kind = tag + elif (w := line_width(line)) is not None: + self.width = _width(w) or self.width + elif (m := _CONFIG_WIDTH.match(line)) is not None: + width = _width(m.group(2)) + if width is not None: + self.config_widths[m.group(1)] = width + + def modal(self, up: str) -> bool: + """G28 / G90 / G91 / G92 (homing, positioning mode, extruder reset); True when `up` was one of them.""" + code = up[:3] + if code not in ("G28", "G90", "G91", "G92") or up[3:4].isdigit(): + return False + if code == "G28": + self.x = self.y = math.nan + elif code == "G92": + words = dict(_WORD.findall(up[3:])) + if "E" in words or not words: + self.e = float(words.get("E", 0.0)) + else: + self.rel_xy = code == "G91" + return True + + def extrusion(self, words: dict[str, str]) -> bool: + """Advance the extruder position by an E word; True when the move extrudes.""" + if "E" not in words: + return False + ev = float(words["E"]) + before = self.e + self.e = before + ev if self.rel_e else ev + return (self.e - before) > 1e-6 + + def target(self, words: dict[str, str]) -> tuple[float, float]: + nx, ny = self.x, self.y + if "X" in words: + nx = float(words["X"]) + (self.x if self.rel_xy else 0.0) + if "Y" in words: + ny = float(words["Y"]) + (self.y if self.rel_xy else 0.0) + return nx, ny + + def moves_in_object(self, g: int, nx: float, ny: float) -> bool: + """An extruding move of the wanted object that goes somewhere (or is an arc).""" + return ( + self.started + and self.owner_ok + and not math.isnan(self.x + self.y + nx + ny) + and (nx != self.x or ny != self.y or g in (2, 3)) + ) + + def move(self, line: str) -> None: + up = line.split(";", 1)[0].strip().upper() + if self.modal(up): + return + m = _MOVE.match(up) + if m is None: + return + g = int(m.group(1)) + words = dict(_WORD.findall(up[m.end() :])) + nx, ny = self.target(words) + if self.extrusion(words) and self.moves_in_object(g, nx, ny): + if self.kind in self.kinds: + self.add(g, words, nx, ny) + elif self.kind == OVERHANG: + self.overhang_mm += math.hypot(nx - self.x, ny - self.y) + self.x, self.y = nx, ny + + def chords( + self, g: int, words: dict[str, str], nx: float, ny: float + ) -> list[tuple[float, float, float, float]]: + """The straight pieces of one move: itself, or the chords of a centre-form arc.""" + if g not in (2, 3): + return [(self.x, self.y, nx, ny)] + if "R" in words: + raise ValueError( + "an arc with an R word (radius form) is not supported: slice with I/J arcs" + ) + ci, cj = float(words.get("I", 0.0)), float(words.get("J", 0.0)) + if math.hypot(ci, cj) > MAX_ARC_RADIUS_MM: + raise ValueError( + f"an arc of radius {math.hypot(ci, cj):.0f} mm is not a print move (limit {MAX_ARC_RADIUS_MM:.0f})" + ) + pieces = Move.arc_move(self.x, self.y, nx, ny, ci, cj, g == 2).chords( + ARC_CHORD_MM + ) + return [ + (c.x0, c.y0, c.x1, c.y1) for c in pieces or [Move(self.x, self.y, nx, ny)] + ] + + def add(self, g: int, words: dict[str, str], nx: float, ny: float) -> None: + cur = self.labels.current + if cur is not None and (cur, self.block) not in self.printed_by: + self.printed_by.append((cur, self.block)) + for x0, y0, x1, y1 in self.chords(g, words, nx, ny): + length = math.hypot(x1 - x0, y1 - y0) + if length > MAX_PIECE_MM: + raise ValueError( + f"a wall piece of {length:.0f} mm is not a print move (limit {MAX_PIECE_MM:.0f})" + ) + self.pieces += 1 + if self.pieces > MAX_PIECES: + raise ValueError( + f"more than {MAX_PIECES} wall pieces: not a print G-code" + ) + self.wall_mm += length + self.seg.append((x0, y0, x1, y1, self.width, self.path_h)) + + def label_changed(self) -> None: + cur = self.labels.current + if cur is not None and cur not in self.objects: + self.objects.append(cur) + self.owner_ok = self.wanted is None or ( + cur is not None and _matches(cur, self.wanted) + ) + + def line(self, raw: str) -> None: + line = raw.strip() + if not line: + return + if line.startswith(_BAMBU_BLOCK): + with contextlib.suppress(ValueError): + self.block = int(line[len(_BAMBU_BLOCK) :]) + if self.labels.feed(line): + self.label_changed() + return + c0 = line[0] + if c0 == ";": + self.comment(line) + elif c0 in "Mm": + head = line.split(";", 1)[0].split()[0].upper() + if head in ("M82", "M83"): + self.rel_e = head == "M83" + elif c0 in "Gg": + self.move(line) + + def fallback_width(self) -> float: + """The footer's line width for the kind asked for (outer-wall kinds: the outer width), else the default.""" + outer = any(k in WALL_FEATURES for k in self.kinds) + keys = ( + _OUTER_KEYS + if outer + else ("inner_wall_line_width", "perimeter_extrusion_width") + ) + return next( + (self.config_widths[k] for k in keys if k in self.config_widths), + DEFAULT_WIDTH_MM, + ) + + +def read_walls( + path: Path | str, + obj: str | None = None, + kinds: tuple[str, ...] = tuple(sorted(WALL_FEATURES)), +) -> Walls: + """The `kinds` lines of every layer of `obj`, in file order. + + `obj` is a project object name (`Plug x2.stl`) or the slicer's full label (`Plug x2.stl id:1 copy 0`); None takes + everything the plate prints, which is only accepted for a plate with at most one labelled object. Raises + `KeyError` when no such wall line exists, `ValueError` when the name matches several copies, the plate has several + objects and none was named, or the file is not a print G-code (limits above) - the messages list the file's labels. + """ + reader = _Reader(obj, kinds) + with open(path, encoding="utf-8", errors="replace") as fh: + for raw in fh: + reader.line(raw) + reader.close_layer() + listing = ", ".join(reader.objects) or "none labelled" + if obj is None and len(reader.objects) > 1: + raise ValueError( + f"the plate holds several objects and none was named: {listing}" + ) + if not reader.layers: + what = f"object '{obj}'" if obj is not None else "the file" + raise KeyError( + f"{what} prints no {'/'.join(kinds)} line; objects in the file: {listing}" + ) + if len(reader.printed_by) > 1: + shown = ", ".join( + f"{label} (block {block})" if block else label + for label, block in reader.printed_by + ) + raise ValueError( + f"object '{obj}' matches several copies: {shown}; " + "name one with its id and copy, e.g. 'name id:1 copy 0' (copies of one name in a Bambu file cannot be told apart: " + "slice them on separate plates)" + ) + fallback = reader.fallback_width() + for layer in reader.layers: + layer.width = np.where(np.isnan(layer.width), fallback, layer.width) + base = min(layer.z - layer.h for layer in reader.layers) + return Walls( + reader.layers, + reader.objects, + max(base, 0.0), + reader.wall_mm, + reader.overhang_mm, + ) diff --git a/tests/e2e/test_assembly_pair.py b/tests/e2e/test_assembly_pair.py new file mode 100644 index 00000000..68956b5e --- /dev/null +++ b/tests/e2e/test_assembly_pair.py @@ -0,0 +1,267 @@ +"""Readable e2e stories for `3d assembly pair`: does a plug screw into a rim, on the model and as printed? + +Every story runs the real `bin/3d` on synthetic threads (tests/pair_fixtures.py): the rim's bore is the plug's thread plus a +0.5 mm radial gap at the aligned phase. No printer, slicer or network. +""" + +from __future__ import annotations + +import json +from pathlib import Path +from typing import Any + +import pair_fixtures as PF +import pytest + +from .workflow_helper import run_cli + +FEATURE = [ + "--a-side", + "inner", + "--b-side", + "outer", + "--pitch", + "3", + "--axis-point", + "0,0,0", +] + + +@pytest.fixture(scope="module") +def parts(tmp_path_factory: pytest.TempPathFactory) -> dict[str, Path]: + tmp = tmp_path_factory.mktemp("pair") + bore, boss = PF.thread_pair(tmp, gap=0.5, amplitude=0.8, pitch=3.0) + tight_dir = tmp / "tight" + tight_dir.mkdir() + tight_bore, tight_boss = PF.thread_pair( + tight_dir, gap=0.05, amplitude=0.8, pitch=3.0 + ) + gcode = PF.write_wall_gcode( + tmp / "plate.gcode", + [ + ("Rim", PF.load(bore), (110.0, 90.0)), + ("Plug", PF.load(boss), (60.0, 40.0)), + ], + ) + return { + "bore": bore, + "boss": boss, + "tight_bore": tight_bore, + "tight_boss": tight_boss, + "gcode": gcode, + } + + +def pair(tmp_path: Path, *args: str) -> Any: + return run_cli(tmp_path, "assembly", "pair", *args, timeout=300) + + +def as_json(res: Any) -> dict[str, Any]: + assert res.stdout.strip(), res.stderr + doc: dict[str, Any] = json.loads(res.stdout) + assert doc["schema"] == "3d.assembly/1" and doc["kind"] == "pair" + return doc + + +def help_entry(lines: list[str], flag: str) -> list[str]: + """The help lines of one options entry (or the `pair` subcommand): its header plus the indented continuation.""" + prefix = " pair " if flag == "pair" else " --" + start = next( + (i for i, ln in enumerate(lines) if ln.startswith(prefix) and flag in ln.split(" Example")[0]), + None, + ) + assert start is not None, flag + block = [lines[start]] + for ln in lines[start + 1 :]: + if not ln.startswith(" ") or not ln.strip(): # the next flag / subcommand / blank line ends the entry + break + block.append(ln) + return block + + +def test_help_documents_pair_and_every_flag_with_an_example(tmp_path: Path) -> None: + res = run_cli(tmp_path, "assembly", "--help") + assert res.returncode == 0 + lines = res.stdout.splitlines() + for flag in ( + "pair", "--scale-a", "--scale-b", "--pose-a", "--pose-b", "--axis", "--axis-point", "--a-side", "--b-side", "--a-z", + "--b-z", "--a-r", "--b-r", "--pitch", "--hand", "--feature-name", "--features", "--min-radial-gap", "--fit", "--shape", + "--nozzle", "--layer", "--min-clearance", "--zone", "--max-overlap-mm3", "--seat-tol", "--screw-to-best", "--gcode-a", + "--gcode-b", "--object-a", "--object-b", "--print-a", "--print-b", "--skip-below", "--max-faces", + ): # fmt: skip + assert any("Example: 3d assembly pair" in ln for ln in help_entry(lines, flag)), flag + + +def test_a_rim_and_plug_that_mate_pass_with_the_radial_gap_of_the_thread( + tmp_path: Path, parts: dict[str, Path] +) -> None: + res = pair( + tmp_path, + str(parts["bore"]), + str(parts["boss"]), + *FEATURE, + "--min-radial-gap", + "0.2", + "--json", + ) + assert res.returncode == 0, res.stderr + doc = as_json(res) + assert doc["ok"] and doc["interpenetration"]["volume_mm3"] == 0.0 + cl = doc["features"][0]["model"]["clearance"] + assert cl["best_phase_min_gap_mm"] == pytest.approx(0.5, abs=0.03) + assert cl["free_play"]["value"] == pytest.approx(0.5 / 0.8 * 3.0, abs=0.15) + text = pair(tmp_path, str(parts["bore"]), str(parts["boss"]), *FEATURE).stdout + assert ( + "feature feature (thread, pitch 3.000 mm): PASS" in text + and "right-handed" in text + ) + assert "interpenetration: PASS" in text + + +def test_a_plug_that_is_too_tight_exits_1_and_says_which_gap( + tmp_path: Path, parts: dict[str, Path] +) -> None: + res = pair( + tmp_path, + str(parts["tight_bore"]), + str(parts["tight_boss"]), + *FEATURE, + "--min-radial-gap", + "0.2", + ) + assert res.returncode == 1 + assert ( + "best-phase radial gap" in res.stdout + and "FAILED: feature:feature" in res.stdout + ) + + +def test_a_plug_posed_in_the_wrong_phase_collides_until_screwed_to_best( + tmp_path: Path, parts: dict[str, Path] +) -> None: + args = [ + str(parts["bore"]), + str(parts["boss"]), + *FEATURE, + "--pose-b", + "{rotate: [0, 0, 180]}", + ] + blocked = pair(tmp_path, *args) + assert blocked.returncode == 1 + assert "interpenetration: FAIL" in blocked.stdout and "move B" in blocked.stdout + screwed = pair(tmp_path, *args, "--screw-to-best", "--json") + assert screwed.returncode == 0, screwed.stdout + screwed.stderr + assert as_json(screwed)["screw_to_best"]["applied_to_b"][ + "rotate_deg" + ] == pytest.approx(-180.0, abs=3.0) + + +def test_negative_window_bounds_are_accepted_as_separate_values( + tmp_path: Path, parts: dict[str, Path] +) -> None: + res = pair( + tmp_path, str(parts["bore"]), str(parts["boss"]), *FEATURE, "--pose-b", "{translate: [0, 0, -3]}", + "--a-z", "-3:6", "--b-z", "-3:6", "--json", + ) # fmt: skip (one pitch down: the same thread phase, so it still screws) + assert res.returncode == 0, res.stderr + assert as_json(res)["features"][0]["windows"]["a"]["z"] == [-3.0, 6.0] + + +def test_the_printed_walls_are_read_from_the_gcode_by_object_name( + tmp_path: Path, parts: dict[str, Path] +) -> None: + res = pair( + tmp_path, str(parts["bore"]), str(parts["boss"]), *FEATURE, "--gcode-a", str(parts["gcode"]), "--object-a", "Rim", + "--gcode-b", str(parts["gcode"]), "--object-b", "plug", "--json", + ) # fmt: skip + assert res.returncode == 0, res.stderr + doc = as_json(res) + assert doc["printed"]["a"]["placement"]["shift_applied_mm"] == pytest.approx( + [110.0, 90.0], abs=0.02 + ) + assert doc["printed"]["b"]["placement"]["shift_applied_mm"] == pytest.approx( + [60.0, 40.0], abs=0.02 + ) + assert doc["features"][0]["as_printed"]["clearance"][ + "best_phase_min_gap_mm" + ] == pytest.approx(0.5, abs=0.04) + text = pair( + tmp_path, + str(parts["bore"]), + str(parts["boss"]), + *FEATURE, + "--gcode-a", + str(parts["gcode"]), + "--object-a", + "Rim", + ).stdout + assert ( + "as printed A:" in text + and "wall vs model: max" in text + and "[as printed (A, rest model)]" in text + ) + + +def test_a_wrong_object_name_exits_2_and_lists_the_objects_in_the_gcode( + tmp_path: Path, parts: dict[str, Path] +) -> None: + res = pair( + tmp_path, + str(parts["bore"]), + str(parts["boss"]), + *FEATURE, + "--gcode-a", + str(parts["gcode"]), + "--object-a", + "Nope", + ) + assert res.returncode == 2 + assert "Rim_id_0_copy_0, Plug_id_0_copy_0" in res.stderr and "--object-a" in res.stderr + + +def test_features_come_from_a_file_and_usage_errors_exit_2( + tmp_path: Path, parts: dict[str, Path] +) -> None: + feats = tmp_path / "features.yaml" + feats.write_text( + "features:\n - {name: rim_thread, pitch: 3, a: {side: inner}, b: {side: outer}, min_gap_mm: 0.2}\n" + ) + ok = pair( + tmp_path, + str(parts["bore"]), + str(parts["boss"]), + "--features", + str(feats), + "--axis-point", + "0,0,0", + "--json", + ) + assert ok.returncode == 0, ok.stderr + assert as_json(ok)["features"][0]["name"] == "rim_thread" + for argv, needle in ( + ([str(parts["bore"])], "needs two parts"), + ([str(parts["bore"]), str(parts["boss"]), "--a-side", "inner"], "--b-side"), + ( + [ + str(parts["bore"]), + str(parts["boss"]), + "--features", + str(feats), + "--pitch", + "3", + ], + "cannot be combined", + ), + ( + [str(parts["bore"]), str(parts["boss"]), *FEATURE, "--a-z", "9:1"], + "LO <= HI", + ), + ([str(parts["bore"]), str(tmp_path / "missing.stl")], "missing.stl"), + ( + [str(parts["bore"]), str(parts["boss"]), "--pose-b", "{rotate: [1, 2"], + "--pose-b", + ), + ): + res = pair(tmp_path, *argv) + assert res.returncode == 2, (argv, res.stdout, res.stderr) + assert needle in res.stderr, (argv, res.stderr) diff --git a/tests/pair_fixtures.py b/tests/pair_fixtures.py new file mode 100644 index 00000000..0f6d92bf --- /dev/null +++ b/tests/pair_fixtures.py @@ -0,0 +1,283 @@ +"""Synthetic threaded and ribbed parts (and the G-code that prints them) for the `3d assembly pair` tests. + +A thread here is a surface of revolution whose radius follows a triangle wave of the helical phase +`u = z/pitch - hand * theta/2pi` (so a perfect single-start thread: one radius per phase), built as an exact watertight +triangle mesh. The numbers a test asserts are read off the construction: a bore profile that is the boss profile plus +`gap` at the same phase has radial gap `gap` at the aligned phase and a free fraction of `gap / amplitude` of a turn. +""" + +from __future__ import annotations + +import math +from collections.abc import Callable +from pathlib import Path + +import numpy as np +import numpy.typing as npt +import shapely +import trimesh +from manifold3d import Manifold + +Radius = Callable[ + [npt.NDArray[np.float64], npt.NDArray[np.float64]], npt.NDArray[np.float64] +] + + +def triangle_thread( + base: float, amplitude: float, pitch: float, hand: int = 1, phase: float = 0.0 +) -> Radius: + """r(theta, z) = base + amplitude * tri(u), u = z/pitch - hand*theta/2pi + phase; tri: 0 at u=0, 1 at u=1/2.""" + + def r( + theta: npt.NDArray[np.float64], z: npt.NDArray[np.float64] + ) -> npt.NDArray[np.float64]: + u = (z / pitch - hand * theta / (2 * math.pi) + phase) % 1.0 + return base + amplitude * (1.0 - np.abs(2.0 * u - 1.0)) + + return r + + +def _grid( + radius: Radius, z0: float, z1: float, n_theta: int, n_z: int +) -> tuple[npt.NDArray[np.float64], npt.NDArray[np.int64]]: + theta = np.linspace(0.0, 2 * math.pi, n_theta, endpoint=False) + z = np.linspace(z0, z1, n_z) + tt, zz = np.meshgrid(theta, z) # rows: z, columns: theta + r = radius(tt, zz) + verts = np.column_stack( + [(r * np.cos(tt)).ravel(), (r * np.sin(tt)).ravel(), zz.ravel()] + ) + idx = np.arange(n_z * n_theta).reshape(n_z, n_theta) + a, b = idx[:-1], idx[1:] + a2, b2 = np.roll(a, -1, axis=1), np.roll(b, -1, axis=1) + faces = np.concatenate( + [ + np.column_stack([a.ravel(), a2.ravel(), b2.ravel()]), + np.column_stack([a.ravel(), b2.ravel(), b.ravel()]), + ] + ) + return verts, faces + + +def threaded_boss( + radius: Radius, + z0: float = 0.0, + z1: float = 9.0, + *, + n_theta: int = 240, + steps_per_pitch: int = 12, + pitch: float = 3.0, +) -> trimesh.Trimesh: + """A solid with a threaded outer surface and flat ends.""" + n_z = max(2, round((z1 - z0) / pitch * steps_per_pitch) + 1) + verts, faces = _grid(radius, z0, z1, n_theta, n_z) + bottom, top = len(verts), len(verts) + 1 + verts = np.vstack([verts, [0, 0, z0], [0, 0, z1]]) + ring0 = np.arange(n_theta) + ring1 = (n_z - 1) * n_theta + ring0 + cap0 = np.column_stack([np.full(n_theta, bottom), np.roll(ring0, -1), ring0]) + cap1 = np.column_stack([np.full(n_theta, top), ring1, np.roll(ring1, -1)]) + mesh = trimesh.Trimesh(verts, np.vstack([faces, cap0, cap1]), process=True) + mesh.fix_normals() + return mesh + + +def threaded_bore( + radius: Radius, + z0: float = 0.0, + z1: float = 9.0, + *, + outer: float = 14.0, + n_theta: int = 240, + steps_per_pitch: int = 12, + pitch: float = 3.0, +) -> trimesh.Trimesh: + """A tube (outer radius `outer`, flat ends) whose bore wall follows `radius`.""" + n_z = max(2, round((z1 - z0) / pitch * steps_per_pitch) + 1) + inner, faces = _grid(radius, z0, z1, n_theta, n_z) + theta = np.linspace(0.0, 2 * math.pi, n_theta, endpoint=False) + z = np.linspace(z0, z1, n_z) + tt, zz = np.meshgrid(theta, z) + shell = np.column_stack( + [(outer * np.cos(tt)).ravel(), (outer * np.sin(tt)).ravel(), zz.ravel()] + ) + verts = np.vstack([inner, shell]) + off = len(inner) + ring0 = np.arange(n_theta) + ring1 = (n_z - 1) * n_theta + ring0 + nxt = np.roll(ring0, -1) + caps = [ + np.column_stack([ring0, nxt + off, ring0 + off]), + np.column_stack([ring0, nxt, nxt + off]), + np.column_stack([ring1, ring1 + off, np.roll(ring1, -1) + off]), + np.column_stack([ring1, np.roll(ring1, -1) + off, np.roll(ring1, -1)]), + ] + mesh = trimesh.Trimesh( + verts, + np.vstack([faces, *caps, _grid_outer(n_theta, n_z, off)]), + process=True, + ) + mesh.fix_normals() + return mesh + + +def _grid_outer(n_theta: int, n_z: int, off: int) -> npt.NDArray[np.int64]: + idx = np.arange(n_z * n_theta).reshape(n_z, n_theta) + off + a, b = idx[:-1], idx[1:] + a2, b2 = np.roll(a, -1, axis=1), np.roll(b, -1, axis=1) + return np.concatenate( + [ + np.column_stack([a.ravel(), a2.ravel(), b2.ravel()]), + np.column_stack([a.ravel(), b2.ravel(), b.ravel()]), + ] + ) + + +def load(path: Path) -> trimesh.Trimesh: + mesh = trimesh.load(path, force="mesh") + assert isinstance(mesh, trimesh.Trimesh) + return mesh + + +def write(mesh: trimesh.Trimesh, path: Path) -> Path: + assert mesh.is_watertight and mesh.volume > 0, "fixture mesh must be a closed solid" + mesh.export(path) + return path + + +def thread_pair( + tmp: Path, + *, + gap: float = 0.5, + amplitude: float = 0.8, + pitch: float = 3.0, + boss_hand: int = 1, + bore_hand: int = 1, + bore_extra_phase: float = 0.0, +) -> tuple[Path, Path]: + """(bore.stl, boss.stl): the bore is the boss's profile plus `gap` at the same phase (aligned = best phase).""" + boss_r = triangle_thread(8.0, amplitude, pitch, boss_hand) + bore_r = triangle_thread(8.0 + gap, amplitude, pitch, bore_hand, bore_extra_phase) + bore = write(threaded_bore(bore_r, pitch=pitch), tmp / "bore.stl") + boss = write(threaded_boss(boss_r, pitch=pitch), tmp / "boss.stl") + return bore, boss + + +def write_stl(path: Path, solid: Manifold) -> Path: + mesh = solid.to_mesh() + trimesh.Trimesh(mesh.vert_properties[:, :3], mesh.tri_verts, process=False).export( + path + ) + return path + + +def ribbed_pair( + tmp: Path, *, n: int = 6, rib_width: float = 2.0, groove_width: float = 2.6 +) -> tuple[Path, Path]: + """(bore.stl, boss.stl): a boss (r 10) with `n` radial ribs out to r 11 in a bore (r 10.8) with `n` grooves out to r 12.""" + height = 8.0 + boss = Manifold.cylinder(height, 10.0, 10.0, 96) + groove_cut = Manifold.cylinder(height + 2, 10.8, 10.8, 96).translate((0, 0, -1)) + bore = Manifold.cylinder(height, 14.0, 14.0, 96) - groove_cut + for k in range(n): + ang = 360.0 * k / n + rib = Manifold.cube((11.0, rib_width, height), False).translate( + (0, -rib_width / 2, 0) + ) + groove = Manifold.cube((12.0, groove_width, height + 2), False).translate( + (0, -groove_width / 2, -1) + ) + boss = boss + rib.rotate((0, 0, ang)) + bore = bore - groove.rotate((0, 0, ang)) + return write_stl(tmp / "rib_bore.stl", bore), write_stl(tmp / "rib_boss.stl", boss) + + +def _layer_polygons(mesh: trimesh.Trimesh, layer_h: float) -> list[list[shapely.Polygon]]: + """The model's section polygons at the mid height of every layer (mesh z = 0 on the bed).""" + from trimesh import intersections + + n_layers = round(float(mesh.bounds[1][2]) / layer_h) + heights = (np.arange(n_layers) + 0.5) * layer_h + lines, _to3d, _ = intersections.mesh_multiplane(mesh, np.zeros(3), np.array([0, 0, 1.0]), heights) + return [list(trimesh.load_path(np.asarray(seg)).polygons_full) if len(seg) else [] for seg in lines] + + +def _object_block( + name: str, polys: list[shapely.Polygon], center: tuple[float, float], grow: float, width: float +) -> list[str]: + """One object's outer-wall lines of one layer, between its EXCLUDE_OBJECT markers.""" + slug = name.replace(" ", "_") + cx, cy = center + out = [f"EXCLUDE_OBJECT_START NAME={slug}_id_0_copy_0", ";TYPE:Outer wall", f";WIDTH:{width}"] + for poly in polys: + for ring in _wall_rings(poly, grow, width): + pts = np.asarray(ring.coords) + out.append(f"G1 X{pts[0][0] + cx:.4f} Y{pts[0][1] + cy:.4f} F3000") + out += [f"G1 X{x + cx:.4f} Y{y + cy:.4f} E0.05" for x, y in pts[1:]] + out.append(f"EXCLUDE_OBJECT_END NAME={slug}_id_0_copy_0") + return out + + +def write_wall_gcode( + path: Path, + objects: list[tuple[str, trimesh.Trimesh, tuple[float, float]]], + *, + layer_h: float = 0.2, + width: float = 0.42, + grow: float = 0.0, +) -> Path: + """An Orca-style G-code holding only the OUTER WALL lines of each `(name, mesh, (cx, cy))` object. + + The mesh is placed with its z = 0 on the bed and its axis at (cx, cy); each layer's wall centre line is the model + section inset by half the line width. `grow` offsets the section outward first (a slicer that prints oversize). + """ + out = ["; generated by tests/pair_fixtures.py", "M83"] + per_object = [_layer_polygons(mesh, layer_h) for _n, mesh, _c in objects] + for name, _m, (cx, cy) in objects: + out.append(f"EXCLUDE_OBJECT_DEFINE NAME={name.replace(' ', '_')}_id_0_copy_0 CENTER={cx},{cy}") + for layer in range(max(len(layers) for layers in per_object)): + out += [";LAYER_CHANGE", f";Z:{(layer + 1) * layer_h:.3f}", f";HEIGHT:{layer_h}"] + for (name, _m, center), layers in zip(objects, per_object, strict=True): + if layer < len(layers): + out += _object_block(name, layers[layer], center, grow, width) + path.write_text("\n".join(out) + "\n") + return path + + +def _wall_rings(poly: object, grow: float, width: float) -> list[shapely.LinearRing]: + assert isinstance(poly, shapely.Polygon) + inset = ( + poly.buffer(grow - width / 2.0, join_style="round") + if grow - width / 2.0 + else poly + ) + geoms = list(inset.geoms) if hasattr(inset, "geoms") else [inset] + rings: list[shapely.LinearRing] = [] + for g in geoms: + if g.is_empty: + continue + rings.append(g.exterior) + rings.extend(g.interiors) + return rings + + +def write_3mf(path: Path, objects: list[tuple[int, trimesh.Trimesh]]) -> Path: + """A minimal 3MF (core spec, one build item per object, identity transform) the stdlib reader accepts.""" + import zipfile + + rows = [] + for oid, mesh in objects: + verts = "".join(f'' for x, y, z in mesh.vertices) + tris = "".join(f'' for a, b, c in mesh.faces) + rows.append( + f'{verts}{tris}' + ) + items = "".join(f'' for oid, _ in objects) + xml = ( + '' + '' + f"{''.join(rows)}{items}" + ) + with zipfile.ZipFile(path, "w") as zf: + zf.writestr("3D/3dmodel.model", xml) + return path diff --git a/tests/test_assembly_pair.py b/tests/test_assembly_pair.py new file mode 100644 index 00000000..69489641 --- /dev/null +++ b/tests/test_assembly_pair.py @@ -0,0 +1,660 @@ +"""`assembly.pair` / `assembly.helical` / `assembly.printed`: do two parts mate? Each test has a teeth case that must FAIL. + +The fixtures (tests/pair_fixtures.py) are exact threads and ribs: the bore profile is the boss profile plus `gap` at the same +phase, so the radial gap at the best phase is `gap` and the free fraction of a turn is `gap / amplitude`. The numbers asserted +below are read off those dimensions, not off a previous run of the code under test. +""" + +from __future__ import annotations + +from pathlib import Path +from typing import Any + +import numpy as np +import pair_fixtures as PF +import pytest +from manifold3d import Manifold +from assembly import helical +from assembly import pair as P +from assembly.scene import pose_matrix +from errors import InputNotFound, InvalidArgument + +GAP = 0.5 +AMPLITUDE = 0.8 +PITCH = 3.0 +FACET_MM = ( + 0.02 # the triangle mesh of the fixture is a faceted surface: radii agree to ~ this +) + + +def thread_feature(**kw: Any) -> P.Feature: + return P.Feature( + "thread", P.FeatureSide("inner"), P.FeatureSide("outer"), pitch_mm=PITCH, **kw + ) + + +def pose(text: str) -> Any: + return P.pose_from_text(text, "pose") + + +def run(a: Path, b: Path, features: list[P.Feature], **kw: Any) -> dict[str, Any]: + pose_b = kw.pop("pose_b", ()) + gcode_a, gcode_b = kw.pop("gcode_a", None), kw.pop("gcode_b", None) + part_a = P.PartInput(a, gcode=gcode_a, gcode_object="bore" if gcode_a else None) + part_b = P.PartInput( + b, pose=pose_b, gcode=gcode_b, gcode_object="boss" if gcode_b else None + ) + return P.run_pair( + P.PairOptions(part_a, part_b, features, axis_point=(0.0, 0.0, 0.0), **kw) + ) + + +@pytest.fixture(scope="module") +def threads(tmp_path_factory: pytest.TempPathFactory) -> tuple[Path, Path]: + return PF.thread_pair( + tmp_path_factory.mktemp("threads"), gap=GAP, amplitude=AMPLITUDE, pitch=PITCH + ) + + +# ---- the radial clearance of a thread ------------------------------------------------------------------------------ + + +def test_thread_radial_gap_free_play_and_handedness_match_the_construction( + threads: tuple[Path, Path], +) -> None: + doc = run(*threads, [thread_feature()]) + feat = doc["features"][0] + model = feat["model"] + cl = model["clearance"] + assert cl["at_pose_min_gap_mm"] == pytest.approx(GAP, abs=FACET_MM) + assert cl["best_phase_min_gap_mm"] == pytest.approx(GAP, abs=FACET_MM) + # slack of the construction: |d| < gap / (2 * amplitude) turns either side -> gap / amplitude of a turn + assert cl["free_play"]["fraction_of_turn"] == pytest.approx( + GAP / AMPLITUDE, abs=0.03 + ) + assert cl["free_play"]["value"] == pytest.approx(GAP / AMPLITUDE * PITCH, abs=0.1) + assert cl["free_play"]["unit"] == "mm axial" + for side in ("a", "b"): + assert model["profiles"][side]["handedness"]["hand"] == "right" + assert model["profiles"][side]["coverage"] == 1.0 + # crest / root: the bore is the boss + gap everywhere + assert cl["male_crest_to_female_root_mm"] == pytest.approx(GAP, abs=FACET_MM) + assert cl["female_crest_to_male_root_mm"] == pytest.approx(GAP, abs=FACET_MM) + # the crest is a sharp ridge the mesh grid samples every pitch/12: its tip is cut by up to a few hundredths + assert model["profiles"]["b"]["r_max_mm"] == pytest.approx(8.0 + AMPLITUDE, abs=0.06) + assert doc["ok"] and doc["interpenetration"]["volume_mm3"] == 0.0 + + +def test_a_bore_window_is_found_without_being_told(threads: tuple[Path, Path]) -> None: + """No radius window on the tube: its own outer wall (r 14) must not be read as part of its bore (r ~ 8.5-9.3).""" + doc = run(*threads, [thread_feature()]) + a = doc["features"][0]["model"]["profiles"]["a"] + assert a["r_max_mm"] < 10.0, a + + +def test_thread_too_tight_for_the_required_gap_fails(tmp_path: Path) -> None: + (tmp_path / "tight").mkdir() + (tmp_path / "ok").mkdir() + bore, boss = PF.thread_pair(tmp_path / "tight", gap=0.05, amplitude=AMPLITUDE, pitch=PITCH) + doc = run(bore, boss, [thread_feature()], min_radial_gap_mm=0.2) + assert not doc["ok"] + feat = doc["features"][0] + assert feat["status"] == "fail" + assert any("best-phase radial gap" in r for r in feat["model"]["reasons"]) + # the same parts with a 0.5 mm gap pass the same requirement + bore2, boss2 = PF.thread_pair(tmp_path / "ok", gap=GAP, pitch=PITCH) + assert run(bore2, boss2, [thread_feature()], min_radial_gap_mm=0.2)["ok"] + + + +def test_a_pose_in_the_wrong_phase_interferes_and_screw_to_best_clears_it( + threads: tuple[Path, Path], +) -> None: + half_turn = pose("{rotate: [0, 0, 180]}") + blocked = run(*threads, [thread_feature()], pose_b=half_turn) + assert ( + blocked["interpenetration"]["volume_mm3"] > 1.0 + ) # the flanks occupy the same space + cl = blocked["features"][0]["model"]["clearance"] + assert cl["at_pose_min_gap_mm"] < 0 < cl["best_phase_min_gap_mm"] + assert any( + "posed phase interferes" in r + for r in blocked["features"][0]["model"]["reasons"] + ) + assert abs(abs(cl["best_phase_move_b"]["rotate_deg"]) - 180.0) < 3.0 + # the advance equivalent to that turn: half a pitch + assert abs(cl["best_phase_move_b"]["axial_mm"]) == pytest.approx( + PITCH / 2, abs=0.05 + ) + + screwed = run(*threads, [thread_feature()], pose_b=half_turn, screw_to_best=True) + assert screwed["interpenetration"]["volume_mm3"] == 0.0 + assert screwed["features"][0]["model"]["clearance"][ + "at_pose_min_gap_mm" + ] == pytest.approx(GAP, abs=FACET_MM) + assert screwed["ok"] + + +def test_screw_to_best_turns_about_the_axis_and_keeps_the_axial_position( + threads: tuple[Path, Path], +) -> None: + doc = run( + *threads, + [thread_feature()], + pose_b=pose("{rotate: [0, 0, 90]}"), + screw_to_best=True, + ) + bbox = doc["parts"]["b"]["bbox_assembled"] + assert bbox["min"][2] == pytest.approx(0.0, abs=1e-3) and bbox["max"][ + 2 + ] == pytest.approx(9.0, abs=1e-3) + # a quarter turn off is neither its own mirror image nor a symmetry of the thread: only the RIGHT turn clears it + assert doc["screw_to_best"]["applied_to_b"]["rotate_deg"] == pytest.approx(-90.0, abs=3.0) + assert doc["interpenetration"]["volume_mm3"] == 0.0 + assert doc["features"][0]["model"]["clearance"]["at_pose_min_gap_mm"] == pytest.approx(GAP, abs=FACET_MM) + assert doc["screw_to_best"]["applied_to_b"]["rotate_deg"] != 0.0 + + +def test_opposite_handedness_cannot_screw_together(tmp_path: Path) -> None: + bore, boss = PF.thread_pair( + tmp_path, gap=GAP, pitch=PITCH, boss_hand=1, bore_hand=-1 + ) + doc = run(bore, boss, [thread_feature()]) + assert not doc["ok"] + reasons = doc["features"][0]["model"]["reasons"] + assert any("handedness mismatch" in r for r in reasons), reasons + assert ( + "clearance" not in doc["features"][0]["model"] + ) # no phase to slide: they never mate + + +def test_a_wrong_pitch_is_reported_not_trusted(threads: tuple[Path, Path]) -> None: + wrong = P.Feature( + "thread", P.FeatureSide("inner"), P.FeatureSide("outer"), pitch_mm=2.0 + ) + doc = run(*threads, [wrong]) + model = doc["features"][0]["model"] + assert model["status"] != "pass", ( + "a pitch that does not describe the thread must not pass silently" + ) + assert any("not decisive" in r or "measured" in r for r in model["reasons"]), model[ + "reasons" + ] + + +def test_a_window_without_surface_is_refused(threads: tuple[Path, Path]) -> None: + feat = P.Feature( + "thread", + P.FeatureSide("inner", z=(50.0, 60.0)), + P.FeatureSide("outer"), + pitch_mm=PITCH, + ) + with pytest.raises(InvalidArgument, match="surface points"): + run(*threads, [feat]) + + +def test_two_bores_do_not_make_a_feature(threads: tuple[Path, Path]) -> None: + feat = P.Feature( + "thread", P.FeatureSide("inner"), P.FeatureSide("inner"), pitch_mm=PITCH + ) + with pytest.raises(InvalidArgument, match="inner"): + run(*threads, [feat]) + + +# ---- ribs and keys ------------------------------------------------------------------------------------------------ + + +def test_ribs_in_wider_grooves_have_a_small_free_turn(tmp_path: Path) -> None: + bore, boss = PF.ribbed_pair(tmp_path, rib_width=2.0, groove_width=2.6) + doc = run( + bore, boss, [P.Feature("keys", P.FeatureSide("inner"), P.FeatureSide("outer"))] + ) + cl = doc["features"][0]["model"]["clearance"] + assert doc["features"][0]["kind"] == "ribs" + assert cl["free_play"]["unit"] == "deg" + assert ( + 2.0 < cl["free_play"]["value"] < 8.0 + ) # 0.6 mm of slack on a radius of ~11 mm is ~3 deg each way + assert cl["best_phase_min_gap_mm"] > 0.5 + assert doc["ok"] + + +def test_ribs_wider_than_their_grooves_never_clear(tmp_path: Path) -> None: + bore, boss = PF.ribbed_pair(tmp_path, rib_width=2.0, groove_width=1.8) + doc = run( + bore, boss, [P.Feature("keys", P.FeatureSide("inner"), P.FeatureSide("outer"))] + ) + assert not doc["ok"] + assert doc["features"][0]["model"]["clearance"]["best_phase_min_gap_mm"] < 0 + assert doc["features"][0]["model"]["clearance"]["free_play"]["value"] == 0.0 + + +def test_rib_clocking_is_recovered_by_screw_to_best(tmp_path: Path) -> None: + bore, boss = PF.ribbed_pair(tmp_path) + feat = P.Feature("keys", P.FeatureSide("inner"), P.FeatureSide("outer")) + off = run( + bore, boss, [feat], pose_b=pose("{rotate: [0, 0, 20]}") + ) # a third of the rib pitch: the ribs hit the bore wall + assert off["interpenetration"]["volume_mm3"] > 1.0 + on = run( + bore, boss, [feat], pose_b=pose("{rotate: [0, 0, 20]}"), screw_to_best=True + ) + assert on["interpenetration"]["volume_mm3"] == 0.0 and on["ok"] + assert on["screw_to_best"]["applied_to_b"]["rotate_deg"] == pytest.approx(-20.0, abs=1.5) + + +# ---- interpenetration and minimum clearance ------------------------------------------------------------------------- + + +def test_overlap_volume_and_minimum_gap_of_two_boxes(tmp_path: Path) -> None: + a = PF.write_stl(tmp_path / "a.stl", Manifold.cube((10, 10, 10), False)) + b = PF.write_stl(tmp_path / "b.stl", Manifold.cube((10, 10, 10), False)) + gap = P.run_pair( + P.PairOptions( + P.PartInput(a), P.PartInput(b, pose=pose("{translate: [12.5, 0, 0]}")) + ) + ) + assert gap["interpenetration"]["volume_mm3"] == 0.0 + assert gap["clearance"]["min_gap_mm"] == pytest.approx(2.5, abs=1e-3) + lap = P.run_pair( + P.PairOptions( + P.PartInput(a), P.PartInput(b, pose=pose("{translate: [9, 0, 0]}")) + ) + ) + assert lap["interpenetration"]["volume_mm3"] == pytest.approx( + 100.0, abs=0.01 + ) # 1 x 10 x 10 + assert lap["interpenetration"]["regions"][0]["depth_mm"] == pytest.approx( + 1.0, abs=0.01 + ) + assert not lap["ok"] + assert lap["clearance"]["min_gap_mm"] == 0.0 + allowed = P.run_pair( + P.PairOptions( + P.PartInput(a), + P.PartInput(b, pose=pose("{translate: [9, 0, 0]}")), + max_overlap_mm3=150.0, + ) + ) + assert allowed["ok"], ( + "a designed interference is allowed only when the caller says how much" + ) + + +def test_a_seat_face_is_a_contact_not_an_interpenetration(tmp_path: Path) -> None: + a = PF.write_stl(tmp_path / "a.stl", Manifold.cube((10, 10, 10), False)) + seated = P.run_pair( + P.PairOptions( + P.PartInput(a), P.PartInput(a, pose=pose("{translate: [9.995, 0, 0]}")) + ) + ) + assert seated["ok"] and seated["interpenetration"]["contacts"]["count"] == 1 + assert seated["interpenetration"]["contacts"]["max_depth_mm"] < 0.02 + sunk = P.run_pair( + P.PairOptions( + P.PartInput(a), P.PartInput(a, pose=pose("{translate: [9.9, 0, 0]}")) + ) + ) + assert not sunk["ok"], "0.1 mm deep is a collision, not a seat" + + +def test_required_minimum_clearance_is_enforced(tmp_path: Path) -> None: + a = PF.write_stl(tmp_path / "a.stl", Manifold.cube((10, 10, 10), False)) + far = P.PairOptions( + P.PartInput(a), + P.PartInput(a, pose=pose("{translate: [10.25, 0, 0]}")), + min_clearance_mm=0.3, + ) + assert P.run_pair(far)["clearance"]["status"] == "fail" + ok = P.PairOptions( + P.PartInput(a), + P.PartInput(a, pose=pose("{translate: [10.4, 0, 0]}")), + min_clearance_mm=0.3, + ) + assert P.run_pair(ok)["ok"] + + +def test_a_zone_measures_the_surfaces_inside_it_only(tmp_path: Path) -> None: + """Two boxes touch on a seat at z=0..10 but are 1 mm apart higher up: the zone above the seat sees the 1 mm.""" + base = Manifold.cube((10, 10, 10), False) + stepped = Manifold.cube((10, 10, 10), False) - Manifold.cube( + (1, 10, 5), False + ).translate((0, 0, 5)) + a = PF.write_stl(tmp_path / "a.stl", base) + b = PF.write_stl(tmp_path / "b.stl", stepped) + pose_b = pose("{translate: [10, 0, 0]}") + whole = P.run_pair(P.PairOptions(P.PartInput(a), P.PartInput(b, pose=pose_b))) + assert whole["clearance"]["min_gap_mm"] == 0.0 + upper = P.run_pair( + P.PairOptions( + P.PartInput(a), P.PartInput(b, pose=pose_b), zone=(0, 0, 6, 12, 10, 10) + ) + ) + assert upper["clearance"]["zone"]["min_gap_mm"] == pytest.approx(1.0, abs=0.05) + + +# ---- the fit table ---------------------------------------------------------------------------------------------------- + + +def test_the_required_radial_gap_comes_from_the_fit_table( + monkeypatch: pytest.MonkeyPatch, tmp_path: Path +) -> None: + monkeypatch.setenv("THREED_FITS_CORRECTIONS", str(tmp_path / "none.yaml")) + monkeypatch.setenv("XDG_CONFIG_HOME", str(tmp_path / "config")) + import fits_sources + + sliding = P.fit_requirement("sliding", "round", None, None, None, None) + press = P.fit_requirement("press", "round", None, None, None, None) + res = fits_sources.resolve("sliding", "round", command="assembly") + assert sliding["required_gap_mm"] == round(res.clearance_mm - res.tol_mm, 4) + assert sliding["required_gap_mm"] > press["required_gap_mm"], ( + "a sliding fit asks for more gap than a press fit" + ) + + +# ---- the as-printed side ---------------------------------------------------------------------------------------------- + + +@pytest.fixture(scope="module") +def printed(tmp_path_factory: pytest.TempPathFactory) -> dict[str, Path]: + tmp = tmp_path_factory.mktemp("printed") + bore, boss = PF.thread_pair(tmp, gap=GAP, amplitude=AMPLITUDE, pitch=PITCH) + bore_m, boss_m = PF.load(bore), PF.load(boss) + faithful = PF.write_wall_gcode( + tmp / "faithful.gcode", + [("bore", bore_m, (110.0, 90.0)), ("boss", boss_m, (60.0, 40.0))], + ) + oversize = PF.write_wall_gcode( + tmp / "oversize.gcode", + [("bore", bore_m, (110.0, 90.0)), ("boss", boss_m, (60.0, 40.0))], + grow=0.15, + ) + return {"bore": bore, "boss": boss, "faithful": faithful, "oversize": oversize} + + +def test_a_faithful_print_reproduces_the_model_and_is_registered_on_its_bed_position( + printed: dict[str, Path], +) -> None: + doc = run( + printed["bore"], + printed["boss"], + [thread_feature()], + gcode_a=printed["faithful"], + gcode_b=printed["faithful"], + ) + info = doc["printed"] + assert info["a"]["placement"]["shift_applied_mm"] == pytest.approx( + [110.0, 90.0], abs=0.02 + ) + assert info["b"]["placement"]["shift_applied_mm"] == pytest.approx( + [60.0, 40.0], abs=0.02 + ) + assert ( + info["a"]["deviation"]["max_mm"] < 0.03 + and info["b"]["deviation"]["max_mm"] < 0.03 + ) + feat = doc["features"][0] + printed_cl = feat["as_printed"]["clearance"] + assert printed_cl["best_phase_min_gap_mm"] == pytest.approx(GAP, abs=0.03) + # the printed radius sits on the model's: the half line width is added back the right way round + for side in ("a", "b"): + diff = feat["printed_vs_model"][side]["radius_printed_minus_model_mm"] + assert abs(diff["mean"]) < 0.02, (side, diff) + + +def test_an_oversize_print_closes_the_gap_the_model_has( + printed: dict[str, Path], +) -> None: + """The slicer prints every wall 0.15 mm into the bore: the model's 0.5 mm gap is 0.35 mm on the plastic.""" + doc = run( + printed["bore"], + printed["boss"], + [ + P.Feature( + "thread", + P.FeatureSide("inner"), + P.FeatureSide("outer"), + pitch_mm=PITCH, + min_gap_mm=0.4, + ) + ], + gcode_a=printed["oversize"], + gcode_b=printed["oversize"], + ) + feat = doc["features"][0] + assert feat["model"]["status"] == "pass", ( + "the model alone clears the 0.4 mm requirement" + ) + printed_cl = feat["as_printed"]["clearance"] + assert printed_cl["best_phase_min_gap_mm"] == pytest.approx( + GAP - 2 * 0.15, abs=0.04 + ) # both walls grow 0.15 + assert feat["as_printed"]["status"] == "fail" and not doc["ok"] + assert doc["printed"]["a"]["deviation"]["max_mm"] == pytest.approx(0.15, abs=0.03) + + +def test_one_printed_part_is_judged_against_the_model_of_the_other( + printed: dict[str, Path], +) -> None: + doc = run( + printed["bore"], + printed["boss"], + [thread_feature()], + gcode_a=printed["oversize"], + ) + feat = doc["features"][0] + assert feat["as_printed_parts"] == ["a"] + assert feat["as_printed"]["clearance"]["best_phase_min_gap_mm"] == pytest.approx( + GAP - 0.15, abs=0.04 + ) + + +def test_an_object_name_that_the_gcode_does_not_have_names_the_ones_it_has( + printed: dict[str, Path], +) -> None: + part = P.PartInput(printed["bore"], gcode=printed["faithful"], gcode_object="nope") + with pytest.raises(InvalidArgument, match="bore"): + P.run_pair( + P.PairOptions( + part, + P.PartInput(printed["boss"]), + [thread_feature()], + axis_point=(0.0, 0.0, 0.0), + ) + ) + + +def test_a_mirrored_model_does_not_register_and_the_deviation_says_so( + printed: dict[str, Path], tmp_path: Path +) -> None: + """A model flipped against its G-code (as when the wrong STL is paired) leaves a large residual, not a quiet pass.""" + mirrored = PF.load(printed["boss"]) + mirrored.apply_transform( + np.diag([1.0, 1.0, -1.0, 1.0]) + ) # z flip: the thread phase runs the other way + mirrored.invert() + path = tmp_path / "boss_flipped.stl" + mirrored.export(path) + doc = run( + printed["bore"], + path, + [ + P.Feature( + "thread", P.FeatureSide("inner"), P.FeatureSide("outer"), pitch_mm=PITCH + ) + ], + gcode_b=printed["faithful"], + ) + assert doc["printed"]["b"]["deviation"]["max_mm"] > 0.1 + + +# ---- helpers -------------------------------------------------------------------------------------------------------------- + + +def test_features_file_round_trip_and_errors(tmp_path: Path) -> None: + good = tmp_path / "f.yaml" + good.write_text( + "features:\n" + " - name: thread\n pitch: 3.0\n hand: right\n min_gap_mm: 0.2\n" + " a: {side: inner, z: [1, 9], r: [0, 12]}\n b: {side: outer, z: [0, 8]}\n" + ) + feats = P.features_from_file(good) + assert ( + feats[0].pitch_mm == 3.0 + and feats[0].a.z == (1.0, 9.0) + and feats[0].b.r is None + and feats[0].min_gap_mm == 0.2 + ) + for body, needle in ( + ("features:\n - {name: x, a: {side: inner}, b: {side: sideways}}\n", "side"), + ( + "features:\n - {name: x, bogus: 1, a: {side: inner}, b: {side: outer}}\n", + "bogus", + ), + ( + "features:\n - {name: x, hand: up, a: {side: inner}, b: {side: outer}}\n", + "hand", + ), + ("features: []\n", "at least one"), + ): + bad = tmp_path / "bad.yaml" + bad.write_text(body) + with pytest.raises(InvalidArgument, match=needle): + P.features_from_file(bad) + with pytest.raises(InputNotFound): + P.features_from_file(tmp_path / "missing.yaml") + + +def test_ranges_and_inline_poses_are_validated() -> None: + assert P.parse_range("-36.78:-28.28", "--a-z") == (-36.78, -28.28) + for bad in ("1", "2:1", "a:b"): + with pytest.raises(InvalidArgument): + P.parse_range(bad, "--a-z") + m = pose_matrix(P.pose_from_text("{scale: 2, translate: [1, 2, 3]}", "--pose-b")) + assert m[:3, 3].tolist() == [1.0, 2.0, 3.0] and m[0, 0] == 2.0 + with pytest.raises(InvalidArgument): + P.pose_from_text("{translate: [1, 2", "--pose-b") + + +def test_the_helical_profile_recovers_a_known_thread_and_its_handedness() -> None: + theta = np.random.default_rng(0).uniform(-np.pi, np.pi, 60000) + z = np.random.default_rng(1).uniform(0, 12, 60000) + r = PF.triangle_thread(10.0, 0.5, 3.0, hand=-1)(theta, z) + pts = np.column_stack([r * np.cos(theta), r * np.sin(theta), z]) + prof, info = helical.pick_handedness(pts, 3.0, "auto") + assert info["hand"] == "left" and info["decided"] + assert prof.coherence_mm < 0.01 and prof.coverage == 1.0 + wrong, _ = helical.pick_handedness(pts, 3.0, "right") + assert wrong.coherence_mm > 10 * prof.coherence_mm, ( + "the wrong handedness must look incoherent" + ) + + +def test_sliding_gap_peaks_at_the_aligned_phase_and_shift_to_motion_has_the_right_sign() -> ( + None +): + n = helical.BINS + u = np.arange(n) / n + male = 8.0 + 0.8 * (1 - np.abs(2 * u - 1)) + female = male + 0.5 + by = helical.sliding(female, male) + assert int(np.argmax(by)) == 0 and by[0] == pytest.approx(0.5) + shifted = np.roll(male, 40) # the male profile moved to higher phase by 40 bins + cl = helical.clearance( + helical.Profile(female, np.ones(n), 0.0, 3.0, 1), + helical.Profile(shifted, np.ones(n), 0.0, 3.0, 1), + ) + assert cl.best_shift_bins == n - 40 # rolling it back by 40 (= -40) realigns it + motion = helical.shift_to_motion(cl.best_shift_bins, 3.0, +1) + assert motion["axial_mm"] == pytest.approx( + -40 / n * 3.0, abs=1e-3 + ) # move it back by 40/360 of a pitch + assert motion["rotate_deg"] == pytest.approx( + +40.0, abs=1e-3 + ) # right-handed: +turn = -advance + + +def test_a_part_printed_upside_down_is_mapped_back_through_its_print_pose( + printed: dict[str, Path], tmp_path: Path +) -> None: + """The boss was printed flipped and turned (rotate 180 about x, 40 about z; lowest point on the bed); `--print-b` says so + and the printed surface still lands on the model's: the bed -> model -> assembly chain is inverted correctly, not + merely identity. (A flip alone would not do: the triangle thread is symmetric under it.)""" + flipped = PF.load(printed["boss"]) + flipped.apply_transform(pose_matrix(P.pose_from_text("{rotate: [180, 0, 40]}", "pose"))) + flipped.apply_translation([0.0, 0.0, -float(flipped.bounds[0][2])]) + gcode = PF.write_wall_gcode(tmp_path / "flipped.gcode", [("boss", flipped, (60.0, 40.0))]) + part_b = P.PartInput( + printed["boss"], + print_pose=P.pose_from_text("{rotate: [180, 0, 40]}", "--print-b"), + gcode=gcode, + gcode_object="boss", + ) + doc = P.run_pair( + P.PairOptions( + P.PartInput(printed["bore"]), part_b, [thread_feature()], axis_point=(0.0, 0.0, 0.0) + ) + ) + feat = doc["features"][0] + assert doc["printed"]["b"]["deviation"]["max_mm"] < 0.03 + diff = feat["printed_vs_model"]["b"]["radius_printed_minus_model_mm"] + assert diff is not None and abs(diff["mean"]) < 0.02 and abs(diff["max"]) < 0.08, diff + assert feat["as_printed"]["clearance"]["best_phase_min_gap_mm"] == pytest.approx(GAP, abs=0.04) + # a print pose that misses the 40 degree turn does not describe these walls: the deviation says so + wrong = P.PartInput( + printed["boss"], + print_pose=P.pose_from_text("{rotate: [180, 0, 0]}", "--print-b"), + gcode=gcode, + gcode_object="boss", + ) + bad = P.run_pair( + P.PairOptions(P.PartInput(printed["bore"]), wrong, [thread_feature()], axis_point=(0.0, 0.0, 0.0)) + ) + assert bad["printed"]["b"]["deviation"]["max_mm"] > 0.1 + + +def test_an_explicit_print_translate_is_used_as_given(printed: dict[str, Path]) -> None: + part_b = P.PartInput( + printed["boss"], + print_pose=P.pose_from_text("{translate: [60, 40, 0]}", "--print-b"), + gcode=printed["faithful"], + gcode_object="boss", + ) + doc = P.run_pair( + P.PairOptions(P.PartInput(printed["bore"]), part_b, [thread_feature()], axis_point=(0.0, 0.0, 0.0)) + ) + placement = doc["printed"]["b"]["placement"] + assert placement["registered"] is False and placement["residual_rms_mm"] is None + assert doc["printed"]["b"]["deviation"]["max_mm"] < 0.03 + # a wrong explicit position is a large deviation, not a silent registration + off = P.PartInput( + printed["boss"], + print_pose=P.pose_from_text("{translate: [65, 40, 0]}", "--print-b"), + gcode=printed["faithful"], + gcode_object="boss", + ) + bad = P.run_pair( + P.PairOptions(P.PartInput(printed["bore"]), off, [thread_feature()], axis_point=(0.0, 0.0, 0.0)) + ) + assert bad["printed"]["b"]["deviation"]["max_mm"] > 1.0 + + +def test_a_3mf_part_is_read_without_lxml_one_object_by_id(threads: tuple[Path, Path], tmp_path: Path) -> None: + bore, boss = (PF.load(p) for p in threads) + plate = PF.write_3mf(tmp_path / "plate.3mf", [(1, bore), (2, boss)]) + path, oid = P.split_part_path(f"{plate}#2") + assert path == plate and oid == 2 + assert P.split_part_path("a#b.stl") == (Path("a#b.stl"), None) + doc = P.run_pair( + P.PairOptions( + P.PartInput(threads[0]), P.PartInput(path, object_id=oid), [thread_feature()], axis_point=(0.0, 0.0, 0.0) + ) + ) + assert doc["features"][0]["model"]["clearance"]["best_phase_min_gap_mm"] == pytest.approx(GAP, abs=FACET_MM) + # two objects and no id: refused, with the ids to pick from + with pytest.raises(InvalidArgument, match="2 objects") as err: + P.run_pair(P.PairOptions(P.PartInput(threads[0]), P.PartInput(plate))) + rendered = err.value.render() + assert "#1:" in rendered and "#2:" in rendered, rendered + with pytest.raises(InvalidArgument, match="#9"): + P.run_pair(P.PairOptions(P.PartInput(threads[0]), P.PartInput(plate, object_id=9))) diff --git a/tests/test_gcode_walls.py b/tests/test_gcode_walls.py new file mode 100644 index 00000000..45cb150d --- /dev/null +++ b/tests/test_gcode_walls.py @@ -0,0 +1,296 @@ +"""`gcode_walls`: the outer-wall lines of ONE object of a sliced G-code, read from the slicer's own labels. + +Each test writes a tiny G-code whose answer is readable off the text; the teeth tests (wrong object, copies, other feature +types, absolute / relative extrusion, G91, G92, arcs, mid-layer heights) fail if the reader attributes a line to the +wrong object, feature or height. +""" + +from __future__ import annotations + +from pathlib import Path + +import numpy as np +import pytest +from gcode_walls import read_walls + +DATA = Path(__file__).parent / "data" + +TWO_OBJECTS = """\ +; outer_wall_line_width = 0.45 +EXCLUDE_OBJECT_DEFINE NAME=Pot_Pepper.stl_id_0_copy_0 CENTER=10,10 +EXCLUDE_OBJECT_DEFINE NAME=Plug_x2.stl_id_1_copy_0 CENTER=50,10 +M83 +G1 X0 Y0 F6000 +;LAYER_CHANGE +;Z:0.2 +;HEIGHT:0.2 +EXCLUDE_OBJECT_START NAME=Pot_Pepper.stl_id_0_copy_0 +;TYPE:Outer wall +;WIDTH:0.42 +G1 X10 Y10 F3000 +G1 X20 Y10 E0.5 +G1 X20 Y20 E0.5 +;TYPE:Inner wall +;WIDTH:0.45 +G1 X15 Y15 E0.5 +EXCLUDE_OBJECT_END NAME=Pot_Pepper.stl_id_0_copy_0 +EXCLUDE_OBJECT_START NAME=Plug_x2.stl_id_1_copy_0 +;TYPE:Outer wall +;WIDTH:0.40 +G1 X50 Y10 F3000 +G1 X60 Y10 E0.5 +EXCLUDE_OBJECT_END NAME=Plug_x2.stl_id_1_copy_0 +;LAYER_CHANGE +;Z:0.4 +;HEIGHT:0.2 +EXCLUDE_OBJECT_START NAME=Pot_Pepper.stl_id_0_copy_0 +;TYPE:Outer wall +;WIDTH:0.42 +G1 X10 Y10 F3000 +G1 X20 Y10 E0.5 +G1 X20 Y10 E0.5 +G1 X30 Y10 E0 +EXCLUDE_OBJECT_END NAME=Pot_Pepper.stl_id_0_copy_0 +""" + + +def write(tmp: Path, text: str, name: str = "t.gcode") -> Path: + p = tmp / name + p.write_text(text) + return p + + +def layer_gcode(body: str, head: str = ";Z:0.2\n;HEIGHT:0.2\n") -> str: + return f";LAYER_CHANGE\n{head};TYPE:Outer wall\n;WIDTH:0.4\n{body}" + + +@pytest.fixture +def gcode(tmp_path: Path) -> Path: + return write(tmp_path, TWO_OBJECTS) + + +def test_only_the_named_objects_outer_wall_is_returned(gcode: Path) -> None: + walls = read_walls( + gcode, "Pot Pepper.stl" + ) # project name with a space matches the EXCLUDE_OBJECT name with `_` + assert [lyr.z for lyr in walls.layers] == [0.2, 0.4] + first = walls.layers[0] + # two outer-wall pieces (10,10 -> 20,10 and 20,10 -> 20,20); the inner wall and the other object are not here + assert first.a.tolist() == [[10.0, 10.0], [20.0, 10.0]] + assert first.b.tolist() == [[20.0, 10.0], [20.0, 20.0]] + assert first.width.tolist() == [0.42, 0.42] + assert first.mid_z == pytest.approx(0.1) + assert walls.objects == ["Pot_Pepper.stl_id_0_copy_0", "Plug_x2.stl_id_1_copy_0"] + + +def test_the_other_object_has_its_own_width_and_pieces(gcode: Path) -> None: + walls = read_walls(gcode, "plug x2.stl") + assert len(walls.layers) == 1 # it prints in the first layer only + assert walls.layers[0].a.tolist() == [[50.0, 10.0]] + assert walls.layers[0].width.tolist() == [0.40] + + +def test_a_non_extruding_or_zero_length_move_is_not_a_wall(gcode: Path) -> None: + second = read_walls(gcode, "Pot_Pepper.stl").layers[1] + # `E0` travel to X30 and the zero-length E move add nothing + assert second.a.tolist() == [[10.0, 10.0]] + assert second.b.tolist() == [[20.0, 10.0]] + + +def test_an_unknown_object_names_the_ones_the_file_has(gcode: Path) -> None: + with pytest.raises( + KeyError, match=r"Pot_Pepper\.stl_id_0_copy_0, Plug_x2\.stl_id_1_copy_0" + ): + read_walls(gcode, "Nothing.stl") + + +def test_a_plate_with_several_objects_needs_a_name(gcode: Path) -> None: + with pytest.raises(ValueError, match="several objects"): + read_walls(gcode) + + +def test_other_feature_kinds_can_be_asked_for(gcode: Path) -> None: + walls = read_walls(gcode, "Pot Pepper.stl", kinds=("Inner wall",)) + assert walls.layers[0].a.tolist() == [[20.0, 20.0]] + assert walls.layers[0].width.tolist() == [0.45] + + +def test_two_copies_of_one_object_are_not_merged(tmp_path: Path) -> None: + text = "M83\n" + "".join( + f";LAYER_CHANGE\n;Z:0.2\n;HEIGHT:0.2\n; printing object Cup.stl id:0 copy {c}\n;TYPE:Outer wall\n;WIDTH:0.4\n" + f"G1 X{10 * c} Y0 F3000\nG1 X{10 * c + 5} Y0 E0.5\n; stop printing object Cup.stl id:0 copy {c}\n" + for c in (0, 1) + ) + p = write(tmp_path, text) + with pytest.raises(ValueError, match="several copies"): + read_walls(p, "Cup.stl") + one = read_walls(p, "Cup.stl id:0 copy 1") + assert one.layers[0].a.tolist() == [[10.0, 0.0]] + + +def test_prusaslicer_comment_labels_and_m486_name_the_object(tmp_path: Path) -> None: + prusa = write( + tmp_path, + "M83\n;LAYER_CHANGE\n;Z:0.2\n;HEIGHT:0.2\n; printing object A id:0 copy 0\n;TYPE:Outer wall\n;WIDTH:0.4\n" + "G1 X1 Y1 F3000\nG1 X5 Y1 E0.5\n; stop printing object A id:0 copy 0\n; printing object B id:1 copy 0\n" + "G1 X30 Y1 F3000\nG1 X35 Y1 E0.5\n; stop printing object B id:1 copy 0\n", + "prusa.gcode", + ) + assert read_walls(prusa, "B").layers[0].a.tolist() == [[30.0, 1.0]] + marlin = write( + tmp_path, + "M83\nM486 S0\n;LAYER_CHANGE\n;Z:0.2\n;HEIGHT:0.2\n;TYPE:Outer wall\n;WIDTH:0.4\nG1 X1 Y1 F3000\nG1 X5 Y1 E0.5\n" + "M486 S1\nG1 X30 Y1 F3000\nG1 X35 Y1 E0.5\nM486 S-1\n", + "marlin.gcode", + ) + assert read_walls(marlin, "object 1").layers[0].a.tolist() == [[30.0, 1.0]] + + +def test_a_bambu_plate_is_read_by_object_name_with_its_z_height_and_layer_height() -> ( + None +): + """`; Z_HEIGHT:` / `; LAYER_HEIGHT:` carry the layer in Bambu files; the name is written `NAME — NAME id:N copy 0`.""" + walls = read_walls(DATA / "bambu_a1mini_path.gcode", "SW") + assert len(walls.layers) == 1 + layer = walls.layers[0] + # the layer marker's step is 0.3, but the object prints at its own path height 0.12 (`; LAYER_HEIGHT: 0.12`) + assert (layer.z, layer.h) == (0.8, 0.12) + assert layer.mid_z == pytest.approx(0.74) + assert layer.width[0] == pytest.approx( + 0.42, abs=1e-3 + ) # `; LINE_WIDTH:` of the outer wall + assert "SW — SW id:23 copy 0" in walls.objects + + +def test_the_first_height_of_a_layer_is_its_thickness_not_the_last( + tmp_path: Path, +) -> None: + """Orca repeats `;HEIGHT:` when a path's height changes (a bridge at the end of the layer): mid z must not move.""" + p = write( + tmp_path, + layer_gcode( + "G1 X1 Y1 F3000\nG1 X5 Y1 E0.5\n;TYPE:Internal bridge\n;HEIGHT:0.4\nG1 X9 Y1 E0.5\n", + ";Z:1.88\n;HEIGHT:0.12\n", + ), + ) + layer = read_walls(p).layers[0] + assert (layer.z, layer.h) == (1.88, 0.12) + assert layer.mid_z == pytest.approx(1.82) + + +def test_absolute_extrusion_draws_only_where_e_grows(tmp_path: Path) -> None: + p = write( + tmp_path, + "M82\nG92 E0\n" + + layer_gcode( + "G1 X1 Y0 F3000\nG1 X2 Y0 E1.0\nG1 X3 Y0 E1.0\nG1 X4 Y0 E0.5\nG92 E0\nG1 X5 Y0 E0.5\n" + ), + ) + # E1.0 draws; the same E again draws nothing; a LOWER E (retraction) draws nothing; after G92 E0 the 0.5 draws + assert read_walls(p).layers[0].b.tolist() == [[2.0, 0.0], [5.0, 0.0]] + + +def test_relative_positioning_g91_moves_from_the_current_point(tmp_path: Path) -> None: + p = write( + tmp_path, + "M83\nG90\n" + + layer_gcode( + "G1 X10 Y10 F3000\nG91\nG1 X5 Y0 E0.5\nG1 X0 Y5 E0.5\nG90\nG1 X0 Y0 E0.5\n" + ), + ) + walls = read_walls(p).layers[0] + assert walls.b.tolist() == [[15.0, 10.0], [15.0, 15.0], [0.0, 0.0]] + + +def test_the_position_is_unknown_after_homing(tmp_path: Path) -> None: + p = write( + tmp_path, + "M83\n" + + layer_gcode("G1 X10 Y10 F3000\nG28\nG1 X20 Y10 E0.5\nG1 X30 Y10 E0.5\n"), + ) + # the first move after G28 starts from an unknown point: it draws nothing, the next one does + assert read_walls(p).layers[0].a.tolist() == [[20.0, 10.0]] + + +def test_arcs_are_followed_as_chords(tmp_path: Path) -> None: + p = write( + tmp_path, + "M83\n" + + layer_gcode( + "G1 X10 Y0 F3000\nG3 X-10 Y0 I-10 J0 E1.0\nG3 X10 Y0 I10 J0 E1.0\n" + ), + ) + pts, _w = read_walls(p).layers[0].points(0.2) + r = np.hypot(pts[:, 0], pts[:, 1]) + assert r.min() == pytest.approx(10.0, abs=0.02) + assert r.max() == pytest.approx(10.0, abs=0.02) + assert len(pts) > 100 # a full circle of radius 10 densified every ~0.2 mm + assert float(np.ptp(pts[:, 0])) == pytest.approx(20.0, abs=0.05) + + +def test_crafted_files_are_refused_not_expanded(tmp_path: Path) -> None: + huge = "G1 X0 Y0 F3000\nG3 X10 Y0 I-3000 J0 E1\n" + with pytest.raises(ValueError, match="not a print move"): + read_walls(write(tmp_path, "M83\n" + layer_gcode(huge))) + long_piece = "G1 X0 Y0 F3000\nG1 X100000000 Y0 E1\n" + with pytest.raises(ValueError, match="not a print move"): + read_walls(write(tmp_path, "M83\n" + layer_gcode(long_piece), "long.gcode")) + radius_form = "G1 X0 Y0 F3000\nG2 X10 Y0 R5 E1\n" + with pytest.raises(ValueError, match="radius form"): + read_walls(write(tmp_path, "M83\n" + layer_gcode(radius_form), "r.gcode")) + + +def test_g2_is_clockwise_and_g3_counter_clockwise(tmp_path: Path) -> None: + cw = read_walls(write(tmp_path, "M83\n" + layer_gcode("G1 X10 Y0 F3000\nG2 X0 Y-10 I-10 J0 E1\n"))).layers[0] + ccw = read_walls(write(tmp_path, "M83\n" + layer_gcode("G1 X10 Y0 F3000\nG3 X0 Y10 I-10 J0 E1\n"), "ccw.gcode")).layers[0] + assert float(cw.b[:, 1].min()) < -9.9 and float(cw.b[:, 1].max()) <= 0.0 # clockwise from (10,0) goes down + assert float(ccw.b[:, 1].max()) > 9.9 and float(ccw.b[:, 1].min()) >= 0.0 + + +def test_prusaslicer_external_perimeter_is_a_wall_and_a_file_without_walls_is_refused(tmp_path: Path) -> None: + prusa = write( + tmp_path, + "M83\n;LAYER_CHANGE\n;Z:0.2\n;HEIGHT:0.2\n;TYPE:External perimeter\n;WIDTH:0.45\nG1 X1 Y1 F3000\nG1 X5 Y1 E0.5\n" + ";TYPE:Overhang wall\nG1 X5 Y5 E0.5\n", + ) + walls = read_walls(prusa) + assert walls.layers[0].b.tolist() == [[5.0, 1.0]] + assert walls.overhang_mm == pytest.approx(4.0) # the left-out overhang wall is counted, not lost + assert read_walls(prusa, kinds=("Overhang wall",)).layers[0].b.tolist() == [[5.0, 5.0]] + with pytest.raises(KeyError, match="prints no"): + read_walls(write(tmp_path, "G1 X1 Y1\n", "empty.gcode")) + + +def test_each_pieces_own_height_gives_the_layer_height_and_a_raft_gives_a_base(tmp_path: Path) -> None: + """The marker says 0.24 (the step above the raft), the object's path says 0.12: mid z and the base follow the path.""" + p = write(tmp_path, "M83\n" + layer_gcode("G1 X1 Y1 F3000\n;HEIGHT:0.12\nG1 X5 Y1 E0.5\n", ";Z:1.34\n;HEIGHT:0.24\n")) + walls = read_walls(p) + assert walls.layers[0].mid_z == pytest.approx(1.28) + assert walls.base_z == pytest.approx(1.22) + + +def test_points_are_spaced_by_the_requested_step(gcode: Path) -> None: + layer = read_walls(gcode, "Pot_Pepper.stl").layers[0] + pts, width = layer.points(0.5) + assert len(pts) == len(width) == 40 # two 10 mm pieces at 0.5 mm + assert np.hypot(*(pts[1] - pts[0])) == pytest.approx(0.5) + + +def test_a_wall_without_a_width_tag_takes_the_footers_line_width( + tmp_path: Path, +) -> None: + p = write( + tmp_path, + "M83\n;LAYER_CHANGE\n;Z:0.2\n;HEIGHT:0.2\n;TYPE:Outer wall\nG1 X1 Y1 F3000\nG1 X5 Y1 E0.2\n" + "; outer_wall_line_width = 0.45\n", + ) + assert read_walls(p).layers[0].width.tolist() == [0.45] + + +def test_a_nonsense_line_width_is_ignored_not_used(tmp_path: Path) -> None: + p = write( + tmp_path, + "M83\n;LAYER_CHANGE\n;Z:0.2\n;HEIGHT:0.2\n;TYPE:Outer wall\nG1 X1 Y1 F3000\nG1 X5 Y1 E0.2\n" + "; outer_wall_line_width = 105%\n", + ) + assert read_walls(p).layers[0].width.tolist() == [0.42] # the default, not 105 mm