function?top（nelxnelyvolfracpenalrmin）;
%?INITIALIZE
nelx=60;
nely=20;
penal=3;
rmin=1.25;
volfrac=0.3;
x（1:nely1:nelx）?=?volfrac;?
loop?=?0;?
change?=?1.;
for?ely=1:nely
????for?elx=1:nelx
????????if?ely<3
????????????passive（elyelx）=1;
????????????x（elyelx）=0.001;
????????else
????????????passive（elyelx）=0;
????????end
????end
end
%?START?ITERATION
while?change?>?0.01??
??loop?=?loop?+?1;
??xold?=?x;
%?FE-ANALYSIS
??[U]=FE（nelxnelyxpenal）;?????????
%?objectIVE?FUNCTION?AND?SENSITIVITY?ANALYSIS
??[KE]?=?lk;
??c?=?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?=?c?+?x（elyelx）^penal*Ue‘*KE*Ue;
??????dc（elyelx）?=?-penal*x（elyelx）^（penal-1）*Ue‘*KE*Ue;
????end
??end
%?FILTERING?OF?SENSITIVITIES
??[dc]???=?check（nelxnelyrminxdc）;????
%?DESIGN?UPDATE?BY?THE?OPTIMALITY?CRITERIA?METHOD
??[x]????=?OC（nelxnelyxvolfracdcpassive）;?
%?PRINT?RESULTS
??change?=?max（max（abs（x-xold）））;
??disp（[‘?It.:?‘?sprintf（‘%4i‘loop）?‘?Obj.:?‘?sprintf（‘%10.4f‘c）?...
???????‘?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?[xnew]=OC（nelxnelyxvolfracdcpassive）??
l1?=?0;?l2?=?100000;?move?=?0.2;
while?（l2-l1?>?1e-4）
??lmid?=?0.5*（l2+l1）;
??
??xnew?=?max（0.001max（x-movemin（1.min（x+movex.*sqrt（-dc./lmid）））））;
??xnew（find（passive））=1;
??if?sum（sum（xnew））?-?volfrac*nelx*nely?>?0;
????l1?=?lmid;
??else
????l2?=?lmid;
??end
end
%%%%%%%%%%?MESH-INDEPENDENCY?FILTER?%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function?[dcn]=check（nelxnelyrminxdc）
dcn=zeros（nelynelx）;
for?i?=?1:nelx
??fo