1 #!/usr/bin/env python3
3 import argparse
5 import matplotlib.pyplot as plt
7 from math import ceil
9 from matplotlib import pylab
10 from mpl_toolkits.mplot3d import Axes3D
12 from linpy import *
14 x, y, z = symbols('x y z')
16 _x, _y, _z = x.asdummy(), y.asdummy(), z.asdummy()
18 def translate(domain, *, dx=0, dy=0, dz=0):
19 domain &= Polyhedron([x - _x + dx, y - _y + dy, z - _z + dz])
20 domain = domain.project([x, y, z])
21 domain = domain.subs({_x: x, _y: y, _z: z})
22 return domain
24 def _menger(domain, size):
25 result = domain
26 result |= translate(domain, dx=0, dy=size, dz=0)
27 result |= translate(domain, dx=0, dy=2*size, dz=0)
28 result |= translate(domain, dx=size, dy=0, dz=0)
29 result |= translate(domain, dx=size, dy=2*size, dz=0)
30 result |= translate(domain, dx=2*size, dy=0, dz=0)
31 result |= translate(domain, dx=2*size, dy=size, dz=0)
32 result |= translate(domain, dx=2*size, dy=2*size, dz=0)
33 result |= translate(domain, dx=0, dy=0, dz=size)
34 result |= translate(domain, dx=0, dy=2*size, dz=size)
35 result |= translate(domain, dx=2*size, dy=0, dz=size)
36 result |= translate(domain, dx=2*size, dy=2*size, dz=size)
37 result |= translate(domain, dx=0, dy=0, dz=2*size)
38 result |= translate(domain, dx=0, dy=size, dz=2*size)
39 result |= translate(domain, dx=0, dy=2*size, dz=2*size)
40 result |= translate(domain, dx=size, dy=0, dz=2*size)
41 result |= translate(domain, dx=size, dy=2*size, dz=2*size)
42 result |= translate(domain, dx=2*size, dy=0, dz=2*size)
43 result |= translate(domain, dx=2*size, dy=size, dz=2*size)
44 result |= translate(domain, dx=2*size, dy=2*size, dz=2*size)
45 return result
47 def menger(domain, count=1, cut=False):
48 size = 1
49 for i in range(count):
50 domain = _menger(domain, size)
51 size *= 3
52 if cut:
53 domain &= Le(x + y + z, ceil(3 * size / 2))
54 return domain
56 if __name__ == '__main__':
57 parser = argparse.ArgumentParser(
58 description='Compute a Menger sponge.')
59 parser.add_argument('-n', '--iterations', type=int, default=1,
60 help='number of iterations (default: 1)')
61 parser.add_argument('-c', '--cut', action='store_true', default=False,
62 help='cut the sponge')
63 args = parser.parse_args()
64 cube = Le(0, x) & Le(x, 1) & Le(0, y) & Le(y, 1) & Le(0, z) & Le(z, 1)
65 fractal = menger(cube, args.iterations, args.cut)
66 fig = plt.figure(facecolor='white')
67 plot = fig.add_subplot(1, 1, 1, projection='3d', aspect='equal')
68 fractal.plot(plot)
69 pylab.show()