# Had somebody implemented the MMA or GCMMA algorithm (Svanberg) in the 99 lines Topology Optimization Matlab code(Sigmund)?

Hello everybody,

I need some help. I have a lot of problems trying to implement the MMA or GCMMA optimizer in the Sigmund's Topology Optimization code (99 lines). Have you ever do that?. Please, is very important for my pregraduate thesis.

Thanks.

Felipe Cárdenas Díaz

Universidad de Concepción - Chile

### SIMP/MMA matlab code (again in 99 lines)

Hi Felipe,

I don't know if you still need this. I just came accross your post when searching for the same thing. Previously there was a web page maintained by Andreas Rietz where the SIMP/MMA code was downloadable (it is mentioned in the Sigmund's paper). I tried to access it but aparently it is not working anymore.

No worries though. I have previously downloaded the code for my Master's thesis and I've found it among my old files. I pasted the whole m file here. I hope it is helpful.

Kazem

%%%% TOPOLOGY OPTIMIZATION CODE BY ANDREAS RIETZ, FEBRUARY 1999 %%%
function TopOriginal(nelx,nely,volfrac,penal);
% INITIALIZE
x(1:nely,1:nelx) = volfrac;
loop = 0;
change = 1.;
% START ITERATION
for loop = 1:80
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
% DESIGN UPDATE BY THE OPTIMALITY CRITERIA METHOD
[x] = MMA(nelx,nely,x,volfrac,dc);
% 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
figure(1); colormap(flipud(gray)); imagesc(x); axis equal; axis tight;
axis off; colorbar; figure(1);
end
%%%%%%%%%% MMA UPDATE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [xnew]=MMA(nelx,nely,x,volfrac,dc)
xlow=0.001; xhigh=1;
L=x-0.1*(xhigh-xlow)*ones(nely,nelx);
high=(x-L).^2.*-dc./(xlow-L).^2;
low=(x-L).^2.*-dc./(xhigh-L).^2;
l2 = min(min(high));
l1 = max(max(low));
for i=1:50
lmid = 0.5*(l2+l1);
xnew=max(xlow,min(xhigh,L+abs(x-L).*sqrt(-dc./lmid)));
if sum(sum(xnew)) - volfrac*nelx*nely > 0;
l2 = lmid;
else
l1 = lmid;
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 = sparse(2*(nely+1)*(nelx+1),1);
for ely = 1:nely
for elx = 1:nelx
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
F(2*(nelx+1)*(nely+1),1) = -1;
fixeddofs = [1:2*(nely+1)];
alldofs = [1:2*(nely+1)*(nelx+1)];
freedofs = setdiff(alldofs,fixeddofs);
% SOLVING
U(freedofs,:) = K(freedofs,freedofs) \ F(freedofs,:);
U(fixeddofs,:)= 0;

ux=reshape(U(1:2:2*(nely+1)*(nelx+1)-1), ...
(nely+1),(nelx+1));
uy=reshape(U(2:2:2*(nely+1)*(nelx+1)), ...
(nely+1),(nelx+1));

figure(2);
quiver(ux,uy);
title('Displacements');

%%%%%%%%%% 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)];
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This Matlab code was written by Andreas Rietz, Department of Mechanics %
% at Linkoping University, Sweden %
% %
% The syntax of the code has been revised since february 1999 to be more %
% similar to the optimization code described in %
% "A 99 line topology optimization code written in Matlab" %
% by Ole Sigmund (To appear in Structural Optimization). %
% The code as well as a postscript version of the paper can be %
% downloaded from the web-site: http://www.topopt.dtu.dk %
% That code includes in addition a mesh independancy filter which makes %
% the result more reliable. %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%