# Amira Abdel-Rahman
# (c) Massachusetts Institute of Technology 2020
function ConTop2D(Macro_struct, Micro_struct, penal, rmin)
## USER-DEFINED LOOP PARAMETERS
maxloop = 200; E0 = 1; Emin = 1e-9; nu = 0.3;
Macro_length = Macro_struct[1]; Macro_width = Macro_struct[2];
Micro_length = Micro_struct[1]; Micro_width = Micro_struct[2];
Macro_nelx = Macro_struct[3]; Macro_nely = Macro_struct[4];
Micro_nelx = Micro_struct[3]; Micro_nely = Micro_struct[4];
Macro_Vol = Macro_struct[5]; Micro_Vol = Micro_struct[5];
Macro_Elex = Macro_length/Macro_nelx; Macro_Eley = Macro_width/Macro_nely;
Macro_nele = Macro_nelx*Macro_nely; Micro_nele = Micro_nelx*Micro_nely;
Macro_ndof = 2*(Macro_nelx+1)*(Macro_nely+1);
# PREPARE FINITE ELEMENT ANALYSIS
[load_x, load_y] = meshgrid(Macro_nelx, Macro_nely/2);
loadnid = load_x*(Macro_nely+1)+(Macro_nely+1-load_y);
F = sparse(2*loadnid[:], 1, -1, 2*(Macro_nelx+1)*(Macro_nely+1),1);
U = zeros(Macro_ndof,1);
[fixed_x, fixed_y] = meshgrid(0, 0:Macro_nely);
fixednid = fixed_x*(Macro_nely+1)+(Macro_nely+1-fixed_y);
fixeddofs = [2*fixednid[:]; 2*fixednid[:]-1];
freedofs = setdiff(1:Macro_ndof,fixeddofs);
nodenrs = reshape(1:(Macro_nely+1)*(Macro_nelx+1),1+Macro_nely,1+Macro_nelx);
edofVec = reshape(2*nodenrs( [1:end-1,1:end-1] )+1,Macro_nele,1);
edofMat = repmat(edofVec,1,8)+repmat( [0 1 2*Macro_nely+[2 3 0 1] -2 -1],Macro_nele,1);
iK = reshape(kron(edofMat,ones(8,1))',64*Macro_nele,1);
jK = reshape(kron(edofMat,ones(1,8))',64*Macro_nele,1);
# PREPARE FILTER
Macro_H,Macro_Hs = filtering2d(Macro_nelx, Macro_nely, Macro_nele, rmin);
Micro_H,Micro_Hs = filtering2d(Micro_nelx, Micro_nely, Micro_nele, rmin);
# INITIALIZE ITERATION
Macro_x = repmat(Macro_Vol,Macro_nely,Macro_nelx);
Micro_x = ones(Micro_nely,Micro_nelx);
for i = 1:Micro_nelx
for j = 1:Micro_nely
if sqrt((i-Micro_nelx/2-0.5)^2+(j-Micro_nely/2-0.5)^2) < min(Micro_nelx,Micro_nely)/3
Micro_x(j,i) = 0;
end
end
end
beta = 1;
Macro_xTilde = Macro_x; Micro_xTilde = Micro_x;
Macro_xPhys = 1-exp(-beta*Macro_xTilde)+Macro_xTilde*exp(-beta);
Micro_xPhys = 1-exp(-beta*Micro_xTilde)+Micro_xTilde*exp(-beta);
loopbeta = 0; loop = 0; Macro_change = 1; Micro_change = 1;
while loop < maxloop || Macro_change > 0.01 || Micro_change > 0.01
loop = loop+1; loopbeta = loopbeta+1;
# FE-ANALYSIS AT TWO SCALES
DH, dDH = EBHM2D(Micro_xPhys, Micro_length, Micro_width, E0, Emin, nu, penal);
Ke = elementMatVec2D(Macro_Elex/2, Macro_Eley/2, DH);
sK = reshape(Ke[:]*(Emin+Macro_xPhys[:]'.^penal*(1-Emin)),64*Macro_nele,1);
K = sparse(iK,jK,sK); K = (K+K')/2;
U[freedofs,:] = K(freedofs,freedofs)\F(freedofs,:);
# OBJECTIVE FUNCTION AND SENSITIVITY ANALYSIS
ce = reshape(sum((U[edofMat]*Ke).*U[edofMat],2),Macro_nely,Macro_nelx);
c = sum(sum((Emin+Macro_xPhys.^penal*(1-Emin)).*ce));
Macro_dc = -penal*(1-Emin)*Macro_xPhys.^(penal-1).*ce;
Macro_dv = ones(Macro_nely, Macro_nelx);
Micro_dc = zeros(Micro_nely, Micro_nelx);
for i = 1:Micro_nele
dDHe = [dDH[1,1][i] dDH[1,2][i] dDH[1,3][i];
dDH[2,1][i] dDH[2,2][i] dDH[2,3][i];
dDH[3,1][i] dDH[3,2][i] dDH[3,3][i]];
dKE = elementMatVec2D(Macro_Elex, Macro_Eley, dDHe);
dce = reshape(sum((U[edofMat]*dKE).*U[edofMat],2),Macro_nely,Macro_nelx);
Micro_dc[i] = -sum(sum((Emin+Macro_xPhys.^penal*(1-Emin)).*dce));
end
Micro_dv = ones(Micro_nely, Micro_nelx);
# FILTERING AND MODIFICATION OF SENSITIVITIES
Macro_dx = beta*exp(-beta*Macro_xTilde)+exp(-beta);
Micro_dx = beta*exp(-beta*Micro_xTilde)+exp(-beta);
Macro_dc[:] = Macro_H*(Macro_dc[:].*Macro_dx[:]./Macro_Hs);
Macro_dv[:] = Macro_H*(Macro_dv[:].*Macro_dx[:]./Macro_Hs);
Micro_dc[:] = Micro_H*(Micro_dc[:].*Micro_dx[:]./Micro_Hs);
Micro_dv[:] = Micro_H*(Micro_dv[:].*Micro_dx[:]./Micro_Hs);
# OPTIMALITY CRITERIA UPDATE MACRO AND MICRO ELELMENT DENSITIES
Macro_x, Macro_xPhys, Macro_change = OC(Macro_x, Macro_dc, Macro_dv, Macro_H, Macro_Hs, Macro_Vol, Macro_nele, 0.2, beta);
Micro_x, Micro_xPhys, Micro_change = OC(Micro_x, Micro_dc, Micro_dv, Micro_H, Micro_Hs, Micro_Vol, Micro_nele, 0.2, beta);
Macro_xPhys = reshape(Macro_xPhys, Macro_nely, Macro_nelx);
Micro_xPhys = reshape(Micro_xPhys, Micro_nely, Micro_nelx);
# PRINT RESULTS
# display(" It.:#5i Obj.:#11.4f Macro_Vol.:#7.3f Micro_Vol.:#7.3f Macro_ch.:#7.3f Micro_ch.:#7.3f\n", loop,c,mean(Macro_xPhys[:] ),mean(Micro_xPhys[:] ), Macro_change, Micro_change);
# colormap(gray); imagesc(1-Macro_xPhys); caxis( [0 1] ); axis equal; axis off; drawnow;
# colormap(gray); imagesc(1-Micro_xPhys); caxis( [0 1] ); axis equal; axis off; drawnow;
## UPDATE HEAVISIDE REGULARIZATION PARAMETER
if beta < 512 && (loopbeta >= 50 || Macro_change <= 0.01 || Micro_change <= 0.01)
beta = 2*beta; loopbeta = 0; Macro_change = 1; Micro_change = 1;
display("Parameter beta increased to $beta.");
end
end
end
## SUB FUNCTION:filtering2D
function filtering2d(nelx, nely, nele, rmin)
iH = ones(nele*(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);
return H,Hs
end
## SUB FUNCTION: EBHM2D
function EBHM2D(den, lx, ly, E0, Emin, nu, penal)
# the initial definitions of the PUC
D0=E0/(1-nu^2)*[1 nu 0; nu 1 0; 0 0 (1-nu)/2]; # the elastic tensor
nely, nelx = size(den);
nele = nelx*nely;
dx = lx/nelx; dy = ly/nely;
Ke = elementMatVec2D(dx/2, dy/2, D0);
Num_node = (1+nely)*(1+nelx);
nodenrs = reshape(1:Num_node,1+nely,1+nelx);
edofVec = reshape(2*nodenrs[1:end-1,1:end-1]+1,nele,1);
edofMat = repmat(edofVec,1,8)+repmat( [0 1 2*nely+[2 3 0 1] -2 -1],nele,1);
# 3D periodic boundary formulation
alldofs = (1:2*(nely+1)*(nelx+1));
n1 = [nodenrs[end,[1,end]],nodenrs[1,[end,1]]];
d1 = reshape( [(2*n1-1);2*n1],1,8);
n3 = [nodenrs[2:end-1,1]',nodenrs[end,2:end-1]];
d3 = reshape( [(2*n3-1);2*n3],1,2*(nelx+nely-2));
n4 = [nodenrs[2:end-1,end]',nodenrs[1,2:end-1]];
d4 = reshape( [(2*n4-1);2*n4],1,2*(nelx+nely-2));
d2 = setdiff(alldofs,[d1,d3,d4] );
e0 = Matrix(1.0I, 3, 3);
ufixed = zeros(8,3);
for j = 1:3
ufixed[3:4,j] = [e0[1,j],e0[3,j]/2;e0[3,j]/2,e0[2,j]]*[lx;0];
ufixed[7:8,j] = [e0[1,j],e0[3,j]/2;e0[3,j]/2,e0[2,j]]*[0;ly];
ufixed[5:6,j] = ufixed[3:4,j]+ufixed[7:8,j];
end
wfixed = [repmat(ufixed[3:4,:],nely-1,1);repmat(ufixed[7:8,:],nelx-1,1)];
# the reduced elastic equilibrium equation to compute the induced displacement field
iK = reshape(kron(edofMat,ones(8,1))',64*nelx*nely,1);
jK = reshape(kron(edofMat,ones(1,8))',64*nelx*nely,1);
sK = reshape(Ke[:]*(Emin+den[:]'.^penal*(1-Emin)),64*nelx*nely,1);
K = sparse(iK,jK,sK); K = (K + K')/2;
Kr = [K[d2,d2],K[d2,d3]+K[d2,d4];K[d3,d2]+K[d4,d2],K[d3,d3]+K[d4,d3]+K[d3,d4]+K[d4,d4]];
U[d1,:]= ufixed;
U[[d2,d3],:] = Kr\(-[K[d2,d1];K[d3,d1]+K[d4,d1]]*ufixed-[K[d2,d4];K[d3,d4]+K[d4,d4]]*wfixed);
U[d4,:] = U[d3,:] + wfixed;
# homogenization to evaluate macroscopic effective properties
DH = zeros(3);
qe = Array{Any,2}(undef, 3, 3);
dDH = Array{Any,2}(undef, 3, 3);
cellVolume = lx*ly;
for i = 1:3
for j = 1:3
U1 = U[:,i]; U2 = U[:,j];
qe[i,j] = reshape(sum((U1[edofMat]*Ke).*U2[edofMat],dims=2),nely,nelx)/cellVolume;
DH[i,j] = sum(sum((Emin+den.^penal*(1-Emin)).*qe[i,j] ));
dDH[i,j] = penal*(1-Emin)*den.^(penal-1).*qe[i,j];
end
end
# disp("--- Homogenized elasticity tensor ---"); disp(DH)
return DH, dDH
end
## SUB FUNCTION: elementMatVec2D
function elementMatVec2D(a, b, DH)
GaussNodes = [-1/sqrt(3); 1/sqrt(3)];
GaussWeigh = [1 1];
L = [1 0 0 0; 0 0 0 1; 0 1 1 0];
Ke = zeros(8,8);
for i = 1:2
for j = 1:2
GN_x = GaussNodes[i];
GN_y = GaussNodes[j];
dN_x = 1/4*[-(1-GN_x) (1-GN_x) (1+GN_x) -(1+GN_x)];
dN_y = 1/4*[-(1-GN_y) -(1+GN_y) (1+GN_y) (1-GN_y)];
J = [dN_x; dN_y]*[ -a a a -a; -b -b b b]';
G = [inv(J) zeros(size(J)); zeros(size(J)) inv(J)];
dN[1,1:2:8] = dN_x;
dN[2,1:2:8] = dN_y;
dN[3,2:2:8] = dN_x;
dN[4,2:2:8] = dN_y;
Be = L*G*dN;
Ke = Ke + GaussWeigh[i]*GaussWeigh[j]*det(J)*Be'*DH*Be;
end
end
return Ke
end
## SUB FUNCTION: OC
function OC(x, dc, dv, H, Hs, volfrac, nele, move, beta)
l1 = 0; l2 = 1e9;
while (l2-l1)/(l1+l2) > 1e-4
lmid = 0.5*(l2+l1);
xnew = max(0,max(x-move,min(1,min(x+move,x.*sqrt(-dc./dv/lmid)))));
xTilde[:] = (H*xnew[:] )./Hs; xPhys = 1-exp(-beta*xTilde)+xTilde*exp(-beta);
if sum(xPhys[:] ) > volfrac*nele
l1 = lmid;
else
l2 = lmid;
end
end
change = max(abs(xnew[:]-x[:] )); x = xnew;
return x, xPhys, change
end
#======================================================================================================================#
# Function ConTop2D: #
# A compact and efficient MATLAB code for Concurrent topology optimization of multiscale composite structures #
# in Matlab. #
# #
# Developed by: Jie Gao, Zhen Luo, Liang Xia and Liang Gao* #
# Email: gaoliang@mail.hust.edu.cn (GabrielJie_Tian@163.com) #
# #
# Main references: #
# #
# (1) Jie Gao, Zhen Luo, Liang Xia, Liang Gao. Concurrent topology optimization of multiscale composite structures #
# in Matlab. Accepted in Structural and multidisciplinary optimization. #
# #
# (2) Xia L, Breitkopf P. Design of materials using topology optimization and energy-based homogenization approach in #
# Matlab. # Structural and multidisciplinary optimization, 2015, 52(6): 1229-1241. #
# #
# ********************************************* Disclaimer ******************************************************* #
# The authors reserve all rights for the programs. The programs may be distributed and used for academic and #
# educational purposes. The authors do not guarantee that the code is free from errors,and they shall not be liable #
# in any event caused by the use of the program. #
#======================================================================================================================#