EEM: Based on the examples in the book "Finite-element method for beams".

By Johan Blaauwendraad http://www.delftacademicpress.nl/f040.php (in Dutch).

This EEM_03 notebook uses the figure 1.1 example in EEM and explains in detail how the FE4_1 skeleton program in the PtFEM package works.

EEM_03 continues where EEM_01 and EEM_02 end. Repeat the initial steps because later on in the explanations we need m.kv, the stiffnes matrix in Skyline format.


In [60]:
using PtFEM, DataFrames, Plots

In [61]:
struct_el = :Rod   # A 1-dimensional structural element that can only handle axial forces.
fin_el = :Line     # The type of finite element used to construct a mesh for the structural element;

In [62]:
l = 1.0       # Total length of structural element [m]
q = 5.0       # Distributed load [N/m]
N = 4        # Number of nodes
els = N - 1   # Number of finite elements
nod = 2       # Number of nodes per finite elements
nodof = 1     # Degrees of freedom for each node, just axial
np_types = 1  # Number of np_types
nip = 1       # Number of numerical integration points
EA = 1.0e5    # Strain stiffness for finite elements;

In [63]:
data = Dict(
  # Rod(nxe, np_types, nip, fin_el(nod, nodof))
  :struc_el => Rod(els, np_types, nip, Line(nod, nodof)),
  :properties => [EA;],
  :x_coords => 0.0:l/els:l,
  :support => [(N, [0])],
  :loaded_nodes => [(i, repeat([0.0], inner=nodof)) for i in 1:N],
  :eq_nodal_forces_and_moments => [(i, repeat([0.0], inner=nodof*nod)) for i in 1:els]
);

In [64]:
for node in 1:els
  data[:eq_nodal_forces_and_moments][node][2][1] = (1/2)*q*l/els
  data[:eq_nodal_forces_and_moments][node][2][2] = (1/2)*q*l/els
end

In [65]:
for node in 1:N-1
  data[:loaded_nodes][node][2][1] += data[:eq_nodal_forces_and_moments][node][2][1]
  data[:loaded_nodes][node+1][2][1] += data[:eq_nodal_forces_and_moments][node][2][2]
end

In [66]:
m = PtFEM.FE4_1(data);


There are 3 equations and the skyline storage is 5.

Here starts the review of all steps in FE4_1

First couple of steps just go over the input data Dictionary and copies the entries to global variables. This is just for demonstration purposes. Inreality all of this is done inside the FE4_1 function.


In [67]:
if :struc_el in keys(data)
    struc_el = data[:struc_el]
end


Out[67]:
PtFEM.Rod(3,1,1,PtFEM.Line(2,1))

In [68]:
ndim = 1
nst = struc_el.np_types;

In [69]:
fin_el = struc_el.fin_el


Out[69]:
PtFEM.Line(2,1)

In [70]:
if typeof(fin_el) == Line # 1D finite element
    (nels, nn) = PtFEM.mesh_size(fin_el, struc_el.nxe)
    elseif typeof(fin_el) == Triangle || typeof(fin_el) == Quadrilateral # 2D finite elements
    (nels, nn) = mesh_size(fin_el, struc_el.nxe, struc_el.nye)
    elseif typeof(fin_el) == Hexahedron  # 3D finite elements
    (nels, nn) = mesh_size(fin_el, struc_el.nxe, struc_el.nye, struc_el.nze)
end


Out[70]:
(3,4)

In [71]:
nodof = fin_el.nodof         # Degrees of freedom per node


Out[71]:
1

In [72]:
ndof = fin_el.nod * nodof    # Degrees of freedom per fin_el, in this case each finite element has 2 nodes


Out[72]:
2

In [73]:
penalty = 1e20               # used to fix fixed_freedoms
if :penalty in keys(data)
    penalty = data[:penalty]
end

In [74]:
if :properties in keys(data)
    prop = zeros(size(data[:properties], 1), size(data[:properties], 2))
    for i in 1:size(data[:properties], 1)
        prop[i, :] = data[:properties][i, :]
    end
