matlab获得蜂窝六边形程序
clear
clc
%%
%下面得到六边形的蜂窝形状。
%这里需要输出一个文件,data=zeros(n,8),其数据格式如下；
%n代表晶界的数量。
%第一到二列：表示晶界第一个端点的横纵坐标，即初始点。
%第三到四列：表示晶界第二个端点的横纵坐标，即被扩展点。
%第五到六列：表示晶界的平移向量。
%第七列：表示晶界的方向。
%第八到九列：表示晶界的旋转正余弦值。
%以三角形的方式向外扩展。
%输入六边形边长；
num=7;
a=10^(-8);
center=[0,0];
point={};
%首先初始化三角形的三条边坐标；
point{1}(1,:)=[a,0,-a/2,a*3^(0.5)/2];
point{1}(2,:)=[-a/2,a*3^(0.5)/2,-a/2,-a*3^(0.5)/2];
point{1}(3,:)=[-a/2,-a*3^(0.5)/2,a,0];
dir1=zeros(3,2);
for i=1:3
dir1(i,:)=((point{1}(i,1:2)-point{1}(i,3:4)))/norm(point{1}(i,1:2)-point{1}(i,3:4));
end
dir2=zeros(3,2);
dir2(1,:)=([a/2,a*3^(0.5)/2]-[-a/2,-a*3^(0.5)/2])/norm([a/2,a*3^(0.5)/2]-[-a/2,-a*3^(0.5)/2]);
dir2(2,:)=[1,0];
dir2(3,:)=([a/2,-a*3^(0.5)/2]-[-a/2,a*3^(0.5)/2])/norm([a/2,-a*3^(0.5)/2]-[-a/2,a*3^(0.5)/2]);
compare=1;
%下面开始扩展,num表示扩展的次数
while(compare<num)
[m,~]=size(point{compare});
point_sum=0;
leng=3*a/2;
for i=1:m
direc=(point{compare}(i,1:2)-point{compare}(i,3:4))/norm((point{compare}(i,1:2)-point{compare}(i,3:4)));
center=(point{compare}(i,1:2)+point{compare}(i,3:4))/2;
%%
judge=zeros(1,3);
for t=1:3
temp1=abs(sum(direc-dir1(t,:)));
temp2=abs(sum(direc+dir1(t,:)));
if(temp1<=1*10^(-15)||temp2<=1*10^(-15))
judge(t)=1;
end
end
if(judge(1)==0&&judge(2)==0&&judge(3)==0)
tt=[direc];
end
%%
if(judge(1))
point_direct=dir2(1,:);
tem_point1=leng*point_direct+center;
tem_point2=-leng*point_direct+center;
%将这两个点的坐标存入point{compare}之中。
point_sum=point_sum+1;
point{compare+1}(point_sum,1:2)=tem_point1;
point{compare+1}(point_sum,3:4)=point{compare}(i,1:2);
point_sum=point_sum+1;
point{compare+1}(point_sum,1:2)=tem_point1;
point{compare+1}(point_sum,3:4)=point{compare}(i,3:4);
point_sum=point_sum+1;
point{compare+1}(point_sum,1:2)=tem_point2;
point{compare+1}(point_sum,3:4)=point{compare}(i,1:2);
point_sum=point_sum+1;
point{compare+1}(point_sum,1:2)=tem_point2;
point{compare+1}(point_sum,3:4)=point{compare}(i,3:4);
end
if(judge(2))
point_direct=dir2(2,:);
tem_point1=leng*point_direct+center;
tem_point2=-leng*point_direct+center;
%将这两个点的坐标存入point{compare}之中。
point_sum=point_sum+1;
point{compare+1}(point_sum,1:2)=tem_point1;
point{compare+1}(point_sum,3:4)=point{compare}(i,1:2);
point_sum=point_sum+1;
point{compare+1}(point_sum,1:2)=tem_point1;
point{compare+1}(point_sum,3:4)=point{compare}(i,3:4);
point_sum=point_sum+1;
point{compare+1}(point_sum,1:2)=tem_point2;
point{compare+1}(point_sum,3:4)=point{compare}(i,1:2);
point_sum=point_sum+1;
point{compare+1}(point_sum,1:2)=tem_point2;
point{compare+1}(point_sum,3:4)=point{compare}(i,3:4);
end
if(judge(3))
point_direct=dir2(3,:);
tem_point1=leng*point_direct+center;
tem_point2=-leng*point_direct+center;
%将这两个点的坐标存入point{compare}之中。
point_sum=point_sum+1;
point{compare+1}(point_sum,1:2)=tem_point1;
point{compare+1}(point_sum,3:4)=point{compare}(i,1:2);
point_sum=point_sum+1;
point{compare+1}(point_sum,1:2)=tem_point1;
point{compare+1}(point_sum,3:4)=point{compare}(i,3:4);
point_sum=point_sum+1;
point{compare+1}(point_sum,1:2)=tem_point2;
point{compare+1}(point_sum,3:4)=point{compare}(i,1:2);
point_sum=point_sum+1;
point{compare+1}(point_sum,1:2)=tem_point2;
point{compare+1}(point_sum,3:4)=point{compare}(i,3:4);
end
end
compare=compare+1;
end
%上面获得了所需点的坐标。
%下面整合该点,并且选出不重复的坐标点
%%
tem1=[];
for i=1:compare
tem1=[tem1;point{i}(:,1:2);point{i}(:,3:4)];
end
tem1=unique(tem1,"rows");
tt=tem1;
[m,~]=size(tem1);
for i=1:m-1
for j=(i+1):m
distance=norm(tem1(i,:)-tem1(j,:));
if(distance<10^(-13))
tem1(j,:)=[10^20,10^20];
end
end
end
pp=0;
points_end=[];
for i=1:m
if(tem1(i,:)~=[10^20,10^20])
pp=pp+1;
points_end(pp,:)=tem1(i,:);
end
end
[m,n]=size(points_end);
for i=1:m
for j=1:n
if(abs(points_end(i,j)-0)<10^(-14))
points_end(i,j)=0;
end
end
end
plot(points_end(:,1),points_end(:,2),"*");
%%
%下面开始将points_end这些点扩展为边；
data=zeros(3*m,9);
%按照顺时针的方向扩展，角度依次为0，120，240度，使用1，2，3来区分。
dir1=[1,0];
dir2=([-a/2,a*3^(0.5)/2])/norm([-a/2,a*3^(0.5)/2]);
dir3=([-a/2,-a*3^(0.5)/2])/norm([-a/2,-a*3^(0.5)/2]);
dir=[dir1;dir2;dir3];
for i=1:m
for j=1:3
data(3*(i-1)+j,7)=j;
data(3*(i-1)+j,1:2)=points_end(i,:);
%接下来扩展下一个点。
data(3*(i-1)+j,3:4)=a*dir(j,:)+points_end(i,:);
end
end
tem1=[data(:,1:2);data(:,3:4)];
plot(tem1(:,1),tem1(:,2),"*");
%%
%下面求得坐标变换规则。
center1=([0,0]+[a,0])/2;
center2=([0,0]+[-a/2,a*3^(0.5)/2])/2;
center3=([0,0]+[-a/2,-a*3^(0.5)/2])/2;
center_point=[center1;center2;center3];
T=[1,0;-1/2,3^(0.5)/2;-1/2,-3^(0.5)/2];%旋转矩阵的正弦和余弦值。
for i=1:m
for j=1:3
center=(data(3*(i-1)+j,1:2)+data(3*(i-1)+j,3:4))/2;
data(3*(i-1)+j,5:6)=center-center_point(j,:);
data(3*(i-1)+j,8:9)=T(j,:);
end
end
for i=1:3*m
for j=1:4
if(abs(data(i,j))<10^(-15))
data(i,j)=0;
end
end
end
%下面我们将其写入到txt文件中
fid=fopen("grain_boundary_message_data.txt","w");
for i=1:3*m
for j=1:9
fprintf(fid,"%g",data(i,j));
fprintf(fid,"%s"," ");
end
fprintf(fid,"%s\n"," ");
end
%%
for i=1:3*m
plot([data(i,1),data(i,3)],[data(i,2),data(i,4)]);
hold on
end
a1=max(max(data(:,1:4)));
a2=min(min(data(:,1:4)));
%%
%这里验证一下坐标变换规则
|