# Amira Abdel-Rahman
# (c) Massachusetts Institute of Technology 2020
#############################################2d Periodic############################################################
# Based on https://link.springer.com/article/10.1007/s00158-015-1294-0
## PERIODIC MATERIAL MICROSTRUCTURE DESIGN
function topX(nelx,nely,volfrac,penal,rmin,ft,optType="bulk")
## MATERIAL PROPERTIES
E0 = 1;
Emin = 1e-9;
nu = 0.3;
## PREPARE FINITE ELEMENT ANALYSIS
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])
nodenrs = reshape(1:(1+nelx)*(1+nely),1+nely,1+nelx)
edofVec = reshape(2*nodenrs[1:end-1,1:end-1].+1,nelx*nely,1)
edofMat = repeat(edofVec,1,8).+repeat([0 1 2*nely.+[2 3 0 1] -2 -1],nelx*nely,1)
iK = convert(Array{Int},reshape(kron(edofMat,ones(8,1))',64*nelx*nely,1))
jK = convert(Array{Int},reshape(kron(edofMat,ones(1,8))',64*nelx*nely,1))
## PREPARE FILTER
iH = ones(convert(Int,nelx*nely*(2*(ceil(rmin)-1)+1)^2),1)
jH = ones(Int,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(vec(iH),vec(jH),vec(sH))
Hs = sum(H,dims=2)
## PERIODIC BOUNDARY CONDITIONS
e0 = Matrix(1.0I, 3, 3);
ufixed = zeros(8,3);
U = zeros(2*(nely+1)*(nelx+1),3);
alldofs = [1:2*(nely+1)*(nelx+1)];
n1 = vcat(nodenrs[end,[1,end]],nodenrs[1,[end,1]]);
d1 = vec(reshape([(2 .* n1 .-1) 2 .*n1]',1,8));
n3 = [vec(nodenrs[2:(end-1),1]');vec(nodenrs[end,2:(end-1)])];
d3 = vec(reshape([(2 .*n3 .-1) 2 .*n3]',1,2*(nelx+nely-2)));
n4 = [vec(nodenrs[2:end-1,end]');vec(nodenrs[1,2:end-1])];
d4 = vec(reshape([(2 .*n4 .-1) 2 .*n4]',1,2*(nelx+nely-2)));
d2 = setdiff(vcat(alldofs...),union(union(d1,d3),d4));
for j = 1:3
ufixed[3:4,j] = [e0[1,j] e0[3,j]/2 ; e0[3,j]/2 e0[2,j]]*[nelx;0];
ufixed[7:8,j] = [e0[1,j] e0[3,j]/2 ; e0[3,j]/2 e0[2,j]]*[0;nely];
ufixed[5:6,j] = ufixed[3:4,j] .+ ufixed[7:8,j];
end
wfixed = [repeat(ufixed[3:4,:],nely-1,1); repeat(ufixed[7:8,:],nelx-1,1)];
## INITIALIZE ITERATION
qe = Array{Any,2}(undef, 3, 3);
Q = zeros(3,3);
dQ = Array{Any,2}(undef, 3, 3);
x = volfrac.*ones(nely,nelx)
for i = 1:nelx
for j = 1:nely
vall=3
if optType=="poisson"
vall=6
end
if sqrt((i-nelx/2-0.5)^2+(j-nely/2-0.5)^2) < min(nelx,nely)/vall
x[j,i] = volfrac/2.0;
end
end
end
xPhys = copy(x);
change = 1;
loop = 0;
xnew=zeros(size(x))
## START ITERATION
while (change > 0.01)
loop = loop +1;
## FE-ANALYSIS
sK = reshape(KE[:]*(Emin .+ xPhys[:]'.^penal*(80 .- Emin)),64*nelx*nely,1);
K = sparse(vec(iK),vec(jK),vec(sK));
K = (K.+K')./2.0;
Kr = vcat(hcat(K[d2,d2] , K[d2,d3]+K[d2,d4]),hcat((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;
## OBJECTIVE FUNCTION AND SENSITIVITY ANALYSIS
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)./(nelx*nely);
Q[i,j] = sum(sum((Emin .+ xPhys.^penal*(E0 .-Emin)).*qe[i,j]));
dQ[i,j] = penal*(E0-Emin)*xPhys.^(penal-1).*qe[i,j];
end
end
if optType=="bulk"
#bulk
c = -(Q[1,1]+Q[2,2]+Q[1,2]+Q[2,1]);
dc = -(dQ[1,1]+dQ[2,2]+dQ[1,2]+dQ[2,1]);
elseif optType=="shear"
#shear
c=-Q[3,3];
dc=-dQ[3,3];
elseif optType=="poisson"
c = Q[1,2]-(0.8^loop)*(Q[1,1]+Q[2,2]);
dc = dQ[1,2]-(0.8^loop)*(dQ[1,1]+dQ[2,2]);
end
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
if optType=="poisson"
l1=0; l2=1e9; move = 0.1;
else
l1 =0; l2 = 1e9; move = 0.2;
end
while (l2-l1 > 1e-9)
lmid = 0.5*(l2+l1);
if optType=="poisson"
xnew = max.(0.0,max.(x.-move,min.(1.0,min.(x.+move,x.*(-dc./dv./lmid)))));
else
xnew = max.(0.0,max.(x.-move,min.(1.0,min.(x.+move,x.*sqrt.(0.0.-dc./dv./lmid)))));
end
if ft == 1
xPhys = copy(xnew);
elseif ft == 2
xPhys[:] = (H*xnew[:])./Hs;
end
if mean(xPhys[:]) > volfrac
l1 = lmid;
else
l2 = lmid;
end
end
change = maximum(abs.(xnew[:].-x[:]))
x = xnew;
## PRINT RESULTS
display(" It:$loop Obj:$c Vol:$(mean(xPhys[:])) ch:$change ")
## PLOT DENSITIES
heatmap(xPhys, aspect_ratio=:equal, legend=false, axis=nothing, foreground_color_subplot=colorant"white",fc=:grays,clims=(0.0, 1.0))
frame(anim)
end
return xnew
end
######################3d Periodic#########################################
#INCORRECT!! FIX LATER
function topX3D(nelx,nely,nelz,volfrac,penal,rmin,ft,optType="bulk")
nel=nelx*nely*nelz
# MATERIAL PROPERTIES
E0 = 1;
Emin = 1e-9;
nu = 0.3;
lx=0.1;ly=0.1;lz=0.1;
vert_cor = [0 lx lx 0 0 lx lx 0;
0 0 ly ly 0 0 ly ly;
0 0 0 0 lz lz lz lz];
cellVolume = lx*ly*lz;
# PREPARE FINITE ELEMENT ANALYSIS
# KE=lk_H8(nu)
# the initial definitions of the PUC
D0 = E0/(1+nu)/(1-2*nu)*
[ 1-nu nu nu 0 0 0 ;
nu 1-nu nu 0 0 0 ;
nu nu 1-nu 0 0 0 ;
0 0 0 (1-2*nu)/2 0 0 ;
0 0 0 0 (1-2*nu)/2 0 ;
0 0 0 0 0 (1-2*nu)/2];
dx = lx/nelx; dy = ly/nely; dz = lz/nelz;
KE = elementMatVec3D(dx/2, dy/2, dz/2, D0);
Num_node = (1+nely)*(1+nelx)*(1+nelz);
nele = nelx*nely*nelz;
nodenrs = reshape(1:Num_node,1+nely,1+nelx,1+nelz);
edofVec = reshape(3*nodenrs[1:end-1,1:end-1,1:end-1] .+1,nelx*nely*nelz,1);
edofMat = repeat(edofVec,1,24).+repeat([0 1 2 3*nely.+[3 4 5 0 1 2] -3 -2 -1 3*(nelx+1)*(nely+1).+[0 1 2 3*nely.+[3 4 5 0 1 2] -3 -2 -1]], nelx*nely*nelz, 1);
# nodenrs = reshape(1:(1+nelx)*(1+nely)*(1+nelz),1+nely,1+nelx,1+nelz)
# edofVec = reshape(3*nodenrs[1:end-1,1:end-1,1:end-1].+1,nelx*nely*nelz,1)
# edofMat = repeat(edofVec,1,24).+repeat([0 1 2 3*nely.+[3 4 5 0 1 2] -3 -2 -1 3*(nely+1)*(nelx+1).+[0 1 2 3*nely.+[3 4 5 0 1 2] -3 -2 -1]],nelx*nely*nelz,1)
iK = reshape(kron(edofMat,ones(24,1))',24*24*nele,1);
jK = reshape(kron(edofMat,ones(1,24))',24*24*nele,1);
## PREPARE FILTER
H,Hs=make_filter3D(nelx,nely,nelz,rmin);
## PERIODIC BOUNDARY CONDITIONS
e0 = Matrix(1.0I, 6, 6);
ufixed = zeros(24,6);
ndof = 3*(nelx+1)*(nely+1)*(nelz+1);
U = zeros(ndof,6);
alldofs = [1:ndof];
# 3D periodic boundary formulation
# the nodes classification
n1 = hcat(nodenrs[end, [1 end], 1], nodenrs[1, [end 1], 1], nodenrs[end, [1 end], end] ,nodenrs[1, [end 1], end]);
d1 = Int.(vec(reshape(vcat(3.0.*n1.-2, 3.0.*n1.-1,3.0.*n1),3*length(n1),1)));
n3 = vcat(vec(reshape(nodenrs[end,1,2:end-1] ,1,length(nodenrs[end,1,2:end-1] ))), # AE
vec(reshape(nodenrs[1, 1, 2:end-1] ,1,length(nodenrs[1, 1, 2:end-1] ))), # DH
vec(reshape(nodenrs[end,2:end-1,1] ,1,length(nodenrs[end,2:end-1,1] ))), # AB
vec(reshape(nodenrs[1, 2:end-1, 1] ,1,length(nodenrs[1, 2:end-1, 1] ))), # DC
vec(reshape(nodenrs[2:end-1, 1, 1] ,1,length(nodenrs[2:end-1, 1, 1] ))), # AD
vec(reshape(nodenrs[2:end-1,1,end] ,1,length(nodenrs[2:end-1,1,end] ))), # EH
vec(reshape(nodenrs[2:end-1, 2:end-1, 1] ,1,length(nodenrs[2:end-1, 2:end-1, 1] ))), # ABCD
vec(reshape(nodenrs[2:end-1, 1, 2:end-1] ,1,length(nodenrs[2:end-1, 1, 2:end-1] ))), # ADHE
vec(reshape(nodenrs[end,2:end-1,2:end-1] ,1,length(nodenrs[end,2:end-1,2:end-1] ))))'; # ABFE
d3 = vec(Int.(reshape( vcat(3.0.*n3.-2 ,3.0.*n3.-1, 3.0.*n3),3*length(n3),1)));
n4 = vcat(vec(reshape(nodenrs[1, end, 2:end-1] ,1,length((nodenrs[1, end, 2:end-1] )))), # CG
vec(reshape(nodenrs[end,end,2:end-1] ,1,length((nodenrs[end,end,2:end-1] )))), # BF
vec(reshape(nodenrs[1, 2:end-1, end] ,1,length((nodenrs[1, 2:end-1, end] )))), # HG
vec(reshape(nodenrs[end,2:end-1,end] ,1,length((nodenrs[end,2:end-1,end] )))), # EF
vec(reshape(nodenrs[2:end-1,end,end] ,1,length((nodenrs[2:end-1,end,end] )))), # FG
vec(reshape(nodenrs[2:end-1, end, 1] ,1,length((nodenrs[2:end-1, end, 1] )))), # BC
vec(reshape(nodenrs[2:end-1,2:end-1,end] ,1,length((nodenrs[2:end-1,2:end-1,end] )))), # EFGH
vec(reshape(nodenrs[2:end-1,end,2:end-1] ,1,length((nodenrs[2:end-1,end,2:end-1] )))), # BCGF
vec(reshape(nodenrs[1, 2:end-1, 2:end-1] ,1,length((nodenrs[1, 2:end-1, 2:end-1] )))))'; # DCGH
d4 = vec(Int.(reshape( vcat(3.0*n4.-2, 3.0*n4.-1 ,3.0*n4),3*length(n4),1)));
n2 = setdiff(nodenrs[:],[n1[:];n3[:];n4[:]] );
d2 = vec(Int.(reshape( [3.0.*n2.-2 3*n2.-1 3.0.*n2],3*length(n2),1)));
for i = 1:6
epsilon = [e0[i,1] e0[i,4]/2 e0[i,6]/2;
e0[i,4]/2 e0[i,2] e0[i,5]/2;
e0[i,6]/2 e0[i,5]/2 e0[i,3]];
ufixed[:,i] = reshape(epsilon*vert_cor,24);
end
# 3D boundary constraint equations
wfixed = [repeat(ufixed[ 7:9,:],length((nodenrs[end,1,2:end-1] )),1); # C
repeat(ufixed[ 4:6,:]-ufixed[10:12,:],length((nodenrs[1, 1, 2:end-1] )),1); # B-D
repeat(ufixed[22:24,:],length((nodenrs[end,2:end-1,1] )),1); # H
repeat(ufixed[13:15,:]-ufixed[10:12,:],length((nodenrs[1, 2:end-1, 1] )),1); # E-D
repeat(ufixed[16:18,:],length((nodenrs[2:end-1, 1, 1] )),1); # F
repeat(ufixed[ 4:6,:]-ufixed[13:15,:],length((nodenrs[2:end-1,1,end] )),1); # B-E
repeat(ufixed[13:15,:],length((nodenrs[2:end-1, 2:end-1, 1] )),1); # E
repeat(ufixed[ 4:6,:],length((nodenrs[2:end-1, 1, 2:end-1] )),1); # B
repeat(ufixed[10:12,:],length((nodenrs[end,2:end-1,2:end-1] )),1)]; # D
## INITIALIZE ITERATION
# homogenization to evaluate macroscopic effective properties
qe = Array{Any,2}(undef, 6, 6);
Q = zeros(6,6);
dQ = Array{Any,2}(undef, 6, 6);
x = volfrac.*ones(nelz,nely,nelx)
for i = 1:nelx
for j = 1:nely
for k = 1:nelz
vall=3
if optType=="poisson"
vall=6
end
if sqrt((i-nelx/2-0.5)^2+(j-nely/2-0.5)^2+(k-nelz/2-0.5)^2) < min(min(nelx,nely),nelz)/vall
x[k,j,i] = volfrac/3.0;
end
end
end
end
xPhys = copy(x);
change = 1;
loop = 0;
xnew=zeros(size(x))
maxIter=50
## START ITERATION
while (change > 0.01 &&loop<maxIter)
loop = loop +1;
## FE-ANALYSIS
sK = reshape(KE[:]*(Emin .+xPhys[:]'.^penal*(1.0.-Emin)),24*24*nele,1);
K = sparse(iK[:], jK[:], sK[:] ); K .= (K.+K')./2;
Kr = vcat(hcat(K[d2,d2] , K[d2,d3]+K[d2,d4]),hcat((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\(-vcat(K[d2,d1], K[d3,d1]+K[d4,d1])*ufixed-vcat(K[d2,d4], K[d3,d4]+K[d4,d4])*wfixed);
U[d4,:] = U[d3,:] + wfixed;
# OBJECTIVE FUNCTION AND SENSITIVITY ANALYSIS
for i = 1:6
for j = 1:6
U1 = U[:,i]; U2 = U[:,j];
qe[i,j] = reshape(sum((U1[edofMat]*KE).*U2[edofMat],dims=2),nelz,nely,nelx);
Q[i,j] = 1.0./cellVolume.*sum(sum(sum((Emin .+ xPhys.^penal*(1 .-Emin)).*qe[i,j])));
dQ[i,j] = 1.0./cellVolume.*(penal*(1-Emin)*xPhys.^(penal-1).*qe[i,j]);
end
end
if optType=="bulk"
#bulk
# c = -(Q[1,1]+Q[2,2]+Q[1,2]+Q[2,1]);
# dc = -(dQ[1,1]+dQ[2,2]+dQ[1,2]+dQ[2,1]);
c = -(Q[1,1]+Q[2,2]+Q[3,3]+Q[1,2]+Q[2,3]+Q[1,3]);
dc = -(dQ[1,1]+dQ[2,2]+dQ[3,3]+dQ[1,2]+dQ[2,3]+dQ[1,3]);
elseif optType=="shear"
#shear
# c=-Q[3,3];
# dc=-dQ[3,3];
c=-(Q[4,4]+Q[5,5]+Q[6,6]);
dc=-(dQ[4,4]+dQ[5,5]+dQ[6,6]);
elseif optType=="poisson"
# c = Q[1,2]-(0.8^loop)*(Q[1,1]+Q[2,2]);
# dc = dQ[1,2]-(0.8^loop)*(dQ[1,1]+dQ[2,2]);
c = Q[1,2]+Q[2,3]+Q[1,3]-(0.8^loop)*(Q[1,1]+Q[2,2]+Q[3,3]);
dc = dQ[1,2]+dQ[2,3]+dQ[1,3]-(0.8^loop)*(dQ[1,1]+dQ[2,2]+dQ[3,3]);
end
dv = ones(nely,nelx,nelz);
## 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
if optType=="poisson"
l1=0; l2=1e9; move = 0.1;
else
l1 =0; l2 = 1e9; move = 0.2;
end
while (l2-l1 > 1e-9)
lmid = 0.5*(l2+l1);
if optType=="poisson"
xnew = max.(0,max.(x.-move,min.(1.0,min.(x.+move,x.*(-dc./dv./lmid)))));
else
xnew =max.(0.0,max.(x.-move,min.(1.0,min.(x.+move,x.*sqrt.(0.0.-dc./dv./lmid)))));
end
if ft == 1
xPhys = copy(xnew);
elseif ft == 2
xPhys[:] = (H*xnew[:])./Hs;
end
if mean(xPhys[:]) > volfrac
l1 = lmid;
else
l2 = lmid;
end
end
change = maximum(abs.(xnew[:].-x[:]))
x = xnew;
## PRINT RESULTS
println(" It:$loop Obj:$c Vol:$(mean(xPhys[:])) ch:$change ")
## PLOT DENSITIES
# heatmap(xPhys, aspect_ratio=:equal, legend=false, axis=nothing, foreground_color_subplot=colorant"white",fc=:grays,clims=(0.0, 1.0))
# volume(xPhys, algorithm = :iso, isorange = 0.2, isovalue = 0.9,colormap=:grays)
# frame(anim)
if (loop%10.0==0)
display(volume(xPhys, algorithm = :iso, isorange = 0.2, isovalue = 0.9,colormap=:grays))
# display(heatmap(1.0.-xPhys[Int(nelz/2),:,:], aspect_ratio=:equal, legend=false, axis=nothing, foreground_color_subplot=colorant"white",fc=:grays))
end
end
display(volume(xnew, algorithm = :iso, isorange = 0.2, isovalue = 0.9,colormap=:grays))
display(Q)
evaluateCH(Q,volfrac)
return xnew
end
##########################with MMA############################################
# 2D
function topXMMA(nelx,nely,volfrac,penal,rmin,ft,optType="bulk",maxEval=200)
## MATERIAL PROPERTIES
E0 = 1;
Emin = 1e-9;
nu = 0.3;
## PREPARE FINITE ELEMENT ANALYSIS
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])
nel=nely*nelx
nodenrs = reshape(1:(1+nelx)*(1+nely),1+nely,1+nelx)
edofVec = reshape(2*nodenrs[1:end-1,1:end-1].+1,nelx*nely,1)
edofMat = repeat(edofVec,1,8).+repeat([0 1 2*nely.+[2 3 0 1] -2 -1],nelx*nely,1)
iK = convert(Array{Int},reshape(kron(edofMat,ones(8,1))',64*nelx*nely,1))
jK = convert(Array{Int},reshape(kron(edofMat,ones(1,8))',64*nelx*nely,1))
## PREPARE FILTER
iH = ones(convert(Int,nelx*nely*(2*(ceil(rmin)-1)+1)^2),1)
jH = ones(Int,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(vec(iH),vec(jH),vec(sH))
Hs = sum(H,dims=2)
## PERIODIC BOUNDARY CONDITIONS
e0 = Matrix(1.0I, 3, 3);
ufixed = zeros(8,3);
U = zeros(2*(nely+1)*(nelx+1),3);
alldofs = [1:2*(nely+1)*(nelx+1)];
n1 = vcat(nodenrs[end,[1,end]],nodenrs[1,[end,1]]);
d1 = vec(reshape([(2 .* n1 .-1) 2 .*n1]',1,8));
n3 = [vec(nodenrs[2:(end-1),1]');vec(nodenrs[end,2:(end-1)])];
d3 = vec(reshape([(2 .*n3 .-1) 2 .*n3]',1,2*(nelx+nely-2)));
n4 = [vec(nodenrs[2:end-1,end]');vec(nodenrs[1,2:end-1])];
d4 = vec(reshape([(2 .*n4 .-1) 2 .*n4]',1,2*(nelx+nely-2)));
d2 = setdiff(vcat(alldofs...),union(union(d1,d3),d4));
for j = 1:3
ufixed[3:4,j] = [e0[1,j] e0[3,j]/2 ; e0[3,j]/2 e0[2,j]]*[nelx;0];
ufixed[7:8,j] = [e0[1,j] e0[3,j]/2 ; e0[3,j]/2 e0[2,j]]*[0;nely];
ufixed[5:6,j] = ufixed[3:4,j] .+ ufixed[7:8,j];
end
wfixed = [repeat(ufixed[3:4,:],nely-1,1); repeat(ufixed[7:8,:],nelx-1,1)];
## INITIALIZE ITERATION
qe = Array{Any,2}(undef, 3, 3);
Q = zeros(3,3);
dQ = Array{Any,2}(undef, 3, 3);
x = volfrac.*ones(nely,nelx)
for i = 1:nelx
for j = 1:nely
vall=3
if optType=="poisson"
vall=6
end
if sqrt((i-nelx/2-0.5)^2+(j-nely/2-0.5)^2) < min(nelx,nely)/vall
x[j,i] = volfrac/2.0;
end
end
end
xPhys = copy(x);
#display(heatmap(xPhys, aspect_ratio=:equal, legend=false, axis=nothing, foreground_color_subplot=colorant"white",fc=:grays))
loop=0
## START ITERATION
function FA(x1::Vector, grad::Vector)
xPhys=reshape(x1,nely,nelx)
#display(heatmap(xPhys, aspect_ratio=:equal, legend=false, axis=nothing, foreground_color_subplot=colorant"white",fc=:grays))
## FE-ANALYSIS
sK = reshape(KE[:]*(Emin .+ xPhys[:]'.^penal*(80 .- Emin)),64*nelx*nely,1);
K = sparse(vec(iK),vec(jK),vec(sK));
K = (K.+K')./2.0;
Kr = vcat(hcat(K[d2,d2] , K[d2,d3]+K[d2,d4]),hcat((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;
## OBJECTIVE FUNCTION AND SENSITIVITY ANALYSIS
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)./(nelx*nely);
Q[i,j] = sum(sum((Emin .+ xPhys.^penal*(E0 .-Emin)).*qe[i,j]));
dQ[i,j] = penal*(E0-Emin)*xPhys.^(penal-1).*qe[i,j];
end
end
c=0
if optType=="bulk"
#bulk
#c = -(Q[1,1]+Q[2,2]+Q[1,2]+Q[2,1]);
#dc = -(dQ[1,1]+dQ[2,2]+dQ[1,2]+dQ[2,1]);
c=-0.5.*(0.5.*(Q[1,1]+Q[2,2])+Q[1,2])
dc=-0.5.*(penal*xPhys.^(penal-1)).*(0.5.*(qe[1,1]+qe[2,2])+qe[1,2])
elseif optType=="shear"
#shear
c=-Q[3,3];
dc=-dQ[3,3];
elseif optType=="poisson"
# c = Q[1,2]-(0.5^loop)*(Q[1,1]+Q[2,2]);
# dc = dQ[1,2]-(0.5^loop)*(dQ[1,1]+dQ[2,2]);
# c = Q[1,2]-(Q[1,1]+Q[2,2]);
# dc = dQ[1,2]-(dQ[1,1]+dQ[2,2]);
# c= (2.0 .*Q[1,2])./(Q[1,1]+Q[2,2])
# dc=2.0 .*(penal*xPhys.^(penal-1)).*(qe[1,2]./(Q[1,1]+Q[2,2]) - Q[1,2] ./ ((Q[1,1]+Q[2,2]).^2) .*(qe[1,1]+qe[2,2]))
c= (5.0 .*Q[1,2])./(Q[1,1]+Q[2,2])
dc= 5.0 .*(penal*xPhys.^(penal-1)).*(qe[1,2]./(Q[1,1]+Q[2,2]) - Q[1,2] ./ ((Q[1,1]+Q[2,2]).^2) .*(qe[1,1]+qe[2,2]))
# c = 2.0*Q[1,2]-(Q[1,1]+Q[2,2]);
# dc = 2.0*dQ[1,2]-(dQ[1,1]+dQ[2,2]);
# c= 2.0*(Q[1,2])./(Q[1,1]+Q[2,2])
# dc=2.0.*((dQ[1,2]./(Q[1,1]+Q[2,2])) .- (Q[1,2] ./ ((Q[1,1]+Q[2,2]).^(2)) .*(dQ[1,1]+dQ[2,2])));
# c= (Q[1,2])./(Q[1,1])
# dc=((dQ[1,2]./(Q[1,1])) .- (Q[1,2] ./ ((Q[1,1]).^(2)) .*(dQ[1,1])))
# c= (Q[1,2])./(Q[1,1]+Q[2,2])
# dc=(Q[1,2]./(dQ[1,1]+dQ[2,2])- dQ[1,2] ./ ((dQ[1,1]+dQ[2,2]).^2) .*(Q[1,1]+Q[2,2]))
end
## FILTERING/MODIFICATION OF SENSITIVITIES
if ft == 1
dc[:] = H*(x[:].*dc[:])./Hs./max.(1e-3,x[:]);
elseif ft == 2
dc[:] = H*(dc[:]./Hs);
end
# display(dc)
# display(qe)
grad[:] .= dc[:]
loop+=1
return c
end
function G(x1::Vector, grad::Vector)
dv = ones(nely,nelx)
if ft == 2
dv[:] = H*(dv[:]./Hs)
end
grad[:] .= dv[:]
return (sum(x1) - volfrac*nel)
end
tol=1e-6
# if optType=="poisson"
# tol=1e-9
# end
opt = Opt(:LD_MMA, nel)
opt.lower_bounds = fill(1e-9,nel)
opt.upper_bounds = fill(1,nel)
opt.xtol_rel = tol
opt.maxeval = maxEval
opt.min_objective = FA
inequality_constraint!(opt, (x,gg) -> G(x,gg), tol)
display(@time (minf,minx,ret) = optimize(opt, xPhys[:]))
numevals = opt.numevals # the number of function evaluations
display("got $minf after $numevals iterations (returned $ret)")
xPhys=reshape(minx,nely,nelx)
display(heatmap(xPhys, aspect_ratio=:equal, legend=false, axis=nothing, foreground_color_subplot=colorant"white",fc=:grays))
display(Q)
return xPhys
end
# 3D
function topXMMA3D(nelx,nely,nelz,volfrac,penal,rmin,ft,optType="bulk",maxEval=200)
nel=nelx*nely*nelz
# MATERIAL PROPERTIES
E0 = 1;
Emin = 1e-9;
nu = 0.3;
lx=0.1;ly=0.1;lz=0.1;
vert_cor = [0 lx lx 0 0 lx lx 0;
0 0 ly ly 0 0 ly ly;
0 0 0 0 lz lz lz lz];
cellVolume = lx*ly*lz;
# PREPARE FINITE ELEMENT ANALYSIS
# KE=lk_H8(nu)
# the initial definitions of the PUC
D0 = E0/(1+nu)/(1-2*nu)*
[ 1-nu nu nu 0 0 0 ;
nu 1-nu nu 0 0 0 ;
nu nu 1-nu 0 0 0 ;
0 0 0 (1-2*nu)/2 0 0 ;
0 0 0 0 (1-2*nu)/2 0 ;
0 0 0 0 0 (1-2*nu)/2];
dx = lx/nelx; dy = ly/nely; dz = lz/nelz;
KE = elementMatVec3D(dx/2, dy/2, dz/2, D0);
Num_node = (1+nely)*(1+nelx)*(1+nelz);
nele = nelx*nely*nelz;
nodenrs = reshape(1:Num_node,1+nely,1+nelx,1+nelz);
edofVec = reshape(3*nodenrs[1:end-1,1:end-1,1:end-1] .+1,nelx*nely*nelz,1);
edofMat = repeat(edofVec,1,24).+repeat([0 1 2 3*nely.+[3 4 5 0 1 2] -3 -2 -1 3*(nelx+1)*(nely+1).+[0 1 2 3*nely.+[3 4 5 0 1 2] -3 -2 -1]], nelx*nely*nelz, 1);
# nodenrs = reshape(1:(1+nelx)*(1+nely)*(1+nelz),1+nely,1+nelx,1+nelz)
# edofVec = reshape(3*nodenrs[1:end-1,1:end-1,1:end-1].+1,nelx*nely*nelz,1)
# edofMat = repeat(edofVec,1,24).+repeat([0 1 2 3*nely.+[3 4 5 0 1 2] -3 -2 -1 3*(nely+1)*(nelx+1).+[0 1 2 3*nely.+[3 4 5 0 1 2] -3 -2 -1]],nelx*nely*nelz,1)
iK = reshape(kron(edofMat,ones(24,1))',24*24*nele,1);
jK = reshape(kron(edofMat,ones(1,24))',24*24*nele,1);
## PREPARE FILTER
H,Hs=make_filter3D(nelx,nely,nelz,rmin);
## PERIODIC BOUNDARY CONDITIONS
e0 = Matrix(1.0I, 6, 6);
ufixed = zeros(24,6);
ndof = 3*(nelx+1)*(nely+1)*(nelz+1);
U = zeros(ndof,6);
alldofs = [1:ndof];
# 3D periodic boundary formulation
# the nodes classification
n1 = hcat(nodenrs[end, [1 end], 1], nodenrs[1, [end 1], 1], nodenrs[end, [1 end], end] ,nodenrs[1, [end 1], end]);
d1 = Int.(vec(reshape(vcat(3.0.*n1.-2, 3.0.*n1.-1,3.0.*n1),3*length(n1),1)));
n3 = vcat(vec(reshape(nodenrs[end,1,2:end-1] ,1,length(nodenrs[end,1,2:end-1] ))), # AE
vec(reshape(nodenrs[1, 1, 2:end-1] ,1,length(nodenrs[1, 1, 2:end-1] ))), # DH
vec(reshape(nodenrs[end,2:end-1,1] ,1,length(nodenrs[end,2:end-1,1] ))), # AB
vec(reshape(nodenrs[1, 2:end-1, 1] ,1,length(nodenrs[1, 2:end-1, 1] ))), # DC
vec(reshape(nodenrs[2:end-1, 1, 1] ,1,length(nodenrs[2:end-1, 1, 1] ))), # AD
vec(reshape(nodenrs[2:end-1,1,end] ,1,length(nodenrs[2:end-1,1,end] ))), # EH
vec(reshape(nodenrs[2:end-1, 2:end-1, 1] ,1,length(nodenrs[2:end-1, 2:end-1, 1] ))), # ABCD
vec(reshape(nodenrs[2:end-1, 1, 2:end-1] ,1,length(nodenrs[2:end-1, 1, 2:end-1] ))), # ADHE
vec(reshape(nodenrs[end,2:end-1,2:end-1] ,1,length(nodenrs[end,2:end-1,2:end-1] ))))'; # ABFE
d3 = vec(Int.(reshape( vcat(3.0.*n3.-2 ,3.0.*n3.-1, 3.0.*n3),3*length(n3),1)));
n4 = vcat(vec(reshape(nodenrs[1, end, 2:end-1] ,1,length((nodenrs[1, end, 2:end-1] )))), # CG
vec(reshape(nodenrs[end,end,2:end-1] ,1,length((nodenrs[end,end,2:end-1] )))), # BF
vec(reshape(nodenrs[1, 2:end-1, end] ,1,length((nodenrs[1, 2:end-1, end] )))), # HG
vec(reshape(nodenrs[end,2:end-1,end] ,1,length((nodenrs[end,2:end-1,end] )))), # EF
vec(reshape(nodenrs[2:end-1,end,end] ,1,length((nodenrs[2:end-1,end,end] )))), # FG
vec(reshape(nodenrs[2:end-1, end, 1] ,1,length((nodenrs[2:end-1, end, 1] )))), # BC
vec(reshape(nodenrs[2:end-1,2:end-1,end] ,1,length((nodenrs[2:end-1,2:end-1,end] )))), # EFGH
vec(reshape(nodenrs[2:end-1,end,2:end-1] ,1,length((nodenrs[2:end-1,end,2:end-1] )))), # BCGF
vec(reshape(nodenrs[1, 2:end-1, 2:end-1] ,1,length((nodenrs[1, 2:end-1, 2:end-1] )))))'; # DCGH
d4 = vec(Int.(reshape( vcat(3.0*n4.-2, 3.0*n4.-1 ,3.0*n4),3*length(n4),1)));
n2 = setdiff(nodenrs[:],[n1[:];n3[:];n4[:]] );
d2 = vec(Int.(reshape( [3.0.*n2.-2 3*n2.-1 3.0.*n2],3*length(n2),1)));
for i = 1:6
epsilon = [e0[i,1] e0[i,4]/2 e0[i,6]/2;
e0[i,4]/2 e0[i,2] e0[i,5]/2;
e0[i,6]/2 e0[i,5]/2 e0[i,3]];
ufixed[:,i] = reshape(epsilon*vert_cor,24);
end
# 3D boundary constraint equations
wfixed = [repeat(ufixed[ 7:9,:],length((nodenrs[end,1,2:end-1] )),1); # C
repeat(ufixed[ 4:6,:]-ufixed[10:12,:],length((nodenrs[1, 1, 2:end-1] )),1); # B-D
repeat(ufixed[22:24,:],length((nodenrs[end,2:end-1,1] )),1); # H
repeat(ufixed[13:15,:]-ufixed[10:12,:],length((nodenrs[1, 2:end-1, 1] )),1); # E-D
repeat(ufixed[16:18,:],length((nodenrs[2:end-1, 1, 1] )),1); # F
repeat(ufixed[ 4:6,:]-ufixed[13:15,:],length((nodenrs[2:end-1,1,end] )),1); # B-E
repeat(ufixed[13:15,:],length((nodenrs[2:end-1, 2:end-1, 1] )),1); # E
repeat(ufixed[ 4:6,:],length((nodenrs[2:end-1, 1, 2:end-1] )),1); # B
repeat(ufixed[10:12,:],length((nodenrs[end,2:end-1,2:end-1] )),1)]; # D
## INITIALIZE ITERATION
# homogenization to evaluate macroscopic effective properties
qe = Array{Any,2}(undef, 6, 6);
Q = zeros(6,6);
dQ = Array{Any,2}(undef, 6, 6);
x = volfrac.*ones(nelz,nely,nelx)
for i = 1:nelx
for j = 1:nely
for k = 1:nelz
vall=3
if optType=="poisson"
vall=6
end
if sqrt((i-nelx/2-0.5)^2+(j-nely/2-0.5)^2+(k-nelz/2-0.5)^2) < min(min(nelx,nely),nelz)/vall
x[k,j,i] = volfrac/3.0;
end
end
end
end
# for i = 1:nelx
# for j = 1:nely
# for k = 1:nelz
# x[k,j,i] = volfrac/2.0;
# end
# end
# end
# xx = ones(nelz,nelx,nelx);
# xx[Int(nelz/2):Int(nelz/2+1),Int(nely/2):Int(nely/2+1),Int(nelx/2):Int(nelx/2+1)] .= Emin;
# beta = 1.0;
# xx = (1.0.-exp.(-beta.*xx).+xx.*exp.(-beta));
# x=xx.*volfrac
# x[x.<Emin].=Emin
xPhys = copy(x);
#display(heatmap(xPhys, aspect_ratio=:equal, legend=false, axis=nothing, foreground_color_subplot=colorant"white",fc=:grays))
loop=0
function FEA(xPhys)
# FE-ANALYSIS
sK = reshape(KE[:]*(Emin .+xPhys[:]'.^penal*(1.0.-Emin)),24*24*nele,1);
K = sparse(iK[:], jK[:], sK[:] ); K .= (K.+K')./2;
Kr = vcat(hcat(K[d2,d2] , K[d2,d3]+K[d2,d4]),hcat((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\(-vcat(K[d2,d1], K[d3,d1]+K[d4,d1])*ufixed-vcat(K[d2,d4], K[d3,d4]+K[d4,d4])*wfixed);
U[d4,:] = U[d3,:] + wfixed;
# OBJECTIVE FUNCTION AND SENSITIVITY ANALYSIS
for i = 1:6
for j = 1:6
U1 = U[:,i]; U2 = U[:,j];
qe[i,j] = reshape(sum((U1[edofMat]*KE).*U2[edofMat],dims=2),nelz,nely,nelx);
Q[i,j] = 1.0./cellVolume.*sum(sum(sum((Emin .+ xPhys.^penal*(1 .-Emin)).*qe[i,j])));
dQ[i,j] = 1.0./cellVolume.*(penal*(1-Emin)*xPhys.^(penal-1).*qe[i,j]);
end
end
#display(Q)
end
# FEA(xx)
# display(Q)
## START ITERATION
function F(x1::Vector, grad::Vector)
xPhys=reshape(x1,nelz,nely,nelx)
#display(heatmap(xPhys, aspect_ratio=:equal, legend=false, axis=nothing, foreground_color_subplot=colorant"white",fc=:grays))
FEA(xPhys)
c=0
if optType=="bulk"
#bulk
#c = -(Q[1,1]+Q[2,2]+Q[1,2]+Q[2,1]);
#dc = -(dQ[1,1]+dQ[2,2]+dQ[1,2]+dQ[2,1]);
c = -(Q[1,1]+Q[2,2]+Q[3,3]+Q[1,2]+Q[2,3]+Q[1,3]);
dc = -(dQ[1,1]+dQ[2,2]+dQ[3,3]+dQ[1,2]+dQ[2,3]+dQ[1,3]);
#c=-0.5.*(0.5.*(Q[1,1]+Q[2,2])+Q[1,2])
#dc=-0.5.*(penal*xPhys.^(penal-1)).*(0.5.*(qe[1,1]+qe[2,2])+qe[1,2])
elseif optType=="shear"
#shear
c=-(Q[4,4]+Q[5,5]+Q[6,6]);
dc=-(dQ[4,4]+dQ[5,5]+dQ[6,6]);
elseif optType=="poisson"
# c = Q[1,2]-(0.8^loop)*(Q[1,1]+Q[2,2]);
# dc = dQ[1,2]-(0.8^loop)*(dQ[1,1]+dQ[2,2]);
# c = (Q[1,2]+Q[2,3]+Q[1,3])-(Q[1,1]+Q[2,2]+Q[3,3]);
# dc = (dQ[1,2]+dQ[2,3]+dQ[1,3])-(dQ[1,1]+dQ[2,2]+dQ[3,3]);
c = (Q[1,2]+Q[2,3]+Q[1,3]) .-(Q[1,1]+Q[2,2]+Q[3,3]);
dc = (dQ[1,2]+dQ[2,3]+dQ[1,3]) .-(dQ[1,1]+dQ[2,2]+dQ[3,3]);
# c= (Q[1,2]+Q[2,3]+Q[1,3])./(Q[1,1]+Q[2,2]+Q[3,3])
# dc=(( (dQ[1,2]+dQ[2,3]+dQ[1,3])./(Q[1,1]+Q[2,2]+Q[3,3])).- ((Q[1,2]+Q[2,3]+Q[1,3]) ./ ((Q[1,1]+Q[2,2]+Q[3,3]).^(2)) .*(dQ[1,1]+dQ[2,2]+dQ[3,3])))
# c= 2.0*(Q[1,2])./(Q[1,1]+Q[2,2])
# dc=2.0.*((dQ[1,2]./(Q[1,1]+Q[2,2])).- (Q[1,2] ./ ((Q[1,1]+Q[2,2]).^(2)) .*(dQ[1,1]+dQ[2,2])))
end
## FILTERING/MODIFICATION OF SENSITIVITIES
if ft == 1
dc[:] = H*(x[:].*dc[:])./Hs./max.(1e-3,x[:]);
elseif ft == 2
dc[:] = H*(dc[:]./Hs);
end
#display(dc)
#display(qe)
println("$(loop) Objective: $(-c) ")
grad[:] .= dc[:]
loop+=1
if (loop%5.0==0)
display(volume(xPhys, algorithm = :iso, isorange = 0.2, isovalue = 0.9,colormap=:grays))
# display(heatmap(1.0.-xPhys[Int(nelz/2),:,:], aspect_ratio=:equal, legend=false, axis=nothing, foreground_color_subplot=colorant"white",fc=:grays))
end
return c
end
function G(x1::Vector, grad::Vector)
dv = ones(nelz,nely,nelx)
if ft == 2
dv[:] = H*(dv[:]./Hs)
end
grad[:] .= dv[:]
return (sum(x1) - volfrac*nel)
end
opt = Opt(:LD_MMA, nel)
opt.lower_bounds = fill(1e-9,nel)
opt.upper_bounds = fill(1,nel)
opt.xtol_rel = 1e-6
opt.maxeval = maxEval
opt.min_objective = F
inequality_constraint!(opt, (x,gg) -> G(x,gg), 1e-6)
display(@time (minf,minx,ret) = optimize(opt, xPhys[:]))
numevals = opt.numevals # the number of function evaluations
display("got $minf after $numevals iterations (returned $ret)")
display(volume(xPhys, algorithm = :iso, isorange = 0.2, isovalue = 0.9,colormap=:grays))
# display(heatmap(1.0.-xPhys[Int(nelz/2),:,:], aspect_ratio=:equal, legend=false, axis=nothing, foreground_color_subplot=colorant"white",fc=:grays))
# display(heatmap(1.0.-xPhys[:,Int(nelz/2),:], aspect_ratio=:equal, legend=false, axis=nothing, foreground_color_subplot=colorant"white",fc=:grays))
# display(heatmap(1.0.-xPhys[:,:,Int(nelz/2)], aspect_ratio=:equal, legend=false, axis=nothing, foreground_color_subplot=colorant"white",fc=:grays))
volume
display(Q)
evaluateCH(Q,volfrac)
return xPhys,Q
end
######################other objective functions###############
function E_MMA(nelx,nely,volfrac,penal,rmin,ft,CStar,optType="bulk",maxEval=200)
## MATERIAL PROPERTIES
E0 = 1;
Emin = 1e-9;
nu = 0.3;
## PREPARE FINITE ELEMENT ANALYSIS
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])
nel=nely*nelx
nodenrs = reshape(1:(1+nelx)*(1+nely),1+nely,1+nelx)
edofVec = reshape(2*nodenrs[1:end-1,1:end-1].+1,nelx*nely,1)
edofMat = repeat(edofVec,1,8).+repeat([0 1 2*nely.+[2 3 0 1] -2 -1],nelx*nely,1)
iK = convert(Array{Int},reshape(kron(edofMat,ones(8,1))',64*nelx*nely,1))
jK = convert(Array{Int},reshape(kron(edofMat,ones(1,8))',64*nelx*nely,1))
## PREPARE FILTER
iH = ones(convert(Int,nelx*nely*(2*(ceil(rmin)-1)+1)^2),1)
jH = ones(Int,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(vec(iH),vec(jH),vec(sH))
Hs = sum(H,dims=2)
## PERIODIC BOUNDARY CONDITIONS
e0 = Matrix(1.0I, 3, 3);
ufixed = zeros(8,3);
U = zeros(2*(nely+1)*(nelx+1),3);
alldofs = [1:2*(nely+1)*(nelx+1)];
n1 = vcat(nodenrs[end,[1,end]],nodenrs[1,[end,1]]);
d1 = vec(reshape([(2 .* n1 .-1) 2 .*n1]',1,8));
n3 = [vec(nodenrs[2:(end-1),1]');vec(nodenrs[end,2:(end-1)])];
d3 = vec(reshape([(2 .*n3 .-1) 2 .*n3]',1,2*(nelx+nely-2)));
n4 = [vec(nodenrs[2:end-1,end]');vec(nodenrs[1,2:end-1])];
d4 = vec(reshape([(2 .*n4 .-1) 2 .*n4]',1,2*(nelx+nely-2)));
d2 = setdiff(vcat(alldofs...),union(union(d1,d3),d4));
for j = 1:3
ufixed[3:4,j] = [e0[1,j] e0[3,j]/2 ; e0[3,j]/2 e0[2,j]]*[nelx;0];
ufixed[7:8,j] = [e0[1,j] e0[3,j]/2 ; e0[3,j]/2 e0[2,j]]*[0;nely];
ufixed[5:6,j] = ufixed[3:4,j] .+ ufixed[7:8,j];
end
wfixed = [repeat(ufixed[3:4,:],nely-1,1); repeat(ufixed[7:8,:],nelx-1,1)];
## INITIALIZE ITERATION
qe = Array{Any,2}(undef, 3, 3);
Q = zeros(3,3);
dQ = Array{Any,2}(undef, 3, 3);
x = volfrac.*ones(nely,nelx)
for i = 1:nelx
for j = 1:nely
vall=3
if optType=="poisson"
vall=6
end
if sqrt((i-nelx/2-0.5)^2+(j-nely/2-0.5)^2) < min(nelx,nely)/vall
x[j,i] = volfrac/2.0;
end
end
end
xPhys = copy(x);
#display(heatmap(xPhys, aspect_ratio=:equal, legend=false, axis=nothing, foreground_color_subplot=colorant"white",fc=:grays))
loop=0
function FE(xPhys)
## FE-ANALYSIS
sK = reshape(KE[:]*(Emin .+ xPhys[:]'.^penal*(80 .- Emin)),64*nelx*nely,1);
K = sparse(vec(iK),vec(jK),vec(sK));
K = (K.+K')./2.0;
Kr = vcat(hcat(K[d2,d2] , K[d2,d3]+K[d2,d4]),hcat((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;
## OBJECTIVE FUNCTION AND SENSITIVITY ANALYSIS
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)./(nelx*nely);
Q[i,j] = sum(sum((Emin .+ xPhys.^penal*(E0 .-Emin)).*qe[i,j]));
dQ[i,j] = penal*(E0-Emin)*xPhys.^(penal-1).*qe[i,j];
end
end
end
# function minimize volume
function F(x1::Vector, grad::Vector)
xPhys=reshape(x1,nely,nelx)
FE(xPhys)
dv = ones(nely,nelx)
if ft == 2
dv[:] = H*(dv[:]./Hs)
end
grad[:] .= dv[:]
#return (sum(x1) - volfrac*nel)
return (sum(x1))
end
function G33(x1::Vector, grad::Vector)
#shear
c=-Q[3,3];
dc=-dQ[3,3];
## FILTERING/MODIFICATION OF SENSITIVITIES
if ft == 1
dc[:] = H*(x[:].*dc[:])./Hs./max.(1e-3,x[:]);
elseif ft == 2
dc[:] = H*(dc[:]./Hs);
end
grad[:] .= dc[:]
return c +CStar[3,3]
end
function G11(x1::Vector, grad::Vector)
#shear
c=-(Q[1,1]);
dc=-(dQ[1,1]);
## FILTERING/MODIFICATION OF SENSITIVITIES
if ft == 1
dc[:] = H*(x[:].*dc[:])./Hs./max.(1e-3,x[:]);
elseif ft == 2
dc[:] = H*(dc[:]./Hs);
end
grad[:] .= dc[:]
return c +CStar[1,1]
end
function G12(x1::Vector, grad::Vector)
#shear
c=-(Q[1,2]);
dc=-(dQ[1,2]);
## FILTERING/MODIFICATION OF SENSITIVITIES
if ft == 1
dc[:] = H*(x[:].*dc[:])./Hs./max.(1e-3,x[:]);
elseif ft == 2
dc[:] = H*(dc[:]./Hs);
end
grad[:] .= dc[:]
return c +CStar[1,2]
end
opt = Opt(:LD_MMA, nel)
opt.lower_bounds = fill(1e-9,nel)
opt.upper_bounds = fill(1,nel)
opt.xtol_rel = 1e-6
opt.maxeval = maxEval
opt.min_objective = F
inequality_constraint!(opt, (x,gg) -> G33(x,gg), 1e-3)
# inequality_constraint!(opt, (x,gg) -> G11(x,gg), 1e-3)
inequality_constraint!(opt, (x,gg) -> G12(x,gg), 1e-6)
display(@time (minf,minx,ret) = optimize(opt, xPhys[:]))
numevals = opt.numevals # the number of function evaluations
display("got $minf after $numevals iterations (returned $ret)")
xPhys=reshape(minx,nely,nelx)
display(heatmap(xPhys, aspect_ratio=:equal, legend=false, axis=nothing, foreground_color_subplot=colorant"white",fc=:grays))
display(Q)
return xPhys
end