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)))