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
  24 from linpy import Le, Polyhedron, symbols
  25
  26
  27 x, y, z = symbols('x y z')
  28
  29 _x, _y, _z = x.asdummy(), y.asdummy(), z.asdummy()
  30
  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
  39 def _menger(domain, size):
  40     result = domain
  41     result |= translate(domain, dx=0, dy=size, dz=0)
  42     result |= translate(domain, dx=0, dy=2*size, dz=0)
  43     result |= translate(domain, dx=size, dy=0, dz=0)
  44     result |= translate(domain, dx=size, dy=2*size, dz=0)
  45     result |= translate(domain, dx=2*size, dy=0, dz=0)
  46     result |= translate(domain, dx=2*size, dy=size, dz=0)
  47     result |= translate(domain, dx=2*size, dy=2*size, dz=0)
  48     result |= translate(domain, dx=0, dy=0, dz=size)
  49     result |= translate(domain, dx=0, dy=2*size, dz=size)
  50     result |= translate(domain, dx=2*size, dy=0, dz=size)
  51     result |= translate(domain, dx=2*size, dy=2*size, dz=size)
  52     result |= translate(domain, dx=0, dy=0, dz=2*size)
  53     result |= translate(domain, dx=0, dy=size, dz=2*size)
  54     result |= translate(domain, dx=0, dy=2*size, dz=2*size)
  55     result |= translate(domain, dx=size, dy=0, dz=2*size)
  56     result |= translate(domain, dx=size, dy=2*size, dz=2*size)
  57     result |= translate(domain, dx=2*size, dy=0, dz=2*size)
  58     result |= translate(domain, dx=2*size, dy=size, dz=2*size)
  59     result |= translate(domain, dx=2*size, dy=2*size, dz=2*size)
  60     return result
  61
  62
  63 def menger(domain, count=1, cut=False):
  64     size = 1
  65     for i in range(count):
  66         domain = _menger(domain, size)
  67         size *= 3
  68     if cut:
  69         domain &= Le(x + y + z, ceil(3 * size / 2))
  70     return domain
  71
  72 if __name__ == '__main__':
  73     parser = argparse.ArgumentParser(
  74         description='Compute a Menger sponge.')
  75     parser.add_argument(
  76         '-n', '--iterations', type=int, default=2,
  77         help='number of iterations (default: 2)')
  78     parser.add_argument(
  79         '-c', '--cut', action='store_true', default=False,
  80         help='cut the sponge')
  81     args = parser.parse_args()
  82     cube = Le(0, x) & Le(x, 1) & Le(0, y) & Le(y, 1) & Le(0, z) & Le(z, 1)
  83     fractal = menger(cube, args.iterations, args.cut)
  84     fig = plt.figure(facecolor='white')
  85     plot = fig.add_subplot(1, 1, 1, projection='3d', aspect='equal')
  86     fractal.plot(plot)
  87     pylab.show()