Commit 18a19ccd authored by Amira Abdel-Rahman's avatar Amira Abdel-Rahman
Browse files

toyota seat

parent 22e34a65
Pipeline #12203 passed with stage
in 23 seconds
This diff is collapsed.
......@@ -63,9 +63,9 @@
11
],
"restrained_degrees_of_freedom": [
false,
true,
false,
true,
true,
true,
true,
true
......@@ -393,7 +393,7 @@
},
"fixedDisplacement": {
"x": 1,
"y": -1,
"y": 1,
"z": 0
},
"angle": {
......@@ -471,9 +471,9 @@
59
],
"restrained_degrees_of_freedom": [
false,
true,
false,
true,
true,
true,
true,
true
......@@ -573,9 +573,9 @@
71
],
"restrained_degrees_of_freedom": [
false,
true,
false,
true,
true,
true,
true,
true
......@@ -1250,7 +1250,7 @@
},
"force": {
"x": 0,
"y": -10,
"y": 0,
"z": 0
},
"displacement": {
......@@ -1301,7 +1301,7 @@
},
"force": {
"x": 0,
"y": 0,
"y": -1,
"z": 0
},
"displacement": {
......@@ -1352,7 +1352,7 @@
},
"force": {
"x": 0,
"y": -10,
"y": 0,
"z": 0
},
"displacement": {
......
This diff is collapsed.
......@@ -809,17 +809,17 @@ function plotTruss(problem,X,scale,threshold=0)
nel=length(X)
E,f,g,idb,ien,ndf,nel,nen,nnp,nsd,xn,len,Ke,Te=problem;
K,F,d,stress,dcomp,g=FEM_truss(problem,X);
p=plot(axis=nothing,ticks=nothing, border=nothing, aspect_ratio=:equal)
p=Plots.plot(axis=nothing,ticks=nothing, border=nothing, aspect_ratio=:equal)
for i in 1:nel
if X[i]>threshold
if threshold>0
p=plot!([xn[1,Int(ien[1,i])],xn[1,Int(ien[2,i])]],
p=Plots.plot!([xn[1,Int(ien[1,i])],xn[1,Int(ien[2,i])]],
[xn[2,Int(ien[1,i])],xn[2,Int(ien[2,i])]],label="",
color=RGB(0.0,0.0,0.0),
linewidth = 3.0)
else
p=plot!([xn[1,Int(ien[1,i])],xn[1,Int(ien[2,i])]],
p=Plots.plot!([xn[1,Int(ien[1,i])],xn[1,Int(ien[2,i])]],
[xn[2,Int(ien[1,i])],xn[2,Int(ien[2,i])]],label="",
color=RGB(mapp(stress[i], minimum(stress[:]), maximum(stress[:]), 1, 0),0.0, mapp(stress[i], minimum(stress[:]), maximum(stress[:]), 0, 1)),
linewidth = X[i]*scale)
......@@ -827,14 +827,14 @@ function plotTruss(problem,X,scale,threshold=0)
# p=annotate!((xn[1,Int(ien[1,i])]+xn[1,Int(ien[2,i])])/2.0, (xn[2,Int(ien[1,i])]+xn[2,Int(ien[2,i])])/2.0, "$(floor(stress[i]/1000))*e3", 8)
end
else
p=plot!([xn[1,Int(ien[1,i])],xn[1,Int(ien[2,i])]],
p=Plots.plot!([xn[1,Int(ien[1,i])],xn[1,Int(ien[2,i])]],
[xn[2,Int(ien[1,i])],xn[2,Int(ien[2,i])]],label="",
color=RGB(1.0,1.0,1.0),
linewidth = 0.0)
end
end
# plot!(axis=nothing,ticks=nothing, border=nothing)
# Plots.plot!(axis=nothing,ticks=nothing, border=nothing)
p
end
# =======================================================================
......@@ -855,23 +855,23 @@ function plotTrussDeformed(problem,X,scale,threshold=0,exageration=100.0)
nel=length(X)
E,f,g,idb,ien,ndf,nel,nen,nnp,nsd,xn,len,Ke,Te=problem;
K,F,d,stress,dcomp,g=FEM_truss(problem,X);
p=plot(axis=nothing,ticks=nothing, border=nothing, aspect_ratio=:equal)
p=Plots.plot(axis=nothing,ticks=nothing, border=nothing, aspect_ratio=:equal)
xnn=zeros(size(xn))
xnn.=xn .+ exageration.*dcomp
for i in 1:nel
if X[i]>threshold
if threshold>0
p=plot!([xn[1,Int(ien[1,i])],xn[1,Int(ien[2,i])]],
p=Plots.plot!([xn[1,Int(ien[1,i])],xn[1,Int(ien[2,i])]],
[xn[2,Int(ien[1,i])],xn[2,Int(ien[2,i])]],label="",
color=RGB(0.0,0.0,0.0),
linewidth = 3.0)
p=plot!([xnn[1,Int(ien[1,i])],xnn[1,Int(ien[2,i])]],
p=Plots.plot!([xnn[1,Int(ien[1,i])],xnn[1,Int(ien[2,i])]],
[xnn[2,Int(ien[1,i])],xnn[2,Int(ien[2,i])]],label="",
color=RGB(0.0,1.0,0.0),
linewidth = 3.0)
else
p=plot!([xn[1,Int(ien[1,i])],xn[1,Int(ien[2,i])]],
p=Plots.plot!([xn[1,Int(ien[1,i])],xn[1,Int(ien[2,i])]],
[xn[2,Int(ien[1,i])],xn[2,Int(ien[2,i])]],label="",
color=RGB(mapp(stress[i], minimum(stress[:]), maximum(stress[:]), 1, 0),0.0, mapp(stress[i], minimum(stress[:]), maximum(stress[:]), 0, 1)),
linewidth = X[i]*scale)
......@@ -879,7 +879,7 @@ function plotTrussDeformed(problem,X,scale,threshold=0,exageration=100.0)
# p=annotate!((xn[1,Int(ien[1,i])]+xn[1,Int(ien[2,i])])/2.0, (xn[2,Int(ien[1,i])]+xn[2,Int(ien[2,i])])/2.0, "$(floor(stress[i]/1000))*e3", 8)
end
else
p=plot!([xn[1,Int(ien[1,i])],xn[1,Int(ien[2,i])]],
p=Plots.plot!([xn[1,Int(ien[1,i])],xn[1,Int(ien[2,i])]],
[xn[2,Int(ien[1,i])],xn[2,Int(ien[2,i])]],label="",
color=RGB(1.0,1.0,1.0),
linewidth = 0.0)
......@@ -894,25 +894,25 @@ function plotTrussDeformed3D(problem,X,scale,threshold=0,exageration=100.0)
nel=length(X)
E,f,g,idb,ien,ndf,nel,nen,nnp,nsd,xn,len,Ke,Te=problem;
K,F,d,stress,dcomp,g=FEM_truss(problem,X);
p=plot(axis=nothing,ticks=nothing, border=nothing, aspect_ratio=:equal)
p=Plots.plot(axis=nothing,ticks=nothing, border=nothing, aspect_ratio=:equal)
xnn=zeros(size(xn))
xnn.=xn .+ exageration.*dcomp
for i in 1:nel
if X[i]>threshold
if threshold>0
p=plot!([xn[1,Int(ien[1,i])],xn[1,Int(ien[2,i])]],
p=Plots.plot!([xn[1,Int(ien[1,i])],xn[1,Int(ien[2,i])]],
[xn[3,Int(ien[1,i])],xn[3,Int(ien[2,i])]],
[xn[2,Int(ien[1,i])],xn[2,Int(ien[2,i])]],label="",
color=RGB(0.0,0.0,0.0),
linewidth = 3.0)
p=plot!([xnn[1,Int(ien[1,i])],xnn[1,Int(ien[2,i])]],
p=Plots.plot!([xnn[1,Int(ien[1,i])],xnn[1,Int(ien[2,i])]],
[xnn[3,Int(ien[1,i])],xnn[3,Int(ien[2,i])]],
[xnn[2,Int(ien[1,i])],xnn[2,Int(ien[2,i])]],label="",
color=RGB(0.0,1.0,0.0),
linewidth = 3.0)
else
p=plot!([xn[1,Int(ien[1,i])],xn[1,Int(ien[2,i])]],
p=Plots.plot!([xn[1,Int(ien[1,i])],xn[1,Int(ien[2,i])]],
[xn[3,Int(ien[1,i])],xn[3,Int(ien[2,i])]],
[xn[2,Int(ien[1,i])],xn[2,Int(ien[2,i])]],label="",
color=RGB(mapp(stress[i], minimum(stress[:]), maximum(stress[:]), 1, 0),0.0, mapp(stress[i], minimum(stress[:]), maximum(stress[:]), 0, 1)),
......@@ -921,7 +921,7 @@ function plotTrussDeformed3D(problem,X,scale,threshold=0,exageration=100.0)
# p=annotate!((xn[1,Int(ien[1,i])]+xn[1,Int(ien[2,i])])/2.0, (xn[2,Int(ien[1,i])]+xn[2,Int(ien[2,i])])/2.0, "$(floor(stress[i]/1000))*e3", 8)
end
else
p=plot!([xn[1,Int(ien[1,i])],xn[1,Int(ien[2,i])]],
p=Plots.plot!([xn[1,Int(ien[1,i])],xn[1,Int(ien[2,i])]],
[xn[3,Int(ien[1,i])],xn[3,Int(ien[2,i])]],
[xn[2,Int(ien[1,i])],xn[2,Int(ien[2,i])]],label="",
color=RGB(1.0,1.0,1.0),
......@@ -929,7 +929,7 @@ function plotTrussDeformed3D(problem,X,scale,threshold=0,exageration=100.0)
end
end
# plot!(axis=nothing,ticks=nothing, border=nothing)
# Plots.plot!(axis=nothing,ticks=nothing, border=nothing)
p
end
......@@ -938,18 +938,18 @@ function plotTruss3D(problem,X,scale,threshold=0)
nel=length(X)
E,f,g,idb,ien,ndf,nel,nen,nnp,nsd,xn,len,Ke,Te=problem;
K,F,d,stress,dcomp,g=FEM_truss(problem,X);
p=plot(axis=nothing,ticks=nothing, border=nothing, aspect_ratio=:equal)
p=Plots.plot(axis=nothing,ticks=nothing, border=nothing, aspect_ratio=:equal)
for i in 1:nel
if X[i]>threshold
if threshold>0
p=plot!([xn[1,Int(ien[1,i])],xn[1,Int(ien[2,i])]],
p=Plots.plot!([xn[1,Int(ien[1,i])],xn[1,Int(ien[2,i])]],
[xn[3,Int(ien[1,i])],xn[3,Int(ien[2,i])]],
[xn[2,Int(ien[1,i])],xn[2,Int(ien[2,i])]],label="",
color=RGB(0.0,0.0,0.0),
linewidth = 3.0)
else
p=plot!([xn[1,Int(ien[1,i])],xn[1,Int(ien[2,i])]],
p=Plots.plot!([xn[1,Int(ien[1,i])],xn[1,Int(ien[2,i])]],
[xn[3,Int(ien[1,i])],xn[3,Int(ien[2,i])]],
[xn[2,Int(ien[1,i])],xn[2,Int(ien[2,i])]],label="",
color=RGB(mapp(stress[i], minimum(stress[:]), maximum(stress[:]), 1, 0),0.0, mapp(stress[i], minimum(stress[:]), maximum(stress[:]), 0, 1)),
......@@ -958,7 +958,7 @@ function plotTruss3D(problem,X,scale,threshold=0)
# p=annotate!((xn[1,Int(ien[1,i])]+xn[1,Int(ien[2,i])])/2.0, (xn[2,Int(ien[1,i])]+xn[2,Int(ien[2,i])])/2.0, "$(floor(stress[i]/1000))*e3", 8)
end
else
p=plot!([xn[1,Int(ien[1,i])],xn[1,Int(ien[2,i])]],
p=Plots.plot!([xn[1,Int(ien[1,i])],xn[1,Int(ien[2,i])]],
[xn[3,Int(ien[1,i])],xn[3,Int(ien[2,i])]],
[xn[2,Int(ien[1,i])],xn[2,Int(ien[2,i])]],label="",
color=RGB(1.0,1.0,1.0),
......@@ -966,7 +966,7 @@ function plotTruss3D(problem,X,scale,threshold=0)
end
end
# plot!(axis=nothing,ticks=nothing, border=nothing)
# Plots.plot!(axis=nothing,ticks=nothing, border=nothing)
p
end
......
......@@ -165,17 +165,178 @@ function getDataFromSetup3D(setup,scale)
return E,f,g,idb,ien,ndf,nel,nen,nnp,nsd,xn,len,Ke,Te,Ls
end
# =======================================================================
function getDataFromSetup2D(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 = 2; # number of spacial dimensions
ndf = 2; # 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:nsd, i] = [(nodes[i]["position"]["x"]/scale) (nodes[i]["position"]["y"]/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)
len[i] = sqrt((x1-x2)^2+(y1-y2)^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
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
#inverter
f[1,i] = nodes[i]["force"]["x"]*P;
f[2,i] = nodes[i]["force"]["y"]*P;
end
Ls=[]
for i in 1:nnp
if (nodes[i]["fixedDisplacement"]["x"]!=0)||(nodes[i]["fixedDisplacement"]["y"]!=0)
L=zeros(ndf,nnp);
L[1,i] = nodes[i]["fixedDisplacement"]["x"];
L[2,i] = nodes[i]["fixedDisplacement"]["y"];
append!(Ls,[L]);
end
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
function getSensitivites(Ke0,dcomp,ien,nel,len,λ)
return E,f,g,idb,ien,ndf,nel,nen,nnp,nsd,xn,len,Ke,Te,Ls
end
# =======================================================================
function getSensitivites(Ke0,dcomp,ien,nel,len,λ,twoD=false)
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])]]
λe=[λ[1,Int(ien[1,i])],λ[2,Int(ien[1,i])],λ[3,Int(ien[1,i])] ,λ[1,Int(ien[2,i])] ,λ[2,Int(ien[2,i])],λ[3,Int(ien[2,i])] ]
if twoD
de=[dcomp[1,Int(ien[1,i])],dcomp[2,Int(ien[1,i])] ,dcomp[1,Int(ien[2,i])] ,dcomp[2,Int(ien[2,i])] ]
λe=[λ[1,Int(ien[1,i])],λ[2,Int(ien[1,i])] ,λ[1,Int(ien[2,i])] ,λ[2,Int(ien[2,i])] ]
else
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])]]
λe=[λ[1,Int(ien[1,i])],λ[2,Int(ien[1,i])],λ[3,Int(ien[1,i])] ,λ[1,Int(ien[2,i])] ,λ[2,Int(ien[2,i])],λ[3,Int(ien[2,i])] ]
end
# println(de)
# println(Ke0[:,:,i])
dKedA=Ke0[:,:,i]
......@@ -187,15 +348,21 @@ end
# =======================================================================
function getSensitivitesSIMP(Ke0,dcomp,ien,nel,len,λ,η,X)
function getSensitivitesSIMP(Ke0,dcomp,ien,nel,len,λ,η,X,twoD=false)
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])]]
λe=[λ[1,Int(ien[1,i])],λ[2,Int(ien[1,i])],λ[3,Int(ien[1,i])] ,λ[1,Int(ien[2,i])] ,λ[2,Int(ien[2,i])],λ[3,Int(ien[2,i])] ]
if twoD
de=[dcomp[1,Int(ien[1,i])],dcomp[2,Int(ien[1,i])] ,dcomp[1,Int(ien[2,i])] ,dcomp[2,Int(ien[2,i])] ]
λe=[λ[1,Int(ien[1,i])],λ[2,Int(ien[1,i])] ,λ[1,Int(ien[2,i])] ,λ[2,Int(ien[2,i])] ]
else
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])]]
λe=[λ[1,Int(ien[1,i])],λ[2,Int(ien[1,i])],λ[3,Int(ien[1,i])] ,λ[1,Int(ien[2,i])] ,λ[2,Int(ien[2,i])],λ[3,Int(ien[2,i])] ]
end
# println(de)
# println(Ke0[:,:,i])
# dKedA=Ke0[:,:,i]
......@@ -206,18 +373,25 @@ function getSensitivitesSIMP(Ke0,dcomp,ien,nel,len,λ,η,X)
return dfdA,dgdA
end
function getAdjoint(K,dcomp,L,free)
# L1=L[:,free]
L1=L[free]
function getAdjoint(K,dcomp,L,free,twoD=false)
if twoD
L1=L[:,free]
else
L1=L[free]
end
λ1=K\L1[:]
λ=zeros(size(dcomp))
# λ[:,free].=reshape(λ1,size(L1))
λ[free].=reshape(λ1,size(L1))
if twoD
λ[:,free].=reshape(λ1,size(L1))
else
λ[free].=reshape(λ1,size(L1))
end
return λ,L1,λ1
end
# =======================================================================
function optimizeCompliantMechanism1(problem,Ls,free,dmax,totalVolFactor,maxeval=500)
function optimizeCompliantMechanism1(problem,Ls,free,dmax,totalVolFactor,maxeval=500,twoD=false)
E,f,g,idb,ien,ndf,nel,nen,nnp,nsd,xn,len,Ke,Te=problem
......@@ -228,9 +402,9 @@ function optimizeCompliantMechanism1(problem,Ls,free,dmax,totalVolFactor,maxeval
function FA(x::Vector, grad::Vector)
K,F,d,stress,dcomp,g=FEM_truss(problem,x);
λ,L1,λ1=getAdjoint(copy(K),copy(dcomp),copy(Ls[1]),copy(free))
λ,L1,λ1=getAdjoint(copy(K),copy(dcomp),copy(Ls[1]),copy(free),twoD)
grad[:] .=getSensitivites(Ke,dcomp,ien,nel,len,copy(λ))[2]
grad[:] .=getSensitivites(Ke,dcomp,ien,nel,len,copy(λ),twoD)[2]
# g=(L1[:]'*d)[1] .- dmax - (λ1'*(K*d -F))[1]
......@@ -265,7 +439,7 @@ function optimizeCompliantMechanism1(problem,Ls,free,dmax,totalVolFactor,maxeval
end
# =======================================================================
function optimizeCompliantMechanism(problem,Ls,free,dmax,amax,maxeval=500,SIMP=false)
function optimizeCompliantMechanism(problem,Ls,free,dmax,amax,maxeval=500,SIMP=false,twoD=false)
E,f,g,idb,ien,ndf,nel,nen,nnp,nsd,xn,len,Ke,Te=problem
......@@ -276,9 +450,9 @@ function optimizeCompliantMechanism(problem,Ls,free,dmax,amax,maxeval=500,SIMP=f
function GA(x::Vector, grad::Vector,num::Int)
K,F,d,stress,dcomp,g=FEM_truss(problem,x);
λ,L1,λ1=getAdjoint(copy(K),copy(dcomp),copy(Ls[num]),copy(free))
λ,L1,λ1=getAdjoint(copy(K),copy(dcomp),copy(Ls[num]),copy(free),twoD)
grad[:] .=getSensitivites(Ke,dcomp,ien,nel,len,copy(λ))[2]
grad[:] .=getSensitivites(Ke,dcomp,ien,nel,len,copy(λ),twoD)[2]
# g=(L1[:]'*d)[1] .- dmax - (λ1'*(K*d -F))[1]
......@@ -292,9 +466,9 @@ function optimizeCompliantMechanism(problem,Ls,free,dmax,amax,maxeval=500,SIMP=f
function GASIMP(x::Vector, grad::Vector,num::Int)
K,F,d,stress,dcomp,g=FEM_truss(problem,x.^η);
λ,L1,λ1=getAdjoint(copy(K),copy(dcomp),copy(Ls[num]),copy(free))
λ,L1,λ1=getAdjoint(copy(K),copy(dcomp),copy(Ls[num]),copy(free),twoD)
grad[:] .=getSensitivitesSIMP(Ke,dcomp,ien,nel,len,copy(λ),η,x)[2]
grad[:] .=getSensitivitesSIMP(Ke,dcomp,ien,nel,len,copy(λ),η,x,twoD)[2]
# g=(L1[:]'*d)[1] .- dmax - (λ1'*(K*d -F))[1]
......@@ -315,6 +489,8 @@ function optimizeCompliantMechanism(problem,Ls,free,dmax,amax,maxeval=500,SIMP=f
return (sum(x.^η .* len ))
end
display(FA(ones(nel)*0.5,ones(nel)*0.5))
display(GA(ones(nel)*0.5,ones(nel)*0.5,1))
opt = Opt(:LD_MMA, nel)
opt.lower_bounds = fill(1e-6,nel)
......
......@@ -54,6 +54,7 @@ include("./problems.jl")
include("./concurrent2D.jl")
include("./concurrent3D.jl")
include("./microstructureDesignUe.jl")
include("./fabrication.jl")
......
......@@ -3,14 +3,13 @@
#####################################FEA KE####################################################################
function lk()
function lk(E=1,nu=0.3)
# 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]);
E=1
nu=0.3
k=[1/2-nu/6,1/8+nu/8,-1/4-nu/12,-1/8+3*nu/8,-1/4+nu/12,-1/8-nu/8,nu/6,1/8-3*nu/8]
KE = E/(1-nu^2)*[ k[0+1] k[1+1] k[2+1] k[3+1] k[4+1] k[5+1] k[6+1] k[7+1];
k[1+1] k[0+1] k[7+1] k[6+1] k[5+1] k[4+1] k[3+1] k[2+1];
......
function MicroTop2DU(θ, Micro_struct,prob, penal,saveItr=5,maxloop = 200,connectivity=false)
# Amira Abdel-Rahman
# (c) Massachusetts Institute of Technology 2021
# use dispAnimation2D to display the macrostructure deformation of resultant microstructure
function MicroTop2DU(Ue,Micro_struct, penal,saveItr=5,maxloop = 200,fab=false)
## USER-DEFINED LOOP PARAMETERS