%%%%?A?SOFT-KILL?BESO?CODE?BY?X.?HUANG?and?Y.M.?XIE?%%%%
%%%%?A?further?review?of?ESO?type?methods?for?topology?optimization?%%%%%
%%%%?SBESO?方法?體積約束下的最小柔度問題%%%%%%%%%%%%
function?sbeso（nelxnelyvolfracerrmin）
nelx=80;
nely=40;
volfrac=0.4;
penal?=?3;??
rmin=2;
er=0.02;
%INITIALIZE
x（1:nely1:nelx）=1.;vol=1.;i=0;change=1.;penal=3.;
%START?iTH?ITERATION
while?change>0.001
????i=i+1;?
????vol=?max（vol*（1-er）volfrac）;
????if?i>1;olddc=dc;
????end
?%FE-ANALYSIS
?[U]=FE（nelxnelyxpenal）;
?%objectIVE?FUNCTION?AND?SENSITIVITY?ANALYSIS
?[KE]=lk;
?c（i）=0.;
?for?ely=1:nely
?????for?elx=1:nelx
?????????n1=（nely+1）*（elx-1）+ely;
?????????n2=（nely+1）*elx+????ely;
?????????Ue=U（[2*n1-1;2*n1;2*n2-1;2*n2;2*n2+1;2*n2+2;2*n1+1;2*n1+2]1）;
?????????c（i）=c（i）+0.5*x（elyelx）^penal*Ue‘*KE*Ue;
?????????dc（elyelx）=0.5*x（elyelx）^（penal-1）*Ue‘*KE*Ue;
?????end
?end
?%FILTERING?OF?SENSITIVITIES
?[dc]=check（nelxnelyrminxdc）;
?%STABILIZATION?OF?EVOLUTIONARY?PROCESS
?if?i>1;
?????dc=（dc+olddc）/2.;
?end
?%BESO?DESIGN?UPDATE
?[x]=ADDDEL（nelxnelyvoldcx）;
?%PRINT?RESULTS
?if?i>10;
?????change=abs（sum（c（i-9:i-5））-sum（c（i-4:i）））/sum（c（i-4:i））;
?end
?disp（[‘?It.:?‘?sprintf（‘%4i‘i）?‘?Obj.:?‘?sprintf（‘%10.4f‘c（i））?...
???????‘?Vol.:?‘?sprintf（‘%6.3f‘sum（sum（x））/（nelx*nely））?...
????????‘?ch.:?‘?sprintf（‘%6.3f‘change?）]）
?%PLOT?DENSITIES
???colormap（gray）;imagesc（-x）;axis?equal;axis?tight;axis?off;pause（1e-6）;
end

%%%%%%%%%%?OPTIMIZATION?CRITERA?UPDATE?%%%%%%%%%%%%%%%%%%%%%%%%%%
function?[x]=ADDDEL（nelxnelyvolfradcx）;
l1=min（min（dc））;l2=max（max（dc））;
while?（（l2-l1）/l2>1.0e-5）
????th=（l1+l2）/2.0;
????x=max（0.001sign（dc-th））;
????if?sum（sum（x））-volfra*（nelx*nely）>0;
????????l1=th;
????else?
????????l2=th;
????end
end

%%%%%%%%%%?MESH-INDEPENDENCY?FILTER?%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function?[dcf]=check（nelxnelyrminxdc）
dcf=zeros（nelynelx）;
for?i=1:nelx
????for?j=1:nely
????????sum=0.0;
????????for?k=max（i-floor（rmin）1）:min（i+floor（rmin）nelx）
????????????for?l=max（j-floor（rmin）1）:min（j+floor（rmin）nely）
????????????????fac=rmin-sqrt（（i-k）^2+（j-1）^2）;
????????????????sum=sum+max（0fac）;
????????????????dcf（ji）=dcf（ji）+max（0fac）*dc（lk）;
????????????end
????????end
????????dcf（ji）=dcf（ji）/sum;
????end
end

%%%%%%%%%%?FE-ANALYSIS?%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function?[U]=FE（nelxnelyxpenal）
[KE]=lk;
K=sparse（2*（nelx+1）*（nely+1）2*（nelx+1）*（nely+1））;
F=sparse（2*（nely+1）*（nelx+1）1）;
U=zeros（2*（nely+1）*（nelx+1）1）;
for?elx=1:nelx
????for?ely=1:nely
????????n1=（nely+1）*（elx-1）+ely;
????????n2=（nely+1）*elx????+ely;
????????edof=[2*n1-1;2*n1;2*n2-1;2*n2;2*n2+1;2*n2+2;2*n1+1;2*n1+2];
????????K（edofedof）=K（edofedof）+x（elyelx）^penal*KE;
????end
end
%?DEFINE?LOADS?AND?SUPPORTS?（Cantilever）
ny=ceil（nely/2）;
????????F（2*（（nely+1）*nelx+（ny+1））1）?=?-1000;
fixeddofs=[1:2*（nely+1）];
alldofs??=[1:2*（nely+1）*（nelx+1）];
freedofs?=setdiff（alldofsfixeddofs）;
%?SOLVIN

?屬性????????????大小?????日期????時間???名稱
-----------?---------??----------?-----??----
?????目錄???????????0??2014-03-27?07:35??基于軟殺死的BESO拓撲優化\
?????文件????????3837??2014-03-26?20:31??基于軟殺死的BESO拓撲優化\sbeso.m