examples/menger.py

   1 #!/usr/bin/env python3
   2
   3 # Plot a Menger sponge.
   4 #
   5 # The construction of a Menger sponge can be described as follows:
   6 #
   7 # 1. Begin with a cube.
   8 # 2. Divide every face of the cube into 9 squares, like a Rubik's Cube. This
   9 #    will sub-divide the cube into 27 smaller cubes.
  10 # 3. Remove the smaller cube in the middle of each face, and remove the smaller
  11 #    cube in the very center of the larger cube, leaving 20 smaller cubes. This
  12 #    is a level-1 Menger sponge (resembling a Void Cube).
  13 # 4. Repeat steps 2 and 3 for each of the remaining smaller cubes, and continue
  14 #    to iterate.
  15
  16 import argparse
  17
  18 import matplotlib.pyplot as plt
  19
  20 from math import ceil
  21
  22 from matplotlib import pylab
  23 from mpl_toolkits.mplot3d import Axes3D
  24
  25 from linpy import *
  26
  27
  28 x, y, z = symbols('x y z')
  29
  30 _x, _y, _z = x.asdummy(), y.asdummy(), z.asdummy()
  31
  32 def translate(domain, *, dx=0, dy=0, dz=0):
  33     domain &= Polyhedron([x - _x + dx, y - _y + dy, z - _z + dz])
  34     domain = domain.project([x, y, z])
  35     domain = domain.subs({_x: x, _y: y, _z: z})
  36     return domain
  37
  38 def _menger(domain, size):
  39     result = domain
  40     result |= translate(domain, dx=0, dy=size, dz=0)
  41     result |= translate(domain, dx=0, dy=2*size, dz=0)
  42     result |= translate(domain, dx=size, dy=0, dz=0)
  43     result |= translate(domain, dx=size, dy=2*size, dz=0)
  44     result |= translate(domain, dx=2*size, dy=0, dz=0)
  45     result |= translate(domain, dx=2*size, dy=size, dz=0)
  46     result |= translate(domain, dx=2*size, dy=2*size, dz=0)
  47     result |= translate(domain, dx=0, dy=0, dz=size)
  48     result |= translate(domain, dx=0, dy=2*size, dz=size)
  49     result |= translate(domain, dx=2*size, dy=0, dz=size)
  50     result |= translate(domain, dx=2*size, dy=2*size, dz=size)
  51     result |= translate(domain, dx=0, dy=0, dz=2*size)
  52     result |= translate(domain, dx=0, dy=size, dz=2*size)
  53     result |= translate(domain, dx=0, dy=2*size, dz=2*size)
  54     result |= translate(domain, dx=size, dy=0, dz=2*size)
  55     result |= translate(domain, dx=size, dy=2*size, dz=2*size)
  56     result |= translate(domain, dx=2*size, dy=0, dz=2*size)
  57     result |= translate(domain, dx=2*size, dy=size, dz=2*size)
  58     result |= translate(domain, dx=2*size, dy=2*size, dz=2*size)
  59     return result
  60
  61 def menger(domain, count=1, cut=False):
  62     size = 1
  63     for i in range(count):
  64         domain = _menger(domain, size)
  65         size *= 3
  66     if cut:
  67         domain &= Le(x + y + z, ceil(3 * size / 2))
  68     return domain
  69
  70 if __name__ == '__main__':
  71     parser = argparse.ArgumentParser(
  72         description='Compute a Menger sponge.')
  73     parser.add_argument('-n', '--iterations', type=int, default=2,
  74         help='number of iterations (default: 2)')
  75     parser.add_argument('-c', '--cut', action='store_true', default=False,
  76         help='cut the sponge')
  77     args = parser.parse_args()
  78     cube = Le(0, x) & Le(x, 1) & Le(0, y) & Le(y, 1) & Le(0, z) & Le(z, 1)
  79     fractal = menger(cube, args.iterations, args.cut)
  80     fig = plt.figure(facecolor='white')
  81     plot = fig.add_subplot(1, 1, 1, projection='3d', aspect='equal')
  82     fractal.plot(plot)
  83     pylab.show()