By Johan Blaauwendraad http://www.delftacademicpress.nl/f040.php (in Dutch).
This Tutorial_02 notebook uses the figure 1.1 example in EEM and explains in detail how the FE4_1 skeleton program in the CSoM package works.
Tutorial_02 continues where Tutorial_01 ends. Repeat the initial steps because later on in the explanations we need m.kv, the stiffnes matrix in Skyline format.
In [1]:
using CSoM, DataFrames, Plots
In [2]:
path = joinpath(Pkg.dir("CSoM"), "src", "Chap04")
# If FE4_1 already loaded, skip the include.
if !isdefined(Main, :FE4_1)
include(path*"/FE4_1.jl")
end
In [3]:
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 [4]:
l = 1.0 # Total length of structural element [m]
q = 5.0 # Distributed load [N/m]
N = 10 # 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 ;
In [5]:
data = Dict(
# Rod(nxe, np_types, nip, fin_el(nod, nodof))
:struc_el => getfield(Main, Symbol(struct_el))(els, 1, 1,
getfield(Main, Symbol(fin_el))(nod, nodof)),
:properties => [1.0e5;],
: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 [6]:
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 [7]:
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 [8]:
m = FE4_1(data);
In [9]:
if :struc_el in keys(data)
struc_el = data[:struc_el]
end
Out[9]:
In [10]:
ndim = 1
nst = struc_el.np_types;
In [11]:
fin_el = struc_el.fin_el
Out[11]:
In [12]:
if typeof(fin_el) == Line # 1D finite element
(nels, nn) = CSoM.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[12]:
In [13]:
nodof = fin_el.nodof # Degrees of freedom per node
Out[13]:
In [14]:
ndof = fin_el.nod * nodof # Degrees of freedom per fin_el, in this case each finite element has 2 nodes
Out[14]:
In [15]:
penalty = 1e20 # used to fix fixed_freedoms
if :penalty in keys(data)
penalty = data[:penalty]
end
In [16]:
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 [17]:
prop
Out[17]:
In [18]:
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[18]:
In [19]:
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[19]:
In [20]:
collect(x_coords)
Out[20]:
In [21]:
#
# 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 [22]:
?CSoM.formnf!
Out[22]:
In [23]:
CSoM.formnf!(nodof, nn, nf)
nf
Out[23]:
Node N (= 10 in this example) is fixed, nodes 1 to 9 represents the degrees of freedom. We need 9 equations to solve for these 9 displacements.
In [24]:
neq = maximum(nf)
Out[24]:
In [25]:
kdiag = zeros(Int64, neq);
In [26]:
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[26]:
In [27]:
for i in 1:nels
num = [i; i+1]
CSoM.num_to_g!(fin_el.nod, nodof, nn, ndof, num, nf, g)
g_g[:, i] = g
CSoM.fkdiag!(ndof, neq, g, kdiag)
end
In [28]:
g_g
Out[28]:
In [29]:
for i in 2:neq
kdiag[i] = kdiag[i] + kdiag[i-1]
end
kdiag
Out[29]:
Above kdiag holds the indices of the diagonal elements of the stiffness matrix in the kv skyline vector.
To illustrate this, we'll borrow kv from the earlier call to FE4_1:
In [30]:
m.kv
Out[30]:
The function fromSkyline(m.kv, m.kdiag) uses the indices in kdiag to reconstruct the (symmetrical) stiffness matrix.
In [31]:
sm = fromSkyline(m.kv, m.kdiag)
Out[31]:
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 [32]:
?AbstractSparseArray
Out[32]:
In [33]:
sparse(sm)
Out[33]:
In [34]:
kv = zeros(kdiag[neq])
gv = zeros(kdiag[neq])
print("There are $(neq) equations,")
println(" and the skyline storage is $(kdiag[neq]).\n")
In [35]:
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[35]:
In [36]:
for i in 1:nels
km = CSoM.rod_km!(km, prop[etype[i], 1], ell[i])
g = g_g[:, i]
CSoM.fsparv!(kv, km, g, kdiag)
end
In [37]:
km
Out[37]:
In [38]:
kv
Out[38]:
In [39]:
# Add radial stress if fin_el is 3d and axisymmetric
if ndim == 3 && struc_el.axisymmetric
nst = 4
end
In [40]:
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[40]:
In [41]:
CSoM.sparin!(kv, kdiag)
In [42]:
kv
Out[42]:
In [43]:
loads[2:end] = CSoM.spabac!(kv, loads[2:end], kdiag)
loads
Out[43]:
In [44]:
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[44]:
In [45]:
loads[1] = 0.0
for i in 1:nels
km = CSoM.rod_km!(km, prop[etype[i], 1], ell[i])
g = g_g[:, i]
eld = loads[g+1]
actions[i, :] = km * eld
end
actions
Out[45]:
In [46]:
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 [47]:
m.displacements
Out[47]:
In [48]:
m.actions
Out[48]: