8d258b13440769cd152d24901657e97cd482259d
[linpy.git] / examples / menger.py
1 #!/usr/bin/env python3
2
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.
15
16 import argparse
17
18 import matplotlib.pyplot as plt
19
20 from math import ceil
21
22 from matplotlib import pylab
23
24 from linpy import Le, Polyhedron, symbols
25
26
27 x, y, z = symbols('x y z')
28
29 _x, _y, _z = x.asdummy(), y.asdummy(), z.asdummy()
30
31
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
37
38
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
61
62
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
71
72 if __name__ == '__main__':
73 parser = argparse.ArgumentParser(
74 description='Compute a Menger sponge.')
75 parser.add_argument(
76 '-n', '--iterations', type=int, default=2,
77 help='number of iterations (default: 2)')
78 parser.add_argument(
79 '-c', '--cut', action='store_true', default=False,
80 help='cut the sponge')
81 args = parser.parse_args()
82 cube = Le(0, x) & Le(x, 1) & Le(0, y) & Le(y, 1) & Le(0, z) & Le(z, 1)
83 fractal = menger(cube, args.iterations, args.cut)
84 fig = plt.figure(facecolor='white')
85 plot = fig.add_subplot(1, 1, 1, projection='3d', aspect='equal')
86 fractal.plot(plot)
87 pylab.show()