diff --git a/src/mechcomp/stl.py b/src/mechcomp/stl.py new file mode 100644 index 0000000..9f6e0aa --- /dev/null +++ b/src/mechcomp/stl.py @@ -0,0 +1,291 @@ +""" +Export a member as a sealed triangle mesh. + +WHAT THIS IS, AND WHAT IT IS NOT + **A member exporter.** ``ROADMAP.md`` §2 names three artifact classes: + members are prismatic and exist, nodes are non-prismatic and do not, panels + are sheet and do not. What is here sweeps a cross-section along +Z and + caps both ends. It will never produce a node, and no amount of extending it + should be attempted -- a node is a different representation, not a harder + prism. + + The module is deliberately named for the artifact rather than the file + format for that reason. ``stl.py`` would invite someone to add nodes to it. + +IDENTIFIERS IN THE OUTPUT ARE PROVENANCE, NOT A CHECKSUM + ``ROADMAP.md`` §7 principle 4: mesh bytes are not reproducible across + toolchain versions. Two exports carrying the same ``build_id`` may differ + byte for byte and both be correct -- a Shapely or GEOS release can move a + triangulation without moving the geometry, which is the F-034 mechanism + applied to tessellation instead of to booleans. + + So the identifiers written into the file say *what produced this*. They do + not say *these bytes*. Anyone diffing two exports of the same design and + concluding the compiler is broken has misread them, and the header says so + in the file itself. + +WHY NO CAD KERNEL + A member is a prism. Caps are a constrained Delaunay triangulation of the + section, walls are a quad strip per ring, and both are decided by the + section's own rings. Nothing here needs a solid modeller, which is why STL + export costs no new dependency while STEP would. + +THE PRECONDITION IS THE SECTION'S, NOT OURS + ``region.is_region_simple`` exists because a section that touches itself at + a point is a valid 2D outline that cannot be tessellated: it measures + perfectly and then fails to extrude into a sealed solid. That is checked + here and refused, rather than emitted as a mesh that opens in a slicer and + prints wrong. +""" + +from __future__ import annotations + +import struct +from typing import Dict, List, NamedTuple, Sequence, Tuple + +Point = Tuple[float, float] +Path = Sequence[Point] +Region = Sequence[Path] +Vertex = Tuple[float, float, float] +Triangle = Tuple[int, int, int] + + +class ExportRefused(Exception): + """ + The section cannot become a sealed solid. + + A refusal rather than an error, in the same sense as ``ProfileRejected``: + the compiler declining to produce something it cannot vouch for is the + compiler working. + """ + + +class Mesh(NamedTuple): + vertices: List[Vertex] + triangles: List[Triangle] + + @property + def volume(self) -> float: + """ + Signed volume by the divergence theorem. + + Positive for outward-facing normals. Comparing this against + ``region_area(section) * length`` checks the winding, the cap + orientation and the sweep length in a single number -- three defects + that are individually hard to see and jointly impossible to miss. + """ + total = 0.0 + v = self.vertices + for i, j, k in self.triangles: + ax, ay, az = v[i] + bx, by, bz = v[j] + cx, cy, cz = v[k] + total += (ax * (by * cz - bz * cy) + - ay * (bx * cz - bz * cx) + + az * (bx * cy - by * cx)) + return total / 6.0 + + +# --------------------------------------------------------------------------- +# Building the mesh +# --------------------------------------------------------------------------- + +def _signed_area(ring: Path) -> float: + total = 0.0 + n = len(ring) + for i in range(n): + x0, y0 = ring[i] + x1, y1 = ring[(i + 1) % n] + total += x0 * y1 - x1 * y0 + return total / 2.0 + + +def _ccw(ring: Path) -> List[Point]: + return list(ring) if _signed_area(ring) >= 0 else list(reversed(ring)) + + +def _cw(ring: Path) -> List[Point]: + return list(ring) if _signed_area(ring) < 0 else list(reversed(ring)) + + +class _Vertices: + """Deduplicating vertex table, keyed on exact coordinates.""" + + def __init__(self) -> None: + self.items: List[Vertex] = [] + self._index: Dict[Vertex, int] = {} + + def add(self, x: float, y: float, z: float) -> int: + key = (x, y, z) + found = self._index.get(key) + if found is None: + found = len(self.items) + self._index[key] = found + self.items.append(key) + return found + + +def prism_mesh(section: Region, length: float) -> Mesh: + """ + Sweep ``section`` from z=0 to z=``length`` and seal both ends. + + Takes a region and a length rather than a ``Result``, which is the whole + reason the signature looks like this. A member, its support, or the support + alone are then three calls with different regions instead of three special + cases inside one function -- and support is a designed part of an artifact + when no printable orientation exists, not something a slicer adds. + + WINDING + ``region_parts`` returns each part as an outer ring clockwise with its + holes counter-clockwise, because ``region_area`` depends on that + convention to make holes subtract. A prism swept along +Z needs the + opposite: traversing the outer boundary counter-clockwise puts the + material on the left and the wall normal on the right, which is + outward. Both are therefore reversed here. The convention is not + assumed from the input -- every ring is forced. + """ + from .geom.region import is_region_simple, region_parts, to_shapely + + if length <= 0: + raise ExportRefused("length must be greater than zero, got %r" % (length,)) + if not section: + raise ExportRefused("the section is empty; there is nothing to sweep") + if not is_region_simple(section): + raise ExportRefused( + "the section touches or crosses itself, so it cannot be sealed " + "into a solid. It is a valid outline and measures correctly; it " + "is not extrudable. See geom.region.is_region_simple.") + + try: + from shapely import constrained_delaunay_triangles + except ImportError as exc: # pragma: no cover + raise ExportRefused( + "constrained_delaunay_triangles requires Shapely 2.1 or newer" + ) from exc + + verts = _Vertices() + tris: List[Triangle] = [] + + # -- caps ------------------------------------------------------------ + # + # The triangulator is given the polygon with its holes, so no triangle + # lands inside a bore. That is asserted in the tests rather than trusted: + # a triangulator that filled the bores would produce a mesh that looks + # right in a slicer and prints solid through the conduit. + tessellation = constrained_delaunay_triangles(to_shapely(section)) + faces = getattr(tessellation, "geoms", [tessellation]) + + cap_count = 0 + for face in faces: + if face.is_empty: + continue + ring = [(x, y) for x, y in list(face.exterior.coords)[:-1]] + if len(ring) != 3: + continue + a, b, c = _ccw(ring) + cap_count += 1 + # Top at z=length: counter-clockwise seen from +Z faces +Z, outward. + tris.append((verts.add(a[0], a[1], length), + verts.add(b[0], b[1], length), + verts.add(c[0], c[1], length))) + # Bottom at z=0: the same triangle reversed faces -Z, also outward. + tris.append((verts.add(c[0], c[1], 0.0), + verts.add(b[0], b[1], 0.0), + verts.add(a[0], a[1], 0.0))) + + if cap_count == 0: + raise ExportRefused( + "the section produced no triangles; it encloses no area") + + # -- walls ----------------------------------------------------------- + for part in region_parts(section): + rings = [_ccw(part[0])] + [_cw(hole) for hole in part[1:]] + for ring in rings: + n = len(ring) + for i in range(n): + x0, y0 = ring[i] + x1, y1 = ring[(i + 1) % n] + b0 = verts.add(x0, y0, 0.0) + b1 = verts.add(x1, y1, 0.0) + t0 = verts.add(x0, y0, length) + t1 = verts.add(x1, y1, length) + tris.append((b0, b1, t1)) + tris.append((b0, t1, t0)) + + return Mesh(verts.items, tris) + + +# --------------------------------------------------------------------------- +# Writing it out +# --------------------------------------------------------------------------- + +HEADER_BYTES = 80 + + +def stl_header(input_id: str = "", build_id: str = "") -> bytes: + """ + The binary format's 80-byte header, carrying provenance. + + A file separated from its design record should still say what it is. Eighty + bytes is not much, so it holds the two identifiers and nothing else -- + about fifty-four characters of the eighty. + + These are NOT a checksum of the bytes that follow. See the module + docstring: mesh bytes are not reproducible across toolchain versions, and + two files with the same ``build_id`` may legitimately differ. + """ + text = "mechcomp" + if input_id: + text += " input=%s" % input_id + if build_id: + text += " build=%s" % build_id + raw = text.encode("ascii", "replace")[:HEADER_BYTES] + return raw + b"\0" * (HEADER_BYTES - len(raw)) + + +def binary_stl(mesh: Mesh, input_id: str = "", build_id: str = "") -> bytes: + """ + The mesh as a binary STL. + + Facet normals are written as zeros. The format has a field for them and + every consumer recomputes from the vertex order anyway; writing a stored + normal that could disagree with the winding would create a second source + of truth for which way a face points. + """ + out = [stl_header(input_id, build_id), + struct.pack(" bytes: + """ + Export whatever ``build()`` returned. + + ``length`` defaults to the report's ``LENGTH_MM``, which is the value the + record's ``VOLUME_MM3`` and ``MASS_G`` were computed against. Sweeping any + other length would make the record describe a different object than the + file shipped beside it, so the default is the only sensible one and an + override exists for callers that know what they are doing. + """ + if length is None: + length = result.report["LENGTH_MM"] + mesh = prism_mesh(result.section, length) + record = getattr(result, "record", None) + return binary_stl(mesh, + record.input_id if record else "", + record.build_id if record else "") + + +def stl_filename(result) -> str: + """A name that identifies the design without needing the record beside it.""" + record = getattr(result, "record", None) + ident = record.input_id if record else "unidentified" + family = result.report.get("FAMILY", "member") + profile = str(result.report.get("PROFILE", "")).replace(" ", "-").lower() + return "mechcomp-%s-%s-%s.stl" % (family, profile, ident) diff --git a/tests/test_stl.py b/tests/test_stl.py new file mode 100644 index 0000000..ba15a82 --- /dev/null +++ b/tests/test_stl.py @@ -0,0 +1,334 @@ +""" +The member mesh: sealed, correctly oriented, and the right size. + +WHY THESE FOUR CHECKS AND NOT A VISUAL ONE + A bad mesh opens in a slicer. That is the whole problem: inverted normals, + a hole filled in, a wall missing or the wrong sweep length all produce a + file that looks plausible and prints wrong, and the failure is discovered + in plastic hours later. + + So each property is asserted arithmetically, and each is chosen to fail + loudly for a different defect: + + edge pairing a missing or duplicated wall, and inconsistent winding + Euler a hole that was filled, or one invented + centroid-in-hole a bore tessellated over -- prints solid through the conduit + volume inverted normals and a wrong length, in one number + + Every negative assertion has a positive control, per ROADMAP principle 11. +""" + +from __future__ import annotations + +import struct + +import pytest + +stl = pytest.importorskip("mechcomp.stl") +profiles = pytest.importorskip("mechcomp.profiles") +region = pytest.importorskip("mechcomp.geom.region") + +LENGTH = 40.0 + + +def every_profile(): + """Each catalogue entry of each family, at its family defaults.""" + for name, family in profiles.FAMILIES.items(): + for profile in sorted(family.catalogue): + yield name, profile + + +def built(family="3x", profile="Y", **params): + return profiles.build(family, profile, params) + + +ALL = list(every_profile()) + + +def ids(pair): + return "%s/%s" % pair + + +# --------------------------------------------------------------------------- +# Sealed and consistently wound +# --------------------------------------------------------------------------- + +@pytest.mark.parametrize("pair", ALL, ids=ids) +def test_every_directed_edge_appears_exactly_once(pair): + """ + Manifoldness and consistent orientation in one assertion. + + In a closed, consistently oriented surface each directed edge (a, b) occurs + once and its reverse (b, a) occurs once. A duplicated directed edge means + two faces wind the same way across a shared edge; a missing reverse means a + hole in the surface. + """ + mesh = stl.prism_mesh(built(*pair).section, LENGTH) + seen = set() + for i, j, k in mesh.triangles: + for edge in ((i, j), (j, k), (k, i)): + assert edge not in seen, "directed edge %r twice" % (edge,) + seen.add(edge) + for a, b in seen: + assert (b, a) in seen, "edge %r has no opposite" % ((a, b),) + + +def test_the_edge_check_can_fail(): + """ + The positive control. Without it, a mesh with no triangles at all would + satisfy the test above vacuously. + """ + mesh = stl.prism_mesh(built().section, LENGTH) + assert len(mesh.triangles) > 0 + broken = stl.Mesh(mesh.vertices, mesh.triangles[:-1]) + seen = set() + for i, j, k in broken.triangles: + for edge in ((i, j), (j, k), (k, i)): + seen.add(edge) + assert any((b, a) not in seen for a, b in seen) + + +# --------------------------------------------------------------------------- +# The right number of holes, still holes +# --------------------------------------------------------------------------- + +@pytest.mark.parametrize("pair", ALL, ids=ids) +def test_euler_characteristic_matches_the_section(pair): + """ + V - E + F = 2P - 2H. + + A prism over a part with H holes is a genus-H handlebody, so its boundary + surface has characteristic 2 - 2H, and P disjoint parts sum. Filling a bore + would raise the left side; inventing a tunnel would lower it. + """ + result = built(*pair) + section = result.section + parts = region.region_parts(section) + p = len(parts) + h = sum(len(part) - 1 for part in parts) + + mesh = stl.prism_mesh(section, LENGTH) + edges = {frozenset((a, b)) + for i, j, k in mesh.triangles + for a, b in ((i, j), (j, k), (k, i))} + chi = len(mesh.vertices) - len(edges) + len(mesh.triangles) + assert chi == 2 * p - 2 * h, ( + "chi=%d for %d part(s) and %d hole(s)" % (chi, p, h)) + + +def test_at_least_one_profile_actually_has_a_hole(): + """ + The positive control for the Euler test and the one below it: if no + section in the catalogue had holes, both would be asserting the easy case + and a triangulator that ignored holes would pass everything. + """ + with_holes = 0 + for pair in ALL: + parts = region.region_parts(built(*pair).section) + with_holes += sum(len(part) - 1 for part in parts) + assert with_holes > 0 + + +# --------------------------------------------------------------------------- +# Nothing tessellated across a bore +# --------------------------------------------------------------------------- + +@pytest.mark.parametrize("pair", ALL, ids=ids) +def test_no_cap_triangle_sits_inside_a_hole(pair): + """ + The check that matters most in practice. + + A triangulator that filled the bores produces a mesh that is sealed, + correctly wound, passes Euler if the holes vanish consistently -- and + prints solid where the strap has to pass. Only the volume test and this one + would notice. + """ + result = built(*pair) + section = result.section + mesh = stl.prism_mesh(section, LENGTH) + + holes = [hole for part in region.region_parts(section) for hole in part[1:]] + if not holes: + pytest.skip("this section has no holes") + + caps = [t for t in mesh.triangles + if len({mesh.vertices[i][2] for i in t}) == 1] + assert caps, "no cap triangles found" + + from shapely.geometry import Polygon + + # Overlapping AREA, not the centroid. A centroid test passes by luck on a + # coarse triangulation: two triangles covering a square with a central bore + # can both have their centroids outside it while jointly paving it over. + # Found by mutation -- ignoring the holes at triangulation time slipped + # through a centroid check and is caught by this one. + hole_polys = [Polygon(h) for h in holes] + for i, j, k in caps: + tri = Polygon([(mesh.vertices[v][0], mesh.vertices[v][1]) + for v in (i, j, k)]) + if tri.area <= 0: + continue + for hp in hole_polys: + assert tri.intersection(hp).area <= 1e-9 * tri.area, ( + "a cap triangle paves over a hole -- the part would print " + "solid where the stock has to pass") + + +# --------------------------------------------------------------------------- +# The right size, the right way out +# --------------------------------------------------------------------------- + +@pytest.mark.parametrize("pair", ALL, ids=ids) +def test_volume_equals_area_times_length(pair): + """ + Catches inverted normals and a wrong sweep length together. + + Compared against ``region_area`` at full precision rather than against the + report's ``SECTION_AREA_MM2``, which is rounded to six significant figures + and would cap this test's resolution at the oracle's rather than the + geometry's. + """ + result = built(*pair) + expected = region.area(result.section) * LENGTH + actual = stl.prism_mesh(result.section, LENGTH).volume + assert actual > 0, "negative volume means the normals face inward" + assert abs(actual - expected) <= 1e-9 * abs(expected) + + +def test_volume_agrees_with_the_report_too(): + """ + The same quantity through the published numbers, at the precision they are + published to. Six significant figures, so the bound is loose on purpose. + """ + result = built() + length = result.report["LENGTH_MM"] + expected = result.report["SECTION_AREA_MM2"] * length + actual = stl.prism_mesh(result.section, length).volume + assert abs(actual - expected) <= 1e-5 * abs(expected) + + +def test_volume_scales_with_length(): + """Positive control: a volume that ignored its length would be constant.""" + section = built().section + assert stl.prism_mesh(section, 80.0).volume == pytest.approx( + 2.0 * stl.prism_mesh(section, 40.0).volume) + + +def test_reversing_a_wall_would_be_caught(): + """ + Positive control for the sign: a mesh with every triangle flipped has + negative volume, so ``actual > 0`` above is a real assertion. + """ + mesh = stl.prism_mesh(built().section, LENGTH) + flipped = stl.Mesh(mesh.vertices, [(k, j, i) for i, j, k in mesh.triangles]) + assert flipped.volume < 0 + + +# --------------------------------------------------------------------------- +# The length that ships +# --------------------------------------------------------------------------- + +@pytest.mark.parametrize("params,expected", [ + ({"preview_length_mm": 37.5}, 37.5), + ({"length_view": "Full Length", "member_length_ft": 10}, 3048.0), +]) +def test_member_stl_sweeps_the_reported_length(params, expected): + """ + The record's VOLUME_MM3 and MASS_G are computed against LENGTH_MM. Sweeping + anything else would make the record describe a different object than the + file beside it. + + NEITHER CASE IS 100 mm, AND THAT IS THE POINT. + An earlier version of this test built at the family defaults, where + ``preview_length_mm`` is 100. A mutation replacing the lookup with a + hardcoded 100.0 passed the whole suite -- the assertion compared a + value against the constant that had replaced it. Both cases here differ + from 100, and the second differs by two orders of magnitude, so no + plausible fixed length satisfies them. + """ + result = built(**params) + assert result.report["LENGTH_MM"] == expected + data = stl.member_stl(result) + count = struct.unpack(" 30 * stl.prism_mesh(preview.section, preview.report["LENGTH_MM"]).volume) + + +# --------------------------------------------------------------------------- +# The file +# --------------------------------------------------------------------------- + +def test_binary_stl_is_exactly_the_right_size(): + mesh = stl.prism_mesh(built().section, LENGTH) + data = stl.binary_stl(mesh) + assert len(data) == 84 + 50 * len(mesh.triangles) + assert struct.unpack("