New example: Menger sponge, to be plotted
[linpy.git] / examples / menger.py
1 #!/usr/bin/env python3
2
3 import argparse
4
5 from pypol import *
6
7 x, y, z = symbols('x y z')
8
9 _x, _y, _z = x.asdummy(), y.asdummy(), z.asdummy()
10
11 def translate(domain, *, dx=0, dy=0, dz=0):
12 domain &= Polyhedron([x - _x + dx, y - _y + dy, z - _z + dz])
13 domain = domain.project([x, y, z])
14 domain = domain.subs({_x: x, _y: y, _z: z})
15 return domain
16
17 def _menger(domain):
18
19 result = domain
20 result |= translate(domain, dx=0, dy=1, dz=0)
21 result |= translate(domain, dx=0, dy=2, dz=0)
22 result |= translate(domain, dx=1, dy=0, dz=0)
23 result |= translate(domain, dx=1, dy=2, dz=0)
24 result |= translate(domain, dx=2, dy=0, dz=0)
25 result |= translate(domain, dx=2, dy=1, dz=0)
26 result |= translate(domain, dx=2, dy=2, dz=0)
27
28 result |= translate(domain, dx=0, dy=0, dz=1)
29 result |= translate(domain, dx=0, dy=2, dz=1)
30 result |= translate(domain, dx=2, dy=0, dz=1)
31 result |= translate(domain, dx=2, dy=2, dz=1)
32
33 result |= translate(domain, dx=0, dy=0, dz=2)
34 result |= translate(domain, dx=0, dy=1, dz=2)
35 result |= translate(domain, dx=0, dy=2, dz=2)
36 result |= translate(domain, dx=1, dy=0, dz=2)
37 result |= translate(domain, dx=1, dy=2, dz=2)
38 result |= translate(domain, dx=2, dy=0, dz=2)
39 result |= translate(domain, dx=2, dy=1, dz=2)
40 result |= translate(domain, dx=2, dy=2, dz=2)
41
42 return result
43
44 def menger(domain, count=1):
45 for i in range(count):
46 domain = _menger(domain)
47 return domain
48
49 if __name__ == '__main__':
50 parser = argparse.ArgumentParser(
51 description='Compute a Menger sponge.')
52 parser.add_argument('-n', '--iterations', type=int, default=1,
53 help='number of iterations (default: 1)')
54 args = parser.parse_args()
55 cube = Le(0, x) & Le(x, 1) & Le(0, y) & Le(y, 1) & Le(0, z) & Le(z, 1)
56 fractal = menger(cube, args.iterations)
57 print('Menger sponge:')
58 print(fractal)
59 print('Number of polyhedra: {}'.format(len(fractal.polyhedra)))