Fix Domain.make_disjoint()
[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 from mpl_toolkits.mplot3d import Axes3D
24
25 from linpy import *
26
27
28 x, y, z = symbols('x y z')
29
30 _x, _y, _z = x.asdummy(), y.asdummy(), z.asdummy()
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 def _menger(domain, size):
39 result = domain
40 result |= translate(domain, dx=0, dy=size, dz=0)
41 result |= translate(domain, dx=0, dy=2*size, dz=0)
42 result |= translate(domain, dx=size, dy=0, dz=0)
43 result |= translate(domain, dx=size, dy=2*size, dz=0)
44 result |= translate(domain, dx=2*size, dy=0, dz=0)
45 result |= translate(domain, dx=2*size, dy=size, dz=0)
46 result |= translate(domain, dx=2*size, dy=2*size, dz=0)
47 result |= translate(domain, dx=0, dy=0, dz=size)
48 result |= translate(domain, dx=0, dy=2*size, dz=size)
49 result |= translate(domain, dx=2*size, dy=0, dz=size)
50 result |= translate(domain, dx=2*size, dy=2*size, dz=size)
51 result |= translate(domain, dx=0, dy=0, dz=2*size)
52 result |= translate(domain, dx=0, dy=size, dz=2*size)
53 result |= translate(domain, dx=0, dy=2*size, dz=2*size)
54 result |= translate(domain, dx=size, dy=0, dz=2*size)
55 result |= translate(domain, dx=size, dy=2*size, dz=2*size)
56 result |= translate(domain, dx=2*size, dy=0, dz=2*size)
57 result |= translate(domain, dx=2*size, dy=size, dz=2*size)
58 result |= translate(domain, dx=2*size, dy=2*size, dz=2*size)
59 return result
60
61 def menger(domain, count=1, cut=False):
62 size = 1
63 for i in range(count):
64 domain = _menger(domain, size)
65 size *= 3
66 if cut:
67 domain &= Le(x + y + z, ceil(3 * size / 2))
68 return domain
69
70 if __name__ == '__main__':
71 parser = argparse.ArgumentParser(
72 description='Compute a Menger sponge.')
73 parser.add_argument('-n', '--iterations', type=int, default=2,
74 help='number of iterations (default: 2)')
75 parser.add_argument('-c', '--cut', action='store_true', default=False,
76 help='cut the sponge')
77 args = parser.parse_args()
78 cube = Le(0, x) & Le(x, 1) & Le(0, y) & Le(y, 1) & Le(0, z) & Le(z, 1)
79 fractal = menger(cube, args.iterations, args.cut)
80 fig = plt.figure(facecolor='white')
81 plot = fig.add_subplot(1, 1, 1, projection='3d', aspect='equal')
82 fractal.plot(plot)
83 pylab.show()