# =======================================================================
function data_FEM(nelx,nely)
# data input for 2D continuum elements
# Author: JV Carstensen, CEE, MIT (JK Guest, Civil Eng, JHU)
# Revised: Aug 22 2017, JVC
# Revised: ADD YOUR NAME HERE
#
# -----------------------------------------------------------------------
# E = Young's modulus
# 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
# iplane = 1 - plane strain, 2 - plane stress
# 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
# nu = Poisson's ratio
# t = thinkness
# xn(nsd,nnp) = nodal coordinates
# =======================================================================
# ---- DIMENSIONS OF THE PROBLEM ----------------------------------------
lx = nelx;
ly = nely;
# nelx = 80;
# nely = 50;
# nelx = 40;
# nely = 26;
nnx = nelx+1; # number of elements in x
nny = nely+1; # number of elements in y
nel = nelx*nely; # number of elements
nnp = nnx*nny; # number of nodal points
# ---- MESH -------------------------------------------------------------
nsd = 2; # number of spacial dimensions
ndf = 2; # number of degrees of freedom per node
nen = 4; # number of element nodes
# ---- NODAL COORDINATES ------------------------------------------------
# xn(i,N) = coordinate i for node N
# N = 1,...,nnp
# i = 1,...,nsd
# -----------------------------------------------------------------------
# ---- ELEMENT CONNECTIVITY ---------------------------------------------
# ien(a,e) = N
# N: global node number - N = 1,...,nnp
# e: element number - e = 1,...,nel
# a: local node number - a = 1,...,nen
# -----------------------------------------------------------------------
# generate mesh
xn,ien = generate_mesh_quad4(nsd,nen,nel,nnp,nnx,nny,lx,ly);
# ---- MATERIAL PROPERTIES ----------------------------------------------
# E = Youngs modulus
# nu = Poissons ratio
# -----------------------------------------------------------------------
E = 1; # Young's modulus
nu = 0.3; # Poisson's ratio
iplane = 2; # iplane = 1 - plane strain, 2 - plane stress
# ---- GEOMETRIC PROPERTIES ---------------------------------------------
# t = element thickness
# -----------------------------------------------------------------------
t = 1;
# ---- 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);
idb[1:2,1:nnx:(nely*nnx+1)] .= 1;
# ---- 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
# -----------------------------------------------------------------------
f = zeros(ndf,nnp);
f[2,Int(round(nely/2)+1)*nnx] = -1.0;
# ---- NUMBER THE EQUATIONS ---------------------------------------------
id,neq = number_eq(idb,ndf,nnp);
# ---- FORM THE ELEMENT STIFFNESS MATRICES ------------------------------
nee = ndf*nen; # number of element equations
Ke0 = zeros(nee,nee,nel);
Imx_nsd = Matrix(1I, nsd, nsd);
zero_nsd = zeros(nsd,nsd);
for i = 1:nel
Ke0[:,:,i] .= Ke_quad4(iplane,E,nu,t,nee,nen,nsd,ndf,ien[:,i],xn,Imx_nsd,zero_nsd);
end
return E,f,g,idb,ien,iplane,ndf,nel,nen,nnp,nsd,nu,t,xn,Ke0
end
# =======================================================================
# =======================================================================
function generate_mesh_quad4(nsd,nen,nel,nnp,nnx,nny,lx,ly)
# function that generates a lattice mesh for quad4 elements
# -----------------------------------------------------------------------
# nsd = number of spacial dimensions
# nen = number of nodes per element
# nel = number of elements
# nnp = number of nodes
# nnx = number of nodes in x
# nny = number of nodes in y
# lx = length of domain in x
# ly = length of domain in y
#
# xn(nsd,nnp) = nodal coordinates
# ien(nen,nel) = element connectivity
#
# =======================================================================
# coordinates
xn = zeros(nsd,nnp);
nelx = nnx-1;
nely = nny-1;
for i = 1:nnx
xn[1,i] = (i-1)*(lx/nelx);
end
for i = 2:nny
loc = (1:nnx) .+ (i-1).*nnx;
xn[1,loc] .= xn[1,1:nnx];
xn[2,loc] .= (i-1)*(ly/nely);
end
# plot(xn(1,:),xn(2,:),'*')
# axis equal
# connectivity
ien=zeros(Int,nen,nel);
for i = 1:nelx
ien[:,i] = [0 1 1+nnx nnx]' .+i;
end
for i = 2:nely
loc = (1:nelx) .+ (i-1)*nelx;
ien[:,loc] .= ien[:,1:nelx] .+(i-1)*nnx;
end
return xn,ien
end
# =======================================================================
function getSensitivity(Ke,de,nel)
dfdA=zeros(nel)
for i in 1:nel
dfdA[i]=(-1.0 .* de[:,i])' * Ke[:,:,i] * de[:,i]
end
return dfdA
end
function getSensitivitySIMP(Ke,de,nel,x,p)
dfdA=zeros(nel)
for i in 1:nel
dfdA[i]=(-1.0 .* de[:,i])' * (p .*(x[i]).^(p-1) .*Ke[:,:,i]) * de[:,i]
end
return dfdA
end
function optimizeDensity(problem,totalVolFactor,maxEval=500)
E,f,g,idb,ien,iplane,ndf,nel,nen,nnp,nsd,snu,t,xn,Ke0=problem
function FA(x::Vector, grad::Vector)
K,F,d,de,Fe=FEM(problem, x .+1e-4);
# display(F)
grad[:] .=getSensitivity(Ke0,de,nel)
Fx=F'*d -d'*K*d +d'*F
# display(Fx)
# display(sum(grad))
return Fx[1]
end
function G(x::Vector, grad::Vector)
grad[:] .=1.0
return (sum(x ) - totalVolFactor*nel)
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-4,nel)
opt.upper_bounds = fill(1,nel)
opt.xtol_rel = 0.001
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).*totalVolFactor))
numevals = opt.numevals # the number of function evaluations
display("got $minf after $numevals iterations (returned $ret)")
return minx
end
function optimizeDensitySIMP(problem,totalVolFactor,maxEval=500,p=3)
E,f,g,idb,ien,iplane,ndf,nel,nen,nnp,nsd,snu,t,xn,Ke0=problem
function FA(x::Vector, grad::Vector)
K,F,d,de,Fe=FEM(problem, ((x).^p .+1e-4));
# display(F)
grad[:] .=getSensitivitySIMP(Ke0,de,nel,x,p)
Fx=F'*d -d'*K*d +d'*F
# display(Fx)
# display(sum(grad))
return Fx[1]
end
function G(x::Vector, grad::Vector)
grad[:] .=1.0
return (sum(x ) - totalVolFactor*nel)
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-8,nel)
opt.upper_bounds = fill(1,nel)
opt.xtol_rel = 0.001
opt.maxeval = maxEval
opt.min_objective = FA
inequality_constraint!(opt, (x,gg) -> G(x,gg), 1e-4)
display(@time (minf,minx,ret) = optimize(opt, ones(nel).*totalVolFactor))
numevals = opt.numevals # the number of function evaluations
display("got $minf after $numevals iterations (returned $ret)")
return minx
end
# =======================================================================