Commit 8d7bca89 authored by Amira Abdel-Rahman's avatar Amira Abdel-Rahman
Browse files

multigrid and 3d optimization

parent 9830a83e
Pipeline #12122 passed with stage
in 26 seconds
This diff is collapsed.
......@@ -109179,7 +109179,7 @@
],
"metadata": {
"kernelspec": {
"display_name": "Julia 1.5.2",
"display_name": "Julia 1.5.3",
"language": "julia",
"name": "julia-1.5"
},
......@@ -109187,7 +109187,7 @@
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.5.2"
"version": "1.5.3"
}
},
"nbformat": 4,
This diff is collapsed.
# Amira Abdel-Rahman
# (c) Massachusetts Institute of Technology 2021
# using Clustering
# using StatsPlots
# using ParallelKMeans
# using Distances
# gr(size=(400,300))
# one microstructure
function ConTop2D(Macro_struct, Micro_struct, penal, rmin,saveItr=5)
## USER-DEFINED LOOP PARAMETERS
......
# Amira Abdel-Rahman
# (c) Massachusetts Institute of Technology 2021
# using Makie
# AbstractPlotting.inline!(true)
# using Clustering
# using StatsPlots
# using ParallelKMeans
# using Distances
# gr(size=(400,300))
# one microstructure
function ConTop3D(Macro_struct, Micro_struct, penal, rmin,saveItr=5)
# USER-DEFINED LOOP PARAMETERS
......@@ -285,12 +294,15 @@ function MultiConTop3D(θ,Macro_struct, Micro_struct, penal, rmin,saveItr=5)
# ce[i]=sum((U[edofMat]*reshape(Kes[:,i],24,24)).*U[edofMat],dims=2)[i]
# end
UU=[]
KKK=[]
for i=1:Macro_nely*Macro_nelx*Macro_nelz
append!(UU,[U[edofMat][i,:]'])
append!(KKK,[reshape(Kes[:,i],24,24)])
end
# UU=[]
# KKK=[]
# for i=1:Macro_nely*Macro_nelx*Macro_nelz
# append!(UU,[U[edofMat][i,:]'])
# append!(KKK,[reshape(Kes[:,i],24,24)])
# end
UU=mapslices(x->[x]'[:], U[edofMat], dims=2)[:]
KKK=mapslices(x->[reshape(x,24,24)], Kes, dims=1)[:]
ce=sum(vcat((UU.*KKK)...).*U[edofMat],dims=2)
ce=reshape(ce,Macro_nely,Macro_nelx,Macro_nelz)
......@@ -437,7 +449,8 @@ function MultiConTop3DU(θ,Macro_struct, Micro_struct,prob, penal,saveItr=5,maxl
Micro_x=[]
for i=1:θ
append!(Micro_x ,[ones(Micro_nely,Micro_nelx,Micro_nelz)]);
Micro_x[i][Int(Micro_nely/2)-5:Int(Micro_nely/2)+5,Int(Micro_nelx/2)-5:Int(Micro_nelx/2)+5,Int(Micro_nelz/2)-5:Int(Micro_nelz/2)+5] .= 0;
# Micro_x[i][Int(Micro_nely/2)-5:Int(Micro_nely/2)+5,Int(Micro_nelx/2)-5:Int(Micro_nelx/2)+5,Int(Micro_nelz/2)-5:Int(Micro_nelz/2)+5] .= 0;
Micro_x[i][Int(Micro_nely/2):Int(Micro_nely/2)+1,Int(Micro_nelx/2):Int(Micro_nelx/2)+1,Int(Micro_nelz/2):Int(Micro_nelz/2)+1] .= 0;
end
beta = 1;
......@@ -508,6 +521,7 @@ function MultiConTop3DU(θ,Macro_struct, Micro_struct,prob, penal,saveItr=5,maxl
U[freedofs] = K[freedofs,freedofs]\Array(F[freedofs]);
# OBJECTIVE FUNCTION AND SENSITIVITY ANALYSIS
# Ke = elementMatVec3D(Macro_Elex/2, Macro_Eley/2, Macro_Elez/2, DH);
# ce = reshape(sum((U[edofMat]*Ke).*U[edofMat],dims=2),Macro_nely,Macro_nelx,Macro_nelz);
#get compliance slow try to make it fast!! ??
......@@ -516,12 +530,16 @@ function MultiConTop3DU(θ,Macro_struct, Micro_struct,prob, penal,saveItr=5,maxl
# ce[i]=sum((U[edofMat]*reshape(Kes[:,i],24,24)).*U[edofMat],dims=2)[i]
# end
UU=[]
KKK=[]
for i=1:Macro_nely*Macro_nelx*Macro_nelz
append!(UU,[U[edofMat][i,:]'])
append!(KKK,[reshape(Kes[:,i],24,24)])
end
# UU=[]
# KKK=[]
# for i=1:Macro_nely*Macro_nelx*Macro_nelz
# append!(UU,[U[edofMat][i,:]'])
# append!(KKK,[reshape(Kes[:,i],24,24)])
# end
UU=mapslices(x->[x]'[:], U[edofMat], dims=2)[:]
KKK=mapslices(x->[reshape(x,24,24)], Kes, dims=1)[:]
ce=sum(vcat((UU.*KKK)...).*U[edofMat],dims=2)
ce=reshape(ce,Macro_nely,Macro_nelx,Macro_nelz)
......@@ -599,7 +617,7 @@ function MultiConTop3DU(θ,Macro_struct, Micro_struct,prob, penal,saveItr=5,maxl
# scene= volume(Macro_xPhys, algorithm = :iso,isorange = 0.3, isovalue = 1.0,colormap=:grays)
display(scene)
save("./img/Macro_xPhys_3_$(Macro_Vol)_$(penal)_$(rmin)_$(loop).png",scene)
save("./img/Macro_xPhys_3_$(Macro_Vol)_$(penal)_$(Macro_rmin)_$(loop).png",scene)
for i=1:θ
temp=copy(Micro_xPhys[i])
......@@ -609,7 +627,7 @@ function MultiConTop3DU(θ,Macro_struct, Micro_struct,prob, penal,saveItr=5,maxl
scene = Scene(resolution = (200, 200))
scene= volume!(temp,colorrange=(0.0, θ), algorithm = :iso,isorange = 0.3, isovalue = i,colormap=cscheme)
display(scene)
save("./img/Micro_xPhys3U_3_$(i)_$(Micro_Vol[i])_$(rmin)_$(loop).png",scene)
save("./img/Micro_xPhys3U_3_$(i)_$(Micro_Vol[i])_$(Micro_rmin)_$(loop).png",scene)
end
# for i=1:θ
......@@ -811,14 +829,20 @@ function CompliantMultiConTop3DU(θ,Macro_struct, Micro_struct, prob,penal,saveI
U2 = U[:,2]
# OBJECTIVE FUNCTION AND SENSITIVITY ANALYSIS
UU1=[]
UU2=[]
KKK=[]
for i=1:Macro_nely*Macro_nelx*Macro_nelz
append!(UU1,[U1[edofMat][i,:]'])
append!(UU2,[U2[edofMat][i,:]'])
append!(KKK,[reshape(Kes[:,i],24,24)])
end
# UU1=[]
# UU2=[]
# KKK=[]
# for i=1:Macro_nely*Macro_nelx*Macro_nelz
# append!(UU1,[U1[edofMat][i,:]'])
# append!(UU2,[U2[edofMat][i,:]'])
# append!(KKK,[reshape(Kes[:,i],24,24)])
# end
UU1=mapslices(x->[x]'[:], U1[edofMat], dims=2)[:]
UU2=mapslices(x->[x]'[:], U2[edofMat], dims=2)[:]
KKK=mapslices(x->[reshape(x,24,24)], Kes, dims=1)[:]
ce=-sum(vcat((UU1.*KKK)...).*U2[edofMat],dims=2)
ce=reshape(ce,Macro_nely,Macro_nelx,Macro_nelz)
......@@ -1056,15 +1080,28 @@ function Uclustering3D(θ,Macro_nely,Macro_nelx,Macro_nelz,Macro_nele,Macro_Vol,
end
# display(heatmap(Macro_xPhys,showaxis = false,xaxis=nothing,yaxis=nothing,legend=nothing,fc=jet,clims=(0.0, 1.0),aspect_ratio=:equal))
scene = Scene(resolution = (400, 400))
scene= volume!(Macro_xPhys, algorithm = :mip,transparency=true,colorrange=(0.0, 1.0),colormap=:lighttest)
display(scene)
save("./img/FM3DU_$(Macro_Vol).png",scene)
##########clustering#########
# X=U[edofMat]'
X=hcat(Macro_xPhys[:],U[edofMat])'
for i=1:3
# display(i)
# XX=reshape(mean(U[edofMat]',dims=1),Macro_nely,Macro_nelx,Macro_nelz)
# XX=reshape(sum(U[edofMat]',dims=1),Macro_nely,Macro_nelx,Macro_nelz)
XX=reshape(U[edofMat]'[i,:],Macro_nely,Macro_nelx,Macro_nelz)
# display(minimum(XX))
# display(maximum(XX))
scene = Scene(resolution = (200, 200))
scene= volume!(XX, algorithm = :mip,transparency=true,colorrange=(minimum(XX), maximum(XX)),colormap=:lighttest)
display(scene)
# save("./img/FM3DU_$(Macro_Vol).png",scene)
end
# cluster X into θ clusters using K-means
# R = kmeans(X, θ; maxiter=200, display=:iter)
R = Clustering.kmeans(X, θ; maxiter=500)
......
......@@ -69,6 +69,53 @@ function lk_H8(nu)
return KE
end
## FUNCTION Ke2D - ELEMENT STIFFNESS MATRIX
function 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]);
return KE
end
## FUNCTION Ke3D - ELEMENT STIFFNESS MATRIX
function Ke3D(nu)
C = [2/9 1/18 1/24 1/36 1/48 5/72 1/3 1/6 1/12];
A11 = [-C[1] -C[3] -C[3] C[2] C[3] C[3]; -C[3] -C[1] -C[3] -C[3] -C[4] -C[5];
-C[3] -C[3] -C[1] -C[3] -C[5] -C[4]; C[2] -C[3] -C[3] -C[1] C[3] C[3];
C[3] -C[4] -C[5] C[3] -C[1] -C[3]; C[3] -C[5] -C[4] C[3] -C[3] -C[1]];
B11 = [C[7] 0 0 0 -C[8] -C[8]; 0 C[7] 0 C[8] 0 0; 0 0 C[7] C[8] 0 0;
0 C[8] C[8] C[7] 0 0; -C[8] 0 0 0 C[7] 0; -C[8] 0 0 0 0 C[7]];
A22 = [-C[1] -C[3] C[3] C[2] C[3] -C[3]; -C[3] -C[1] C[3] -C[3] -C[4] C[5];
C[3] C[3] -C[1] C[3] C[5] -C[4]; C[2] -C[3] C[3] -C[1] C[3] -C[3];
C[3] -C[4] C[5] C[3] -C[1] C[3]; -C[3] C[5] -C[4] -C[3] C[3] -C[1]];
B22 = [C[7] 0 0 0 -C[8] C[8]; 0 C[7] 0 C[8] 0 0; 0 0 C[7] -C[8] 0 0;
0 C[8] -C[8] C[7] 0 0; -C[8] 0 0 0 C[7] 0; C[8] 0 0 0 0 C[7]];
A12 = [C[6] C[3] C[5] -C[4] -C[3] -C[5]; C[3] C[6] C[5] C[3] C[2] C[3];
-C[5] -C[5] C[4] -C[5] -C[3] -C[4]; -C[4] C[3] C[5] C[6] -C[3] -C[5];
-C[3] C[2] C[3] -C[3] C[6] C[5]; C[5] -C[3] -C[4] C[5] -C[5] C[4]];
B12 = [-C[9] 0 -C[9] 0 C[8] 0; 0 -C[9] -C[9] -C[8] 0 -C[8]; C[9] C[9] -C[9] 0 C[8] 0;
0 -C[8] 0 -C[9] 0 C[9]; C[8] 0 -C[8] 0 -C[9] -C[9]; 0 C[8] 0 -C[9] C[9] -C[9]];
A13 = [-C[4] -C[5] -C[3] C[6] C[5] C[3]; -C[5] -C[4] -C[3] -C[5] C[4] -C[5];
C[3] C[3] C[2] C[3] C[5] C[6]; C[6] -C[5] -C[3] -C[4] C[5] C[3];
C[5] C[4] -C[5] C[5] -C[4] -C[3]; -C[3] C[5] C[6] -C[3] C[3] C[2]];
B13 = [0 0 C[8] -C[9] -C[9] 0; 0 0 C[8] C[9] -C[9] C[9]; -C[8] -C[8] 0 0 -C[9] -C[9];
-C[9] C[9] 0 0 0 -C[8]; -C[9] -C[9] C[9] 0 0 C[8]; 0 -C[9] -C[9] C[8] -C[8] 0];
A14 = [C[2] C[5] C[5] C[4] -C[5] -C[5]; C[5] C[2] C[5] C[5] C[6] C[3];
C[5] C[5] C[2] C[5] C[3] C[6]; C[4] C[5] C[5] C[2] -C[5] -C[5];
-C[5] C[6] C[3] -C[5] C[2] C[5]; -C[5] C[3] C[6] -C[5] C[5] C[2]];
B14 = [-C[9] 0 0 -C[9] C[9] C[9]; 0 -C[9] 0 -C[9] -C[9] 0; 0 0 -C[9] -C[9] 0 -C[9];
-C[9] -C[9] -C[9] -C[9] 0 0; C[9] -C[9] 0 0 -C[9] 0; C[9] 0 -C[9] 0 0 -C[9]];
A23 = [C[2] C[5] -C[5] C[4] -C[5] C[5]; C[5] C[2] -C[5] C[5] C[6] -C[3];
-C[5] -C[5] C[2] -C[5] -C[3] C[6]; C[4] C[5] -C[5] C[2] -C[5] C[5];
-C[5] C[6] -C[3] -C[5] C[2] -C[5]; C[5] -C[3] C[6] C[5] -C[5] C[2]];
B23 = [-C[9] 0 0 -C[9] C[9] -C[9]; 0 -C[9] 0 -C[9] -C[9] 0; 0 0 -C[9] C[9] 0 -C[9];
-C[9] -C[9] C[9] -C[9] 0 0; C[9] -C[9] 0 0 -C[9] 0; -C[9] 0 -C[9] 0 0 -C[9]];
KE = 1/(1+nu)/(2*nu-1)*([A11 A12 A13 A14; A12' A22 A23 A13'; A13' A23' A22 A12'; A14' A13 A12 A11] +
nu*[B11 B12 B13 B14; B12' B22 B23 B13'; B13' B23' B22 B12'; B14' B13 B12 B11]);
return KE
end
## SUB FUNCTION: elementMatVec2D
function elementMatVec2D(a, b, DH)
GaussNodes = [-1/sqrt(3); 1/sqrt(3)];
......
......@@ -94,7 +94,7 @@ function topologyOptimization(nelx,nely,volfrac,rmin,penal)
change = maximum(abs.(xnew[:].-x[:]))
x = xnew
m=mean(xPhys[:])
display(" It:$loop Obj:$c Vol:$m ch:$change ")
println(" It:$loop Obj:$(round.(c,digits=2)) Vol:$(round.(m,digits=2)) ch:$(round.(change,digits=2)) ")
if loop<20||mod(loop,10)==0
xxx=1 .- clamp.(xPhys,0,1)
display(Plots.heatmap(xxx,xaxis=nothing,showaxis = false,yaxis=nothing,legend=nothing,fc=:grays,clims=(0.0, 1.0),aspect_ratio=:equal))
......
# Amira Abdel-Rahman
# (c) Massachusetts Institute of Technology 2021
using LinearAlgebra
#### 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)
......@@ -9,7 +11,7 @@ function top2dmgcg(nelx,nely,volfrac,penal,rmin,ft,nl,cgtol,cgmax)
E0 = 1; Emin = 1e-9;
nu = 0.3;
## PREPARE FINITE ELEMENT ANALYSIS
KE =lk();
KE =Ke2D(nu);
# Prepare fine grid
nelem = nelx*nely;
nx = nelx+1; ny = nely+1;
......@@ -30,22 +32,30 @@ function top2dmgcg(nelx,nely,volfrac,penal,rmin,ft,nl,cgtol,cgmax)
for l = 1:nl-1
Pu[l,1] = Prepcoarse(nely/2^(l-1),nelx/2^(l-1));
end
# Define loads and supports (cantilever)
# # Define loads and supports (cantilever)
F = sparse([Int.(2*(nelx*(nely+1)+nely/2+1))],[1],[-1],2*(nely+1)*(nelx+1),1);
U = zeros(Int.(2*(nely+1)*(nelx+1)),1);
fixeddofs = 1:2*(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))
# Null space elimination of supports
N = ones(ndof);
N[fixeddofs] .= 0;
# Null = spdiags(N,0,ndof,ndof);
Null = sparse(diagm(N));
# Null1 = sparse(diagm(N));
Null = spdiagm(ndof,ndof,0 => N);
# display(sum(Null1-Null))
H ,Hs=make_filter(nelx,nely,rmin);
x = fill(volfrac,nely,nelx);
xPhys = x;
loop = 0;
change = 1;
## START ITERATION
while change > 1e-3 && loop < 500
while change > 1e-3 && loop < 100
loop = loop+1;
## FE-ANALYSIS
K = Array{Any}(undef, nl);
......@@ -53,12 +63,13 @@ function top2dmgcg(nelx,nely,volfrac,penal,rmin,ft,nl,cgtol,cgmax)
K[1,1] = sparse(vec(iK),vec(jK),vec(sK));
K[1,1] = Null'*K[1,1]*Null - (Null .-sparse(Matrix(1.0I, ndof, ndof)));
K[1,1] = Null'*K[1,1]*Null - (Null .-sparse(1.0I, ndof, ndof));
for l = 1:nl-1
K[l+1,1] = Pu[l,1]'*(K[l,1]*Pu[l,1]);
end
Lfac = cholesky(Hermitian(K[nl,1])).L; Ufac = Lfac';
# Lfac = cholesky(Hermitian(K[nl,1])).L; Ufac = Lfac';
Lfac = factorize(Matrix(Hermitian(K[nl,1]))).L ;Ufac = Lfac';
cgiters,cgres,U = mgcg(K,F,U,Lfac,Ufac,Pu,nl,1,cgtol,cgmax);
## OBJECTIVE FUNCTION AND SENSITIVITY ANALYSIS
......@@ -69,7 +80,7 @@ function top2dmgcg(nelx,nely,volfrac,penal,rmin,ft,nl,cgtol,cgmax)
## FILTERING/MODIFICATION OF SENSITIVITIES
if ft == 1
dc[:] = H*(x[:].*dc[:])./Hs./max(1e-3,x[:]);
dc[:] = H*(x[:].*dc[:])./Hs./max.(1e-3,x[:]);
elseif ft == 2
dc[:] = H*(dc[:]./Hs);
dv[:] = H*(dv[:]./Hs);
......@@ -96,7 +107,7 @@ function top2dmgcg(nelx,nely,volfrac,penal,rmin,ft,nl,cgtol,cgmax)
xPhys[:] = (H*xnew[:])./Hs;
end
## PRINT RESULTS
println(" Iter.: $loop Comp.: $c Vol.: $(mean(xPhys[:])) Chg.: $change CGres: $cgres CGits: $cgiters Penal: $penal");
println(" It: $loop Comp.: $(round.(c,digits=2)) Vol.: $(round.((mean(xPhys[:])),digits=2)) Chg.: $(round.(change,digits=2)) CGres: $(round.(cgres,digits=2)) CGits: $(round.(cgiters,digits=2)) Penal: $penal");
## PLOT DENSITIES
if loop<20||mod(loop,10)==0
xxx=1 .- clamp.(xPhys,0,1)
......@@ -108,6 +119,144 @@ function top2dmgcg(nelx,nely,volfrac,penal,rmin,ft,nl,cgtol,cgmax)
end
return xPhys,anim;
end
# To reproduce example #2 from the article, run:
# top3dmgcg(48,24,24,0.12,3.0,sqrt[3],1,4,1e-10,100)
function top3dmgcg(nelx,nely,nelz,volfrac,penal,rmin,ft,nl,cgtol,cgmax)
## MATERIAL PROPERTIES
E0 = 1;
Emin = 1e-9;
nu = 0.3;
## PREPARE FINITE ELEMENT ANALYSIS
KE = Ke3D(nu);
# Prepare fine grid
nelem = nelx*nely*nelz;
nx = nelx+1; ny = nely+1; nz = nelz+1;
ndof = 3*nx*ny*nz;
nodenrs = reshape(1:ny*nz*nx,ny,nz,nx)
edofVec = reshape(3*nodenrs[1:end-1,1:end-1,1:end-1].+1,nelem,1)
edofMat = repeat(edofVec,1,24).+repeat([0 1 2 3*nely.+[3 4 5 0 1 2] -3 -2 -1 3*(nely+1)*(nelz+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*nelem,1))
jK = convert(Array{Int},reshape(kron(edofMat,ones(1,24))',24*24*nelem,1))
# Prologation operators
# Pu = cell(nl-1,1);
Pu=Array{Any}(undef, nl-1);
for l = 1:nl-1
Pu[l,1] = prepcoarse3(nelz/2^(l-1),nely/2^(l-1),nelx/2^(l-1));
end
# Define loads and supports (cantilever)
F = sparse(Int.(3 .*nodenrs[1:nely+1,1,nelx+1][:]),fill(1,length(nodenrs[1:nely+1,1,nelx+1])),-sin.((0:nely)./nely.*pi),Int(ndof),1); # Sine load, bottom right
U = zeros(ndof,1);
fixeddofs = 1:3*(nely+1)*(nelz+1);
freedofs = setdiff(1:ndof,fixeddofs);
# Null space elimination of supports
N = ones(ndof); N[fixeddofs] .= 0;
# Null = spdiags(N,0,ndof,ndof);
# Null = sparse(diagm(N));
Null = spdiagm(ndof,ndof,0 => N);
# ## PREPARE FILTER
# H,Hs = make_filter3D(nelx, nely, nelz, rmin);
iH = ones(Int(nelx*nely*nelz*(2*(ceil(rmin)-1)+1)^3));
jH = ones(size(iH));
sH = zeros(size(iH));
k = 0;
for i1 = 1:nelx
for k1 = 1:nelz
for j1 = 1:nely
e1 = (i1-1)*nely*nelz + (k1-1)*nely + j1;
for i2 = max(i1-(ceil(rmin)-1),1):min(i1+(ceil(rmin)-1),nelx)
for k2 = max(k1-(ceil(rmin)-1),1):min(k1+(ceil(rmin)-1),nelz)
for j2 = max(j1-(ceil(rmin)-1),1):min(j1+(ceil(rmin)-1),nely)
e2 = (i2-1)*nely*nelz + (k2-1)*nely + j2;
k = k + 1;
iH[k] = e1;
jH[k] = e2;
sH[k] = max(0,rmin-sqrt((i1-i2)^2+(j1-j2)^2+(k1-k2)^2));
end
end
end
end
end
end
H = sparse(vec(iH),vec(jH),vec(sH));
Hs = sum(H,dims=2);
## INITIALIZE ITERATION
x = volfrac.*ones(nely,nelz,nelx);
xPhys = x;
loop = 0;
change = 1;
## START ITERATION
while change > 1e-2 && loop < 50
loop = loop+1;
#tic
#sK = reshape(KE(:)*(Emin+xPhys(:)'.^penal*(E0-Emin)),24*24*nelem,1);
#K = sparse(iK,jK,sK); K = (K+K')/2;
#U(freedofs,:) = K(freedofs,freedofs)\F(freedofs,:);
#toc
## FE-ANALYSIS
# K = cell(nl,1);
K = Array{Any}(undef, nl);
sK = reshape(KE[:]*(Emin.+xPhys[:]'.^penal*(E0-Emin)),576*nelem,1)
K[1,1] = sparse(vec(iK),vec(jK),vec(sK));
K[1,1] = Null'*K[1,1]*Null - (Null .-sparse(1.0I, ndof, ndof));
for l = 1:nl-1
K[l+1,1] = Pu[l,1]'*(K[l,1]*Pu[l,1]);
end
# Lfac = cholesky(Hermitian(K[nl,1])).L; Ufac = Lfac';
Lfac = factorize(Matrix(Hermitian(K[nl,1]))).L ;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],dims=2),nely,nelz,nelx);
c = sum(sum((Emin .+xPhys.^penal*(E0-Emin)).*ce));
dc = -penal*(E0-Emin)*xPhys.^(penal-1).*ce;
dv = ones(nely,nelz,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;xnew=0;
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.((0.0.-dc)./dv./lmid)))))
gt = g + sum((dv[:].*(xnew[:]-x[:])));
if gt > 0
l1 = lmid;
else
l2 = lmid;
end
end
change = maximum(abs.(xnew[:].-x[:]));
x = xnew;
## FILTERING OF DESIGN VARIABLES
if ft == 1
xPhys = xnew;
elseif ft == 2
xPhys[:] = (H*xnew[:])./Hs;
end
if loop%5==0
scene= volume(xPhys, algorithm = :iso, isorange = 0.3, isovalue = 1.0,colormap=:grays)
display(scene)
end
## PRINT RESULTS
println(" It: $loop Comp.: $(round.(c,digits=2)) Vol.: $(round.((mean(xPhys[:])),digits=2)) Chg.: $(round.(change,digits=2)) CGres: $(round.(cgres,digits=2)) CGits: $(round.(cgiters,digits=2)) Penal: $penal");
end
return xPhys
end
## FUNCTION mgcg - MULTIGRID PRECONDITIONED CONJUGATE GRADIENTS
function mgcg(A,b,u,Lfac,Ufac,Pu,nl,nswp,tol,maxiter)
r = b - A[1,1]*u;
......@@ -119,39 +268,42 @@ function mgcg(A,b,u,Lfac,Ufac,Pu,nl,nswp,tol,maxiter)
invD[l,1]= 1.0 ./diag(A[l,1],0);
end
z = VCycle(A,r,Lfac,Ufac,Pu,1,nl,invD,omega,nswp);
p = z;
rho = r'*z;
rho_p= rho;
#z = VCycle(A,r,Lfac,Ufac,Pu,1,nl,invD,omega,nswp);
#p = z;
#rho = r'*z;
rho_p= 0;relres=0;p=0;
ii=1
relres=0
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
for i = 1:maxiter
if ii ==1 || relres >= tol
# display("mgcg $i")
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;
ii=ii+1
end
ii=ii+1
#if relres < tol || i>=maxiter
#break
#end
end
return ii,relres,u
end
## FUNCTION VCycle - COARSE GRID CORRECTION
function VCycle(A,r,Lfac,Ufac,Pu,l,nl,invD,omega,nswp)
# display("VCycle $l")
z = 0*r;
z = smthdmpjac(z,A[l,1],r,invD[l,1],omega,nswp);
Az = A[l,1]*z;
......@@ -178,7 +330,7 @@ end
function Prepcoarse(ney,nex)
# Assemble state variable prolongation
maxnum =