end

In [75]:
prop


Out[75]:
1×1 Array{Float64,2}:
 100000.0

In [76]:
nf = ones(Int64, nodof, nn)
if :support in keys(data)
    for i in 1:size(data[:support], 1)
        nf[:, data[:support][i][1]] = data[:support][i][2]
    end
end
(size(nf), nf)


Out[76]:
((1,4),
[1 1 1 0])

In [77]:
x_coords = zeros(nn)
if :x_coords in keys(data)
    x_coords = data[:x_coords]
end
  
y_coords = zeros(nn)
if :y_coords in keys(data)
    y_coords = data[:y_coords]
end
  
z_coords = zeros(nn)
if :z_coords in keys(data)
    z_coords = data[:z_coords]
end

etype = ones(Int64, nels)
if :etype in keys(data)
    etype = data[:etype]
end
x_coords


Out[77]:
0.0:0.3333333333333333:1.0

In [78]:
collect(x_coords)


Out[78]:
4-element Array{Float64,1}:
 0.0     
 0.333333
 0.666667
 1.0     

In [79]:
#
# Initialize dynamic arrays stored in the FEM object
#

points = zeros(struc_el.nip, ndim)
g = zeros(Int64, ndof)
g_coord = zeros(ndim,nn)
fun = zeros(fin_el.nod)
coord = zeros(fin_el.nod, ndim)
gamma = zeros(nels)
jac = zeros(ndim, ndim)
g_num = zeros(Int64, fin_el.nod, nels)
der = zeros(ndim, fin_el.nod)
deriv = zeros(ndim, fin_el.nod)
bee = zeros(nst,ndof)
km = zeros(ndof, ndof)
mm = zeros(ndof, ndof)
gm = zeros(ndof, ndof)
kg = zeros(ndof, ndof)
eld = zeros(ndof)
weights = zeros(struc_el.nip)
g_g = zeros(Int64, ndof, nels)
num = zeros(Int64, fin_el.nod)
actions = zeros(nels, ndof)
displacements = zeros(size(nf, 1), ndim)
gc = ones(ndim, ndim)
dee = zeros(nst,nst)
sigma = zeros(nst)
axial = zeros(nels);

Ok, most arrays have been initialized. Time to start the real work. First determine the global numbering:


In [80]:
?PtFEM.formnf!


Out[80]:

formnf!

Returns nodal freedom numbering array nf

Function

formnf!(nodof::Int64, nn::Int64, nf::Matrix{Int64})

Arguments

* nodof::Int64       : Number of degrees of freedom for each node
* nn::Int64          : Number of nodes in mesh
* nf::Array{Int64,2} : Nodal freedom matrix (updated)

In [81]:
PtFEM.formnf!(nodof, nn, nf)
nf


Out[81]:
1×4 Array{Int64,2}:
 1  2  3  0

Node N (= 4 in this example) is fixed, nodes 1 to 3 represents the degrees of freedom. We need 9 equations to solve for these 3 displacements.


In [82]:
neq = maximum(nf)


Out[82]:
3

In [83]:
kdiag = zeros(Int64, neq);

In [84]:
ell = zeros(nels) # Used to hold element length
if :x_coords in keys(data)
    for i in 1:length(data[:x_coords])-1
        ell[i] = data[:x_coords][i+1] - data[:x_coords][i]
    end
end
ell


Out[84]:
3-element Array{Float64,1}:
 0.333333
 0.333333
 0.333333

In [85]:
for i in 1:nels
    num = [i; i+1]
    PtFEM.num_to_g!(fin_el.nod, nodof, nn, ndof, num, nf, g)
    g_g[:, i] = g
    PtFEM.fkdiag!(ndof, neq, g, kdiag)
end

Matrix g_g now holds the globally numbered degrees of freedom for the Line elements. Note that the 2nd node of element 3 is fixed.


