1 #!/usr/bin/env python3
3 import argparse
5 from pypol import *
7 x, y, z = symbols('x y z')
9 _x, _y, _z = x.asdummy(), y.asdummy(), z.asdummy()
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
17 def _menger(domain):
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)
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)
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)
42 return result
44 def menger(domain, count=1):
45 for i in range(count):
46 domain = _menger(domain)
47 return domain
49 if __name__ == '__main__':
50 parser = argparse.ArgumentParser(
51 description='Compute a Menger sponge.')