%%%%%?A?SOFT-KILL?BESO?CODE?BY?X.?HUANG?and?Y.M.?Xie??%%%%%
function?[x]?=?softbeso（nelxnelyvolfracerrmin）
%?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）;
%?STABLIZATION?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?
%%%%%%%%%%?OPTIMALITY?CRITERIA?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