%%%%?A?99?LINE?TOPOLOGY?OPTIMIZATION?CODE?BY?OLE?SIGMUND?OCTOBER?1999?%%%?
function?top（nelxnelyvolfracpenalrmin）;%nelx和nely分別是水平和垂直方向上的元素?cái)?shù)量volfrac容積率%penal懲罰因子一般為3.0，rmin過(guò)濾半徑
nelx=60;nely=30;volfrac=0.5;penal=3;rmin=1.2;
%?INITIALIZE?
x（1:nely1:nelx）?=?volfrac;%均勻分布材料
loop?=?0;?
change?=?1.;?
%?START?ITERATION?
while?change?>?0.01?
loop?=?loop?+?1;?
xold?=?x;?
%?FE-ANALYSIS?
[U]=FE（nelxnelyxpenal）;%有限元子程序?
%?objectIVE?FUNCTION?AND?SENSITIVITY?ANALYSIS?目的函數(shù)與靈敏度分析
[KE]?=?lk;%單元?jiǎng)偠染仃囎映绦蛑徽{(diào)用一次，因?yàn)閷?duì)固體來(lái)說(shuō)所有元素的單元?jiǎng)偠染仃囅嗤?br/>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;?%全局柔度最小，剛度最大，目標(biāo)函數(shù)
?????dc（elyelx）?=?-penal*x（elyelx）^（penal-1）*?Ue‘*KE*Ue;%目標(biāo)函數(shù)的靈敏度
???end?
end?
%?FILTERING?OF?SENSITIVITIES?網(wǎng)格獨(dú)立過(guò)濾器
[dc]?=?check（nelxnelyrminxdc）;
%?DESIGN?UPDATE?BY?THE?OPTIMALITY?CRITERIA?METHOD?最優(yōu)標(biāo)準(zhǔn)優(yōu)化器
[x]?=?OC（nelxnelyxvolfracdc）;?%%定義優(yōu)化準(zhǔn)則函數(shù)
%?PRINT?RESULTS?
%?A=sum（sum（x））
change?=?max（max（abs（x-xold）））;%如果設(shè)計(jì)變量的變化（第30行中確定的變化）小于1％，則終止主循環(huán)，abs表示取絕對(duì)值
disp（[‘It.:‘?sprintf（‘%4i‘loop）???‘Obj.:‘??sprintf（‘%10.4f‘c）?...?
‘?Vol.:?‘?sprintf（‘%6.3f‘sum（sum（x））/?（nelx*nely））?...?
‘?ch.:?‘?sprintf（‘%6.3f‘change?）]）?%sum（x）表示[x]的列總和sum（sum（x））/?（nelx*nely）表示體積分?jǐn)?shù)，It是迭代次數(shù)，c是全局剛度，ch是設(shè)計(jì)變量的變化
%?PLOT?DENSITIES最后結(jié)果的密度分布圖
colormap（gray）;?imagesc（-x）;?axis?equal;?axis?tight;?axis?off;pause（1e-6）;?
end?
%%%%%%%%%%?OPTIMALITY?CRITERIA?UPDATE?%%%%%%%%%?基于優(yōu)化程序的最優(yōu)化準(zhǔn)則oc法
function?[xnew]=OC（nelxnelyxvolfracdc）?
l1?=?0;?l2?=?100000;?move?=?0.2;?
while?（l2-l1?>?1e-4）%?bi-sectioning算法求出滿(mǎn)足體積約束的拉格朗日乘數(shù)數(shù)值（41-49）小于l1，大于l2的這個(gè)范圍，限定拉格朗日乘數(shù)的間隔大小會(huì)一直減半直到他的大小小于收斂判別準(zhǔn)則
lmid?=?0.5*（l2+l1）;?
%?A=x+move
xnew?=?max（0.001??max（x-move??min（1.?min（x+movex.*sqrt（-dc./lmid）））））;?%開(kāi)根號(hào)?1/2次方代表數(shù)值阻尼系數(shù)η?x.*sqrt（-dc./lmid）為論文中的XeBeη
if?sum（sum（xnew））?-?volfrac*nelx*nely?>?0?
l1?=?lmid;?
else
l2?=?lmid;?
end?
end?
%%%%%%%%%%?MESH-INDEPENDENCY?FILTER?%%%%%%%%%%%?網(wǎng)格獨(dú)立性過(guò)濾技術(shù)
function?[dcn]=check（nelxnelyrminxdc）?
dcn=zeros（nelynelx）;
for?i?=?1:nelx?
????for?j?=?1:nely?
?????sum=