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
23 from mpl_toolkits.mplot3d import Axes3D
25 from linpy import Le, Polyhedron, symbols
28 x, y, z = symbols('x y z')
30 _x, _y, _z = x.asdummy(), y.asdummy(), z.asdummy()
33 def translate(domain, *, dx=0, dy=0, dz=0):
34 domain &= Polyhedron([x - _x + dx, y - _y + dy, z - _z + dz])
35 domain = domain.project([x, y, z])
36 domain = domain.subs({_x: x, _y: y, _z: z})
37 return domain
40 def _menger(domain, size):
41 result = domain
42 result |= translate(domain, dx=0, dy=size, dz=0)
43 result |= translate(domain, dx=0, dy=2*size, dz=0)
44 result |= translate(domain, dx=size, dy=0, dz=0)
45 result |= translate(domain, dx=size, dy=2*size, dz=0)
46 result |= translate(domain, dx=2*size, dy=0, dz=0)
47 result |= translate(domain, dx=2*size, dy=size, dz=0)
48 result |= translate(domain, dx=2*size, dy=2*size, dz=0)
49 result |= translate(domain, dx=0, dy=0, dz=size)
50 result |= translate(domain, dx=0, dy=2*size, dz=size)
51 result |= translate(domain, dx=2*size, dy=0, dz=size)
52 result |= translate(domain, dx=2*size, dy=2*size, dz=size)
53 result |= translate(domain, dx=0, dy=0, dz=2*size)
54 result |= translate(domain, dx=0, dy=size, dz=2*size)
55 result |= translate(domain, dx=0, dy=2*size, dz=2*size)
56 result |= translate(domain, dx=size, dy=0, dz=2*size)
57 result |= translate(domain, dx=size, dy=2*size, dz=2*size)
58 result |= translate(domain, dx=2*size, dy=0, dz=2*size)
59 result |= translate(domain, dx=2*size, dy=size, dz=2*size)
60 result |= translate(domain, dx=2*size, dy=2*size, dz=2*size)
61 return result
64 def menger(domain, count=1, cut=False):
65 size = 1
66 for i in range(count):
67 domain = _menger(domain, size)
68 size *= 3
69 if cut:
70 domain &= Le(x + y + z, ceil(3 * size / 2))
71 return domain
73 if __name__ == '__main__':
74 parser = argparse.ArgumentParser(
75 description='Compute a Menger sponge.')