Making a fractal in R3, a Sierpinski pyramid

First, we define some useful routines. Tetra takes a list of the verticies of a tetrahedron, and returns a list of its faces (4 triangles).

Tetra := proc(verts) 
       local a,b,c,d; 
       a := verts[1]; 
       b := verts[2]; 
       c := verts[3];
       d := verts[4];
       RETURN( [[a,b,c], [a,b,d], [a,c,d], [b,c,d]]);
end:
Midpoint returns the midpoint of the segment joining two points.

Midpoint := proc(a,b) 
    RETURN( [(a[1]+b[1])/2, (a[2]+b[2])/2, (a[3]+b[3])/2]); 
  end:
Now for the guts.
Refine takes a tetrahedron (a list of 4 points) and returns the 4 smaller tetrahedra defined by those points and the midpoints of each edge.

Refine := proc (tetra) 
      local a,b,c,d; 
      a:=tetra[1]; 
      b:=tetra[2]; 
      c:=tetra[3]; 
      d:=tetra[4]; 
       RETURN( [a, Midpoint(a,b), Midpoint(a,c), Midpoint(a,d)],  
               [b, Midpoint(a,b), Midpoint(b,c), Midpoint(b,d)], 
               [c, Midpoint(c,b), Midpoint(a,c), Midpoint(c,d)], 
               [d, Midpoint(d,b), Midpoint(d,c), Midpoint(a,d)]); 
      end:
RefineAll applies Refine to each element of a list of tetrahedra. Note that it uses map. we could have written it with seq, but this is shorter and has the same effect.

RefineAll := proc(tlist) 
    RETURN(map(Refine,tlist)); 
  end:  
Recurse is a general purpose routine that computes the n-fold composition of f applied to the second argument. that is, f(f(f(f(...f(input)...))))

Recurse := proc(f, input, n) 
    if ( n <= 0) then 
       RETURN(input); 
    else 
       RETURN(Recurse(f, f(input), n-1)); 
    fi; 
  end:

Now we can put it all together. FirstT is our initial tetrahedron, equilateral of side length 2. I know the height looks funny, but thats it. Really.

FirstT := [[-1,-1,0],[0,sqrt(3)-1,0],[1,-1,0],[0,sqrt(3)/3-1,2*sqrt(15)/5]]:
Here we go. Here's the 3rd stage. (Notice that we have to map Tetra onto each element of the list that comes out of the recursion, so that we will have triangles instead of the verticies of tetrahedra.

with(plots):

setoptions3d(axes=none, scaling=constrained,shading=XY,style=patch,orientation=[60,48]);

polygonplot3d(map(Tetra,[Refine(FirstT)]),title=`Sierpinski Pyramid, level 1`); 

Things get a little difficult to see when you start refining, although if you make the picture in Maple and can spin it, it helps a lot. Below is level 4, and if you click on the image, you get level 5.

polygonplot3d(map(Tetra,Recurse(RefineAll,[FirstT],4)),title=`Sierpinski Pyramid, level 4`);
What would you suspect the dimension of the limiting object is? (remember we started with solid pyramids.) If you do the calculation, you find it is log 4/log2, that is, its 2! so we have a fractal with an integer dimension (its just the "wrong" one, from a topological perspective).