1 #!/usr/bin/env python3
3 # Plot a Menger sponge.
4 #
5 # The construction of a Menger sponge can be described as follows:
6 #
7 # 1. Begin with a cube.
8 # 2. Divide every face of the cube into 9 squares, like a Rubik's Cube. This
9 # will sub-divide the cube into 27 smaller cubes.
10 # 3. Remove the smaller cube in the middle of each face, and remove the smaller
11 # cube in the very center of the larger cube, leaving 20 smaller cubes. This
12 # is a level-1 Menger sponge (resembling a Void Cube).
13 # 4. Repeat steps 2 and 3 for each of the remaining smaller cubes, and continue
14 # to iterate.
16 import argparse
18 import matplotlib.pyplot as plt
20 from math import ceil
22 from matplotlib import pylab
24 from linpy import Le, Polyhedron, symbols
27 x, y, z = symbols('x y z')
29 _x, _y, _z = x.asdummy(), y.asdummy(), z.asdummy()
32 def translate(domain, *, dx=0, dy=0, dz=0):
33 domain &= Polyhedron([x - _x + dx, y - _y + dy, z - _z + dz])
34 domain = domain.project([x, y, z])
35 domain = domain.subs({_x: x, _y: y, _z: z})
36 return domain
39 def _menger(domain, size):
40 result = domain
41 result |= translate(domain, dx=0, dy=size, dz=0)
42 result |= translate(domain, dx=0, dy=2*size, dz=0)
43 result |= translate(domain, dx=size, dy=0, dz=0)
44 result |= translate(domain, dx=size, dy=2*size, dz=0)
45 result |= translate(domain, dx=2*size, dy=0, dz=0)
46 result |= translate(domain, dx=2*size, dy=size, dz=0)
47 result |= translate(domain, dx=2*size, dy=2*size, dz=0)
48 result |= translate(domain, dx=0, dy=0, dz=size)
49 result |= translate(domain, dx=0, dy=2*size, dz=size)
50 result |= translate(domain, dx=2*size, dy=0, dz=size)
51 result |= translate(domain, dx=2*size, dy=2*size, dz=size)
52 result |= translate(domain, dx=0, dy=0, dz=2*size)
53 result |= translate(domain, dx=0, dy=size, dz=2*size)
54 result |= translate(domain, dx=0, dy=2*size, dz=2*size)
55 result |= translate(domain, dx=size, dy=0, dz=2*size)
56 result |= translate(domain, dx=size, dy=2*size, dz=2*size)
57 result |= translate(domain, dx=2*size, dy=0, dz=2*size)
58 result |= translate(domain, dx=2*size, dy=size, dz=2*size)
59 result |= translate(domain, dx=2*size, dy=2*size, dz=2*size)
60 return result
63 def menger(domain, count=1, cut=False):
64 size = 1
65 for i in range(count):
66 domain = _menger(domain, size)
67 size *= 3
68 if cut:
69 domain &= Le(x + y + z, ceil(3 * size / 2))
70 return domain
72 if __name__ == '__main__':
73 parser = argparse.ArgumentParser(
74 description='Compute a Menger sponge.')