064219ed7fdf045156cb945656dc213412b156f2
[linpy.git] / examples / menger.py
1 #!/usr/bin/env python3
2 #
3 # Copyright 2014 MINES ParisTech
4 #
5 # This file is part of LinPy.
6 #
7 # LinPy is free software: you can redistribute it and/or modify
8 # it under the terms of the GNU General Public License as published by
9 # the Free Software Foundation, either version 3 of the License, or
10 # (at your option) any later version.
11 #
12 # LinPy is distributed in the hope that it will be useful,
13 # but WITHOUT ANY WARRANTY; without even the implied warranty of
14 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 # GNU General Public License for more details.
16 #
17 # You should have received a copy of the GNU General Public License
18 # along with LinPy. If not, see <http://www.gnu.org/licenses/>.
19
20 import argparse
21
22 import matplotlib.pyplot as plt
23
24 from math import ceil
25
26 from matplotlib import pylab
27 from mpl_toolkits.mplot3d import Axes3D
28
29 from linpy import *
30
31 x, y, z = symbols('x y z')
32
33 _x, _y, _z = x.asdummy(), y.asdummy(), z.asdummy()
34
35 def translate(domain, *, dx=0, dy=0, dz=0):
36 domain &= Polyhedron([x - _x + dx, y - _y + dy, z - _z + dz])
37 domain = domain.project([x, y, z])
38 domain = domain.subs({_x: x, _y: y, _z: z})
39 return domain
40
41 def _menger(domain, size):
42 result = domain
43 result |= translate(domain, dx=0, dy=size, dz=0)
44 result |= translate(domain, dx=0, dy=2*size, dz=0)
45 result |= translate(domain, dx=size, dy=0, dz=0)
46 result |= translate(domain, dx=size, dy=2*size, dz=0)
47 result |= translate(domain, dx=2*size, dy=0, dz=0)
48 result |= translate(domain, dx=2*size, dy=size, dz=0)
49 result |= translate(domain, dx=2*size, dy=2*size, dz=0)
50 result |= translate(domain, dx=0, dy=0, dz=size)
51 result |= translate(domain, dx=0, dy=2*size, dz=size)
52 result |= translate(domain, dx=2*size, dy=0, dz=size)
53 result |= translate(domain, dx=2*size, dy=2*size, dz=size)
54 result |= translate(domain, dx=0, dy=0, dz=2*size)
55 result |= translate(domain, dx=0, dy=size, dz=2*size)
56 result |= translate(domain, dx=0, dy=2*size, dz=2*size)
57 result |= translate(domain, dx=size, dy=0, dz=2*size)
58 result |= translate(domain, dx=size, dy=2*size, dz=2*size)
59 result |= translate(domain, dx=2*size, dy=0, dz=2*size)
60 result |= translate(domain, dx=2*size, dy=size, dz=2*size)
61 result |= translate(domain, dx=2*size, dy=2*size, dz=2*size)
62 return result
63
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
72
73 if __name__ == '__main__':
74 parser = argparse.ArgumentParser(
75 description='Compute a Menger sponge.')
76 parser.add_argument('-n', '--iterations', type=int, default=1,
77 help='number of iterations (default: 1)')
78 parser.add_argument('-c', '--cut', action='store_true', default=False,
79 help='cut the sponge')
80 args = parser.parse_args()
81 cube = Le(0, x) & Le(x, 1) & Le(0, y) & Le(y, 1) & Le(0, z) & Le(z, 1)
82 fractal = menger(cube, args.iterations, args.cut)
83 fig = plt.figure(facecolor='white')
84 plot = fig.add_subplot(1, 1, 1, projection='3d', aspect='equal')
85 fractal.plot(plot)
86 pylab.show()