% This script plots up a zero-mean-curvature surface with the IWP topology % % We make use of the level set approximation to the minimal surface. % The equations were taken from http://www.msri.org/publications/sgp/jim/geom/level/library/triper/mainc.html % This is the work of James Hoffman at the MSRI. - jim@msri.org % Other references are, % a) Paul Harper's thesis (1996) Princeton. % b) Piotr Garstecki and Robert Holyst papers. % % % BRIEF INTRO ON THE SURFACE CURVATURE % Consider a point at the top of a hemi-sphere (radius R). % Standing on top of the hemisphere, no matter what direction you walk, you go down. % If we walk off in the x-direction, or the y-direction, the surface curves away a rate 1/R. % The mean curvature of the point on top of the hemisphere is 1/R+1/R = 2/R % % Consider a point on the top of an infinite hemi-cylinder (radius R) (axis along x-axis). % Along the "top" of the hemisphere (x-axis) the surface does not fall at all so it has zero curvature. % Across the hemisphere (y-direction) the surface falls away at a rate 1/R. % Thus, the mean curvature is 0 + 1/R = 1/R. % % Consider now a saddle. Moving on the x-axis (along the spine of the horse), it rises up in either direction at a rate 1/Rx. % However, moving in the y-direction (the girth of the horse) it falls in either direction at a rate 1/Ry. % In the x-axis the curvature is -1/Rx, and y-axis the curvature is 1/Ry. % The mean curvature is 1/Ry - 1/Rx. % If Rx=Ry, the mean curvature is zero. % A surface need not be flat to have zero mean curvature. % % BRIEF INTRO TO THE IWP % The IWP is a surface with the following properties % 1. It is triply periodic (ie. the surface repeats in the x, y and z directions) with space-group Im3m. % 2. The surface does not intersect itself at any point. % 3. At every point on the surface, the mean curvature is zero. % This is not a unique definition of the IWP, but it's a good start. % % The IWP divides space into two separate, continuous networks. % % NETWORK ONE - 4-fold connectors % These are located at % 0, 1/2, 1/2 (connected to 0, 3/2, 1/2 ; 0 -1/2, 1/2 ; 0 1/2 3/2 ; 0 1/2 -1/2) % 1/2, 0 , 1/2 % 1/2, 1/2, 0 % % NETWORK TWO - 8-fold connectors % 1/2,1/2,1/2 connected to all corners. % % BRIEF DESCRIPTION OF LEVEL SETS % Level sets are a convenient way to describe smooth surfaces. % Consider a function V(x,y,z). % The equation V(x,y,z)=C defines a 2D surface. % By picking an appropriate function V(x,y,z) and value C, one can define a 2D surface with the desired properties. % % To plot the surface, in Matlab, we set up a grid in 3D space and compute V(x,y,z) on that grid. % Then, using the function isosurface, we find the surfaces on which V(x,y,z)=Constant. % Finally, we plot those. % % For the IWP, the function V(x,y,z) = cos(x).*cos(y)+cos(y).*cos(z)+cos(z).*cos(x) = 0 % defines a surface with zero mean curvature (almost) that is periodic on 0< x,y,z < 2*pi. % Furthermore, the surface has the topology of the IWP. % % V(x,y,z)=1.5 picks out just the network intersections for network one. % V(x,y,z)=-1.5 picks out just the network intersections for network two. % V(x,y,z)=C (-1.5< C< 1.5) picks out a surface that surrounds one of the two networks. % % If we want to plot the ZMC surface we plot V(x,y,z)=0 % If we want to plot a surface enclosing network one we plot V(x,y,z)=C (0<=C<1.5) % If we want to plot a surface enclosing network two we plot V(x,y,z)=-C. % % % GEST April 2004 %Set up a 3-D grid on which to evaluate functions and establish the level surface. Ngrid = 40; % Make a grid in x*y*z that is Ngrid*Ngrid*Ngrid [x, y ,z] = meshgrid(1:Ngrid,1:Ngrid,1:Ngrid); % Define our level set function on the grid Nunit_cells = 1; % Number of unit cells on the grid. If 1, plot just a unit cell. If 2, plot 2*2*2 unit cells. x2 = Nunit_cells* 2*pi*(x-0.5)/Ngrid; y2 = Nunit_cells*2*pi*(y-0.5)/Ngrid; z2 = Nunit_cells*2*pi*(z-0.5)/Ngrid; V = cos(x2).*cos(y2)+cos(y2).*cos(z2)+cos(z2).*cos(x2) - 0.2*(cos(2*x2)+cos(2*y2)+cos(2*z2)); % this version of V works better for our purposes. % Now we just have to pick out surfaces for which V(x,y,z) is a constant. % Figure One - The ZMC surface. figure(1);clf; C = 0.0; handle1=patch(isosurface(x,y,z,V,C)); % Pick out the zero mean curvature surface view(20,10); % Set the viewing angle daspect([1,1,1]) ; % Fix the x,y,z aspect ratio camlight headlight % set the lighting lighting gouraud axis on set(handle1,'FaceColor',[1 0.7 0.3]); % Change the surface colour to gold set(handle1,'EdgeColor','none'); % Remove outlining of the edge of each facet. % Figure Two - Plot a particular pair of interlocking networks. figure(2); clf hold on C1 = -0.8; % Pick out the 4-fold network C2 = 0.3; % Pick out the 8-fold network handle2=patch(isosurface(x,y,z,V,C1)); % Network 1 handle2caps = patch(isocaps(x,y,z,V,C1,'enclose','below')); % The surface is going to be cut at the edges of the unit cell. We want to cap these off with isocaps so the viewer is not thrown off. handle3=patch(isosurface(x,y,z,V,C2)); % Network 2 handle3caps = patch(isocaps(x,y,z,V,C2,'enclose','above')); % Set the viewing options view(20,10); % Angle to view from initially. daspect([1,1,1]); % Aspect Ratio camlight % Lighting lighting gouraud axis off % Make the first network red. set(handle2,'FaceColor',[1 0.0 0.0]); set(handle2,'EdgeColor','none'); set(handle2caps,'FaceColor',[1.0 0.0 0.0]); set(handle2caps,'EdgeColor','none'); % Make the second network blue. set(handle3,'FaceColor',[0.0 0.0 1.0]); set(handle3,'EdgeColor','none'); set(handle3caps,'FaceColor',[0.0 0.0 1.0]); set(handle3caps,'EdgeColor','none'); % Compute the volume enclosed. Volume1 = sum(sum(sum( ceil((C1-V)/4.001)))) / (Ngrid*Ngrid*Ngrid); Volume2 = sum(sum(sum( ceil((V-C2)/4.001)) )) / (Ngrid*Ngrid*Ngrid); fprintf(1,'\n The volume enclosed by the four-fold network in Figure Two is %f percent.', Volume1*100); fprintf(1,'\n The volume enclosed by the eight-fold network in Figure Two is %f percent.',Volume2*100);