%%%%?A?99?LINE?TOPOLOGY?OPTIMIZATION?CODE?BY?OLE?SIGMUND?JANUARY?2000?%%%
	%%%%?CODE?MODIFIED?FOR?INCREASED?SPEED?September?2002?BY?OLE?SIGMUND?%%%
	%%%%?一個由?OLE?SIGMUND編寫的99行拓撲優化代碼，2000年1月?%%%
	%%%%?為加速而修改的代碼，2002年9月由OLE?SIGMUND編寫?%%
	function?topsun（nelxnelyvolfracpenalrmin）?
	%?INITIALIZE
	x（1:nely1:nelx）?=?volfrac;??
%?????for?ely=1:nely
%????????for?elx=1:nelx
%?????????????if?sqrt（（ely-nely）^2+（elx）^2）>（nely+1）
%??????????????????passive（elyelx）=1;
%??????????????????x（elyelx）=0.001;
%?????????????else
%??????????????????passive（elyelx）=0;
%?????????????end
%????????end
%?????end
	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
	[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（[n1;n2;n2+1;n1+1]1）;
%?????Ue1?=?U（[n1;n2;n2+1;n1+1]1）;
%?????Ue2?=?U（[n1;n2;n2+1;n1+1]2）;
	c?=?c?+（0.001+0.999*x（elyelx）^penal）*Ue‘*KE*Ue;
	dc（elyelx）?=?-0.999*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（nelxnelyxvolfracdc）;?
	%?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（nelxnelyxvolfracdc）?
	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））=0.001;
	if?sum（sum（xnew））?-?volfrac*nelx*nely?>?0;
	l1?=?lmid;
	else
	l2?=?lmid;
	end
	end
	%%%%%%%%%%?MESH-INDEPENDENCY?FILTER?%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	function?[dcn]=check（nelxnelyrminxdc）
	dcn=zeros（nelyne