晶界研究程序

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)));
%%
%这里验证一下坐标变换规则