%%%% 2D TOPOLOGY OPTIMIZATION CODE, MGCG ANALYSIS %%%% % Example: top2dmgcg(360,240,0.2,3.0,2.5,2,4,1e-10,100) function top2dmgcg(nelx,nely,volfrac,penal,rmin,ft,nl,cgtol,cgmax) %% MATERIAL PROPERTIES E0 = 1; Emin = 1e-9; nu = 0.3; %% PREPARE FINITE ELEMENT ANALYSIS KE = Ke2D(nu); % Prepare fine grid nelem = nelx*nely; nx = nelx+1; ny = nely+1; ndof = 2*nx*ny; nodenrs(1:ny,1:nx) = reshape(1:ny*nx,ny,nx); edofVec(1:nelem,1) = reshape(2*nodenrs(1:ny-1,1:nx-1)+1,nelem,1); edofMat = repmat(edofVec,1,8)+repmat([0 1 2*nely+[2 3 0 1] -2 -1],nelx*nely,1); iK = reshape(kron(edofMat,ones(8,1))',64*nelx*nely,1); jK = reshape(kron(edofMat,ones(1,8))',64*nelx*nely,1); % Prologation operators Pu = cell(nl-1,1); for l = 1:nl-1 [Pu{l,1}] = Prepcoarse(nely/2^(l-1),nelx/2^(l-1)); end % Define loads and supports (cantilever) F = sparse(2*(nelx*(nely+1)+nely/2+1),1,-1,2*(nely+1)*(nelx+1),1); U = zeros(2*(nely+1)*(nelx+1),1); fixeddofs = 1:2*(nely+1); % Null space elimination of supports N = ones(ndof,1); N(fixeddofs) = 0; Null = spdiags(N,0,ndof,ndof); %% PREPARE FILTER iH = ones(nelx*nely*(2*(ceil(rmin)-1)+1)^2,1); jH = ones(size(iH)); sH = zeros(size(iH)); k = 0; for i1 = 1:nelx for j1 = 1:nely e1 = (i1-1)*nely+j1; for i2 = max(i1-(ceil(rmin)-1),1):min(i1+(ceil(rmin)-1),nelx) for j2 = max(j1-(ceil(rmin)-1),1):min(j1+(ceil(rmin)-1),nely) e2 = (i2-1)*nely+j2; k = k+1; iH(k) = e1; jH(k) = e2; sH(k) = max(0,rmin-sqrt((i1-i2)^2+(j1-j2)^2)); end end end end H = sparse(iH,jH,sH); Hs = sum(H,2); %% INITIALIZE ITERATION x = repmat(volfrac,nely,nelx); xPhys = x; loop = 0; change = 1; %% START ITERATION while change > 1e-3 && loop < 500 loop = loop+1; %% FE-ANALYSIS K = cell(nl,1); sK = reshape(KE(:)*(Emin+xPhys(:)'.^penal*(E0-Emin)),64*nelem,1); K{1,1} = sparse(iK,jK,sK); K{1,1} = Null'*K{1,1}*Null - (Null-speye(ndof,ndof)); for l = 1:nl-1 K{l+1,1} = Pu{l,1}'*(K{l,1}*Pu{l,1}); end Lfac = chol(K{nl,1},'lower'); Ufac = Lfac'; [cgiters,cgres,U] = mgcg(K,F,U,Lfac,Ufac,Pu,nl,1,cgtol,cgmax); %% OBJECTIVE FUNCTION AND SENSITIVITY ANALYSIS ce = reshape(sum((U(edofMat)*KE).*U(edofMat),2),nely,nelx); c = sum(sum((Emin+xPhys.^penal*(E0-Emin)).*ce)); dc = -penal*(E0-Emin)*xPhys.^(penal-1).*ce; dv = ones(nely,nelx); %% FILTERING/MODIFICATION OF SENSITIVITIES if ft == 1 dc(:) = H*(x(:).*dc(:))./Hs./max(1e-3,x(:)); elseif ft == 2 dc(:) = H*(dc(:)./Hs); dv(:) = H*(dv(:)./Hs); end %% OPTIMALITY CRITERIA UPDATE OF DESIGN VARIABLES AND PHYSICAL DENSITIES g = mean(xPhys(:)) - volfrac; l1 = 0; l2 = 1e9; move = 0.2; while (l2-l1)/(l1+l2) > 1e-6 lmid = 0.5*(l2+l1); xnew = max(0,max(x-move,min(1,min(x+move,x.*sqrt(-dc./dv/lmid))))); gt = g + sum((dv(:).*(xnew(:)-x(:)))); if gt > 0, l1 = lmid; else l2 = lmid; end end change = max(abs(xnew(:)-x(:))); x = xnew; %% FILTERING OF DESIGN VARIABLES if ft == 1, xPhys = xnew; elseif ft == 2, xPhys(:) = (H*xnew(:))./Hs; end %% PRINT RESULTS fprintf(' Iter.: %4i Comp.: %6.3e Vol.: %6.3e Chg.: %4.2e CGres: %4.2e CGits: %3i Penal: %3.1f \n',... loop,c,mean(xPhys(:)),change,cgres,cgiters,penal); %% PLOT DENSITIES figure(1); clf; imagesc(1-xPhys); caxis([0 1]); axis equal; axis off; drawnow; end end %% FUNCTION mgcg - MULTIGRID PRECONDITIONED CONJUGATE GRADIENTS function [i,relres,u] = mgcg(A,b,u,Lfac,Ufac,Pu,nl,nswp,tol,maxiter) r = b - A{1,1}*u; res0 = norm(b); % Jacobi smoother omega = 0.8; invD = cell(nl-1,1); for l = 1:nl-1 invD{l,1}= 1./spdiags(A{l,1},0); end for i = 1:1e6 z = VCycle(A,r,Lfac,Ufac,Pu,1,nl,invD,omega,nswp); rho = r'*z; if i == 1 p = z; else beta = rho/rho_p; p = beta*p + z; end q = A{1,1}*p; dpr = p'*q; alpha = rho/dpr; u = u + alpha*p; r = r - alpha*q; rho_p = rho; relres = norm(r)/res0; if relres < tol || i>=maxiter break end end end %% FUNCTION VCycle - COARSE GRID CORRECTION function z = VCycle(A,r,Lfac,Ufac,Pu,l,nl,invD,omega,nswp) z = 0*r; z = smthdmpjac(z,A{l,1},r,invD{l,1},omega,nswp); Az = A{l,1}*z; d = r - Az; dh2 = Pu{l,1}'*d; if (nl == l+1) vh2 = Ufac \ (Lfac \ dh2); else vh2 = VCycle(A,dh2,Lfac,Ufac,Pu,l+1,nl,invD,omega,nswp); end v = Pu{l,1}*vh2; z = z + v; z = smthdmpjac(z,A{l,1},r,invD{l,1},omega,nswp); end %% FUNCTIODN smthdmpjac - DAMPED JACOBI SMOOTHER function [u] = smthdmpjac(u,A,b,invD,omega,nswp) for i = 1:nswp u = u - omega*invD.*(A*u) + omega*invD.*b; end end %% FUNCTION Ke2D - ELEMENT STIFFNESS MATRIX function KE = Ke2D(nu) A11 = [12 3 -6 -3; 3 12 3 0; -6 3 12 -3; -3 0 -3 12]; A12 = [-6 -3 0 3; -3 -6 -3 -6; 0 -3 -6 3; 3 -6 3 -6]; B11 = [-4 3 -2 9; 3 -4 -9 4; -2 -9 -4 -3; 9 4 -3 -4]; B12 = [ 2 -3 4 -9; -3 2 9 -2; 4 9 2 3; -9 -2 3 2]; KE = 1/(1-nu^2)/24*([A11 A12;A12' A11]+nu*[B11 B12;B12' B11]); end %% FUNCTION prepcoarse - PREPARE MG PROLONGATION OPERATOR function [Pu] = Prepcoarse(ney,nex) % Assemble state variable prolongation maxnum = nex*ney*20; iP = zeros(maxnum,1); jP = zeros(maxnum,1); sP = zeros(maxnum,1); nexc = nex/2; neyc = ney/2; % weights for fixed distances to neighbors on a structured grid vals = [1,0.5,0.25]; cc = 0; for nx = 1:nexc+1 for ny = 1:neyc+1 col = (nx-1)*(neyc+1)+ny; % Coordinate on fine grid nx1 = nx*2 - 1; ny1 = ny*2 - 1; % Loop over fine nodes within the rectangular domain for k = max(nx1-1,1):min(nx1+1,nex+1) for l = max(ny1-1,1):min(ny1+1,ney+1) row = (k-1)*(ney+1)+l; % Based on squared dist assign weights: 1.0 0.5 0.25 ind = 1+((nx1-k)^2+(ny1-l)^2); cc = cc+1; iP(cc) = 2*row-1; jP(cc) = 2*col-1; sP(cc) = vals(ind); cc = cc+1; iP(cc) = 2*row; jP(cc) = 2*col; sP(cc) = vals(ind); end end end end % Assemble matrices Pu = sparse(iP(1:cc),jP(1:cc),sP(1:cc)); end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % This code was written by Niels Aage and Boyan S. Lazarov, % % Department of Mechanical Engineering, Technical University of Denmark, % % and Oded Amir, Faculty of Civil & Environmental Engineering, % % Technion - Israel Institute of Technology. % % % % 3D code available upon request by email % % Contact: naage@mek.dtu.dk, bsl@mek.dtu.dk, odedamir@cv.technion.ac.il % % % % Details are discussed in the paper % % "On multigrid-CG for efficient topology optimization" accepted October % % 2013 to Structural and Multidisciplinary Optimization. % % % % Disclaimer: % % The authors reserve all rights but do not guarantee that the code is % % free from errors. Furthermore, the authors shall not be liable in any % % event caused by the use of the program. % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%