In [86]:
g_g


Out[86]:
2×3 Array{Int64,2}:
 1  2  3
 2  3  0

Compute the locations in the skyline vector kv which form the diagonal elements of the full stiffness matrix.


In [87]:
for i in 2:neq
    kdiag[i] = kdiag[i] + kdiag[i-1]
end
kdiag


Out[87]:
3-element Array{Int64,1}:
 1
 3
 5

Above kdiag holds the indices of the diagonal elements of the stiffness matrix in the kv skyline vector.


In [88]:
kv = zeros(kdiag[neq])
gv = zeros(kdiag[neq])
print("There are $(neq) equations,")
println(" and the skyline storage is $(kdiag[neq]).\n")


There are 3 equations, and the skyline storage is 5.


In [89]:
loads = zeros(neq+1)
if :loaded_nodes in keys(data)
    for i in 1:size(data[:loaded_nodes], 1)
        loads[nf[:, data[:loaded_nodes][i][1]]+1] = 
            data[:loaded_nodes][i][2]
    end
end
nf


Out[89]:
1×4 Array{Int64,2}:
 1  2  3  0

In [90]:
for i in 1:nels
    km = PtFEM.rod_km!(km, prop[etype[i], 1], ell[i])
    g = g_g[:, i]
    PtFEM.fsparv!(kv, km, g, kdiag)
end

In [91]:
km


Out[91]:
2×2 Array{Float64,2}:
  300000.0  -300000.0
 -300000.0   300000.0

In [92]:
kv


Out[92]:
5-element Array{Float64,1}:
  300000.0
 -300000.0
  600000.0
 -300000.0
  600000.0

In [93]:
3EA/l


Out[93]:
300000.0

The function fromSkyline(kv, kdiag) uses the indices in kdiag to reconstruct the (symmetrical) stiffness matrix.


In [94]:
sm = fromSkyline(kv, kdiag)    # Compare eq 1.48 in EEM.
sm


Out[94]:
3×3 Array{Float64,2}:
  300000.0  -300000.0        0.0
 -300000.0   600000.0  -300000.0
       0.0  -300000.0   600000.0

Julia has it's own sparse matrix representation. At some point the intenstion is to replace the skyline format with the Julia SparseArrays representation.


In [95]:
?AbstractSparseArray


search: AbstractSparseArray AbstractSparseMatrix AbstractSparseVector

Out[95]:

No documentation found.

Summary:

abstract AbstractSparseArray{Tv,Ti,N} <: AbstractArray{Tv,N}

Subtypes:

Base.SparseArrays.CHOLMOD.Sparse{Tv<:Union{Complex{Float64},Float64}}
SparseMatrixCSC{Tv,Ti<:Integer}
SparseVector{Tv,Ti<:Integer}

In [96]:
sparse(sm)


Out[96]:
3×3 sparse matrix with 7 Float64 nonzero entries:
	[1, 1]  =  300000.0
	[2, 1]  =  -300000.0
	[1, 2]  =  -300000.0
	[2, 2]  =  600000.0
	[3, 2]  =  -300000.0
	[2, 3]  =  -300000.0
	[3, 3]  =  600000.0

In [97]:
# Add radial stress if fin_el is 3d and axisymmetric
if ndim == 3 && struc_el.axisymmetric
    nst = 4
end

In [98]:
fixed_freedoms = 0
if :fixed_freedoms in keys(data)
    fixed_freedoms = size(data[:fixed_freedoms], 1)
end
no = zeros(Int64, fixed_freedoms)
node = zeros(Int64, fixed_freedoms)
sense = zeros(Int64, fixed_freedoms)
value = zeros(Float64, fixed_freedoms)
if :fixed_freedoms in keys(data) && fixed_freedoms > 0
    for i in 1:fixed_freedoms
      node[i] = data[:fixed_freedoms][i][1]
      sense[i] = data[:fixed_freedoms][i][2]
      no[i] = nf[sense[i], node[i]]
      value[i] = data[:fixed_freedoms][i][3]
    end
    kv[kdiag[no]] = kv[kdiag[no]] + penalty
    loads[no+1] = kv[kdiag[no]] .* value
