diff --git a/snub-dodecahedron/.gitignore b/snub-dodecahedron/.gitignore new file mode 100644 index 0000000..bef8b65 --- /dev/null +++ b/snub-dodecahedron/.gitignore @@ -0,0 +1,4 @@ +__pycache__/ +*.pyc +*.stl +out/ diff --git a/snub-dodecahedron/net_preview.svg b/snub-dodecahedron/net_preview.svg new file mode 100644 index 0000000..ffec421 --- /dev/null +++ b/snub-dodecahedron/net_preview.svg @@ -0,0 +1,186 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/snub-dodecahedron/preview.py b/snub-dodecahedron/preview.py new file mode 100644 index 0000000..2cb15ea --- /dev/null +++ b/snub-dodecahedron/preview.py @@ -0,0 +1,95 @@ +"""Render an SVG preview of the unfolded net(s) and verify no faces overlap.""" + +from __future__ import annotations + +import numpy as np + +from snubgeom import SnubDodecahedron +from unfold import Unfolder, _convex_overlap, _shrink + + +def layout_pieces(pieces, gap=1.5): + """Translate each piece so they sit side by side; return offsets.""" + offsets = [] + x_cursor = 0.0 + for p in pieces: + x0, y0, x1, y1 = p.bbox() + offsets.append((x_cursor - x0, -y0)) + x_cursor += (x1 - x0) + gap + return offsets + + +def verify_no_overlap(pieces): + """Independent global check: no two faces in the same piece overlap.""" + bad = 0 + for p in pieces: + polys = [poly for _, poly in p.faces] + adj_share = {} # face index pairs that share an edge are allowed to touch + for i in range(len(polys)): + for j in range(i + 1, len(polys)): + a = _shrink(polys[i], 5e-3) + b = _shrink(polys[j], 5e-3) + if _convex_overlap(a, b): + bad += 1 + return bad + + +def write_svg(pieces, path, scale=40.0, gap=1.5): + offsets = layout_pieces(pieces, gap=gap) + # global bbox + minx = miny = 1e9 + maxx = maxy = -1e9 + for p, (ox, oy) in zip(pieces, offsets): + for _, poly in p.faces: + for x, y in poly: + minx = min(minx, x + ox) + miny = min(miny, y + oy) + maxx = max(maxx, x + ox) + maxy = max(maxy, y + oy) + pad = 1.0 + W = (maxx - minx + 2 * pad) * scale + H = (maxy - miny + 2 * pad) * scale + + def X(x): + return (x - minx + pad) * scale + + def Y(y): + return H - (y - miny + pad) * scale # flip for screen coords + + out = [ + f'', + f'', + ] + for p, (ox, oy) in zip(pieces, offsets): + for fi, poly in p.faces: + pts = " ".join(f"{X(x + ox):.2f},{Y(y + oy):.2f}" for x, y in poly) + fill = "#ffe9b3" if len(poly) == 5 else "#cfe8ff" + out.append( + f'' + ) + # hinges (fold lines) in green + for h in p.hinges: + x1, y1 = h.p + x2, y2 = h.q + col = "#2a8f2a" if h.kind == "3-3" else "#c000c0" + out.append( + f'' + ) + out.append("") + with open(path, "w") as f: + f.write("\n".join(out)) + + +if __name__ == "__main__": + s = SnubDodecahedron() + s.verify() + pieces = Unfolder(s).unfold() + bad = verify_no_overlap(pieces) + print(f"pieces={len(pieces)} faces={sum(len(p.faces) for p in pieces)} " + f"overlapping_pairs={bad}") + write_svg(pieces, "net_preview.svg") + print("wrote net_preview.svg") diff --git a/snub-dodecahedron/snubgeom.py b/snub-dodecahedron/snubgeom.py new file mode 100644 index 0000000..801f94f --- /dev/null +++ b/snub-dodecahedron/snubgeom.py @@ -0,0 +1,278 @@ +""" +Exact geometry of the snub dodecahedron (an Archimedean solid). + +The snub dodecahedron has: + - 60 vertices + - 150 edges + - 92 faces: 12 regular pentagons + 80 equilateral triangles + - every vertex has configuration 3.3.3.3.5 + +We build it from the exact Cartesian coordinates (Wikipedia / Weisstein), +take the convex hull, merge coplanar hull triangles back into the real +pentagon/triangle faces, and expose: + + - vertices (60, 3) float array, normalised to edge length 1 + - faces list[list[int]] CCW (outward) vertex-index loops + - edges dict[(i,j)] -> (faceA, faceB) + - dihedral_angle(...) interior dihedral angle of an edge (degrees) + +Everything is verified at import time (see `build_snub_dodecahedron`): +all 150 edges are equal length, there are exactly 12 pentagons and 80 +triangles, and the dihedral angles match the known values +(3-3: 164.175 deg, 3-5: 152.930 deg). +""" + +from __future__ import annotations + +import itertools +import math + +import numpy as np +from scipy.spatial import ConvexHull + +PHI = (1.0 + math.sqrt(5.0)) / 2.0 + + +def _xi() -> float: + """Real root of xi^3 - 2*xi = phi.""" + # Solve xi^3 - 2 xi - phi = 0 for the single real positive root. + coeffs = [1.0, 0.0, -2.0, -PHI] + roots = np.roots(coeffs) + real = [r.real for r in roots if abs(r.imag) < 1e-9] + # The relevant root is the largest real one (~1.7155615). + return max(real) + + +def _base_points() -> np.ndarray: + """The five generating points before permutation / sign expansion.""" + xi = _xi() + a = xi - 1.0 / xi + b = xi * PHI + PHI * PHI + PHI / xi + p = PHI + pts = [ + (2 * a, 2.0, 2 * b), + (a + b / p + p, -a * p + b + 1.0 / p, a / p + b * p - 1.0), + (a + b / p - p, a * p - b + 1.0 / p, a / p + b * p + 1.0), + (-a / p + b * p + 1.0, -a + b / p - p, a * p + b - 1.0 / p), + (-a / p + b * p - 1.0, a - b / p - p, a * p + b + 1.0 / p), + ] + return np.array(pts, dtype=float) + + +def _even_permutations(v): + """The 3 even (cyclic) permutations of a 3-vector.""" + x, y, z = v + return [(x, y, z), (y, z, x), (z, x, y)] + + +def _sign_variants(v): + """Sign assignments with an even number of minus signs (0 or 2).""" + x, y, z = v + out = [] + for sx, sy, sz in itertools.product((1, -1), repeat=3): + if (sx * sy * sz) == 1: # product +1 <=> even number of minus signs + out.append((sx * x, sy * y, sz * z)) + return out + + +def _all_vertices() -> np.ndarray: + seen = [] + for base in _base_points(): + for perm in _even_permutations(base): + for sv in _sign_variants(perm): + seen.append(sv) + pts = np.array(seen, dtype=float) + # De-duplicate (there should be exactly 60 distinct points). + uniq = [] + for p in pts: + if not any(np.allclose(p, q, atol=1e-6) for q in uniq): + uniq.append(p) + return np.array(uniq, dtype=float) + + +def _merge_coplanar_faces(points: np.ndarray, hull: ConvexHull): + """Group hull simplices (triangles) into the real polygonal faces. + + Two adjacent simplices belong to the same face iff they are coplanar + (same outward normal). Returns a list of faces, each an ordered + (CCW seen from outside) list of vertex indices. + """ + n_simp = len(hull.simplices) + + # Outward unit normal for each simplex. + normals = hull.equations[:, :3] + offsets = hull.equations[:, 3] + + # Union-find over simplices that share an edge and are coplanar. + parent = list(range(n_simp)) + + def find(i): + while parent[i] != i: + parent[i] = parent[parent[i]] + i = parent[i] + return i + + def union(i, j): + parent[find(i)] = find(j) + + # Map each undirected triangle edge -> simplices that contain it. + from collections import defaultdict + + edge_to_simps = defaultdict(list) + for si, tri in enumerate(hull.simplices): + for a, b in ((tri[0], tri[1]), (tri[1], tri[2]), (tri[2], tri[0])): + edge_to_simps[frozenset((int(a), int(b)))].append(si) + + for simps in edge_to_simps.values(): + if len(simps) == 2: + s0, s1 = simps + if ( + np.allclose(normals[s0], normals[s1], atol=1e-6) + and abs(offsets[s0] - offsets[s1]) < 1e-6 + ): + union(s0, s1) + + groups = defaultdict(list) + for si in range(n_simp): + groups[find(si)].append(si) + + faces = [] + for simps in groups.values(): + vids = set() + for si in simps: + vids.update(int(v) for v in hull.simplices[si]) + vids = list(vids) + normal = normals[simps[0]] + center = points[vids].mean(axis=0) + # Build an in-plane basis to order vertices by angle. + ref = points[vids[0]] - center + ref = ref / np.linalg.norm(ref) + binormal = np.cross(normal, ref) + + def angle(vi): + d = points[vi] - center + return math.atan2(float(d @ binormal), float(d @ ref)) + + vids.sort(key=angle) + # Ensure CCW with respect to the OUTWARD normal. + ordered = np.array([points[v] for v in vids]) + area_vec = np.zeros(3) + for i in range(len(ordered)): + area_vec += np.cross(ordered[i], ordered[(i + 1) % len(ordered)]) + if area_vec @ normal < 0: + vids.reverse() + faces.append(vids) + return faces + + +class SnubDodecahedron: + def __init__(self): + pts = _all_vertices() + assert len(pts) == 60, f"expected 60 vertices, got {len(pts)}" + + hull = ConvexHull(pts) + # Normalise so that edge length == 1. + faces = _merge_coplanar_faces(pts, hull) + + # Measure a representative edge length and rescale. + i, j = faces[0][0], faces[0][1] + edge_len = np.linalg.norm(pts[i] - pts[j]) + pts = pts / edge_len + + self.vertices = pts + self.faces = faces + + # Build edge -> faces map. + from collections import defaultdict + + edge_faces = defaultdict(list) + for fi, face in enumerate(faces): + n = len(face) + for k in range(n): + a, b = face[k], face[(k + 1) % n] + edge_faces[frozenset((a, b))].append(fi) + self.edge_faces = edge_faces + + def face_normal(self, fi): + face = self.faces[fi] + v = self.vertices[face] + c = v.mean(axis=0) + n = np.cross(v[1] - v[0], v[2] - v[0]) + n = n / np.linalg.norm(n) + # Make it point outward (away from origin, which is the centroid). + if n @ c < 0: + n = -n + return n + + def dihedral_angle(self, fa, fb): + """Interior dihedral angle between two faces (degrees).""" + na, nb = self.face_normal(fa), self.face_normal(fb) + ang_between_normals = math.degrees( + math.acos(max(-1.0, min(1.0, float(na @ nb)))) + ) + return 180.0 - ang_between_normals + + # ---- verification ------------------------------------------------- + def verify(self): + nv = len(self.vertices) + nf = len(self.faces) + ne = len(self.edge_faces) + tris = sum(1 for f in self.faces if len(f) == 3) + pents = sum(1 for f in self.faces if len(f) == 5) + assert nv == 60, nv + assert nf == 92, nf + assert ne == 150, ne + assert tris == 80, tris + assert pents == 12, pents + + # All edges equal length. + lengths = [] + for e in self.edge_faces: + a, b = tuple(e) + lengths.append(np.linalg.norm(self.vertices[a] - self.vertices[b])) + lengths = np.array(lengths) + assert lengths.std() < 1e-6, lengths.std() + assert abs(lengths.mean() - 1.0) < 1e-6, lengths.mean() + + # Every edge shared by exactly 2 faces. + assert all(len(v) == 2 for v in self.edge_faces.values()) + + # Dihedral angles. + d33 = [] + d35 = [] + for e, (fa, fb) in self.edge_faces.items(): + la, lb = len(self.faces[fa]), len(self.faces[fb]) + d = self.dihedral_angle(fa, fb) + if la == 3 and lb == 3: + d33.append(d) + else: + d35.append(d) + # There are no pentagon-pentagon edges. + assert len(d35) == 60, len(d35) # 12 pentagons * 5 edges + assert len(d33) == 90, len(d33) # remaining + m33, m35 = float(np.mean(d33)), float(np.mean(d35)) + assert abs(m33 - 164.1750) < 0.01, m33 + assert abs(m35 - 152.9299) < 0.01, m35 + return { + "vertices": nv, + "faces": nf, + "edges": ne, + "triangles": tris, + "pentagons": pents, + "edge_len": float(lengths.mean()), + "dihedral_3_3": m33, + "dihedral_3_5": m35, + } + + +def build_snub_dodecahedron() -> SnubDodecahedron: + s = SnubDodecahedron() + s.verify() + return s + + +if __name__ == "__main__": + s = build_snub_dodecahedron() + import json + + print(json.dumps(s.verify(), indent=2)) diff --git a/snub-dodecahedron/unfold.py b/snub-dodecahedron/unfold.py new file mode 100644 index 0000000..778cbee --- /dev/null +++ b/snub-dodecahedron/unfold.py @@ -0,0 +1,330 @@ +""" +Unfold the snub dodecahedron into one or more flat, non-overlapping nets. + +We grow each net greedily over the face-adjacency graph. A face is attached +to a growing net by rotating it (in 3D) about the shared edge until it is +coplanar with the rest of the net, then projecting to 2D. If attaching a +face would make it overlap a face already in the net, we leave it for a +later net. The result is a "spanning forest": a small number of connected +flat pieces that together tile the whole surface and can be glued up. + +Each piece records: + - faces: [(face_index, [(x, y), ...]), ...] 2D polygons + - hinges: [Hinge(p, q, dihedral, kind), ...] fold lines (grooved) + - boundary: [Boundary(p, q, dihedral, inward, kind)] cut/glue edges (chamfered) +""" + +from __future__ import annotations + +import math +from dataclasses import dataclass, field + +import numpy as np + +from snubgeom import SnubDodecahedron + + +@dataclass +class Hinge: + p: tuple # 2D endpoint + q: tuple + dihedral: float + kind: str # "3-3" or "3-5" + + +@dataclass +class Boundary: + p: tuple + q: tuple + dihedral: float + inward: tuple # unit 2D vector pointing into the panel (where to chamfer) + kind: str + + +@dataclass +class Piece: + faces: list = field(default_factory=list) + hinges: list = field(default_factory=list) + boundary: list = field(default_factory=list) + + def bbox(self): + xs, ys = [], [] + for _, poly in self.faces: + for x, y in poly: + xs.append(x) + ys.append(y) + return min(xs), min(ys), max(xs), max(ys) + + +# --------------------------------------------------------------------------- +# small linear-algebra helpers (homogeneous 4x4) +# --------------------------------------------------------------------------- +def _skew(v): + x, y, z = v + return np.array([[0, -z, y], [z, 0, -x], [-y, x, 0]], dtype=float) + + +def _rot_align(a, b): + """3x3 rotation taking unit vector a onto unit vector b.""" + a = a / np.linalg.norm(a) + b = b / np.linalg.norm(b) + v = np.cross(a, b) + c = float(a @ b) + s = np.linalg.norm(v) + if s < 1e-12: + if c > 0: + return np.eye(3) + # antiparallel: 180 deg about any axis perpendicular to a + perp = np.array([1.0, 0.0, 0.0]) + if abs(a[0]) > 0.9: + perp = np.array([0.0, 1.0, 0.0]) + axis = np.cross(a, perp) + axis = axis / np.linalg.norm(axis) + K = _skew(axis) + return np.eye(3) + 2 * (K @ K) + K = _skew(v) + return np.eye(3) + K + K @ K * ((1 - c) / (s * s)) + + +def _mat4(R3=None, t=None): + M = np.eye(4) + if R3 is not None: + M[:3, :3] = R3 + if t is not None: + M[:3, 3] = t + return M + + +def _rot_about_line(P, d, angle): + """4x4: rotate by `angle` about the line through point P with direction d.""" + d = d / np.linalg.norm(d) + K = _skew(d) + R3 = np.eye(3) + math.sin(angle) * K + (1 - math.cos(angle)) * (K @ K) + T1 = _mat4(t=-np.asarray(P, float)) + T2 = _mat4(t=np.asarray(P, float)) + return T2 @ _mat4(R3=R3) @ T1 + + +def _apply(M, v): + h = np.array([v[0], v[1], v[2], 1.0]) + return (M @ h)[:3] + + +# --------------------------------------------------------------------------- +# 2D geometry helpers +# --------------------------------------------------------------------------- +def _side(a, b, p): + """Signed area sign of point p relative to directed line a->b.""" + return (b[0] - a[0]) * (p[1] - a[1]) - (b[1] - a[1]) * (p[0] - a[0]) + + +def _poly_centroid(poly): + arr = np.array(poly) + return arr.mean(axis=0) + + +def _shrink(poly, factor): + c = _poly_centroid(poly) + return [tuple(c + (np.array(p) - c) * (1 - factor)) for p in poly] + + +def _convex_overlap(poly_a, poly_b): + """SAT overlap test for two convex polygons (True if they overlap).""" + for poly in (poly_a, poly_b): + n = len(poly) + for i in range(n): + x1, y1 = poly[i] + x2, y2 = poly[(i + 1) % n] + ax, ay = -(y2 - y1), (x2 - x1) # axis = edge normal + amin = min((ax * px + ay * py) for px, py in poly_a) + amax = max((ax * px + ay * py) for px, py in poly_a) + bmin = min((ax * px + ay * py) for px, py in poly_b) + bmax = max((ax * px + ay * py) for px, py in poly_b) + if amax <= bmin or bmax <= amin: + return False # separating axis found + return True + + +# --------------------------------------------------------------------------- +# unfolding +# --------------------------------------------------------------------------- +class Unfolder: + def __init__(self, solid: SnubDodecahedron, overlap_margin=2e-3): + self.s = solid + self.margin = overlap_margin + # face adjacency: fi -> list of (neighbor_fi, (va, vb)) + self.adj = {fi: [] for fi in range(len(solid.faces))} + for e, fs in solid.edge_faces.items(): + a, b = tuple(e) + fa, fb = fs + self.adj[fa].append((fb, (a, b))) + self.adj[fb].append((fa, (a, b))) + + def _flatten_root(self, fi): + """4x4 mapping world -> plane(z=0) for the root face fi.""" + n = self.s.face_normal(fi) + R = _rot_align(n, np.array([0.0, 0.0, 1.0])) + M = _mat4(R3=R) + # translate so face sits on z=0 + v0 = _apply(M, self.s.vertices[self.s.faces[fi][0]]) + return _mat4(t=np.array([0.0, 0.0, -v0[2]])) @ M + + def _poly2d(self, M, fi): + return [tuple(_apply(M, self.s.vertices[v])[:2]) for v in self.s.faces[fi]] + + def _child_matrix(self, parent_M, fa, fb, edge): + """Matrix that places face fb coplanar with the net built around fa.""" + va, vb = edge + P = self.s.vertices[va] + Q = self.s.vertices[vb] + dih = self.s.dihedral_angle(fa, fb) + fold = math.radians(180.0 - dih) + + edge2 = ( + tuple(_apply(parent_M, P)[:2]), + tuple(_apply(parent_M, Q)[:2]), + ) + parent_c = _poly_centroid(self._poly2d(parent_M, fa)) + ps = _side(edge2[0], edge2[1], parent_c) + + best = None + for ang in (fold, -fold): + U = _rot_about_line(P, Q - P, ang) + M = parent_M @ U + poly = self._poly2d(M, fb) + cc = _poly_centroid(poly) + cs = _side(edge2[0], edge2[1], cc) + # planarity check: all z near 0 + zmax = max(abs(_apply(M, self.s.vertices[v])[2]) for v in self.s.faces[fb]) + if zmax < 1e-6 and ps * cs < 0: + best = M + break + if best is None: + # fall back to the better of the two even if marginal + best = M + return best, edge2, dih + + def unfold(self, root_face=None): + n_faces = len(self.s.faces) + # prefer to start pieces on pentagons (nice clusters of 5 triangles) + order = sorted(range(n_faces), key=lambda f: (len(self.s.faces[f]) != 5, f)) + if root_face is not None: + order = [root_face] + [f for f in order if f != root_face] + + unplaced = set(range(n_faces)) + placed_M = {} # fi -> matrix (global, but each piece has own frame) + placed_poly = {} # fi -> 2D polygon + piece_of = {} # fi -> piece index + hinge_records = [] # (piece, edge2, dihedral) + pieces = [] + + for seed in order: + if seed not in unplaced: + continue + pidx = len(pieces) + piece_faces = {} + M = self._flatten_root(seed) + placed_M[seed] = M + poly = self._poly2d(M, seed) + placed_poly[seed] = poly + piece_faces[seed] = poly + piece_of[seed] = pidx + unplaced.discard(seed) + + # grow piece maximally with repeated passes + changed = True + while changed: + changed = False + for f in list(piece_faces.keys()): + for nb, edge in self.adj[f]: + if nb not in unplaced: + continue + M_child, edge2, dih = self._child_matrix( + placed_M[f], f, nb, edge + ) + cand = self._poly2d(M_child, nb) + cand_s = _shrink(cand, self.margin) + overlap = False + for of, opoly in piece_faces.items(): + if of == f: + continue + if _convex_overlap(cand_s, _shrink(opoly, self.margin)): + overlap = True + break + if overlap: + continue + placed_M[nb] = M_child + placed_poly[nb] = cand + piece_faces[nb] = cand + piece_of[nb] = pidx + unplaced.discard(nb) + hinge_records.append((pidx, f, nb, edge2, dih)) + changed = True + + pieces.append(piece_faces) + + # Assemble Piece objects with hinges + boundary edges. + result = [] + hinge_by_piece = {} + hinge_edge_set = {} # piece -> set of frozenset(face pair) that are hinges + for pidx, fa, fb, edge2, dih in hinge_records: + hinge_by_piece.setdefault(pidx, []).append((edge2, dih, fa, fb)) + hinge_edge_set.setdefault(pidx, set()).add(frozenset((fa, fb))) + + for pidx, piece_faces in enumerate(pieces): + piece = Piece() + for fi, poly in piece_faces.items(): + piece.faces.append((fi, poly)) + + for edge2, dih, fa, fb in hinge_by_piece.get(pidx, []): + kind = "3-3" if (len(self.s.faces[fa]) == 3 and len(self.s.faces[fb]) == 3) else "3-5" + piece.hinges.append(Hinge(edge2[0], edge2[1], dih, kind)) + + # boundary edges: every face edge whose partner is NOT a hinge in this piece + hset = hinge_edge_set.get(pidx, set()) + seen_boundary = set() + for fi, poly in piece_faces.items(): + face = self.s.faces[fi] + n = len(face) + for k in range(n): + va, vb = face[k], face[(k + 1) % n] + nb = None + for cand_nb, e in self.adj[fi]: + if frozenset(e) == frozenset((va, vb)): + nb = cand_nb + break + if nb is not None and frozenset((fi, nb)) in hset: + continue # this edge is an internal hinge + key = (fi, frozenset((va, vb))) + if key in seen_boundary: + continue + seen_boundary.add(key) + p2 = poly[k] + q2 = poly[(k + 1) % n] + dih = self.s.dihedral_angle(fi, nb) if nb is not None else 180.0 + c = _poly_centroid(poly) + mid = ((p2[0] + q2[0]) / 2, (p2[1] + q2[1]) / 2) + inward = np.array([c[0] - mid[0], c[1] - mid[1]]) + inward = inward / (np.linalg.norm(inward) + 1e-12) + kind = "3-3" if (nb is not None and len(self.s.faces[fi]) == 3 and len(self.s.faces[nb]) == 3) else "3-5" + piece.boundary.append( + Boundary(p2, q2, dih, tuple(inward), kind) + ) + result.append(piece) + + return result + + +if __name__ == "__main__": + s = SnubDodecahedron() + s.verify() + pieces = Unfolder(s).unfold() + total_faces = sum(len(p.faces) for p in pieces) + print(f"pieces: {len(pieces)}, total faces: {total_faces}") + for i, p in enumerate(pieces): + x0, y0, x1, y1 = p.bbox() + print( + f" piece {i}: {len(p.faces):2d} faces, " + f"{len(p.hinges):2d} hinges, {len(p.boundary):2d} boundary, " + f"size {x1 - x0:.2f} x {y1 - y0:.2f}" + )