New example: Menger sponge, to be plotted

diff --git a/examples/menger.py b/examples/menger.py
new file mode 100755 (executable)
index 0000000..ee48f5d
--- /dev/null
+++ b/examples/menger.py
@@ -0,0 +1,59 @@
+#!/usr/bin/env python3
+
+import argparse
+
+from pypol import *
+
+x, y, z = symbols('x y z')
+
+_x, _y, _z = x.asdummy(), y.asdummy(), z.asdummy()
+
+def translate(domain, *, dx=0, dy=0, dz=0):
+    domain &= Polyhedron([x - _x + dx, y - _y + dy, z - _z + dz])
+    domain = domain.project([x, y, z])
+    domain = domain.subs({_x: x, _y: y, _z: z})
+    return domain
+
+def _menger(domain):
+
+    result = domain
+    result |= translate(domain, dx=0, dy=1, dz=0)
+    result |= translate(domain, dx=0, dy=2, dz=0)
+    result |= translate(domain, dx=1, dy=0, dz=0)
+    result |= translate(domain, dx=1, dy=2, dz=0)
+    result |= translate(domain, dx=2, dy=0, dz=0)
+    result |= translate(domain, dx=2, dy=1, dz=0)
+    result |= translate(domain, dx=2, dy=2, dz=0)
+
+    result |= translate(domain, dx=0, dy=0, dz=1)
+    result |= translate(domain, dx=0, dy=2, dz=1)
+    result |= translate(domain, dx=2, dy=0, dz=1)
+    result |= translate(domain, dx=2, dy=2, dz=1)
+
+    result |= translate(domain, dx=0, dy=0, dz=2)
+    result |= translate(domain, dx=0, dy=1, dz=2)
+    result |= translate(domain, dx=0, dy=2, dz=2)
+    result |= translate(domain, dx=1, dy=0, dz=2)
+    result |= translate(domain, dx=1, dy=2, dz=2)
+    result |= translate(domain, dx=2, dy=0, dz=2)
+    result |= translate(domain, dx=2, dy=1, dz=2)
+    result |= translate(domain, dx=2, dy=2, dz=2)
+
+    return result
+
+def menger(domain, count=1):
+    for i in range(count):
+        domain = _menger(domain)
+    return domain
+
+if __name__ == '__main__':
+    parser = argparse.ArgumentParser(
+        description='Compute a Menger sponge.')
+    parser.add_argument('-n', '--iterations', type=int, default=1,
+        help='number of iterations (default: 1)')
+    args = parser.parse_args()
+    cube = Le(0, x) & Le(x, 1) & Le(0, y) & Le(y, 1) & Le(0, z) & Le(z, 1)
+    fractal = menger(cube, args.iterations)
+    print('Menger sponge:')
+    print(fractal)
+    print('Number of polyhedra: {}'.format(len(fractal.polyhedra)))