Topology Optimization - Multiple Loads with a circle

I'm reading Ole Sigmund's article called 'A 99 line toplogy optimization code written in MATLAB'. I'm following the instruction to modify the topology matlab code to recreate the shape optimized with a radius.

But my results ended up as follows:

Did I not substitute the lines correctly ? I've reread the article and didn't think I'm missing anything but if anyone can point out what I missed or not put code on the right line that'd be great.

Answer 1

[Torsten]

That's how I understand what changes have to be made to the code (although the free upper left zone isn't properly resolved).

nelx = 45;
nely = 30;
volfrac = 0.5;
penal = 3.0;
rmin = 1.5;
top(nelx,nely,volfrac,penal,rmin)

 It.:    1 Obj.:   412.2715 Vol.:  0.500 ch.:  0.200
 It.:    2 Obj.:   151.6200 Vol.:  0.500 ch.:  0.200
 It.:    3 Obj.:    87.9676 Vol.:  0.500 ch.:  0.200
 It.:    4 Obj.:    68.3247 Vol.:  0.500 ch.:  0.200
 It.:    5 Obj.:    62.7476 Vol.:  0.500 ch.:  0.200
 It.:    6 Obj.:    60.3618 Vol.:  0.500 ch.:  0.177
 It.:    7 Obj.:    59.1843 Vol.:  0.500 ch.:  0.137
 It.:    8 Obj.:    58.5941 Vol.:  0.500 ch.:  0.100
 It.:    9 Obj.:    58.2952 Vol.:  0.500 ch.:  0.079
 It.:   10 Obj.:    58.1486 Vol.:  0.500 ch.:  0.072
 It.:   11 Obj.:    58.0727 Vol.:  0.500 ch.:  0.077
 It.:   12 Obj.:    58.0302 Vol.:  0.500 ch.:  0.076
 It.:   13 Obj.:    57.9873 Vol.:  0.500 ch.:  0.071
 It.:   14 Obj.:    57.9625 Vol.:  0.500 ch.:  0.060
 It.:   15 Obj.:    57.9408 Vol.:  0.500 ch.:  0.052
 It.:   16 Obj.:    57.9246 Vol.:  0.500 ch.:  0.046
 It.:   17 Obj.:    57.9121 Vol.:  0.500 ch.:  0.040
 It.:   18 Obj.:    57.8949 Vol.:  0.500 ch.:  0.033
 It.:   19 Obj.:    57.8900 Vol.:  0.500 ch.:  0.027
 It.:   20 Obj.:    57.8876 Vol.:  0.500 ch.:  0.022
 It.:   21 Obj.:    57.8866 Vol.:  0.500 ch.:  0.018
 It.:   22 Obj.:    57.8857 Vol.:  0.500 ch.:  0.014
 It.:   23 Obj.:    57.8763 Vol.:  0.500 ch.:  0.012
 It.:   24 Obj.:    57.8825 Vol.:  0.500 ch.:  0.010

%%%% A 99 LINE TOPOLOGY OPTIMIZATION CODE BY OLE SIGMUND, JANUARY 2000 %%%
%%%% CODE MODIFIED FOR INCREASED SPEED, September 2002, BY OLE SIGMUND %%%
function top(nelx,nely,volfrac,penal,rmin);
% INITIALIZE
x(1:nely,1:nelx) = volfrac; 
passive = zeros(nely,nelx);
for ely = 1:nely
  for elx = 1:nelx
    if sqrt((ely-nely/2.)^2+(elx-nelx/3.)^2) < nely/3.
      passive(ely,elx) = 1;
      x(ely,elx) = 0.001;
    else
      passive(ely,elx) = 0;
    end
  end
end
loop = 0; 
change = 1.;
% START ITERATION
while change > 0.01  
  loop = loop + 1;
  xold = x;
% FE-ANALYSIS
  [U]=FE(nelx,nely,x,penal);         
% 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(ely,elx)^penal*Ue'*KE*Ue;
      dc(ely,elx) = -penal*x(ely,elx)^(penal-1)*Ue'*KE*Ue;
    end
  end
% FILTERING OF SENSITIVITIES
  [dc]   = check(nelx,nely,rmin,x,dc);    
% DESIGN UPDATE BY THE OPTIMALITY CRITERIA METHOD
  [x]    = OC(nelx,nely,x,volfrac,dc,passive); 
