DEV Community

Black
Black

Posted on

Penrose tiling in Python: generate a 3D-printable aperiodic ornament

A Penrose tiling never repeats and has fivefold symmetry, which is why it looks so odd on a 3D-printed ornament. You can generate one in about 40 lines of Python with manifold3d. The trick is Robinson triangle subdivision: start with a wheel of 10 triangles and split each one by the golden ratio, again and again.

import math, numpy as np
phi = (1 + 5 ** .5) / 2

def subdivide(tris):
    out = []
    for col, A, B, C in tris:          # points are complex numbers
        if col == 0:                   # thin (red) triangle
            P = A + (B - A) / phi
            out += [(0, C, P, B), (1, P, C, A)]
        else:                          # thick (blue) triangle
            Q = B + (A - B) / phi
            R = B + (C - B) / phi
            out += [(1, R, C, A), (1, Q, R, B), (0, R, Q, A)]
    return out

Rd = 29.0
tris = []
for i in range(10):
    B = Rd * 1.3 * np.exp(1j * (2 * i - 1) * math.pi / 10)
    C = Rd * 1.3 * np.exp(1j * (2 * i + 1) * math.pi / 10)
    if i % 2 == 0:
        B, C = C, B
    tris.append((0, 0j, B, C))
for _ in range(4):
    tris = subdivide(tris)
Enter fullscreen mode Exit fullscreen mode

After 4 rounds you have 890 triangles. The rhombus edges of the tiling are the A-B and C-A sides of every triangle, so draw those and skip the rest. Deduplicate shared edges, then thicken each one into a thin rectangle with CrossSection:

from manifold3d import Manifold, CrossSection
import manifold3d as M3

def seg(a, b, w=0.62):
    d = b - a
    sq = CrossSection.square([abs(d), 2 * w], True)
    return sq.rotate(math.degrees(math.atan2(d.imag, d.real))).translate(
        [(a.real + b.real) / 2, (a.imag + b.imag) / 2])

edges, segs = set(), []
for col, A, B, C in tris:
    for p, q in ((A, B), (C, A)):
        k = tuple(sorted(((round(p.real, 3), round(p.imag, 3)),
                          (round(q.real, 3), round(q.imag, 3)))))
        if k not in edges:
            edges.add(k)
            segs.append(seg(p, q))

disc = CrossSection.circle(Rd, 128)
pat = CrossSection.batch_boolean(segs, M3.OpType.Add) ^ disc
rim = disc - CrossSection.circle(Rd - 2.0, 128)
part = (pat + rim).extrude(2.4)      # flat 2.4 mm plate, prints in about 15 minutes
Enter fullscreen mode Exit fullscreen mode

Gotchas

  • Round the endpoints before hashing edges, or floating-point noise leaves duplicate edges and slivers.
  • Use 4 subdivisions for a 58 mm disc. At 5 the tiles are 3 mm wide and clog a 0.4 mm nozzle.
  • Check CrossSection.decompose() for loose islands before you extrude. Keep only the largest piece so nothing floats.
  • Line width of 1.2 mm (3 nozzle widths) prints cleanly with no supports.

If you would rather skip the code, the finished ornament is on Cults3D for $1.27.

Top comments (0)