24小时热门版块排行榜

>论坛更新日志 (2975)
>虫友互识 (718)
>文献求助 (156)
>休闲灌水 (133)
>导师招生 (132)
>考博 (123)
>教师之家 (79)
>博后之家 (71)
>硕博家园 (71)
>找工作 (71)
>论文投稿 (68)
>论文道贺祈福 (64)
>基金申请 (61)
>考研 (54)
>公派出国 (49)
>招聘信息布告栏 (31)

返回列表

当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖

爱笑的眼睛na

金虫 (初入文坛)

应助: 0 (幼儿园)
金币: 1118.1
帖子: 17
在线: 80.7小时
虫号: 2027185
注册: 2012-09-24
性别: MM
专业: 煤炭地下开采

[求助] 随机生长法模拟多孔介质的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

回复此楼