examples/menger.py

   1 #!/usr/bin/env python3
   2
   3 import argparse
   4
   5 import matplotlib.pyplot as plt
   6
   7 from math import ceil
   8
   9 from matplotlib import pylab
  10 from mpl_toolkits.mplot3d import Axes3D
  11
  12 from linpy import *
  13
  14 x, y, z = symbols('x y z')
  15
  16 _x, _y, _z = x.asdummy(), y.asdummy(), z.asdummy()
  17
  18 def translate(domain, *, dx=0, dy=0, dz=0):
  19     domain &= Polyhedron([x - _x + dx, y - _y + dy, z - _z + dz])
  20     domain = domain.project([x, y, z])
  21     domain = domain.subs({_x: x, _y: y, _z: z})
  22     return domain
  23
  24 def _menger(domain, size):
  25     result = domain
  26     result |= translate(domain, dx=0, dy=size, dz=0)
  27     result |= translate(domain, dx=0, dy=2*size, dz=0)
  28     result |= translate(domain, dx=size, dy=0, dz=0)
  29     result |= translate(domain, dx=size, dy=2*size, dz=0)
  30     result |= translate(domain, dx=2*size, dy=0, dz=0)
  31     result |= translate(domain, dx=2*size, dy=size, dz=0)
  32     result |= translate(domain, dx=2*size, dy=2*size, dz=0)
  33     result |= translate(domain, dx=0, dy=0, dz=size)
  34     result |= translate(domain, dx=0, dy=2*size, dz=size)
  35     result |= translate(domain, dx=2*size, dy=0, dz=size)
  36     result |= translate(domain, dx=2*size, dy=2*size, dz=size)
  37     result |= translate(domain, dx=0, dy=0, dz=2*size)
  38     result |= translate(domain, dx=0, dy=size, dz=2*size)
  39     result |= translate(domain, dx=0, dy=2*size, dz=2*size)
  40     result |= translate(domain, dx=size, dy=0, dz=2*size)
  41     result |= translate(domain, dx=size, dy=2*size, dz=2*size)
  42     result |= translate(domain, dx=2*size, dy=0, dz=2*size)
  43     result |= translate(domain, dx=2*size, dy=size, dz=2*size)
  44     result |= translate(domain, dx=2*size, dy=2*size, dz=2*size)
  45     return result
  46
  47 def menger(domain, count=1, cut=False):
  48     size = 1
  49     for i in range(count):
  50         domain = _menger(domain, size)
  51         size *= 3
  52     if cut:
  53         domain &= Le(x + y + z, ceil(3 * size / 2))
  54     return domain
  55
  56 if __name__ == '__main__':
  57     parser = argparse.ArgumentParser(
  58         description='Compute a Menger sponge.')
  59     parser.add_argument('-n', '--iterations', type=int, default=1,
  60         help='number of iterations (default: 1)')
  61     parser.add_argument('-c', '--cut', action='store_true', default=False,
  62         help='cut the sponge')
  63     args = parser.parse_args()
  64     cube = Le(0, x) & Le(x, 1) & Le(0, y) & Le(y, 1) & Le(0, z) & Le(z, 1)
  65     fractal = menger(cube, args.iterations, args.cut)
  66     fig = plt.figure(facecolor='white')
  67     plot = fig.add_subplot(1, 1, 1, projection='3d', aspect='equal')
  68     fractal.plot(plot)
  69     pylab.show()