X-Git-Url: https://scm.cri.ensmp.fr/git/linpy.git/blobdiff_plain/8b80ed2263c6e0cf7d6589f36eb3338eaeaaa7e3..51e97eade63b2f4c7b500feb503436cc4a886e59:/examples/menger.py

diff --git a/examples/menger.py b/examples/menger.py
index ee48f5d..b6ab90d 100755
--- a/examples/menger.py
+++ b/examples/menger.py
@@ -1,7 +1,31 @@
 #!/usr/bin/env python3
 
+"""
+    This file is part of Linpy.
+
+    Linpy is free software: you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    Linpy is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with Linpy.  If not, see <http://www.gnu.org/licenses/>.
+"""
+
 import argparse
 
+import matplotlib.pyplot as plt
+
+from math import ceil
+
+from matplotlib import pylab
+from mpl_toolkits.mplot3d import Axes3D
+
 from pypol import *
 
 x, y, z = symbols('x y z')
@@ -14,36 +38,36 @@ def translate(domain, *, dx=0, dy=0, dz=0):
     domain = domain.subs({_x: x, _y: y, _z: z})
     return domain
 
-def _menger(domain):
-
+def _menger(domain, size):
     result = domain
-    result |= translate(domain, dx=0, dy=1, dz=0)
-    result |= translate(domain, dx=0, dy=2, dz=0)
-    result |= translate(domain, dx=1, dy=0, dz=0)
-    result |= translate(domain, dx=1, dy=2, dz=0)
-    result |= translate(domain, dx=2, dy=0, dz=0)
-    result |= translate(domain, dx=2, dy=1, dz=0)
-    result |= translate(domain, dx=2, dy=2, dz=0)
-
-    result |= translate(domain, dx=0, dy=0, dz=1)
-    result |= translate(domain, dx=0, dy=2, dz=1)
-    result |= translate(domain, dx=2, dy=0, dz=1)
-    result |= translate(domain, dx=2, dy=2, dz=1)
-
-    result |= translate(domain, dx=0, dy=0, dz=2)
-    result |= translate(domain, dx=0, dy=1, dz=2)
-    result |= translate(domain, dx=0, dy=2, dz=2)
-    result |= translate(domain, dx=1, dy=0, dz=2)
-    result |= translate(domain, dx=1, dy=2, dz=2)
-    result |= translate(domain, dx=2, dy=0, dz=2)
-    result |= translate(domain, dx=2, dy=1, dz=2)
-    result |= translate(domain, dx=2, dy=2, dz=2)
-
+    result |= translate(domain, dx=0, dy=size, dz=0)
+    result |= translate(domain, dx=0, dy=2*size, dz=0)
+    result |= translate(domain, dx=size, dy=0, dz=0)
+    result |= translate(domain, dx=size, dy=2*size, dz=0)
+    result |= translate(domain, dx=2*size, dy=0, dz=0)
+    result |= translate(domain, dx=2*size, dy=size, dz=0)
+    result |= translate(domain, dx=2*size, dy=2*size, dz=0)
+    result |= translate(domain, dx=0, dy=0, dz=size)
+    result |= translate(domain, dx=0, dy=2*size, dz=size)
+    result |= translate(domain, dx=2*size, dy=0, dz=size)
+    result |= translate(domain, dx=2*size, dy=2*size, dz=size)
+    result |= translate(domain, dx=0, dy=0, dz=2*size)
+    result |= translate(domain, dx=0, dy=size, dz=2*size)
+    result |= translate(domain, dx=0, dy=2*size, dz=2*size)
+    result |= translate(domain, dx=size, dy=0, dz=2*size)
+    result |= translate(domain, dx=size, dy=2*size, dz=2*size)
+    result |= translate(domain, dx=2*size, dy=0, dz=2*size)
+    result |= translate(domain, dx=2*size, dy=size, dz=2*size)
+    result |= translate(domain, dx=2*size, dy=2*size, dz=2*size)
     return result
 
-def menger(domain, count=1):
+def menger(domain, count=1, cut=False):
+    size = 1
     for i in range(count):
-        domain = _menger(domain)
+        domain = _menger(domain, size)
+        size *= 3
+    if cut:
+        domain &= Le(x + y + z, ceil(3 * size / 2))
     return domain
 
 if __name__ == '__main__':
@@ -51,9 +75,14 @@ if __name__ == '__main__':
         description='Compute a Menger sponge.')
     parser.add_argument('-n', '--iterations', type=int, default=1,
         help='number of iterations (default: 1)')
+    parser.add_argument('-c', '--cut', action='store_true', default=False,
+        help='cut the sponge')
     args = parser.parse_args()
     cube = Le(0, x) & Le(x, 1) & Le(0, y) & Le(y, 1) & Le(0, z) & Le(z, 1)
-    fractal = menger(cube, args.iterations)
-    print('Menger sponge:')
-    print(fractal)
-    print('Number of polyhedra: {}'.format(len(fractal.polyhedra)))
+    fractal = menger(cube, args.iterations, args.cut)
+    fig = plt.figure(facecolor='white')
+    plot = fig.add_subplot(1, 1, 1, projection='3d', aspect='equal')
+    fractal.plot(plot)
+    pylab.show()
+    
+# Copyright 2014 MINES ParisTech