% This script plots up a zero-mean-curvature surface with the double-diamond 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 DOUBLE DIAMOND % The DD 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 Pn3m. % 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 DD but it's a good start. % % The DD divides space into two separate, conintuous networks. % Each network consists of a 4-fold connectors, just as in diamond. % NETWORK ONE % Standard Diamond (FCC lattice with units at % 0,0,0 and 1/4, 1/4, 1/4 % % NETWORK TWO % Second Diamond (FCC Lattice at) % 1/2,1/2,1/2 and 3/4,3/4,3/4 % % 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 DD, the function V(x,y,z) = sin(x2).*sin(y2).*sin(z2)+sin(x2).*cos(y2).*cos(z2)+cos(x2).*sin(y2).*cos(z2)+cos(x2).*cos(y2).*sin(z2) = 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 DD. % % V(x,y,z)=1.414 picks out just the network intersections for network one. % V(x,y,z)=-1.414 picks out just the network intersections for network two. % V(x,y,z)=C (-1.414< C< 1.414) 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.414) % 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 = 20; % 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 +pi/2; y2 = Nunit_cells*2*pi*(y-0.5)/Ngrid+pi/2; z2 = Nunit_cells*2*pi*(z-0.5)/Ngrid+pi/2; V = sin(x2).*sin(y2).*sin(z2)+sin(x2).*cos(y2).*cos(z2)+cos(x2).*sin(y2).*cos(z2)+cos(x2).*cos(y2).*sin(z2); % 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; handle1=patch(isosurface(x,y,z,V,C)); % Pick out the zero mean curvature surface view(30,30); % Set the viewing angle daspect([1,1,1]) ; % Fix the x,y,z aspect ratio camlight headlight % set the lighting lighting gouraud axis off 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 C = 0.7; % Pick out a surface close to the networks handle2=patch(isosurface(x,y,z,V,C)); % Network 1 handle2caps = patch(isocaps(x,y,z,V,C,'enclose','above')); % 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,-C)); % Network 2 handle3caps = patch(isocaps(x,y,z,V,-C,'enclose','below')); % Set the viewing options view(30,30); % 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. Volume = sum(sum(sum( ceil((V-C)/3.001)))) / (Ngrid*Ngrid*Ngrid); fprintf(1,'\n The volume enclosed by each network in Figure Two is %f percent.', Volume*100); % Plot on the Unit cell axes plot3(Ngrid/Nunit_cells*[0 1], [0 0], [0 0 ],'g'); plot3([0 0], Ngrid/Nunit_cells*[0 1], [0 0], 'g'); plot3([0 0], [0 0], Ngrid/Nunit_cells*[0 1], 'g');