1 #!/usr/bin/env python3
3 import argparse
5 import matplotlib.pyplot as plt
7 from math import ceil
9 from matplotlib import pylab
10 from mpl_toolkits.mplot3d import Axes3D
12 from linpy import *
14 x, y, z = symbols('x y z')
16 _x, _y, _z = x.asdummy(), y.asdummy(), z.asdummy()
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
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
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
56 if __name__ == '__main__':
57 parser = argparse.ArgumentParser(
58 description='Compute a Menger sponge.')