| 查看: 1663 | 回复: 3 | ||
爱笑的眼睛na金虫 (初入文坛)
|
[求助]
随机生长法模拟多孔介质的matlab二维代码改写三维，结果出来还是二维 已有2人参与
|
|
本人根据http://muchong.com/html/201111/3749191.html随机生长法模拟多孔介质的二维代码编写了一个三维代码，但结果出来还是二维的，请求各位大神查一下错误在哪儿，并给出解决的方法，谢谢啦！ 代码如下： %======设定初值============= clc;clear; rand('state',0); set(gcf,'DoubleBuffer','on');%消除振动 Cd=0.01;%生长核分布概率数 Vr=0.01;%生长相的体积分数 global kn kn=500; p=[1/8 1/8 1/8 1/8 1/4 1/4 1/4 1/4];%方向增长概率数 A=rand(kn,kn);%初始网格kn×kn A(find(A<=Cd))=0;%黑色，生长核 A(find(A>Cd))=1;%白色 m=0;%迭代步骤 N=kn*kn*kn;%格点数 [x,y,z]=find(A==0);%生长核脚标 [y1,z] = ind2sub([kn,kn],y); Ixy1z = sub2ind([kn,kn,kn],x,y1,z); gn=length(Ixy1z);%生长核计数 Vs=gn/N %体积分数 %=========================== while Vs<0.001 xW=x-1;xE=x+1; xW(xW==0)=1;% 对近邻进行边界条件处理 xE(xE==kn+1)=kn;% 对近邻进行边界条件处理 y1N=y1-1;y1S=y1+1; y1N(y1N==0)=1;% 对近邻进行边界条件处理 y1S(y1S==kn+1)=kn;% 对近邻进行边界条件处理 zB=z-1;zT=z+1; zB(zB==0)=1;% 对近邻进行边界条件处理 zT(zT==kn+1)=kn;% 对近邻进行边界条件处理 InE=sub2ind([kn,kn],xE,z); % 右方转化脚标为索引 InB=sub2ind([kn,kn],x,zB); % 下方转化脚标为索引 InW=sub2ind([kn,kn],xW,z); % 左方转化脚标为索引 InT=sub2ind([kn,kn],x,zT); % 上方转化脚标为索引 InS=sub2ind([kn,kn],x,y1N); % 下方转化脚标为索引 InN=sub2ind([kn,kn],x,y1S); % 上方转化脚标为索引 InET=sub2ind([kn,kn],xE,zT); % 右上方转化脚标为索引 InEB=sub2ind([kn,kn],xE,zB); % 右下方转化脚标为索引 InWB=sub2ind([kn,kn],xW,zB); % 左下方转化脚标为索引 InWT=sub2ind([kn,kn],xW,zT); % 左上方转化脚标为索引 InEN=sub2ind([kn,kn],xE,y1N); % 右上方转化脚标为索引 InES=sub2ind([kn,kn],xE,y1S); % 右下方转化脚标为索引 InWS=sub2ind([kn,kn],xW,y1S); % 左下方转化脚标为索引 InWN=sub2ind([kn,kn],xW,y1N); % 左上方转化脚标为索引 InTS=sub2ind([kn,kn],y1S,zT); % 右上方转化脚标为索引 InBS=sub2ind([kn,kn],y1S,zB); % 右下方转化脚标为索引 InBN=sub2ind([kn,kn],y1N,zB); % 左下方转化脚标为索引 InTN=sub2ind([kn,kn],y1N,zT); % 左上方转化脚标为索引 %==========做图======== CInE=setdiff(InE,Ixy1z); En=rand(1,length(CInE));[i,j,k]=find(En<=p(5));pE=sub2ind([1,length(CInE)],i,j,k); A(CInE(pE))=0; [x,y1,z]=find(A==0); Ixy1z=sub2ind([kn,kn,kn],x,y1,z); CInT= setdiff(InT,Ixy1z);Tn=rand(1,length(CInT));[i,j,k]=find(Tn<=p(6));pT=sub2ind([1,length(CInT)],i,j,k); A(CInT(pT))=0; [x,y1,z]=find(A==0); Ixy1z=sub2ind([kn,kn,kn],x,y1,z); CInW=setdiff(InW,Ixy1z);Wn=rand(1,length(CInW));[i,j,k]=find(Wn<=p(7));pW=sub2ind([1,length(CInW)],i,j,k); A(CInW(pW))=0; [x,y1,z]=find(A==0); Ixy1z=sub2ind([kn,kn,kn],x,y1,z); CInB=setdiff(InB,Ixy1z);Bn=rand(1,length(CInB));[i,j,k]=find(Bn<=p(8));pB=sub2ind([1,length(CInB)],i,j,k); A(CInB(pB))=0; [x,y1,z]=find(A==0); Ixy1z=sub2ind([kn,kn,kn],x,y1,z); CInEB=setdiff(InEB,Ixy1z);EBn=rand(1,length(CInEB));[i,j,k]=find(EBn<=p(8));pEB=sub2ind([1,length(CInEB)],i,j,k); A(CInEB(pEB))=0; [x,y1,z]=find(A==0); Ixy1z=sub2ind([kn,kn,kn],x,y1,z); CInET=setdiff(InET,Ixy1z);ETn=rand(1,length(CInET));[i,j,k]=find(ETn<=p(8));pET=sub2ind([1,length(CInET)],i,j,k); A(CInET(pET))=0; [x,y1,z]=find(A==0); Ixy1z=sub2ind([kn,kn,kn],x,y1,z); CInWT=setdiff(InWT,Ixy1z);WTn=rand(1,length(CInWT));[i,j,k]=find(WTn<=p(8));pWT=sub2ind([1,length(CInWT)],i,j,k); A(CInWT(pWT))=0; [x,y1,z]=find(A==0); Ixy1z=sub2ind([kn,kn,kn],x,y1,z); CInWB=setdiff(InWB,Ixy1z);WBn=rand(1,length(CInWB));[i,j,k]=find(WBn<=p(8));pWB=sub2ind([1,length(CInWB)],i,j,k); A(CInWB(pWB))=0; [x,y1,z]=find(A==0); Ixy1z=sub2ind([kn,kn,kn],x,y1,z); CInS= setdiff(InS,Ixy1z);Sn=rand(1,length(CInS));[i,j,k]=find(Sn<=p(5));pS=sub2ind([1,length(CInS)],i,j,k); A(CInS(pS))=0; [x,y1,z]=find(A==0); Ixy1z=sub2ind([kn,kn,kn],x,y1,z); CInN=setdiff(InN,Ixy1z);Nn=rand(1,length(CInN));[i,j,k]=find(Nn<=p(8));pN=sub2ind([1,length(CInN)],i,j,k); A(CInN(pN))=0; [x,y1,z]=find(A==0); Ixy1z=sub2ind([kn,kn,kn],x,y1,z); CInEN=setdiff(InEN,Ixy1z);ENn=rand(1,length(CInEN));[i,j,k]=find(ENn<=p(8));pEN=sub2ind([1,length(CInEN)],i,j,k); A(CInEN(pEN))=0; [x,y1,z]=find(A==0); Ixy1z=sub2ind([kn,kn,kn],x,y1,z); CInES=setdiff(InES,Ixy1z);ESn=rand(1,length(CInES));[i,j,k]=find(ESn<=p(8));pES=sub2ind([1,length(CInES)],i,j,k); A(CInES(pES))=0; [x,y1,z]=find(A==0); Ixy1z=sub2ind([kn,kn,kn],x,y1,z); CInWS=setdiff(InWS,Ixy1z);WSn=rand(1,length(CInWS));[i,j,k]=find(WSn<=p(8));pWS=sub2ind([1,length(CInWS)],i,j,k); A(CInWS(pWS))=0; [x,y1,z]=find(A==0); Ixy1z=sub2ind([kn,kn,kn],x,y1,z); CInWN=setdiff(InWN,Ixy1z);WNn=rand(1,length(CInWN));[i,j,k]=find(WNn<=p(8));pWN=sub2ind([1,length(CInWN)],i,j,k); A(CInWN(pWN))=0; [x,y1,z]=find(A==0); Ixy1z=sub2ind([kn,kn,kn],x,y1,z); CInTS=setdiff(InTS,Ixy1z); TSn=rand(1,length(CInTS));[i,j,k]=find(TSn<=p(8));pTS=sub2ind([1,length(CInTS)],i,j,k); A(CInTS (pTS))=0; [x,y1,z]=find(A==0); Ixy1z=sub2ind([kn,kn,kn],x,y1,z); CInBS=setdiff(InBS,Ixy1z);BSn=rand(1,length(CInBS));[i,j,k]=find(BSn<=p(8));pBS=sub2ind([1,length(CInBS)],i,j,k); A(CInBS(pBS))=0; [x,y1,z]=find(A==0); Ixy1z=sub2ind([kn,kn,kn],x,y1,z); CInBN=setdiff(InBN,Ixy1z); BNn=rand(1,length(CInBN));[i,j,k]=find(BNn<=p(8));pBN=sub2ind([1,length(CInBN)],i, j,k); A(CInBN (pBN))=0; [x,y1,z]=find(A==0); Ixy1z=sub2ind([kn,kn,kn],x,y1,z); CInTN=setdiff(InTN,Ixy1z); TNn=rand(1,length(CInTN));[i,j,k]=find(TNn<=p(8));pTN=sub2ind([1,length(CInTN)],i, j,k); A(CInTN (pTN))=0; [x,y1,z]=find(A==0); Ixy1z=sub2ind([kn,kn,kn],x,y1,z); [x,y1,z]=find(A==0); Ixy1z=sub2ind([kn,kn,kn],x,y1,z); gn=length(Ixy1z); Vs=gn/N %体积分数 m=m+1 pause(0.1) imshow(A,'InitialMagnification','fit')%更新图像 title('time=','Fontsize',14,'Fontname','Times new roman'); % 更新时间参数 end |
回复此楼» 猜你喜欢
售SCI一区T0P文章，我:8.O.5.5.1.O.5.4,科目齐全,可+急
已经有3人回复
-
评审有感
已经有27人回复
-
论文撤稿了
已经有7人回复
-
售SCI一区T0P文章，我:8.O.5.5.1.O.5.4,科目齐全,可+急
已经有4人回复
-
售SCI一区T0P文章，我:8.O.5.5.1.O.5.4,科目齐全,可+急
已经有3人回复
-
售SCI一区T0P文章，我:8.O.5.5.1.O.5.4,科目齐全,可+急
已经有3人回复
-
售SCI一区T0P文章，我:8.O.5.5.1.O.5.4,科目齐全,可+急
已经有5人回复
-
售SCI一区T0P文章，我:8.O.5.5.1.O.5.4,科目齐全,可+急
已经有4人回复
-
河北省自然科学基金
已经有9人回复
-
售SCI一区T0P文章，我:8.O.5.5.1.O.5.4,科目齐全,可+急
已经有8人回复
2楼2016-06-22 19:51:30
3楼2016-06-22 19:53:48
4楼2018-08-06 17:03:37
5