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);
In [67]:
if :struc_el in keys(data)
struc_el = data[:struc_el]
end
Out[67]:
In [68]:
ndim = 1
nst = struc_el.np_types;
In [69]:
fin_el = struc_el.fin_el
Out[69]:
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]:
In [71]:
nodof = fin_el.nodof # Degrees of freedom per node
Out[71]:
In [72]:
ndof = fin_el.nod * nodof # Degrees of freedom per fin_el, in this case each finite element has 2 nodes
Out[72]:
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]:
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]:
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]:
In [78]:
collect(x_coords)
Out[78]:
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]:
In [81]:
PtFEM.formnf!(nodof, nn, nf)
nf
Out[81]:
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]:
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]:
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
In [86]:
g_g
Out[86]:
In [87]:
for i in 2:neq
kdiag[i] = kdiag[i] + kdiag[i-1]
end
kdiag
Out[87]:
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")
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]:
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]:
In [92]:
kv
Out[92]:
In [93]:
3EA/l
Out[93]:
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]:
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
Out[95]:
In [96]:
sparse(sm)
Out[96]:
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]:
In [99]:
PtFEM.sparin!(kv, kdiag)
In [100]:
kv
Out[100]:
In [101]:
chol(sm)
Out[101]:
In [102]:
loads
Out[102]:
In [103]:
loads[2:end] = PtFEM.spabac!(kv, loads[2:end], kdiag)
loads
Out[103]:
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]:
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]:
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]:
In [108]:
m.actions
Out[108]:
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]
);
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]: