# Amira Abdel-Rahman
# (c) Massachusetts Institute of Technology 2020
function getDataFromSetup3D(setup,scale)
# -----------------------------------------------------------------------
# A(nel,1) = cross-sectional area of elements
# E(nel,1) = Young's modulus of elements
# f(ndf,nnp) = prescribed nodal forces
# g(ndf,nnp) = prescribed nodal displacements
# idb(ndf,nnp) = 1 if the degree of freedom is prescribed, 0 otherwise
# ien(nen,nel) = element connectivities
# ndf = number of degrees of freedom per node
# nel = number of elements
# nen = number of element equations
# nnp = number of nodes
# nsd = number of spacial dimensions
# xn(nsd,nnp) = nodal coordinates
# =======================================================================
nodes=setup["nodes"]
edges=setup["edges"]
# ---- MESH -------------------------------------------------------------
nsd = 3; # number of spacial dimensions
ndf = 3; # number of degrees of freedom per node
nen = 2; # number of element nodes
nel = length(edges); # number of elements
nnp = length(nodes); # number of nodal points
# ---- NODAL COORDINATES ------------------------------------------------
# xn(i,N) = coordinate i for node N
# N = 1,...,nnp
# i = 1,...,nsd
# -----------------------------------------------------------------------
xn = zeros(nsd,nnp);
for i in 1:nnp
xn[1:3, i] = [(nodes[i]["position"]["x"]/scale) (nodes[i]["position"]["y"]/scale) (nodes[i]["position"]["z"]/scale)]';
end
# ---- NODAL COORDINATES ------------------------------------------------
# ien(a,e) = N
# N: global node number - N = 1,...,nnp
# e: element number - e = 1,...,nel
# a: local node number - a = 1,...,nen
# -----------------------------------------------------------------------
ien = zeros(nen,nel);
for i in 1:nel
ien[1:2,i] = [(edges[i]["source"]+1) (edges[i]["target"]+1)]' ;
end
len=zeros(nel);
for i in 1:nel
x1=(nodes[(edges[i]["source"]+1)]["position"]["x"]/scale)
x2=(nodes[(edges[i]["target"]+1)]["position"]["x"]/scale)
y1=(nodes[(edges[i]["source"]+1)]["position"]["y"]/scale)
y2=(nodes[(edges[i]["target"]+1)]["position"]["y"]/scale)
z1=(nodes[(edges[i]["source"]+1)]["position"]["z"]/scale)
z2=(nodes[(edges[i]["target"]+1)]["position"]["z"]/scale)
len[i] = sqrt((x1-x2)^2+(y1-y2)^2+(z1-z2)^2);
end
# ---- MATERIAL PROPERTIES ----------------------------------------------
# E(e) = E_mat
# e: element number - e = 1,...,nel
# -----------------------------------------------------------------------
E_mat = 29000*1000; # Young's modulus # lbf/in^2
E = E_mat*ones(nel,1);
# todo change to make it parameter
# ---- GEOMETRIC PROPERTIES ---------------------------------------------
# A(e) = A_bar
# e: element number - e = 1,...,nel
# -----------------------------------------------------------------------
#A = ones(nel);
#A[:] .= 400; # mm^2
# ---- BOUNDARY CONDITIONS ----------------------------------------------
# prescribed displacement flags (essential boundary condition)
#
# idb(i,N) = 1 if the degree of freedom i of the node N is prescribed
# = 0 otherwise
#
# 1) initialize idb to 0
# 2) enter the flag for prescribed displacement boundary conditions
# -----------------------------------------------------------------------
idb = zeros(ndf,nnp);
for i in 1:nnp
if nodes[i]["restrained_degrees_of_freedom"][1]
idb[1,i]=1;
end
if nodes[i]["restrained_degrees_of_freedom"][2]
idb[2,i]=1;
end
if nodes[i]["restrained_degrees_of_freedom"][3]
idb[3,i]=1;
end
end
# ---- BOUNDARY CONDITIONS: PRESCRIBED NODAL DISPLACEMENTS --------------
# g(i,N) = prescribed displacement magnitude
# N = 1,...,nnp
# i = 1,...,nsd
#
# 1) initialize g to 0
# 2) enter the values
# -----------------------------------------------------------------------
g = zeros(ndf,nnp);
# ---- BOUNDARY CONDITIONS: PRESCRIBED NODAL FORCES ---------------------
# f(i,N) = prescribed force magnitude
# N = 1,...,nnp
# i = 1,...,nsd
#
# 1) initialize f to 0
# 2) enter the values
# -----------------------------------------------------------------------
P = 20.0* 1000; # lbf
f = zeros(ndf,nnp);
for i in 1:nnp
f[1,i] = nodes[i]["force"]["x"]*P;
f[2,i] = nodes[i]["force"]["y"]*P;
f[3,i] = nodes[i]["force"]["z"]*P;
end
A=ones(nel)
# ---- NUMBER THE EQUATIONS ---------------------------------------------
# line 380
id,neq = number_eq(idb,ndf,nnp);
# ---- FORM THE ELEMENT STIFFNESS MATRICES ------------------------------
# line 324
nee = ndf*nen; # number of element equations
Ke = zeros(nee,nee,nel);
Te = zeros(nen*1,nen*nsd,nel); # *1 is specific to truss
for i = 1:nel
Ke[:,:,i],Te[:,:,i] = Ke_truss(A[i],E[i],ien[:,i],nee,nsd,xn);
end
return E,f,g,idb,ien,ndf,nel,nen,nnp,nsd,xn,len,Ke,Te
end
# =======================================================================
function getSensitivites(Ke0,dcomp,ien,nel,len)
dfdA=zeros(nel)
dgdA=zeros(nel)
for i in 1:nel
de=[dcomp[1,Int(ien[1,i])],dcomp[2,Int(ien[1,i])],dcomp[3,Int(ien[1,i])] ,dcomp[1,Int(ien[2,i])] ,dcomp[2,Int(ien[2,i])] ,dcomp[3,Int(ien[2,i])]]
# println(de)
# println(Ke0[:,:,i])
dKedA=Ke0[:,:,i]
dfdA[i]=(-1.0 .*de)'*dKedA*de
dgdA[i]=len[i]
end
return dfdA,dgdA
end
# =======================================================================
function getSensitivitesFinite(nel,epsilon,len)
dfdA=zeros(nel)
dgdA=zeros(nel)
A1 = ones(nel)
A1[:] .= 0.25;
Ke,K,F,d,ien,xn,stress,dcomp,g=FEM_truss(data_truss2,A1);
for i in 1:nel
A[:] .= 0.25; # mm^2
A[i]=0.25+epsilon
Ke,Kee,Fe,de,ien,xn,stress,dcomp,g=FEM_truss(data_truss2,A);
Fx=F'*d -d'*K*d +d'*F
Fxe=Fe'*de - de'*Kee*de + de'*Fe
dfdA[i]=(Fxe[1] - Fx[1])/epsilon
dgdA[i]=((sum(A .* len ) - 480) - (sum(A1 .* len ) - 480))/epsilon
end
return dfdA,dgdA
end
# =======================================================================
function optimizeTruss(problem,totalVolFactor,maxeval=500)
E,f,g,idb,ien,ndf,nel,nen,nnp,nsd,xn,len,Ke,Te=problem
nel=length(len)
function FA(x::Vector, grad::Vector)
K,F,d,stress,dcomp,g=FEM_truss(problem,x);
# display(F)
grad[:] .=getSensitivites(Ke,dcomp,ien,nel,len)[1]
Fx=F'*d -d'*K*d +d'*F
return Fx[1]
end
function FASIMP(x::Vector, grad::Vector)
print("hena")
K,F,d,stress,dcomp,g=FEM_frame(problem,x.^η);
# display(F)
grad[:] .=getSensitivitesSIMP(Ke0,dcomp,ien,nel,len,η,x)[1]
Fx=F'*d -d'*K*d +d'*F
return Fx[1]
end
function G(x::Vector, grad::Vector)
grad[:] .=len[:]
return (sum(x .* len ) - totalVolFactor)
end
FA(ones(nel)*totalVolFactor, fill(0.25,nel))
G(ones(nel)*totalVolFactor, fill(0.25,nel))
opt = Opt(:LD_MMA, nel)
opt.lower_bounds = fill(1e-3,nel)
opt.xtol_rel = 0.00001
opt.maxeval = maxeval
opt.min_objective = FA
inequality_constraint!(opt, (x,gg) -> G(x,gg), 1e-6)
display(@time (minf,minx,ret) = optimize(opt, ones(nel)))
numevals = opt.numevals # the number of function evaluations
display("got $minf after $numevals iterations (returned $ret)")
return minx
end
# =======================================================================