end
loads


Out[98]:
4-element Array{Float64,1}:
 0.833333
 0.833333
 1.66667 
 1.66667 

In [99]:
PtFEM.sparin!(kv, kdiag)

In [100]:
kv


Out[100]:
5-element Array{Float64,1}:
  547.723
 -547.723
  547.723
 -547.723
  547.723

In [101]:
chol(sm)


Out[101]:
3×3 UpperTriangular{Float64,Array{Float64,2}}:
 547.723  -547.723     0.0  
    ⋅      547.723  -547.723
    ⋅         ⋅      547.723

In [102]:
loads


Out[102]:
4-element Array{Float64,1}:
 0.833333
 0.833333
 1.66667 
 1.66667 

In [103]:
loads[2:end] = PtFEM.spabac!(kv, loads[2:end], kdiag)
loads


Out[103]:
4-element Array{Float64,1}:
 0.833333  
 2.5e-5    
 2.22222e-5
 1.38889e-5

In [104]:
displacements = zeros(size(nf))
for i in 1:size(displacements, 1)
    for j in 1:size(displacements, 2)
      if nf[i, j] > 0
        displacements[i,j] = loads[nf[i, j]+1]
      end
    end
end
displacements = displacements'


Out[104]:
4×1 Array{Float64,2}:
 2.5e-5    
 2.22222e-5
 1.38889e-5
 0.0       

In [105]:
loads[1] = 0.0
for i in 1:nels
    km = PtFEM.rod_km!(km, prop[etype[i], 1], ell[i])
    g = g_g[:, i]
    eld = loads[g+1]
    actions[i, :] = km * eld
end
actions


Out[105]:
3×2 Array{Float64,2}:
 0.833333  -0.833333
 2.5       -2.5     
 4.16667   -4.16667 

In [106]:
m=FEM(struc_el, fin_el, ndim, nels, nst, ndof, nn, nodof, neq, penalty,
    etype, g, g_g, g_num, kdiag, nf, no, node, num, sense, actions, 
    bee, coord, gamma, dee, der, deriv, displacements, eld, fun, gc,
    g_coord, jac, km, mm, gm, kv, gv, loads, points, prop, sigma, value,
    weights, x_coords, y_coords, z_coords, axial);

In [107]:
m.displacements


Out[107]:
4×1 Array{Float64,2}:
 2.5e-5    
 2.22222e-5
 1.38889e-5
 0.0       

In [108]:
m.actions


Out[108]:
3×2 Array{Float64,2}:
 0.833333  -0.833333
 2.5       -2.5     
 4.16667   -4.16667 

In [109]:
dis_df = DataFrame(
    x_translation = m.displacements[:, 1],
)
fm_df = DataFrame(
    normal_force_1 = m.actions[:, 1],
    normal_force_2 = m.actions[:, 2],
    normal_force_1_corrected = m.actions[:, 1],
    normal_force_2_corrected = m.actions[:, 2]
);

Correct element forces and moments for equivalent nodal forces and moments introduced for loading between nodes


In [110]:
if :eq_nodal_forces_and_moments in keys(data)
    eqfm = data[:eq_nodal_forces_and_moments]
    k = data[:struc_el].fin_el.nod * data[:struc_el].fin_el.nodof
    for t in eqfm
      for i in 1:k
        fm_df[t[1], i+2] = round(fm_df[t[1], i] - t[2][i], 5)
      end
    end
end

In [111]:
fm_df


Out[111]:
normal_force_1normal_force_2normal_force_1_correctednormal_force_2_corrected
10.8333333333333321-0.8333333333333321-0.0-1.66667
22.5-2.51.66667-3.33333
34.166666666666664-4.1666666666666643.33333-5.0