VertexsMove.zip_二维不规则多边形晶粒构建；_六边晶粒_晶粒

% 通过移动六边形格子顶点得到不规格晶粒 clc;clear; hnum = 20; % x方向格子数 vnum = 20; % y方向格子数 LatticeVerticalLen = 1.0; % 格子对边距离 GrainMeanRadius = 5.0; %cm GrainRadiusStd = 2; %cm % 晶粒中心点权重 rng(1,'twister'); GrainRadius = random('norm',GrainMeanRadius,GrainRadiusStd,[vnum*hnum,1]); % 晶粒尺寸服从高斯分布 GrainRadius( find(GrainRadius<1) ) = 1; %晶粒尺寸下限 GrainRadius( find(GrainRadius>10) ) = 10; %晶粒尺寸上限 %hist(GrainRadius); GrainW = GrainRadius.^(-1); % 相应权重 %GrainW = ones(size(GrainRadius)); % mean = 1; % 格子平均权重 % std = 0; % 格子权重标准差 % GrainW = random('norm',mean,std,[vnum*hnum,1]); % 生成服从正太分布的权重系数 % GrainW(find(GrainW<0)) = 0.1; % 权重系数须大于0 % 构建顶点 % ---------------------------------------------------------------------- % 中心在原点的六面体格子顶点坐标 const1 = LatticeVerticalLen/sqrt(3.0); const2 = const1/2.0; const3 = LatticeVerticalLen/2.0; p1 = [const1 0]; p2 = [const2 -const3]; p3 = [-const2 -const3]; p4 = -p1; p5 = -p2; p6 = -p3; % -------六面体格子------- % 5 * 6 % * * % 4 o 1 % * * % 3 * 2 % ------------------------ point1 = p3-[0,LatticeVerticalLen/2.0]; point2 = p4-[0,LatticeVerticalLen/2.0]; Points1 = []; for j=1:vnum point1 = point1 + [0,LatticeVerticalLen]; point2 = point2 + [0,LatticeVerticalLen]; Points1 = [Points1;point1;point2]; end Points1 = [Points1;p5+[0,(vnum-0.5)*LatticeVerticalLen]]; Points2(:,1) = Points1(:,1) + (const1+const2); Points2(:,2) = Points1(:,2) + const3; Vertexs=Points1; for i=2:hnum if mod(i,2)==0 Pointsnow = []; Pointsnow = [Points2(:,1)+(i-2)*(const1+const2),Points2(:,2)]; Pointsnow = [[p2(1)+(i-2)*(const1+const2),p2(2)+LatticeVerticalLen/2.0];Pointsnow]; Vertexs = [Vertexs;Pointsnow]; else Pointsnow = []; Pointsnow = [Points1(:,1)+(i-1)*(const1+const2),Points1(:,2)]; Pointsnow = [Pointsnow;[p6(1)+(i-2)*(const1+const2),p6(2)+(vnum)*LatticeVerticalLen]]; Vertexs = [Vertexs;Pointsnow]; end end if mod(hnum,2) == 0 Pointsnow = []; Pointsnow = [Points1(:,1)+hnum*(const1+const2),Points1(:,2)]; Pointsnow = [Pointsnow(2:end,:);[p6(1)+(hnum-1)*(const1+const2),p6(2)+(vnum)*LatticeVerticalLen]]; Vertexs = [Vertexs;Pointsnow]; else Pointsnow = []; Pointsnow = [Points2(:,1)+(hnum-1)*(const1+const2),Points2(:,2)]; Pointsnow = [[p2(1)+(hnum-1)*(const1+const2),p2(2)+LatticeVerticalLen/2.0];Pointsnow(1:end-1,:)]; Vertexs = [Vertexs;Pointsnow]; end Vertexs(:,2) = Vertexs(:,2) - LatticeVerticalLen/2.0; % ------------------------------------------------------------------------- % 格子中心 % --------------------------------------------------------------- Grain=zeros(hnum*vnum,2); for i= 1:hnum for j=1:vnum grainno = (i-1)*vnum + j; if mod(i,2)==0 %偶数 xloc = (i-1)*(const1+const2); yloc = (j-0.5)*LatticeVerticalLen; else xloc = (i-1)*(const1+const2); yloc = (j-1)*LatticeVerticalLen; end Grain(grainno,:)=[xloc,yloc]; end end % ---------------------------------------------------------- % 格子对应的顶点序号 % ----------------------------------------------------- p123456 = [p1;p2;p3;p4;p5;p6]; GrainVertexs=zeros(size(Grain,1),6); for grainno = 1:size(Grain,1) graincenter = Grain(grainno,:); pp(:,1) = p123456(:,1) + graincenter(1); pp(:,2) = p123456(:,2) + graincenter(2); for pnum = 1:6 pointnow = pp(pnum,:); Pointsnow1 = Vertexs(:,1)-pointnow(1); Pointsnow2 = Vertexs(:,2)-pointnow(2); Pointsnow = [Pointsnow1,Pointsnow2]; PointsnowSq = sum(Pointsnow.^2,2); I = find(PointsnowSq<1e-2); if length(I) == 1 GrainVertexs(grainno,pnum) = I; else error(''); end end end % -------------------------------------------------------- % 格子中心位置在格子内随机放置 GrainNow = Grain; for grainno = 1:size(Grain,1) radius = rand*LatticeVerticalLen/2.0; theta = rand*2*pi; dis = [radius*cos(theta),radius*sin(theta)]; GrainNow(grainno,:) = Grain(grainno,:)+dis; end % 各顶点移动至其连接的格子中心点构成的三角形的加权形心位置 % ------------------------------------------------------- % 确定每个顶点所连接的格子 VertexGrains= cell(size(Vertexs,1),1); for i = 1:size(Vertexs,1) I1 = find(GrainVertexs(:,1)==i); I2 = find(GrainVertexs(:,2)==i); I3 = find(GrainVertexs(:,3)==i); I4 = find(GrainVertexs(:,4)==i); I5 = find(GrainVertexs(:,5)==i); I6 = find(GrainVertexs(:,6)==i); I7 = union(I1,I2); I8 = union(I7,I3); I9 = union(I8,I4); I10 = union(I9,I5); I11 = union(I10,I6); VertexGrains{i,:}=I11'; end VertexsNow = Vertexs; for i = 1:size(Vertexs,1) neighborgrainsindex = VertexGrains{i,:}; neighborgrainnum = length(neighborgrainsindex); switch neighborgrainnum case 1 case 2 grainindex1 = neighborgrainsindex(1); grainindex2 = neighborgrainsindex(2); VertexsNow(i,:) = ( GrainNow(grainindex1,:)*GrainW(grainindex1) + ... GrainNow(grainindex2,:)*GrainW(grainindex2) ) ... /(GrainW(grainindex1) + GrainW(grainindex2)); case 3 grainindex1 = neighborgrainsindex(1); grainindex2 = neighborgrainsindex(2); grainindex3 = neighborgrainsindex(3); VertexsNow(i,:) = ( GrainNow(grainindex1,:)*GrainW(grainindex1) + ... GrainNow(grainindex2,:)*GrainW(grainindex2) + ... GrainNow(grainindex3,:)*GrainW(grainindex3) ) /... (GrainW(grainindex1) + GrainW(grainindex2) + GrainW(grainindex3) ); otherwise error(''); end end % ------------------------------------------------------------------------ % 显示初始晶格 % for grainno = 1:size(Grain,1) % grainvertexsindex = GrainVertexs(grainno,:); % poly = Vertexs(grainvertexsindex,:); % p = patch(poly(:,1),poly(:,2),rand(1,3)); % set(p,'EdgeColor','b','FaceColor','none','Marker','.','LineStyle',':'); % hold on; % end % axis equal; % axis off; % 显示初始格子中心 % plot(Grain(:,1),Grain(:,2),'b*'); % hold on; % 显示变形后晶格 hold off; for grainno = 1:size(Grain,1) grainvertexsindex = GrainVertexs(grainno,:); poly = VertexsNow(grainvertexsindex,:); p = patch(poly(:,1),poly(:,2),rand(1,3)); set(p,'EdgeColor','r','FaceColor','none','Marker','.','LineWidth',2); hold on; end %显示当前格子中心 % plot(GrainNow(:,1),GrainNow(:,2),'ro','MarkerFaceColor','r','MarkerSize',6); % hold on; % 显示晶粒半径 % for i=1:size(GrainNow,1) % str = num2str(GrainRadius(i)); % if length(str)>4 % str = str(1:4); % end % text(GrainNow(i,1),GrainNow(i,2),str,'FontSize',6); % end % 显示权重系数 % for i=1:size(GrainNow,1) % str = num2str(GrainW(i)); % if length(str)>4 % str = str(1:4); % end % text(GrainNow(i,1),GrainNow(i,2),str,'FontSize',6); % end axis equal; axis off; print('-djpeg','-r300','MoveVertexs') % 还可以研究通过移动顶点得到的晶粒尺寸分布规律与设想的分布规律是否相符 % 边的创建，切割 % 晶粒尺寸分布规律 GrainAreasNow = zeros(size(Grain,1),1); for grainno = 1:size(Grain,1) grainvertexs = GrainVertexs(grainno,:); grainvertexsc = [grainvertexs,grainvertexs(1)]; vertexs = VertexsNow(grainvertexsc,:); GrainAreasNow(grainno,1) = polyarea(vertexs(:,1),vertexs(:,2)); end GrainAreasNow = GrainAreasNow/(LatticeVerticalLen^2/2.0*sqrt(3))*(pi*GrainMeanRadius^2); % 顶点移动后的晶粒面积 GrainAreas = pi*GrainR