% 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 
end
%%%%%%%%%% OPTIMALITY CRITERIA UPDATE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [xnew]=OC(nelx,nely,x,volfrac,dc,passive)  
l1 = 0; l2 = 100000; move = 0.2;
while (l2-l1 > 1e-4)
  lmid = 0.5*(l2+l1);
  xnew = max(0.001,max(x-move,min(1.,min(x+move,x.*sqrt(-dc./lmid)))));
  xnew(find(passive)) = 0.001;
  if sum(sum(xnew)) - volfrac*nelx*nely > 0;
    l1 = lmid;
  else
    l2 = lmid;
  end
end
end
%%%%%%%%%% MESH-INDEPENDENCY FILTER %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [dcn]=check(nelx,nely,rmin,x,dc)
dcn=zeros(nely,nelx);
for i = 1:nelx
  for j = 1:nely
    sum=0.0; 
    for k = max(i-floor(rmin),1):min(i+floor(rmin),nelx)
      for l = max(j-floor(rmin),1):min(j+floor(rmin),nely)
        fac = rmin-sqrt((i-k)^2+(j-l)^2);
        sum = sum+max(0,fac);
        dcn(j,i) = dcn(j,i) + max(0,fac)*x(l,k)*dc(l,k);
      end
    end
    dcn(j,i) = dcn(j,i)/(x(j,i)*sum);
  end
end
end
%%%%%%%%%% FE-ANALYSIS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [U]=FE(nelx,nely,x,penal)
[KE] = lk; 
K = sparse(2*(nelx+1)*(nely+1), 2*(nelx+1)*(nely+1));
F = sparse(2*(nely+1)*(nelx+1),1); U = zeros(2*(nely+1)*(nelx+1),1);
for elx = 1:nelx
  for ely = 1:nely
    n1 = (nely+1)*(elx-1)+ely; 
    n2 = (nely+1)* elx   +ely;
    edof = [2*n1-1; 2*n1; 2*n2-1; 2*n2; 2*n2+1; 2*n2+2; 2*n1+1; 2*n1+2];
    K(edof,edof) = K(edof,edof) + x(ely,elx)^penal*KE;
  end
end
% DEFINE LOADS AND SUPPORTS (HALF MBB-BEAM)
F(2,1) = -1;
fixeddofs   = union([1:2:2*(nely+1)],[2*(nelx+1)*(nely+1)]);
alldofs     = [1:2*(nely+1)*(nelx+1)];
freedofs    = setdiff(alldofs,fixeddofs);
% SOLVING
U(freedofs,:) = K(freedofs,freedofs) \ F(freedofs,:);      
U(fixeddofs,:)= 0;
end
%%%%%%%%%% ELEMENT STIFFNESS MATRIX %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [KE]=lk
E = 1.; 
nu = 0.3;
k=[ 1/2-nu/6   1/8+nu/8 -1/4-nu/12 -1/8+3*nu/8 ... 
   -1/4+nu/12 -1/8-nu/8  nu/6       1/8-3*nu/8];
KE = E/(1-nu^2)*[ k(1) k(2) k(3) k(4) k(5) k(6) k(7) k(8)
                  k(2) k(1) k(8) k(7) k(6) k(5) k(4) k(3)
                  k(3) k(8) k(1) k(6) k(7) k(4) k(5) k(2)
                  k(4) k(7) k(6) k(1) k(8) k(3) k(2) k(5)
                  k(5) k(6) k(7) k(8) k(1) k(2) k(3) k(4)
                  k(6) k(5) k(4) k(3) k(2) k(1) k(8) k(7)
                  k(7) k(4) k(5) k(2) k(3) k(8) k(1) k(6)
                  k(8) k(3) k(2) k(5) k(4) k(7) k(6) k(1)];
end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This Matlab code was written by Ole Sigmund, Department of Solid         %
% Mechanics, Technical University of Denmark, DK-2800 Lyngby, Denmark.     %
% Please sent your comments to the author: [e-mail]              %
%                                                                          %
% The code is intended for educational purposes and theoretical details    %
% are discussed in the paper                                               %
% "A 99 line topology optimization code written in Matlab"                 %
% by Ole Sigmund (2001), Structural and Multidisciplinary Optimization,    %
% Vol 21, pp. 120--127.                                                    %
%                                                                          %
% The code as well as a postscript version of the paper can be             %
% downloaded from the web-site: http://www.topopt.dtu.dk                   %
%                                                                          %
% Disclaimer:                                                              %
% The author reserves all rights but does not guaranty that the code is    %
% free from errors. Furthermore, he shall not be liable in any event       %
% caused by the use of the program.                                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%