function FEM(data_FEM,ρ)
# FEM solver for 2D elements
# Author: JV Carstensen, CEE, MIT (JK Guest, Civil Eng, JHU)
# Revised: Oct 26 2017, JVC
# Revised: Amira Abdel-Rahman
# =======================================================================
# ---- READ IN DATA -----------------------------------------------------
# read in external file
E,f,dis,idb,ien,iplane,ndf,nel,nen,nnp,nsd,snu,t,xn,Ke0 = data_FEM;
g=zeros(size(dis))
g.=dis
# ---- NUMBER THE EQUATIONS ---------------------------------------------
id,neq = number_eq(idb,ndf,nnp);
# ---- FORM THE ELEMENT STIFFNESS MATRICES ------------------------------
nee = ndf*nen; # number of element equations
Ke = zeros(nee,nee,nel);
Imx_nsd = Matrix(1I, nsd, nsd);
zero_nsd = zeros(nsd,nsd);
for i = 1:nel
Ke[:,:,i] .= (ρ[i]).*Ke_quad4(iplane,E,snu,t,nee,nen,nsd,ndf,ien[:,i],xn,Imx_nsd,zero_nsd);
end
# ---- PERFORM GLOBAL TO LOCAL MAPPING ----------------------------------
LM = zeros(nee,nel);
for i = 1:nel
LM[:,i] = get_local_id(id,ien[:,i],ndf,nee,nen);
end
# ---- IF THERE IS FREE DEGREES OF FREEDOM, THEN SOLVE THE EQUILIBRIUM --
if (neq > 0)
# get global force vector (line 275)
F = globalF(f,g,id,ien,Ke,LM,ndf,nee,nel,nen,neq,nnp);
# solve Kd - F = 0 (line 521)
d,K = solveEQ(F,LM,Ke,nee,nel,neq);
end
# display(Ke)
de,Fe=post_processing(d,E,g,id,ien,Ke,ndf,nee,nel,nen,nnp)
return K,F,d,de,Fe
end
# =======================================================================
# =======================================================================
function add_d2dcomp(dcomp,d,id,ndf,nnp)
# function that adds the displacements of the free degrees of freedom to
# the nodal displacements
# -----------------------------------------------------------------------
# dcomp(ndf,nnp) = nodal displacements
# d(neq,1) = displacement at free degrees of freedom
# id(ndf,nnp) = equation numbers of degrees of freedom
# ndf = number of degrees of freedom per node
# nnp = number of nodal points
#------------------------------------------------------------------------
# loop over nodes and degrees of freedom
for n=1:nnp
for i=1:ndf
# if it is a free dof then add the global displacement
if (id[i,n]>0)
dcomp[i,n] = dcomp[i,n]+d[Int(id[i,n])];
end
end
end
return dcomp
end
# =======================================================================
# =======================================================================
function add_loads_to_force(F,f,id,ndf,nnp)
# function that adds nodal forces to the global force vector
# -----------------------------------------------------------------------
# F(neq,1) = global force vector
# f(ndf,nnp) = prescribed nodal forces
# id(ndf,nnp) = equation numbers of degrees of freedom
# ndf = number of degrees of freedom per node
# nnp = number of nodal points
#------------------------------------------------------------------------
# loop over nodes and degrees of freedom
for n = 1:nnp
for i = 1:ndf
# get the global equation number
M = Int(id[i,n]);
# if free degree of freedom, then add nodal load to global force
# vector
if (M > 0)
F[M] = F[M]+f[i,n];
end
end
end
return F
end
# =======================================================================
# =======================================================================
function addforce(F,Fe,LM,nee)
# function that adds element forces to the global force vector
# -----------------------------------------------------------------------
# F(neq,1) = global force vector
# Fe(nee,1) = element force vector
# LM(nee,nel) = global to local map for the element
# nee = number of element equations
# =======================================================================
# loop over rows of Fe
for i = 1:nee
# get the global equation number for local equation i
M = Int(LM[i]);
# if free dof (eqn number > 0) add to F vector
if (M > 0)
F[M]=F[M]+Fe[i];
end
end
return F
end
# =======================================================================
# =======================================================================
function addstiff(K,Ke,LM,nee)
# function that adds the element stiffness matrices to the global
# stiffness matrix
# -----------------------------------------------------------------------
# K(neq,neq) = global stiffness matrix
# Ke(nee,nee,1) = element stiffness matrix
# LM(nee,nel) = global to local map for the element
# nee = number of element equations
# =======================================================================
# loop over rows of Ke
for i = 1:nee
# loop over columns of Ke
for j = 1:nee
Mr = Int(LM[i]);
Mc = Int(LM[j]);
if (Mr > 0 && Mc > 0)
# if equation #'s are non-zero add element contribution to the
# stiffness matrix
K[Mr,Mc] = K[Mr,Mc]+Ke[i,j];
end
end
end
return K
end
# =======================================================================
# =======================================================================
function B_2d_elastic(dNdx,dNdy,nen,ndf)
# Computes the strain-displacement matrix for an elastic 2d problem
# -----------------------------------------------------------------------
# B(3,nen*ndf) = strain displacement matrix at current gauss point
# dNdx = derivative of shape function at current gauss point
# dNdy = derivative of shape function at current gauss point
# nen = number of nodes per element
# ndf = number of degrees of freedom per node
#
# =======================================================================
# Initialize
B = zeros(3,nen*ndf);
# Compute B
for i = 1:nen
B[:,(i-1)*ndf+1:ndf*i] = [ dNdx[i] 0
0 dNdy[i]
dNdy[i] dNdx[i] ];
end
# =======================================================================
return B
end
# =======================================================================
# =======================================================================
function D_2d(E,snu,iplane,nstr)
# Computes the elastic constitutive matrix for a 2d problem
# -----------------------------------------------------------------------
# E = Young's modulus
# snu = Poisson's ratio
# iplane = 1 - plane strain, 2 - plane stress
# nstr = number of independent stress components
#
# D(nstr,nstr) = elastic constitutive matrix for element e
#
# =======================================================================
# initialize and define D
D = zeros(Int(nstr),Int(nstr));
if (iplane == 2) # plane stress
coeff = E/(1-(snu*snu));
D[1,1] = coeff;
D[2,2] = D[1,1];
D[1,2] = coeff*snu;
D[2,1] = D[1,2];
D[3,3] = coeff*(1-snu)/2;
elseif (iplane == 1) # plane strain
coeff = E/((1+snu)*(1-2*snu));
D[1,1] = coeff*(1-snu);
D[2,2] = D[1,1];
D[1,2] = coeff*snu;
D[2,1] = D[1,2];
D[3,3] = coeff*(1-2*snu)/2;
end
return D
end
# =======================================================================
# =======================================================================
function get_de_from_dcomp(dcomp,ien,ndf,nen)
# extracts element displacement vector from complete displacement vector
#------------------------------------------------------------------------
# dcomp(ndf,nnp) = nodal displacements
# ien(nen,1) = element connectivity
# ndf = number of degrees of freedom per node
# nen = number of element equations
#
# de(nen,1) = element displacements
#
#------------------------------------------------------------------------
de = zeros((nen-1)*ndf+ndf,1);
# loop over number of element nodes
for i = 1:nen
# loop over number of degrees of freedom per node
for j = 1:ndf
# get the local element number and place displacement in de
leq = (i-1)*ndf+j;
de[leq] = dcomp[j,ien[i]];
end
end
return de
end
# =======================================================================
# =======================================================================
function gausspoints_quad4(ndf,nsd)
# function that defines integration points and weigths for quad4 elements
#------------------------------------------------------------------------
# ndf = number of degrees of freedom per node
# nsd = number of spacial dimensions
#
# Nint = number of gauss points in x and y
# point = location of gauss points in x and y
# weight = weigthing of gauss points for numerical integration
# nstre = number of independent stress components
# nstr1 = total number of stress components
# ngp = number of gauss points per element
#
#------------------------------------------------------------------------
Nint=zeros(Int,2)
Nint[1] = 2; # number of gauss points in x
Nint[2] = 2; # number of gauss points in y
# total number of gauss points
ngp = Nint[1]*Nint[2];
# Initialize
point = zeros(Nint[2],Nint[1]);
weight = zeros(Nint[2],Nint[1]);
weight .= point;
# get gauss point locations and weights for integration
point[1,1] = (-0.577350269189626);
point[1,2] = (-0.577350269189626);
point[2,1] = ( 0.577350269189626);
point[2,2] = ( 0.577350269189626);
weight .= weight .+1;
# define number of individual stress/strain components
nstre = nsd*(ndf+1)/2;
# define number of stress components in each element (4 in 2d)
nstr1 = 4;
return Nint,point,weight,nstre,nstr1,ngp
end
# =======================================================================
# =======================================================================
function get_local_id(id,ien,ndf,nee,nen)
# functions that performs the global to local mapping of equation numbers
# -----------------------------------------------------------------------
# id(ndf,nnp) = equation numbers of degrees of freedom
# ien(nen,1) = element connectivity
# ndf = number of degrees of freedom per node
# nee = number of element equations
# nen = number of element equations
#
# LM(nee,1) = global to local map for element
# =======================================================================
# initialize global-local mapping matrix
LM = zeros(nee,1);
# initialize local equation number counter
k = 0;
# loop over element nodes
for i = 1:nen
# loop over degrees of freedom at each node
for j = 1:ndf
# update counter and prescribe global equation number
k = k+1;
LM[k] = id[j,ien[i]];
end
end
return LM
end
# =======================================================================
# =======================================================================
function globalF(f,g,id,ien,Ke,LM,ndf,nee,nel,nen,neq,nnp)
# function that assembles the global load vector
# -----------------------------------------------------------------------
# id(ndf,nnp) = equation numbers of degrees of freedom
# f(ndf,nnp) = prescribed nodal forces
# g(ndf,nnp) = prescribed nodal displacements
# ien(nen,nel) = element connectivities
# Ke(nee,nee,nel) = element stiffness matrices
# ndf = number of degrees of freedom per node
# nee = number of element equations
# nel = number of elements
# nen = number of element equations
# neq = number of equations
# nnp = number of nodal points
#
# F(neq,1) = global force vector
# =======================================================================
# initialize
F = zeros(neq,1);
# Insert applied loads into F
F = add_loads_to_force(F,f,id,ndf,nnp);
# Compute forces from applied displacements (ds~=0) and insert into F
Fse = zeros(nee,nel);
# loop over elements
for i = 1:nel
# get dse for current element
dse = get_de_from_dcomp(g,ien[:,i],ndf,nen);
# compute element force
Fse[:,i] = -Ke[:,:,i]*dse;
# assemble elem force into global force vector
F = addforce(F,Fse[:,i],LM[:,i],nee);
end
return F
end
# =======================================================================
# =======================================================================
function jacobian_2d(dNdr,dNds,xn,nen,Imx,Jaco)
# computes the Jacobian, its determinate and inverse at current gauss
# point in element e
# -----------------------------------------------------------------------
# dNdr = derivative of shape function at current gauss point
# dNds = derivative of shape function at current gauss point
# xn(nsd,nen) = nodal coordinates for element e
# nen = number of nodes per element
# Imx = eye(nsd)
# Jaco = zeros(nsd)
#
# Jaco = Jacobian matrix for current gauss point
# detJ = determinant of the Jacobian for current gauss point
# InvJ = inverse of the Jacobian matrix for current gauss point
# =======================================================================
# Add to Jacobian
for i = 1:nen
#for j = 1:nsd
Jaco[1,1] = Jaco[1,1]+dNdr[i]*xn[1,i];
Jaco[1,2] = Jaco[1,2]+dNdr[i]*xn[2,i];
Jaco[2,1] = Jaco[2,1]+dNds[i]*xn[1,i];
Jaco[2,2] = Jaco[2,2]+dNds[i]*xn[2,i];
#end
end
# Find its determinant
detJ = det(Jaco);
# Find its inverse
invJ = Imx/Jaco;
return Jaco,detJ,invJ
end
# =======================================================================
# =======================================================================
function ke_elem_quad4(Nint,point,weight,ien, xn,ndf,nen,nee,t,D,Imx_nsd,zero_nsd)
# loops over element gauss points and computes element stiffness matrix
# -----------------------------------------------------------------------
# Nint = number of gauss points in x and y
# point = location of gauss points
# weight = weigthing of gauss points in numerical integration
# ien(nen,nel) = element connectivity
# xn(nsd,nnp) = nodal coordinates
# nen = number of nodes per element
# nee = number of element equations
# t = element thickness
# D = elastic constititive matrix for the element
# Imx = eye(nsd)
# Jaco = zeros(nsd)
#
# kee(nee,nee,1) = element stiffness matrix for quad 4 element
# =======================================================================
# =======================================================================
# Initialize
kee = zeros(nee,nee);
coeff = 0;
for i = 1:Nint[1]
ri = point[i,1];
wi = weight[i,1];
for j = 1:Nint[2]
sj = point[j,2];
wj = weight[j,2];
sN,dNdx,dNdy,detJ,invJ = shape_quad4(ri,sj,xn,ien,nen,Imx_nsd,zero_nsd);
B = B_2d_elastic(dNdx,dNdy,nen,ndf);
coeff = t*wi*wj*detJ;
# Compute kee = kee + transpose(B)*D*B*coeff
kee .= kee .+(B')*D*B*coeff;
end
end
return kee
end
# =======================================================================
# =======================================================================
function Ke_quad4(iplane,E,snu,t,nee,nen,nsd,ndf,ien,xn, Imx_nsd,zero_nsd)
# computes element stiffness matrix
# -----------------------------------------------------------------------
# iplane = 1 - plane strain, 2 - plane stress
# E = Young's modulus
# snu = Poisson's ratio
# t = element thickness
# nee = number of element equations
# nen = number of nodes per element
# nsd = number of spacial dimmensions
# ndf = number of degrees of freedom
# weight = weigthing of gauss points in numerical integration
# ien(nen,nel) = element connectivity
# xn(nsd,nnp) = nodal coordinates
# Imx = eye(nsd)
# Jaco = zeros(nsd)
#
# kee(nee,nee,1) = element stiffness matrix for quad 4 element
# =======================================================================
# Define location and weights of integration points
Nint,point,weight,nstre,nstr1,ngp = gausspoints_quad4(ndf,nsd);
# get constitutive matrix
nstr = nsd*(nsd+1)/2;
D = D_2d(E,snu,iplane,nstr);
# Compute element matrix
kee = ke_elem_quad4(Nint,point,weight,ien,xn,ndf,nen,nee,t,D,Imx_nsd,zero_nsd);
return kee
end
# =======================================================================
# =======================================================================
function number_eq(idb,ndf,nnp)
# function that numbers the unknown degrees of freedom (equations)
# -----------------------------------------------------------------------
# idb(ndf,nnp) = 1 if the degree of freedom is prescribed, 0 otherwise
# ndf = number of degrees of freedom per node
# nnp = number of nodal points
#
# id(ndf,nnp) = equation numbers of degrees of freedom
# neq = number of equations (tot number of degrees of freedom)
# =======================================================================
# initialize id and neq
id = zeros(ndf,nnp);
neq = 0;
# loop over nodes
for n = 1:nnp
# loop over degrees of freedom
for i = 1:ndf
if idb[i,n] == 0
# udate # of equations
neq = neq + 1;
# if no prescribed displacement at dof i of node n
# => give an equation # different from 0
id[i,n] = neq;
end
end
end
return id,neq
end
# =======================================================================
# =======================================================================
function solveEQ(F,LM,Ke,nee,nel,neq)
# function that solves the equilibrium condition
# -----------------------------------------------------------------------
# F(neq,1) = global force vector
# LM(nee,nel) = global to local maps
# Ke(nee,nee,nel) = element stiffness matrices
# nee = number of element equations
# nel = number of elements
# neq = number of equations
#
# d(neq,1) = displacements at free degrees of freedom
# =======================================================================
# assemble global stiffness matrix
# K = zeros(neq,neq); # Use 'sparse' for more efficient memory usage
K=spzeros(neq,neq)
for i = 1:nel
K = addstiff(K,Ke[:,:,i],LM[:,i],nee);
end
# solve the equlibrium
d = K\F;
return d,K
end
# =======================================================================
# =======================================================================
function shape_quad4(r,s,xn,ien,nen, Imx_nsd,zero_nsd)
# computes the shape functions for quad4 element e at gauss point (r,s)
# -----------------------------------------------------------------------
# r = gauss point coordinate in parent domain
# s = gauss point coordinate in parent domain
# xn(nsd,nnp) = nodal coordinates
# ien(nen,1) = element connectivity for element e
# nen = number of nodes per element
# Imx = eye(nsd)
# Jaco = zeros(nsd)
#
# sN(nen,1) = shape functions at current gauss point
# dNdx(nen,1) = shape function derivatives at current gauss point in
# (x,y)
# dNdy(nen,1) = shape function derivatives at current gauss point in
# (x,y)
# detJ = determinant of the Jacobian for current gauss point
# InvJ = inverse of the Jacobian matrix for current gauss point
#
# =======================================================================
# Initialize
dNdx = zeros(nen,1);
dNdy = zeros(nen,1);
dNdy .= dNdx;
# compute shape functions at(r,s)
sN,dNdr,dNds = shapefunct_quad4(r,s);
# compute Jacobian, its deterninant and inverse
Jaco,detJ,invJ = jacobian_2d(dNdr,dNds,xn[:,ien],nen,Imx_nsd,zero_nsd);
# compute the derivatives of the shape functions
for i = 1:nen
dNdx[i] = invJ[1,1]*dNdr[i]+invJ[1,2]*dNds[i];
dNdy[i] = invJ[2,1]*dNdr[i]+invJ[2,2]*dNds[i];
end
return sN,dNdx,dNdy,detJ,invJ
end
# =======================================================================
# =======================================================================
function shapefunct_quad4(r,s)
# evaluates the shape functions and their derivatives at point (s,r)
# -----------------------------------------------------------------------
# r = gauss point coordinate in parent domain
# s = gauss point coordinate in parent domain
#
# sN(nen,1) = shape functions at current gauss point
# dNdr(nen,1) = shape function derivatives at current gauss point in
# (r,s)
# dNds(nen,1) = shape function derivatives at current gauss point in
# (r,s)
#
# =======================================================================
# Initialize
sN = zeros(4,1);
dNdr = zeros(4,1);
dNds = zeros(4,1);
dNdr .= sN;
dNds .= sN;
# Shape functions
sN[1] = (1-r)*(1-s);
sN[2] = (1+r)*(1-s);
sN[3] = (1+r)*(1+s);
sN[4] = (1-r)*(1+s);
sN .= 0.25*sN;
# Derivatives
dNdr[1] = -(1-s);
dNdr[2] = (1-s);
dNdr[3] = (1+s);
dNdr[4] = -(1+s);
dNdr .= 0.25*dNdr;
dNds[1] = -(1-r);
dNds[2] = -(1+r);
dNds[3] = (1+r);
dNds[4] = (1-r);
dNds .= 0.25*dNds;
return sN,dNdr,dNds
end
# =======================================================================
# =======================================================================
function post_processing(d,E,g,id,ien,Ke,ndf,nee,nel,nen,nnp)
# function that performs post processing for truss elements
# -----------------------------------------------------------------------
# A(nel,1) = cross-sectional area of elements
# d(neq,1) = displacements at free degrees of freedom
# E(nel,1) = Young's modulus of elements
# g(ndf,nnp) = prescribed nodal displacements
# id(ndf,nnp) = equation numbers of degrees of freedom
# ien(nen,nel) = element connectivities
# Ke(nee,nee,nel) = element stiffness matrices
# ndf = number of degrees of freedom per node
# nee = number of element equations
# nel = number of elements
# nen = number of element equations
# nnp = number of nodes
# Te(nee,nee,nel) = element transformation matrices
#
# dcomp(ndf,nnp) = nodal displacements
# axial(nel,1) = axial element forces
# stress(nel,1) = element stresses
# strain(nel,1) = element strains
# Fe(nee,nel) = element forces
# =======================================================================
# get the total displacement of the structure in matrix form dcomp(nsd,nnp)
dcomp = add_d2dcomp(g,d,id,ndf,nnp);
# initalize evaluation of global element forces Fe, local element forces
# fe, axial forces, element stresses and strains
Fe = zeros(nee,nel);
de = zeros(nee,nel);
fe = zeros(1*nen,nel); # element local force vector
axial = zeros(nel,1);
stress = zeros(nel,1);
strain = zeros(nel,1); # element axial, stress, strain
# loop over elements
for i=1:nel
# get the element displacaments
de[:,i] = get_de_from_dcomp(dcomp,ien[:,i],ndf,nen);
# compute the element forces
Fe[:,i] = Ke[:,:,i]*de[:,i];
end
# display("Fe")
# display(Fe)
return de,Fe
end