双向渐进结构拓扑优化BESO

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