Sierpinski-Menger Sponge Code

by Robert M. Dickau

Developer Support Group
Wolfram Research, Inc.
100 Trade Center Drive
Champaign, IL 61820
email: robertd@wri.com

The parameter "iters," below, is the number of iterations used to create the sponge. Execution with iters=2 can be completed in a reasonable amount of time; execution with iters=3 will use nearly twelve megabytes of memory, and take about ten minutes to compute and twenty minutes to render using a NeXT kernel connected to a Macintosh PowerBook 180 front end.

iters = 3. (* change iters to 2 if you're short on time or RAM;
			  if anyone runs it with iters=4, I'd like to see
			  the result. *);

side = 3. ^ iters ;

cubmat (* cuboid-matrix, that is *) =

        Table[
        
        If[i==side + 1. || j==side + 1. || k==side + 1.,
                
                (* Pad the table's edges with zeroes; if you want
                   to see the complement of the sponge, transpose
                   the 0. and 1. directly below. *)
                   
                0., 1.],
                
                {i,1.,side + 1.},{j,1.,side+1.},{k,1.,side+1.}];

Do[ If[ 
        
        (Mod[Round[i/3.^n + 0.5],3]==2 &&
                (Mod[Round[j/3.^n + 0.5],3]==2 ||
                        Mod[Round[k/3.^n + 0.5],3]==2)) ||
                        
        (Mod[Round[j/3.^n + 0.5],3]==2 &&
                (Mod[Round[i/3.^n + 0.5],3]==2 ||
                        Mod[Round[k/3.^n + 0.5],3]==2)) ||
                        
        (Mod[Round[k/3.^n + 0.5],3]==2 &&
                (Mod[Round[i/3.^n + 0.5],3]==2 ||
                        Mod[Round[j/3.^n + 0.5],3]==2)),
         
         (* then--taking advantage of eightfold symmetry--... *)
         
        (cubmat[[i,j,k]]=0.;
         cubmat[[side+1-i,j,k]]=0.;
          cubmat[[i, side+1-j,k]]=0.;
          cubmat[[i,j,side+1-k]]=0.;
         cubmat[[side+1-i, side+1-j,k]]=0.;
          cubmat[[side+1-i,j,side+1-k]]=0.;
          cubmat[[i, side+1-j, side+1-k]]=0.;
         cubmat[[side+1-i,side+1-j,side+1-k]]=0.;)
         
          (* ...no cuboid goes there *)],
         
                {i,(side+1)/2},{j,(side+1)/2},{k,(side+1)/2},
                	{n,0.,iters-1.}]

faces = {};

              (* Instead of using the Cuboid graphics primitive,
                 we show only the polygons visible from
                 viewpoints in the default octant. *)

Do[ 
        If[ cubmat[[i,j,k]]==1. && cubmat[[i,j,k+1.]]==0.
        
        		(* That is, if a face belongs at {i,j,k}
        		   and there's nothing hiding it, add the
        		   appropriate polygon to the list. *),
        
                AppendTo[ faces,
                
                (* cuboid tops... *)
                {{i,j,k+1.},{i,j+1.,k+1.},
                        {i+1.,j+1.,k+1.},{i+1.,j,k+1.}}
                        
                        ] ],
                        
        {i,1.,side},{j,1.,side},{k,1.,side} ];

(* Since the figure looks the same regardless of which axis
   is vertical, the polygon-corner list "faces" is computed
   only for the tops of the cuboids, then rotated twice to get
   lists of sides and fronts. *)
   
faces = Join[ faces (*tops*),
        (*sides *)	Map[ RotateLeft[#,2]&,faces,{2}],
        (*fronts*)	Map[ 
        	RotateLeft[#,1]*{1,-1,1}+{0,side+2,0}&,
        		faces,{2}]
        			];

Show[Graphics3D[ {EdgeForm[], Map[ Polygon, faces ]}],
                Boxed->False]

28 January 1994