拓扑优化99行程序

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