-
Amira Abdel-Rahman authoredAmira Abdel-Rahman authored
multimaterial.jl 36.46 KiB
# Amira Abdel-Rahman
# (c) Massachusetts Institute of Technology 2020
#############################################multimaterial############################################################
#Based on: https://link.springer.com/article/10.1007/s00158-013-0999-1
function multitop(nx,ny,tol_out,tol_f,iter_max_in,iter_max_out,p,q,e,v,rf)
anim=Animation()
alpha = zeros(Float64,nx*ny,p)
for i = 1:p
alpha[:,i] .= v[i]
end
# MAKE FILTER
H,Hs = make_filter(nx,ny,rf)
obj=0
change_out = 2*tol_out; iter_out = 0;
while (iter_out < iter_max_out) && (change_out > tol_out)
alpha_old = copy(alpha)
for a = 1:p
for b = a+1:p
obj,alpha = bi_top(a,b,nx,ny,p,v,e,q,alpha,H,Hs,iter_max_in);
end
end
iter_out = iter_out + 1;
change_out = norm(alpha[:].-alpha_old[:],Inf);
display("Iter:$iter_out Obj.:$obj change:$change_out");
# UPDATE FILTER
if (change_out < tol_f) && (rf>3)
tol_f = 0.99*tol_f; rf = 0.99*rf;
H,Hs = make_filter(nx,ny,rf);
end
# SCREEN OUT TEMPORAL TOPOLOGY EVERY 5 ITERATIONS
if mod(iter_out,5)==0
I = make_bitmap(p,nx,ny,alpha);
display(RGB.(I[:,:,1],I[:,:,2],I[:,:,3]))
I = permutedims(I, [3, 1, 2])
img = colorview(RGB, I)
heatmap(img,xaxis=nothing,yaxis=nothing,aspect_ratio=:equal)
frame(anim)
end
end
return anim,alpha
end
function bi_top(a,b,nx,ny,p,v,e,q,alpha_old,H,Hs,maxIter)
alpha = copy(alpha_old)
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]);
nelx=nx
nely=ny
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))
# DEFINE LOADS AND SUPPORTS (HALF MBB-BEAM)
F = sparse([2],[1],[-1.0],2*(nely+1)*(nelx+1),1)
U = zeros(2*(nely+1)*(nelx+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)
o=0;alpha_a=0
# inner iteration
for i =0: (maxIter-1)
# FE-ANALYSIS
E = e[1]*alpha[:,1].^q
for phase = 2:p
E = E + e[phase]*alpha[:,phase].^q;
end
sK = reshape(KE[:]*E[:]',64*nelx*nely,1)
K = sparse(vec(iK),vec(jK),vec(sK)); K = (K+K')/2
@timed U[freedofs] = K[freedofs,freedofs] \ Array(F[freedofs])
# Objective function and sensitivity analysis
ce = sum((U[edofMat]*KE).*U[edofMat],dims=2)
o = sum(sum(E.*ce))
# FILTERING OF SENSITIVITIES
dc = 0.0 .-(q .*(e[a].-e[b]).*alpha[:,a].^(q-1)).*ce;
dc = H*(alpha[:,a].*dc)./Hs./max.(1e-3,alpha[:,a])
dc = min.(dc,0.0)
# UPDATE LOWER AND UPPER BOUNDS OF DESIGN VARIABLES
move = 0.2;
r = ones(nx*ny,1)
for k = 1:p
#if (k ~= a) && (k ~= b)
if (k != a) && (k != b)
# display("k $k a $a b $b")
r = r .- alpha[:,k]
end
end
l = max.(0,alpha[:,a] .-move)
u = min.(r,alpha[:,a] .+move)
# OPTIMALITY CRITERIA UPDATE OF DESIGN VARIABLES
l1 = 0; l2 = 1e9
while (l2-l1)/(l1+l2) > 1e-3
lmid = 0.5*(l2+l1)
alpha_a = max.(l,min.(u,alpha[:,a].*sqrt.(-dc./lmid)));
if sum(alpha_a) > nx*ny*v[a]
l1 = lmid
else
l2 = lmid
end
end
alpha[:,a] = alpha_a
alpha[:,b] = r .-alpha_a
end
return o,alpha
end
########
function multitop_compliant(nx,ny,tol_out,tol_f,iter_max_in,iter_max_out,p,q,e,v,rf,Load,Support,Spring,DOut)
anim=Animation()
alpha = zeros(Float64,nx*ny,p)
for i = 1:p
alpha[:,i] .= v[i]
end
# MAKE FILTER
H,Hs = make_filter(nx,ny,rf)
obj=0
change_out = 2*tol_out; iter_out = 0; while (iter_out < iter_max_out) && (change_out > tol_out)
alpha_old = copy(alpha)
for a = 1:p
for b = a+1:p
obj,alpha = bi_top_compliant(a,b,nx,ny,p,v,e,q,alpha,H,Hs,iter_max_in);
end
end
iter_out = iter_out + 1;
change_out = norm(alpha[:].-alpha_old[:],Inf);
display("Iter:$iter_out Obj.:$obj change:$change_out");
# UPDATE FILTER
if (change_out < tol_f) && (rf>3)
tol_f = 0.99*tol_f; rf = 0.99*rf;
H,Hs = make_filter(nx,ny,rf);
end
# SCREEN OUT TEMPORAL TOPOLOGY EVERY 5 ITERATIONS
if mod(iter_out,5)==0
I = make_bitmap_compliant(p,nx,ny,alpha);
display(RGB.(I[:,:,1],I[:,:,2],I[:,:,3]))
I = permutedims(I, [3, 1, 2])
img = colorview(RGB, I)
heatmap(img,xaxis=nothing,yaxis=nothing,aspect_ratio=:equal)
frame(anim)
end
end
return anim,alpha
end
function bi_top_compliant(a,b,nx,ny,p,v,e,q,alpha_old,H,Hs,maxIter,Load,Support,Spring,DOut)
alpha = copy(alpha_old)
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]);
nelx=nx
nely=ny
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))
# # DEFINE LOADS AND SUPPORTS (HALF MBB-BEAM)
# F = sparse([2],[1],[-1.0],2*(nely+1)*(nelx+1),1)
# U = zeros(2*(nely+1)*(nelx+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)
# DEFINE LOADS AND SUPPORTS
F = sparse(2 .*Load[:,1] .-2 .+ Load[:,2],ones(size(Load,1)),Load[:,3],2*(nely+1)*(nelx+1),2)
DofDOut = 2 * DOut[1] - 2 +DOut[2]; #only one
F=Array(F)
F[DofDOut,2]=-1
fixeddofs = 2 .*Support[:,1] .-2 .+ Support[:,2]
NSpring = size(Spring,1);
s = sparse(2 .*Spring[:,1] .-2 .+ Spring[:,2],ones(size(Spring,1)),Spring[:,3],2*(nely+1)*(nelx+1),2)
m=Array(s)[:,1]
S= sparse(diagm(m))
U = zeros(2*(nely+1)*(nelx+1),2)
alldofs = 1:(2*(nely+1)*(nelx+1))
freedofs = setdiff(alldofs,fixeddofs)
o=0;alpha_a=0
# inner iteration
for i =0: (maxIter-1)
# FE-ANALYSIS
E = e[1]*alpha[:,1].^q
for phase = 2:p
E = E + e[phase]*alpha[:,phase].^q;
end
sK = reshape(KE[:]*E[:]',64*nelx*nely,1)
K = sparse(vec(iK),vec(jK),vec(sK));
K[din,din] = K[din,din] + 0.1;
K[dout,dout] = K[dout,dout] + 0.1;
K = (K+K')/2
# @timed U[freedofs] = K[freedofs,freedofs] \ Array(F[freedofs])
@timed U[freedofs,:] = K[freedofs,freedofs] \ Array(F[freedofs,:])
U1=U[:,1]
U2=U[:,2]
# Objective function and sensitivity analysis
# ce = sum((U[edofMat]*KE).*U[edofMat],dims=2)
ce =-sum((U1[edofMat]*KE).*U2[edofMat],dims=2)
# o = sum(sum(E.*ce))
o=U[DofDOut,1]
# FILTERING OF SENSITIVITIES
dc = 0.0 .-(q .*(e[a].-e[b]).*alpha[:,a].^(q-1)).*ce;
dc = H*(alpha[:,a].*dc)./Hs./max.(1e-3,alpha[:,a])
dc = min.(dc,0.0)
# UPDATE LOWER AND UPPER BOUNDS OF DESIGN VARIABLES
move=0.1;
r = ones(nx*ny,1)
for k = 1:p
#if (k ~= a) && (k ~= b)
if (k != a) && (k != b)
# display("k $k a $a b $b")
r = r .- alpha[:,k]
end
end
l = max.(0,alpha[:,a] .-move)
u = min.(r,alpha[:,a] .+move)
# OPTIMALITY CRITERIA UPDATE OF DESIGN VARIABLES
l1 = 0; l2 = 1e9
while (l2-l1)/(l2+l1) > 1e-4 && l2>1e-40
# while (l2-l1)/(l1+l2) > 1e-3
lmid = 0.5*(l2+l1)
alpha_a = max.(l,min.(u,alpha[:,a].*((max.(1e-10,-dc./lmid)).^0.3 )));
if sum(alpha_a) > nx*ny*v[a]
l1 = lmid
else
l2 = lmid
end
end
alpha[:,a] = alpha_a
alpha[:,b] = r .-alpha_a
end
return o,alpha
end
########
function multitop3D(nx,ny,nz,tol_out,tol_f,iter_max_in,iter_max_out,p,q,e,v,rf)
anim=Animation()
alpha = zeros(Float64,nx*ny*nz,p)
for i = 1:p
alpha[:,i] .= v[i]
end
# MAKE FILTER
H,Hs = make_filter3D(nx,ny,nz,rf)
obj=0
change_out = 2*tol_out; iter_out = 0;
while (iter_out < iter_max_out) && (change_out > tol_out)
alpha_old = copy(alpha)
for a = 1:p
for b = a+1:p
obj,alpha = bi_top3D(a,b,nx,ny,nz,p,v,e,q,alpha,H,Hs,iter_max_in);
end
end
iter_out = iter_out + 1;
change_out = norm(alpha[:].-alpha_old[:],Inf);
display("Iter:$iter_out Obj.:$obj change:$change_out");
# UPDATE FILTER
if (change_out < tol_f) && (rf>3)
tol_f = 0.99*tol_f; rf = 0.99*rf;
H,Hs = make_filter3D(nx,ny,nz,rf)
end
# SCREEN OUT TEMPORAL TOPOLOGY EVERY 5 ITERATIONS
if mod(iter_out,5)==0
I = make_bitmap_3d(p,nx,ny,nz,alpha)
threshold=0.05
display(topplotmulti3d(nx,ny,nz,I,threshold))
frame(anim)
end
end
return anim,alpha
end
function bi_top3D(a,b,nx,ny,nz,p,v,e,q,alpha_old,H,Hs,maxIter)
alpha = copy(alpha_old)
nu = 0.3;
# PREPARE FINITE ELEMENT ANALYSIS
KE=lk_H8(nu)
nelx=nx
nely=ny
nelz=nz
nele = nelx*nely*nelz
ndof = 3*(nelx+1)*(nely+1)*(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 = convert(Array{Int},reshape(kron(edofMat,ones(24,1))',24*24*nele,1))
jK = convert(Array{Int},reshape(kron(edofMat,ones(1,24))',24*24*nele,1))
# USER-DEFINED LOAD DOFs
m= Matlab.meshgrid(nelx:nelx, 0:0, 0:nelz)
il=m[1]
jl=m[2]
kl=m[3]
loadnid = kl.*(nelx+1).*(nely+1).+il.*(nely+1).+(nely+1 .-jl)
loaddof = 3 .*loadnid[:] .- 1;
# USER-DEFINED SUPPORT FIXED DOFs m= Matlab.meshgrid(0:0,0:nely,0:nelz)# Coordinates
iif=m[1]
jf=m[2]
kf=m[3]
fixednid = kf.*(nelx+1).*(nely+1) .+iif .*(nely .+1) .+(nely .+1 .-jf) # Node IDs
fixeddof = [3 .*fixednid[:]; 3 .*fixednid[:].-1; 3 .*fixednid[:].-2] # DOFs
# DEFINE LOADS AND SUPPORTS (HALF MBB-BEAM)
F= sparse(loaddof,fill(1,length(loaddof)),fill(-1,length(loaddof)),ndof,1);
U = zeros(ndof,1)
alldofs = 1:ndof
freedofs = setdiff(1:ndof,fixeddof)
o=0;alpha_a=0
# inner iteration
for i =0: (maxIter-1)
# FE-ANALYSIS
E = e[1]*alpha[:,1].^q
for phase = 2:p
E = E + e[phase]*alpha[:,phase].^q;
end
sK = reshape(KE[:]*E[:]',24*24*nelx*nely*nelz,1)
K = sparse(vec(iK),vec(jK),vec(sK)); K = (K+K')/2
@timed U[freedofs] = K[freedofs,freedofs] \ Array(F[freedofs])
# Objective function and sensitivity analysis
ce = sum((U[edofMat]*KE).*U[edofMat],dims=2)
o = sum(sum(E.*ce))
# FILTERING OF SENSITIVITIES
dc = 0.0 .-(q .*(e[a].-e[b]).*alpha[:,a].^(q-1)).*ce;
dc = H*(alpha[:,a].*dc)./Hs./max.(1e-3,alpha[:,a])
dc = min.(dc,0.0)
# UPDATE LOWER AND UPPER BOUNDS OF DESIGN VARIABLES
move = 0.2;
r = ones(nx*ny*nz,1)
for k = 1:p
#if (k ~= a) && (k ~= b)
if (k != a) && (k != b)
# display("k $k a $a b $b")
r = r .- alpha[:,k]
end
end
l = max.(0,alpha[:,a] .-move)
u = min.(r,alpha[:,a] .+move)
# OPTIMALITY CRITERIA UPDATE OF DESIGN VARIABLES
l1 = 0; l2 = 1e9
while (l2-l1)/(l1+l2) > 1e-3
lmid = 0.5*(l2+l1)
alpha_a = max.(l,min.(u,alpha[:,a].*sqrt.(-dc./lmid)));
if sum(alpha_a) > nx*ny*v[a]
l1 = lmid
else
l2 = lmid
end
end
alpha[:,a] = alpha_a
alpha[:,b] = r .-alpha_a
end
return o,alpha
end
########
function multitop3D_compliant(nx,ny,nz,tol_out,tol_f,iter_max_in,iter_max_out,p,q,e,v,rf,getProblem=inverter)
anim=Animation()
alpha = zeros(Float64,nx*ny*nz,p)
for i = 1:p
alpha[:,i] .= v[i]
end
# MAKE FILTER
H,Hs = make_filter3D(nx,ny,nz,rf)
obj=0
change_out = 2*tol_out; iter_out = 0;
while (iter_out < iter_max_out) && (change_out > tol_out)
alpha_old = copy(alpha)
for a = 1:p
for b = a+1:p
obj,alpha = bi_top3D_compliant(a,b,nx,ny,nz,p,v,e,q,alpha,H,Hs,iter_max_in,getProblem);
end
end
iter_out = iter_out + 1;
change_out = norm(alpha[:].-alpha_old[:],Inf);
display("Iter:$iter_out Obj.:$obj change:$change_out");
# UPDATE FILTER
if (change_out < tol_f) && (rf>3)
tol_f = 0.99*tol_f; rf = 0.99*rf;
H,Hs = make_filter3D(nx,ny,nz,rf)
end
# SCREEN OUT TEMPORAL TOPOLOGY EVERY 5 ITERATIONS
if mod(iter_out,5)==0
I = make_bitmap_3d(p,nx,ny,nz,alpha)
threshold=0.05
display(topplotmulti3d(nx,ny,nz,I,threshold))
frame(anim)
end
end
return anim,alpha
end
function bi_top3D_compliant(a,b,nx,ny,nz,p,v,e,q,alpha_old,H,Hs,maxIter,getProblem=inverter)
alpha = copy(alpha_old)
nu = 0.3;
# PREPARE FINITE ELEMENT ANALYSIS
KE=lk_H8(nu)
nelx=nx
nely=ny
nelz=nz
nele = nelx*nely*nelz
ndof = 3*(nelx+1)*(nely+1)*(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 = convert(Array{Int},reshape(kron(edofMat,ones(24,1))',24*24*nele,1))
jK = convert(Array{Int},reshape(kron(edofMat,ones(1,24))',24*24*nele,1))
# # USER-DEFINED LOAD DOFs
# il = [0 nelx]; jl = [nely nely]; kl = [0 0];
# loadnid = kl .*(nelx+1) .*(nely+1) .+il .*(nely .+1) .+(nely+1 .-jl);
# loaddof = 3 .*loadnid[:] .- 2;
# din = loaddof[1]; dout = loaddof[2];
# F = sparse(loaddof,[1;2],[1;-1],ndof,2);
# # USER-DEFINED SUPPORT FIXED DOFs
# # Top face
# m = Matlab.meshgrid(0:nelx,0:nelz)
# iif=m[1] # kf=m[2]
# fixednid = kf .*(nelx+1) .*(nely+1) .+iif .*(nely+1) .+1;
# fixeddof_t = 3 .*fixednid[:] .-1;
# # Side face
# m = Matlab.meshgrid(0:nelx,0:nely)
# iif=m[1]
# jf=m[2]
# fixednid = iif .*(nely+1) .+(nely+1 .-jf);
# fixeddof_s = 3*fixednid[:];
# # Two pins
# iif = [0;0]; jf = [0;1]; kf = [nelz;nelz];
# fixednid = kf .*(nelx+1) .*(nely .+1) .+iif .*(nely .+1) .+(nely+1 .-jf)
# fixeddof_p = [3 .*fixednid[:]; 3 .*fixednid[:].-1; 3 .*fixednid[:].-2] # DOFs
# # Fixed DOFs
# fixeddof = union(fixeddof_t,fixeddof_s);
# fixeddof = union(fixeddof, fixeddof_p);
# fixeddof = sort(fixeddof);
# # Displacement vector
# U = zeros(ndof,2);
# freedofs = setdiff(1:ndof,fixeddof)
U,F,freedofs,din,dout=getProblem(nelx,nely,nelz)
o=0;alpha_a=0
# inner iteration
for i =0: (maxIter-1)
# FE-ANALYSIS
E = e[1]*alpha[:,1].^q
for phase = 2:p
E = E + e[phase]*alpha[:,phase].^q;
end
sK = reshape(KE[:]*E[:]',24*24*nelx*nely*nelz,1)
K = sparse(vec(iK),vec(jK),vec(sK));
K[din,din] = K[din,din] + 0.1;
K[dout,dout] = K[dout,dout] + 0.1;
K = (K+K')/2;
# @timed U[freedofs] = K[freedofs,freedofs] \ Array(F[freedofs])
@timed U[freedofs,:] = K[freedofs,freedofs] \ Array(F[freedofs,:])
U1 = U[:,1]
U2 = U[:,2]
# Objective function and sensitivity analysis
ce = -sum((U1[edofMat]*KE).*U2[edofMat],dims=2)
# o = sum(sum(E.*ce))
o = U[dout,1]
# FILTERING OF SENSITIVITIES
dc = 0.0 .-(q .*(e[a].-e[b]).*alpha[:,a].^(q-1)).*ce;
dc = H*(alpha[:,a].*dc)./Hs./max.(1e-3,alpha[:,a])
dc = min.(dc,0.0)
# UPDATE LOWER AND UPPER BOUNDS OF DESIGN VARIABLES
move = 0.1;
r = ones(nx*ny*nz,1)
for k = 1:p
#if (k ~= a) && (k ~= b)
if (k != a) && (k != b)
# display("k $k a $a b $b")
r = r .- alpha[:,k]
end
end
l = max.(0,alpha[:,a] .-move)
u = min.(r,alpha[:,a] .+move) # OPTIMALITY CRITERIA UPDATE OF DESIGN VARIABLES
l1 = 0; l2 = 1e9
# while (l2-l1)/(l1+l2) > 1e-3
while (l2-l1)/(l1+l2) > 1e-4 && l2 > 1e-40
lmid = 0.5*(l2+l1)
# alpha_a = max.(l,min.(u,alpha[:,a].*sqrt.(-dc./lmid)));
alpha_a = max.(l,min.(u,alpha[:,a].*((max.(1e-10,-dc./lmid)).^0.3 )));
if sum(alpha_a) > nx*ny*v[a]
l1 = lmid
else
l2 = lmid
end
end
alpha[:,a] = alpha_a
alpha[:,b] = r .-alpha_a
end
return o,alpha
end
########
function multitop3DKE(nx,ny,nz,tol_out,tol_f,iter_max_in,iter_max_out,p,q,KEs,Es,v,rf)
anim=Animation()
alpha = zeros(Float64,nx*ny*nz,p)
for i = 1:p
alpha[:,i] .= v[i]
end
# MAKE FILTER
H,Hs = make_filter3D(nx,ny,nz,rf)
obj=0
change_out = 2*tol_out; iter_out = 0;
while (iter_out < iter_max_out) && (change_out > tol_out)
alpha_old = copy(alpha)
for a = 1:p
for b = a+1:p
obj,alpha = bi_top3DKE(a,b,nx,ny,nz,p,v,KEs,Es,q,alpha,H,Hs,iter_max_in);
end
end
iter_out = iter_out + 1;
change_out = norm(alpha[:].-alpha_old[:],Inf);
display("Iter:$iter_out Obj.:$obj change:$change_out");
# UPDATE FILTER
if (change_out < tol_f) && (rf>3)
tol_f = 0.99*tol_f; rf = 0.99*rf;
H,Hs = make_filter3D(nx,ny,nz,rf)
end
# SCREEN OUT TEMPORAL TOPOLOGY EVERY 5 ITERATIONS
if mod(iter_out,5)==0
I = make_bitmap_3d(p,nx,ny,nz,alpha)
threshold=0.05*2.0
display(topplotmulti3d(nx,ny,nz,I,threshold))
frame(anim)
end
end
return anim,alpha
end
function bi_top3DKE(a,b,nx,ny,nz,p,v,KEs,Es,q,alpha_old,H,Hs,maxIter)
alpha = copy(alpha_old)
nu = 0.3;
# PREPARE FINITE ELEMENT ANALYSIS
# KE=lk_H8(nu)
nelx=nx
nely=ny
nelz=nz
nele = nelx*nely*nelz
ndof = 3*(nelx+1)*(nely+1)*(nelz+1) nnx = nelx+1; nny = nely+1; nnz = 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 = convert(Array{Int},reshape(kron(edofMat,ones(24,1))',24*24*nele,1))
jK = convert(Array{Int},reshape(kron(edofMat,ones(1,24))',24*24*nele,1))
# USER-DEFINED LOAD DOFs
m= Matlab.meshgrid(nelx:nelx, 0:0, 0:nelz)
il=m[1]
jl=m[2]
kl=m[3]
loadnid = kl.*(nelx+1).*(nely+1).+il.*(nely+1).+(nely+1 .-jl)
loaddof = 3 .*loadnid[:] .- 1;
# USER-DEFINED SUPPORT FIXED DOFs
m= Matlab.meshgrid(0:0,0:nely,0:nelz)# Coordinates
iif=m[1]
jf=m[2]
kf=m[3]
fixednid = kf.*(nelx+1).*(nely+1) .+iif .*(nely .+1) .+(nely .+1 .-jf) # Node IDs
fixeddof = [3 .*fixednid[:]; 3 .*fixednid[:].-1; 3 .*fixednid[:].-2] # DOFs
# DEFINE LOADS AND SUPPORTS (HALF MBB-BEAM)
F= sparse(loaddof,fill(1,length(loaddof)),fill(-1,length(loaddof)),ndof,1);
U = zeros(ndof,1)
alldofs = 1:ndof
freedofs = setdiff(1:ndof,fixeddof)
o=0;alpha_a=0
# inner iteration
for i =0: (maxIter-1)
# FE-ANALYSIS
# KE[:]*E[:]'
KEE = KEs[1][:]*alpha[:,1][:]'.^q
for phase = 2:p
KEE = KEE + KEs[phase][:]*alpha[:,phase][:]'.^q;
end
EE = mean(Es[1])*alpha[:,1].^q
for phase = 2:p
EE = EE + mean(Es[phase])*alpha[:,phase].^q;
end
sK = reshape(KEE[:],24*24*nelx*nely*nelz,1)
K = sparse(vec(iK),vec(jK),vec(sK)); K = (K+K')/2
@timed U[freedofs] = K[freedofs,freedofs] \ Array(F[freedofs])
KE=reshape(sum(KEE,dims=2),24,24)
# dc=zeros(nelx* nely* nelz)
# for i=1:Int(nelx* nely* nelz)
# KE=(q .*(KEs[a].-KEs[b]).*alpha[i,a].^(q-1))
# # ce[i]=sum((U[edofMat][i,:]'*KE).*U[edofMat][i,:]')
# dc[i]=-sum((U[edofMat][i,:]'*KE).*U[edofMat][i,:]')
# end
ce=zeros(nelx* nely* nelz)
for i=1:Int(nelx* nely* nelz)
KE = KEs[1].*alpha[i,1]
for phase = 2:p
KE = KE + KEs[phase].*alpha[i,phase]
end ce[i]=sum((U[edofMat][i,:]'*KE).*U[edofMat][i,:]')
end
# display(size(KEE))
# OBJECTIVE FUNCTION AND SENSITIVITY ANALYSIS
# cc = 0.0;
# dc = zeros(nelx, nely, nelz);
# for elz = 1:nelz
# for ely = 1:nely
# for elx = 1:nelx
# # KE = get_KE(truss, x, Es, vs, Gs, elx, ely, elz, false);
# # DKE = get_KE(truss, x, Es, vs, Gs, elx, ely, elz, true);
# # KE=reshape(KEEE[:,get_num(nelx, nely, nelz, elx, ely, elz)] ,24,24)
# # display(get_num1(nelx, nely, nelz, elx, ely, elz))
# # display(get_num(nnx, nny, nnz, elx, ely, elz))
# # KE = (KEs[a].*reshape(alpha[:,a][:],nelx, nely, nelz)[elx, ely, elz])-(KEs[b].*reshape(alpha[:,b][:],nelx, nely, nelz)[elx, ely, elz])
# # KE = KEs[1].*reshape(alpha[:,1][:],nelx, nely, nelz)[elx, ely, elz]
# # for phase = 2:p
# # KE = KE + KEs[phase].*reshape(alpha[:,phase][:],nelx, nely, nelz)[elx, ely, elz];
# # end
# DKE=KE
# dofs1 = get_elem_dofs(nnx, nny, nnz, elx, ely, elz);
# dofs1 = get_elem_dofs(nelx, nely, nelz, elx, ely, elz);
# # display(dofs1)
# Ue = U[dofs1];
# cc = cc + Ue'*KE*Ue;
# dc[elx,ely,elz] = Ue'*DKE*Ue;
# end
# end
# end
# display(size(dc))
# display(c)
# ce=dc[:]
# o=cc
# display(ce)
# display(sum(o))
# dc=dc[:]
# Objective function and sensitivity analysis
# ce= sum((U[edofMat]*KE).*U[edofMat],dims=2)
o = sum(sum(EE.*ce))
# o = sum(sum(-dc))
# display(ce)
# display(sum(o1))
# display(size(ce))
# display(size(EE))
# FILTERING OF SENSITIVITIES
dc = 0.0 .-(q .*mean(Es[a].-Es[b]).*alpha[:,a].^(q-1)).*ce;
dc = H*(alpha[:,a].*dc)./Hs./max.(1e-3,alpha[:,a])
dc = min.(dc,0.0)
# UPDATE LOWER AND UPPER BOUNDS OF DESIGN VARIABLES
move = 0.2;
r = ones(nx*ny*nz,1)
for k = 1:p #if (k ~= a) && (k ~= b)
if (k != a) && (k != b)
# display("k $k a $a b $b")
r = r .- alpha[:,k]
end
end
l = max.(0,alpha[:,a] .-move)
u = min.(r,alpha[:,a] .+move)
# OPTIMALITY CRITERIA UPDATE OF DESIGN VARIABLES
l1 = 0; l2 = 1e9
while (l2-l1)/(l1+l2) > 1e-3
lmid = 0.5*(l2+l1)
alpha_a = max.(l,min.(u,alpha[:,a].*sqrt.(-dc./lmid)));
if sum(alpha_a) > nx*ny*v[a]
l1 = lmid
else
l2 = lmid
end
end
alpha[:,a] = alpha_a
alpha[:,b] = r .-alpha_a
end
return o,alpha
end
function bi_top3DKEOldKindofWorking(a,b,nx,ny,nz,p,v,KEs,Es,q,alpha_old,H,Hs,maxIter)
alpha = copy(alpha_old)
nu = 0.3;
# PREPARE FINITE ELEMENT ANALYSIS
# KE=lk_H8(nu)
nelx=nx
nely=ny
nelz=nz
nele = nelx*nely*nelz
ndof = 3*(nelx+1)*(nely+1)*(nelz+1)
nnx = nelx+1; nny = nely+1; nnz = 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 = convert(Array{Int},reshape(kron(edofMat,ones(24,1))',24*24*nele,1))
jK = convert(Array{Int},reshape(kron(edofMat,ones(1,24))',24*24*nele,1))
# USER-DEFINED LOAD DOFs
m= Matlab.meshgrid(nelx:nelx, 0:0, 0:nelz)
il=m[1]
jl=m[2]
kl=m[3]
loadnid = kl.*(nelx+1).*(nely+1).+il.*(nely+1).+(nely+1 .-jl)
loaddof = 3 .*loadnid[:] .- 1;
# USER-DEFINED SUPPORT FIXED DOFs
m= Matlab.meshgrid(0:0,0:nely,0:nelz)# Coordinates
iif=m[1]
jf=m[2]
kf=m[3]
fixednid = kf.*(nelx+1).*(nely+1) .+iif .*(nely .+1) .+(nely .+1 .-jf) # Node IDs
fixeddof = [3 .*fixednid[:]; 3 .*fixednid[:].-1; 3 .*fixednid[:].-2] # DOFs
# DEFINE LOADS AND SUPPORTS (HALF MBB-BEAM)
F= sparse(loaddof,fill(1,length(loaddof)),fill(-1,length(loaddof)),ndof,1);
U = zeros(ndof,1)
alldofs = 1:ndof
freedofs = setdiff(1:ndof,fixeddof)
o=0;alpha_a=0
# inner iteration
for i =0: (maxIter-1)
# FE-ANALYSIS
# KE[:]*E[:]'
KEE = KEs[1][:]*alpha[:,1][:]'.^q
for phase = 2:p
KEE = KEE + KEs[phase][:]*alpha[:,phase][:]'.^q;
end
EE = mean(Es[1])*alpha[:,1].^q
for phase = 2:p
EE = EE + mean(Es[phase])*alpha[:,phase].^q;
end
sK = reshape(KEE[:],24*24*nelx*nely*nelz,1)
K = sparse(vec(iK),vec(jK),vec(sK)); K = (K+K')/2
@timed U[freedofs] = K[freedofs,freedofs] \ Array(F[freedofs])
KE=reshape(sum(KEE,dims=2),24,24)
ce=zeros(nelx* nely* nelz)
# display(size(U[edofMat]))
for i=1:Int(nelx* nely* nelz)
KE = KEs[1].*alpha[i,1]
for phase = 2:p
KE = KE + KEs[phase].*alpha[i,phase]
end
# ce[i]=sum((U[edofMat][i,:]'*reshape(KE[:,i],24,24)).*U[edofMat][i,:]')
ce[i]=sum((U[edofMat][i,:]'*KE).*U[edofMat][i,:]')
end
# display(size(KEE))
# OBJECTIVE FUNCTION AND SENSITIVITY ANALYSIS
# cc = 0.0;
# dc = zeros(nelx, nely, nelz);
# for elz = 1:nelz
# for ely = 1:nely
# for elx = 1:nelx
# # KE = get_KE(truss, x, Es, vs, Gs, elx, ely, elz, false);
# # DKE = get_KE(truss, x, Es, vs, Gs, elx, ely, elz, true);
# # KE=reshape(KEEE[:,get_num(nelx, nely, nelz, elx, ely, elz)] ,24,24)
# # display(get_num1(nelx, nely, nelz, elx, ely, elz))
# # display(get_num(nnx, nny, nnz, elx, ely, elz))
# # KE = (KEs[a].*reshape(alpha[:,a][:],nelx, nely, nelz)[elx, ely, elz])-(KEs[b].*reshape(alpha[:,b][:],nelx, nely, nelz)[elx, ely, elz])
# # KE = KEs[1].*reshape(alpha[:,1][:],nelx, nely, nelz)[elx, ely, elz]
# # for phase = 2:p
# # KE = KE + KEs[phase].*reshape(alpha[:,phase][:],nelx, nely, nelz)[elx, ely, elz];
# # end
# DKE=KE
# dofs1 = get_elem_dofs(nnx, nny, nnz, elx, ely, elz);
# dofs1 = get_elem_dofs(nelx, nely, nelz, elx, ely, elz);
# # display(dofs1)
# Ue = U[dofs1];
# cc = cc + Ue'*KE*Ue;
# dc[elx,ely,elz] = Ue'*DKE*Ue;
# end
# end
# end
# display(size(dc)) # display(c)
# ce=dc[:]
# o=cc
# display(ce)
# display(sum(o))
# dc=dc[:]
# Objective function and sensitivity analysis
# ce= sum((U[edofMat]*KE).*U[edofMat],dims=2)
o = sum(sum(EE.*ce))
# display(ce)
# display(sum(o1))
# display(size(ce))
# display(size(EE))
# FILTERING OF SENSITIVITIES
dc = 0.0 .-(q .*mean(Es[a].-Es[b]).*alpha[:,a].^(q-1)).*ce;
dc = H*(alpha[:,a].*dc)./Hs./max.(1e-3,alpha[:,a])
dc = min.(dc,0.0)
# UPDATE LOWER AND UPPER BOUNDS OF DESIGN VARIABLES
move = 0.2;
r = ones(nx*ny*nz,1)
for k = 1:p
#if (k ~= a) && (k ~= b)
if (k != a) && (k != b)
# display("k $k a $a b $b")
r = r .- alpha[:,k]
end
end
l = max.(0,alpha[:,a] .-move)
u = min.(r,alpha[:,a] .+move)
# OPTIMALITY CRITERIA UPDATE OF DESIGN VARIABLES
l1 = 0; l2 = 1e9
while (l2-l1)/(l1+l2) > 1e-3
lmid = 0.5*(l2+l1)
alpha_a = max.(l,min.(u,alpha[:,a].*sqrt.(-dc./lmid)));
if sum(alpha_a) > nx*ny*v[a]
l1 = lmid
else
l2 = lmid
end
end
alpha[:,a] = alpha_a
alpha[:,b] = r .-alpha_a
end
return o,alpha
end
####
########
function multitop3D_compliantKE(nx,ny,nz,tol_out,tol_f,iter_max_in,iter_max_out,p,q,KEs,Es,v,rf,getProblem=inverter)
anim=Animation()
alpha = zeros(Float64,nx*ny*nz,p)
for i = 1:p
alpha[:,i] .= v[i]
end
# MAKE FILTER H,Hs = make_filter3D(nx,ny,nz,rf)
obj=0
change_out = 2*tol_out; iter_out = 0;
while (iter_out < iter_max_out) && (change_out > tol_out)
alpha_old = copy(alpha)
for a = 1:p
for b = a+1:p
obj,alpha = bi_top3D_compliantKE(a,b,nx,ny,nz,p,v,KEs,Es,q,alpha,H,Hs,iter_max_in,getProblem);
end
end
iter_out = iter_out + 1;
change_out = norm(alpha[:].-alpha_old[:],Inf);
display("Iter:$iter_out Obj.:$obj change:$change_out");
# UPDATE FILTER
if (change_out < tol_f) && (rf>3)
tol_f = 0.99*tol_f; rf = 0.99*rf;
H,Hs = make_filter3D(nx,ny,nz,rf)
end
# SCREEN OUT TEMPORAL TOPOLOGY EVERY 5 ITERATIONS
if mod(iter_out,5)==0
I = make_bitmap_3d(p,nx,ny,nz,alpha)
threshold=0.05
display(topplotmulti3d(nx,ny,nz,I,threshold))
frame(anim)
end
end
return anim,alpha
end
function bi_top3D_compliantKE(a,b,nx,ny,nz,p,v,KEs,Es,q,alpha_old,H,Hs,maxIter,getProblem=inverter)
alpha = copy(alpha_old)
nu = 0.3;
# PREPARE FINITE ELEMENT ANALYSIS
# KE=lk_H8(nu)
nelx=nx
nely=ny
nelz=nz
nele = nelx*nely*nelz
ndof = 3*(nelx+1)*(nely+1)*(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 = convert(Array{Int},reshape(kron(edofMat,ones(24,1))',24*24*nele,1))
jK = convert(Array{Int},reshape(kron(edofMat,ones(1,24))',24*24*nele,1))
# # USER-DEFINED LOAD DOFs
# il = [0 nelx]; jl = [nely nely]; kl = [0 0];
# loadnid = kl .*(nelx+1) .*(nely+1) .+il .*(nely .+1) .+(nely+1 .-jl);
# loaddof = 3 .*loadnid[:] .- 2;
# din = loaddof[1]; dout = loaddof[2];
# F = sparse(loaddof,[1;2],[1;-1],ndof,2);
# # USER-DEFINED SUPPORT FIXED DOFs
# # Top face
# m = Matlab.meshgrid(0:nelx,0:nelz)
# iif=m[1]
# kf=m[2]
# fixednid = kf .*(nelx+1) .*(nely+1) .+iif .*(nely+1) .+1;
# fixeddof_t = 3 .*fixednid[:] .-1;
# # Side face
# m = Matlab.meshgrid(0:nelx,0:nely)
# iif=m[1]
# jf=m[2]
# fixednid = iif .*(nely+1) .+(nely+1 .-jf); # fixeddof_s = 3*fixednid[:];
# # Two pins
# iif = [0;0]; jf = [0;1]; kf = [nelz;nelz];
# fixednid = kf .*(nelx+1) .*(nely .+1) .+iif .*(nely .+1) .+(nely+1 .-jf)
# fixeddof_p = [3 .*fixednid[:]; 3 .*fixednid[:].-1; 3 .*fixednid[:].-2] # DOFs
# # Fixed DOFs
# fixeddof = union(fixeddof_t,fixeddof_s);
# fixeddof = union(fixeddof, fixeddof_p);
# fixeddof = sort(fixeddof);
# # Displacement vector
# U = zeros(ndof,2);
# freedofs = setdiff(1:ndof,fixeddof)
U,F,freedofs,din,dout=getProblem(nelx,nely,nelz)
o=0;alpha_a=0
# inner iteration
for i =0: (maxIter-1)
# FE-ANALYSIS
# E = e[1]*alpha[:,1].^q
# for phase = 2:p
# E = E + e[phase]*alpha[:,phase].^q;
# end
KEE = KEs[1][:]*alpha[:,1][:]'.^q
for phase = 2:p
KEE = KEE + KEs[phase][:]*alpha[:,phase][:]'.^q
end
# EE = mean(Es[1])*alpha[:,1].^q
# for phase = 2:p
# EE = EE + mean(Es[phase])*alpha[:,phase].^q;
# end
sK = reshape(KEE[:],24*24*nelx*nely*nelz,1)
# sK = reshape(KE[:]*E[:]',24*24*nelx*nely*nelz,1)
K = sparse(vec(iK),vec(jK),vec(sK));
K[din,din] = K[din,din] + 0.1;
K[dout,dout] = K[dout,dout] + 0.1;
K = (K+K')/2;
# @timed U[freedofs] = K[freedofs,freedofs] \ Array(F[freedofs])
@timed U[freedofs,:] = K[freedofs,freedofs] \ Array(F[freedofs,:])
U1 = U[:,1]
U2 = U[:,2]
KE=reshape(sum(KEE,dims=2),24,24)
# ce=zeros(nelx* nely* nelz)
# for i=1:Int(nelx* nely* nelz)
# KE = KEs[1].*alpha[i,1]
# for phase = 2:p
# KE = KE + KEs[phase].*alpha[i,phase]
# end
# ce[i]=-sum((U1[edofMat][i,:]'*KE).*U2[edofMat][i,:]')
# end
ce=zeros(nelx* nely* nelz)
for i=1:Int(nelx* nely* nelz)
KE = KEs[1].*alpha[i,1]
for phase = 2:p
KE = KE + KEs[phase].*alpha[i,phase]
end
ce[i]=-sum((U1[edofMat][i,:]'*KE).*U2[edofMat][i,:]') end
# Objective function and sensitivity analysis
# ce = -sum((U1[edofMat]*KE).*U2[edofMat],dims=2)
# o = sum(sum(E.*ce))
o = U[dout,1]
# FILTERING OF SENSITIVITIES
# dc = 0.0 .-(q .*(e[a].-e[b]).*alpha[:,a].^(q-1)).*ce;
dc = 0.0 .-(q .*mean(Es[a].-Es[b]).*alpha[:,a].^(q-1)).*ce;
dc = H*(alpha[:,a].*dc)./Hs./max.(1e-3,alpha[:,a])
dc = min.(dc,0.0)
# UPDATE LOWER AND UPPER BOUNDS OF DESIGN VARIABLES
move = 0.1;
r = ones(nx*ny*nz,1)
for k = 1:p
#if (k ~= a) && (k ~= b)
if (k != a) && (k != b)
# display("k $k a $a b $b")
r = r .- alpha[:,k]
end
end
l = max.(0,alpha[:,a] .-move)
u = min.(r,alpha[:,a] .+move)
# OPTIMALITY CRITERIA UPDATE OF DESIGN VARIABLES
l1 = 0; l2 = 1e9
# while (l2-l1)/(l1+l2) > 1e-3
while (l2-l1)/(l1+l2) > 1e-4 && l2 > 1e-40
lmid = 0.5*(l2+l1)
# alpha_a = max.(l,min.(u,alpha[:,a].*sqrt.(-dc./lmid)));
alpha_a = max.(l,min.(u,alpha[:,a].*((max.(1e-10,-dc./lmid)).^0.3 )));
if sum(alpha_a) > nx*ny*v[a]
l1 = lmid
else
l2 = lmid
end
end
alpha[:,a] = alpha_a
alpha[:,b] = r .-alpha_a
end
return o,alpha
end