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.

top_circle(45,30,0.5,3.0,1.5)

 It.:    1 Obj.:   370.8613 Vol.:  0.500 ch.:  0.200
 It.:    2 Obj.:   136.5481 Vol.:  0.500 ch.:  0.200
 It.:    3 Obj.:    79.9951 Vol.:  0.500 ch.:  0.200
 It.:    4 Obj.:    61.5256 Vol.:  0.500 ch.:  0.200
 It.:    5 Obj.:    55.0482 Vol.:  0.500 ch.:  0.200
 It.:    6 Obj.:    51.4752 Vol.:  0.500 ch.:  0.200
 It.:    7 Obj.:    49.3565 Vol.:  0.500 ch.:  0.200
 It.:    8 Obj.:    48.1120 Vol.:  0.500 ch.:  0.200
 It.:    9 Obj.:    47.3086 Vol.:  0.500 ch.:  0.200
 It.:   10 Obj.:    46.7398 Vol.:  0.500 ch.:  0.200
 It.:   11 Obj.:    46.3488 Vol.:  0.500 ch.:  0.200
 It.:   12 Obj.:    46.0484 Vol.:  0.500 ch.:  0.200
 It.:   13 Obj.:    45.8206 Vol.:  0.500 ch.:  0.190
 It.:   14 Obj.:    45.6221 Vol.:  0.500 ch.:  0.177
 It.:   15 Obj.:    45.4580 Vol.:  0.500 ch.:  0.148
 It.:   16 Obj.:    45.3244 Vol.:  0.500 ch.:  0.142
 It.:   17 Obj.:    45.1887 Vol.:  0.500 ch.:  0.119
 It.:   18 Obj.:    45.0904 Vol.:  0.500 ch.:  0.112
 It.:   19 Obj.:    44.9828 Vol.:  0.500 ch.:  0.106
 It.:   20 Obj.:    44.8901 Vol.:  0.500 ch.:  0.094
 It.:   21 Obj.:    44.8093 Vol.:  0.500 ch.:  0.090
 It.:   22 Obj.:    44.7224 Vol.:  0.500 ch.:  0.091
 It.:   23 Obj.:    44.6303 Vol.:  0.500 ch.:  0.094
 It.:   24 Obj.:    44.5481 Vol.:  0.500 ch.:  0.092
 It.:   25 Obj.:    44.4638 Vol.:  0.500 ch.:  0.086
 It.:   26 Obj.:    44.3994 Vol.:  0.500 ch.:  0.086
 It.:   27 Obj.:    44.3300 Vol.:  0.500 ch.:  0.097
 It.:   28 Obj.:    44.2618 Vol.:  0.500 ch.:  0.101
 It.:   29 Obj.:    44.2031 Vol.:  0.500 ch.:  0.103
 It.:   30 Obj.:    44.1347 Vol.:  0.500 ch.:  0.118
 It.:   31 Obj.:    44.0238 Vol.:  0.500 ch.:  0.150
 It.:   32 Obj.:    43.8487 Vol.:  0.500 ch.:  0.169
 It.:   33 Obj.:    43.6123 Vol.:  0.500 ch.:  0.135
 It.:   34 Obj.:    43.4040 Vol.:  0.500 ch.:  0.090
 It.:   35 Obj.:    43.3324 Vol.:  0.500 ch.:  0.085
 It.:   36 Obj.:    43.3033 Vol.:  0.500 ch.:  0.078
 It.:   37 Obj.:    43.2664 Vol.:  0.500 ch.:  0.083
 It.:   38 Obj.:    43.2299 Vol.:  0.500 ch.:  0.093
 It.:   39 Obj.:    43.1722 Vol.:  0.500 ch.:  0.098
 It.:   40 Obj.:    43.1052 Vol.:  0.500 ch.:  0.105
 It.:   41 Obj.:    43.0321 Vol.:  0.500 ch.:  0.103
 It.:   42 Obj.:    42.9657 Vol.:  0.500 ch.:  0.093
 It.:   43 Obj.:    42.8893 Vol.:  0.500 ch.:  0.097
 It.:   44 Obj.:    42.8261 Vol.:  0.500 ch.:  0.093
 It.:   45 Obj.:    42.7773 Vol.:  0.500 ch.:  0.091
 It.:   46 Obj.:    42.7308 Vol.:  0.500 ch.:  0.093
 It.:   47 Obj.:    42.6651 Vol.:  0.500 ch.:  0.101
 It.:   48 Obj.:    42.6019 Vol.:  0.500 ch.:  0.116
 It.:   49 Obj.:    42.5175 Vol.:  0.500 ch.:  0.134
 It.:   50 Obj.:    42.3993 Vol.:  0.500 ch.:  0.140
 It.:   51 Obj.:    42.2894 Vol.:  0.500 ch.:  0.078
 It.:   52 Obj.:    42.2153 Vol.:  0.500 ch.:  0.070
 It.:   53 Obj.:    42.1817 Vol.:  0.500 ch.:  0.063
 It.:   54 Obj.:    42.1584 Vol.:  0.500 ch.:  0.064
 It.:   55 Obj.:    42.1575 Vol.:  0.500 ch.:  0.067
 It.:   56 Obj.:    42.1394 Vol.:  0.500 ch.:  0.076
 It.:   57 Obj.:    42.1184 Vol.:  0.500 ch.:  0.092
 It.:   58 Obj.:    42.0698 Vol.:  0.500 ch.:  0.119
 It.:   59 Obj.:    41.9871 Vol.:  0.500 ch.:  0.146
 It.:   60 Obj.:    41.8348 Vol.:  0.500 ch.:  0.148
 It.:   61 Obj.:    41.6270 Vol.:  0.500 ch.:  0.093
 It.:   62 Obj.:    41.5037 Vol.:  0.500 ch.:  0.075
 It.:   63 Obj.:    41.4652 Vol.:  0.500 ch.:  0.069
 It.:   64 Obj.:    41.4498 Vol.:  0.500 ch.:  0.063
 It.:   65 Obj.:    41.4404 Vol.:  0.500 ch.:  0.057
 It.:   66 Obj.:    41.4307 Vol.:  0.500 ch.:  0.059
 It.:   67 Obj.:    41.4355 Vol.:  0.500 ch.:  0.061
 It.:   68 Obj.:    41.4252 Vol.:  0.500 ch.:  0.060
 It.:   69 Obj.:    41.4249 Vol.:  0.500 ch.:  0.056
 It.:   70 Obj.:    41.4245 Vol.:  0.500 ch.:  0.053
 It.:   71 Obj.:    41.4326 Vol.:  0.500 ch.:  0.056
 It.:   72 Obj.:    41.4337 Vol.:  0.500 ch.:  0.052
 It.:   73 Obj.:    41.4443 Vol.:  0.500 ch.:  0.048
 It.:   74 Obj.:    41.4451 Vol.:  0.500 ch.:  0.046
 It.:   75 Obj.:    41.4482 Vol.:  0.500 ch.:  0.053
 It.:   76 Obj.:    41.4498 Vol.:  0.500 ch.:  0.061
 It.:   77 Obj.:    41.4545 Vol.:  0.500 ch.:  0.070
 It.:   78 Obj.:    41.4501 Vol.:  0.500 ch.:  0.079
 It.:   79 Obj.:    41.4458 Vol.:  0.500 ch.:  0.083
 It.:   80 Obj.:    41.4275 Vol.:  0.500 ch.:  0.080
 It.:   81 Obj.:    41.4170 Vol.:  0.500 ch.:  0.099
 It.:   82 Obj.:    41.3895 Vol.:  0.500 ch.:  0.123
 It.:   83 Obj.:    41.3262 Vol.:  0.500 ch.:  0.143
 It.:   84 Obj.:    41.2241 Vol.:  0.500 ch.:  0.140
 It.:   85 Obj.:    41.1081 Vol.:  0.500 ch.:  0.112
 It.:   86 Obj.:    41.0167 Vol.:  0.500 ch.:  0.052
 It.:   87 Obj.:    40.9928 Vol.:  0.500 ch.:  0.038
 It.:   88 Obj.:    40.9961 Vol.:  0.500 ch.:  0.036
 It.:   89 Obj.:    41.0010 Vol.:  0.500 ch.:  0.037
 It.:   90 Obj.:    41.0024 Vol.:  0.500 ch.:  0.036
 It.:   91 Obj.:    41.0041 Vol.:  0.500 ch.:  0.034
 It.:   92 Obj.:    41.0054 Vol.:  0.500 ch.:  0.035
 It.:   93 Obj.:    40.9973 Vol.:  0.500 ch.:  0.035
 It.:   94 Obj.:    40.9992 Vol.:  0.500 ch.:  0.035
 It.:   95 Obj.:    41.0016 Vol.:  0.500 ch.:  0.034
 It.:   96 Obj.:    41.0062 Vol.:  0.500 ch.:  0.034
 It.:   97 Obj.:    41.0023 Vol.:  0.500 ch.:  0.032
 It.:   98 Obj.:    41.0080 Vol.:  0.500 ch.:  0.031
 It.:   99 Obj.:    41.0042 Vol.:  0.500 ch.:  0.028
 It.:  100 Obj.:    41.0108 Vol.:  0.500 ch.:  0.028
 It.:  101 Obj.:    41.0065 Vol.:  0.500 ch.:  0.025
 It.:  102 Obj.:    41.0127 Vol.:  0.500 ch.:  0.024
 It.:  103 Obj.:    41.0089 Vol.:  0.500 ch.:  0.023
 It.:  104 Obj.:    41.0069 Vol.:  0.500 ch.:  0.022
 It.:  105 Obj.:    41.0121 Vol.:  0.500 ch.:  0.021
 It.:  106 Obj.:    41.0120 Vol.:  0.500 ch.:  0.022
 It.:  107 Obj.:    41.0111 Vol.:  0.500 ch.:  0.021
 It.:  108 Obj.:    41.0197 Vol.:  0.500 ch.:  0.021
 It.:  109 Obj.:    41.0207 Vol.:  0.500 ch.:  0.022
 It.:  110 Obj.:    41.0220 Vol.:  0.500 ch.:  0.022
 It.:  111 Obj.:    41.0215 Vol.:  0.500 ch.:  0.023
 It.:  112 Obj.:    41.0214 Vol.:  0.500 ch.:  0.023
 It.:  113 Obj.:    41.0208 Vol.:  0.500 ch.:  0.023
 It.:  114 Obj.:    41.0320 Vol.:  0.500 ch.:  0.023
 It.:  115 Obj.:    41.0249 Vol.:  0.500 ch.:  0.022
 It.:  116 Obj.:    41.0380 Vol.:  0.500 ch.:  0.023
 It.:  117 Obj.:    41.0407 Vol.:  0.500 ch.:  0.023
 It.:  118 Obj.:    41.0351 Vol.:  0.500 ch.:  0.021
 It.:  119 Obj.:    41.0445 Vol.:  0.500 ch.:  0.021
 It.:  120 Obj.:    41.0461 Vol.:  0.500 ch.:  0.020
 It.:  121 Obj.:    41.0493 Vol.:  0.500 ch.:  0.020
 It.:  122 Obj.:    41.0503 Vol.:  0.500 ch.:  0.019
 It.:  123 Obj.:    41.0514 Vol.:  0.500 ch.:  0.018
 It.:  124 Obj.:    41.0516 Vol.:  0.500 ch.:  0.018
 It.:  125 Obj.:    41.0517 Vol.:  0.500 ch.:  0.018
 It.:  126 Obj.:    41.0506 Vol.:  0.500 ch.:  0.017
 It.:  127 Obj.:    41.0597 Vol.:  0.500 ch.:  0.017
 It.:  128 Obj.:    41.0609 Vol.:  0.500 ch.:  0.017
 It.:  129 Obj.:    41.0522 Vol.:  0.500 ch.:  0.017
 It.:  130 Obj.:    41.0592 Vol.:  0.500 ch.:  0.017
 It.:  131 Obj.:    41.0578 Vol.:  0.500 ch.:  0.017
 It.:  132 Obj.:    41.0590 Vol.:  0.500 ch.:  0.018
 It.:  133 Obj.:    41.0575 Vol.:  0.500 ch.:  0.017
 It.:  134 Obj.:    41.0564 Vol.:  0.500 ch.:  0.017
 It.:  135 Obj.:    41.0548 Vol.:  0.500 ch.:  0.016
 It.:  136 Obj.:    41.0638 Vol.:  0.500 ch.:  0.016
 It.:  137 Obj.:    41.0642 Vol.:  0.500 ch.:  0.016
 It.:  138 Obj.:    41.0556 Vol.:  0.500 ch.:  0.016
 It.:  139 Obj.:    41.0627 Vol.:  0.500 ch.:  0.016
 It.:  140 Obj.:    41.0619 Vol.:  0.500 ch.:  0.016
 It.:  141 Obj.:    41.0651 Vol.:  0.500 ch.:  0.018
 It.:  142 Obj.:    41.0676 Vol.:  0.500 ch.:  0.018
 It.:  143 Obj.:    41.0700 Vol.:  0.500 ch.:  0.018
 It.:  144 Obj.:    41.0712 Vol.:  0.500 ch.:  0.018
 It.:  145 Obj.:    41.0727 Vol.:  0.500 ch.:  0.018
 It.:  146 Obj.:    41.0739 Vol.:  0.500 ch.:  0.017
 It.:  147 Obj.:    41.0749 Vol.:  0.500 ch.:  0.017
 It.:  148 Obj.:    41.0754 Vol.:  0.500 ch.:  0.017
 It.:  149 Obj.:    41.0762 Vol.:  0.500 ch.:  0.017
 It.:  150 Obj.:    41.0773 Vol.:  0.500 ch.:  0.016
 It.:  151 Obj.:    41.0791 Vol.:  0.500 ch.:  0.017
 It.:  152 Obj.:    41.0798 Vol.:  0.500 ch.:  0.016
 It.:  153 Obj.:    41.0901 Vol.:  0.500 ch.:  0.017
 It.:  154 Obj.:    41.0830 Vol.:  0.500 ch.:  0.015
 It.:  155 Obj.:    41.0931 Vol.:  0.500 ch.:  0.016
 It.:  156 Obj.:    41.0841 Vol.:  0.500 ch.:  0.014
 It.:  157 Obj.:    41.0946 Vol.:  0.500 ch.:  0.015
 It.:  158 Obj.:    41.0952 Vol.:  0.500 ch.:  0.015
 It.:  159 Obj.:    41.0881 Vol.:  0.500 ch.:  0.014
 It.:  160 Obj.:    41.0956 Vol.:  0.500 ch.:  0.014
 It.:  161 Obj.:    41.0960 Vol.:  0.500 ch.:  0.014
 It.:  162 Obj.:    41.0972 Vol.:  0.500 ch.:  0.014
 It.:  163 Obj.:    41.0966 Vol.:  0.500 ch.:  0.014
 It.:  164 Obj.:    41.0960 Vol.:  0.500 ch.:  0.013
 It.:  165 Obj.:    41.0945 Vol.:  0.500 ch.:  0.012
 It.:  166 Obj.:    41.1027 Vol.:  0.500 ch.:  0.014
 It.:  167 Obj.:    41.0932 Vol.:  0.500 ch.:  0.013
 It.:  168 Obj.:    41.1023 Vol.:  0.500 ch.:  0.013
 It.:  169 Obj.:    41.1024 Vol.:  0.500 ch.:  0.014
 It.:  170 Obj.:    41.0939 Vol.:  0.500 ch.:  0.013
 It.:  171 Obj.:    41.1006 Vol.:  0.500 ch.:  0.013
 It.:  172 Obj.:    41.1001 Vol.:  0.500 ch.:  0.013
 It.:  173 Obj.:    41.1009 Vol.:  0.500 ch.:  0.013
 It.:  174 Obj.:    41.1001 Vol.:  0.500 ch.:  0.013
 It.:  175 Obj.:    41.0994 Vol.:  0.500 ch.:  0.013
 It.:  176 Obj.:    41.0998 Vol.:  0.500 ch.:  0.013
 It.:  177 Obj.:    41.1105 Vol.:  0.500 ch.:  0.014
 It.:  178 Obj.:    41.1029 Vol.:  0.500 ch.:  0.013
 It.:  179 Obj.:    41.1138 Vol.:  0.500 ch.:  0.013
 It.:  180 Obj.:    41.1046 Vol.:  0.500 ch.:  0.012
 It.:  181 Obj.:    41.1151 Vol.:  0.500 ch.:  0.012
 It.:  182 Obj.:    41.1155 Vol.:  0.500 ch.:  0.012
 It.:  183 Obj.:    41.1084 Vol.:  0.500 ch.:  0.011
 It.:  184 Obj.:    41.1165 Vol.:  0.500 ch.:  0.011
 It.:  185 Obj.:    41.1174 Vol.:  0.500 ch.:  0.011
 It.:  186 Obj.:    41.1200 Vol.:  0.500 ch.:  0.010
 It.:  187 Obj.:    41.1208 Vol.:  0.500 ch.:  0.010
 It.:  188 Obj.:    41.1215 Vol.:  0.500 ch.:  0.010

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