examples/menger.py

   1 #!/usr/bin/env python3
   2 #
   3 # Copyright 2014 MINES ParisTech
   4 #
   5 # This file is part of LinPy.
   6 #
   7 # LinPy is free software: you can redistribute it and/or modify
   8 # it under the terms of the GNU General Public License as published by
   9 # the Free Software Foundation, either version 3 of the License, or
  10 # (at your option) any later version.
  11 #
  12 # LinPy is distributed in the hope that it will be useful,
  13 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  14 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15 # GNU General Public License for more details.
  16 #
  17 # You should have received a copy of the GNU General Public License
  18 # along with LinPy.  If not, see <http://www.gnu.org/licenses/>.
  19
  20 import argparse
  21
  22 import matplotlib.pyplot as plt
  23
  24 from math import ceil
  25
  26 from matplotlib import pylab
  27 from mpl_toolkits.mplot3d import Axes3D
  28
  29 from linpy import *
  30
  31 x, y, z = symbols('x y z')
  32
  33 _x, _y, _z = x.asdummy(), y.asdummy(), z.asdummy()
  34
  35 def translate(domain, *, dx=0, dy=0, dz=0):
  36     domain &= Polyhedron([x - _x + dx, y - _y + dy, z - _z + dz])
  37     domain = domain.project([x, y, z])
  38     domain = domain.subs({_x: x, _y: y, _z: z})
  39     return domain
  40
  41 def _menger(domain, size):
  42     result = domain
  43     result |= translate(domain, dx=0, dy=size, dz=0)
  44     result |= translate(domain, dx=0, dy=2*size, dz=0)
  45     result |= translate(domain, dx=size, dy=0, dz=0)
  46     result |= translate(domain, dx=size, dy=2*size, dz=0)
  47     result |= translate(domain, dx=2*size, dy=0, dz=0)
  48     result |= translate(domain, dx=2*size, dy=size, dz=0)
  49     result |= translate(domain, dx=2*size, dy=2*size, dz=0)
  50     result |= translate(domain, dx=0, dy=0, dz=size)
  51     result |= translate(domain, dx=0, dy=2*size, dz=size)
  52     result |= translate(domain, dx=2*size, dy=0, dz=size)
  53     result |= translate(domain, dx=2*size, dy=2*size, dz=size)
  54     result |= translate(domain, dx=0, dy=0, dz=2*size)
  55     result |= translate(domain, dx=0, dy=size, dz=2*size)
  56     result |= translate(domain, dx=0, dy=2*size, dz=2*size)
  57     result |= translate(domain, dx=size, dy=0, dz=2*size)
  58     result |= translate(domain, dx=size, dy=2*size, dz=2*size)
  59     result |= translate(domain, dx=2*size, dy=0, dz=2*size)
  60     result |= translate(domain, dx=2*size, dy=size, dz=2*size)
  61     result |= translate(domain, dx=2*size, dy=2*size, dz=2*size)
  62     return result
  63
  64 def menger(domain, count=1, cut=False):
  65     size = 1
  66     for i in range(count):
  67         domain = _menger(domain, size)
  68         size *= 3
  69     if cut:
  70         domain &= Le(x + y + z, ceil(3 * size / 2))
  71     return domain
  72
  73 if __name__ == '__main__':
  74     parser = argparse.ArgumentParser(
  75         description='Compute a Menger sponge.')
  76     parser.add_argument('-n', '--iterations', type=int, default=1,
  77         help='number of iterations (default: 1)')
  78     parser.add_argument('-c', '--cut', action='store_true', default=False,
  79         help='cut the sponge')
  80     args = parser.parse_args()
  81     cube = Le(0, x) & Le(x, 1) & Le(0, y) & Le(y, 1) & Le(0, z) & Le(z, 1)
  82     fractal = menger(cube, args.iterations, args.cut)
  83     fig = plt.figure(facecolor='white')
  84     plot = fig.add_subplot(1, 1, 1, projection='3d', aspect='equal')
  85     fractal.plot(plot)
  86     pylab.show()