传热结构拓扑优化

	%%%% 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