基础教程Matlab模拟六边形曲面细分

R=1; % expected hexadron circumcircle radius a0=sqrt(3)*R; b0=3*R; % in rectangle a0 x b0 there are 2 points and 2 hexadrons expected % rectangle a0 x b0 is minimal possible periodic rectangle in honeycombs % structure na=3; nb=2; % use biger periodic area a=na*a0; b=nb*b0; % % in rectangle a x b there are na*nb*2 points and na*nb*2 hexadrons expected N=2*(na*nb); % , number of points in one of 3*3 rectangles N_ex=N*(3*3); % number of points in 3*3 rectangles itm=400; % number of iteration trm=50; % number of tries cr=1-0.015; % contraction ration mde=sqrt(a*b/N); % estimation for mean distance between points rs=mde; % start search zone size r=[a*rand(1,N); b*rand(1,N)]; % start from some random distribution r_ex=zeros(2,N_ex); % extend to 3 x 3 rectangles periodically: lc=1; for xc=-1:1 for yc=-1:1 i1=1+(lc-1)*N; i2=lc*N; r00=[xc*a; yc*b]; r_ex(1,i1:i2)=r(1,:)+r00(1); r_ex(2,i1:i2)=r(2,:)+r00(2); lc=lc+1; end end figure('units','normalized','position',[0.3 0.3 0.26 0.4]); axes; mxab=max([a b]); [vx,vy]=voronoi(r_ex(1,:),r_ex(2,:)); plot(vx,vy,'b-'); hold on; plot(r_ex(1,:),r_ex(2,:),'.k'); axis equal; xlim([-0.1*a 1.1*a]); ylim([-0.1*b 1.1*b]); drawnow; rr_min=r; Ls_ex_min=1e10; vx_min=vx; vy_min=vy; clear F; fc=1; F1=getframe(gcf); [im,map] = rgb2ind(F1.cdata,256,'nodither'); im(1,1,1,20) = 0; for it=1:itm % iteration count for tr=1:trm % random tireis counter rr=r+rs*(rand(2,N)-0.5); % current random position rr(1,:)=mod(rr(1,:),a); % return to [0; a] range rr(2,:)=mod(rr(2,:),b); % return to [0; b] range % extend to 3 x 3 rectangles periodically: lc=1; for xc=-1:1 for yc=-1:1 i1=1+(lc-1)*N; i2=lc*N; r00=[xc*a; yc*b]; r_ex(1,i1:i2)=rr(1,:)+r00(1); r_ex(2,i1:i2)=rr(2,:)+r00(2); lc=lc+1; end end % minimal distance: mds=1e10; for n1=1:N_ex-1 for n2=n1+1:N_ex dr=r_ex(:,n1)-r_ex(:,n2); dr2=sum(dr.^2); mdst=sqrt(dr2); if mdst<mds mds=mdst; end end end [vx,vy]=voronoi(r_ex(1,:),r_ex(2,:)); Lvx=length(vx); Ls=0; % sum of all edges length inside a x b rectangle for vc=1:Lvx if (vx(1,vc)==vx(2,vc))&&(vy(1,vc)==vy(2,vc)) % it is a point, not segment % do nothing else rseg=[vx(1,vc), vx(2,vc); vy(1,vc), vy(2,vc)]; L=segment_rectange(a,b,rseg); Ls=Ls+L; end end Ls_ex=Ls+0.2*a*b/mds; if Ls_ex<Ls_ex_min Ls_ex_min=Ls_ex; rr_min=rr; vx_min=vx; vy_min=vy; end end r=rr_min; % next iteration starts from curent minimum found rs=rs*cr; % contract search region cla; title(['iteration ' num2str(it) '; rs=' num2str(rs,'%1.5f')]); plot(vx_min,vy_min,'b-'); hold on; lc=1; for xc=-1:1 for yc=-1:1 i1=1+(lc-1)*N; i2=lc*N; r00=[xc*a; yc*b]; r_ex(1,i1:i2)=r(1,:)+r00(1); r_ex(2,i1:i2)=r(2,:)+r00(2); lc=lc+1; end end plot(r_ex(1,:),r_ex(2,:),'.k'); axis equal; xlim([-0.1*a 1.1*a]); ylim([-0.1*b 1.1*b]); drawnow; if mod(it,10)==0 F1=getframe(gcf); im(:,:,1,fc) = rgb2ind(F1.cdata,map,'nodither'); fc=fc+1; end end clc; disp('movie making finished'); imwrite(im,map,'hon_cb.gif','DelayTime',0,'LoopCount',inf);