By Johan Blaauwendraad http://www.delftacademicpress.nl/f040.php (in Dutch).
This tutorial notebook uses the figure 1.1 example in EEM and explanations how this can implemented using the CSoM package.
In [1]:
#Pkg.update() # Only needed if Julia is not up to date.
In [2]:
using CSoM, DataFrames, Plots
The Julia DataFrames package is used to display the final results. The Julai Plots package is used to plot the results, in this notebook using the GR backend. See https://juliaplots.github.io
The input data, e.g. for the example defined in p4.1.1.jl. In this tutorial notebook the contents of p4.1.1 is re-created. The p4.1.1.jl Julia file itself can be found in the "examples/4 Static Equilibrium" subdirectory of the CSoM.jl package.
A program skeleton as described in the PtFEM book, e.g. Program 4.1 is FE4_1.jl in subdirectory "src/Chap04" and can be included accordingly (see below).
The skeleton programs call the appropriate Julia functions in subdirectory "src/CSoM".
In [3]:
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
Note: FE4_1.jl, the skeleton program used in this example, is an almost straight translation from Fortran in PtFEM to Julia in CSoM.jl (with the exception of the program input and output sections). This is true for all FEx_x.jl files.
In [4]:
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 [5]:
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 [6]:
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 [7]:
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 [8]:
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 [9]:
data
Out[9]:
Notice that in this 'flagpole' example the x axis goes from top to bottom and the distributed load is compressive. The clamped ('encastré') end node is node N.
In [10]:
m = FE4_1(data);
In [11]:
@time m = FE4_1(data);
In [12]:
m.displacements
Out[12]:
In [13]:
m.actions
Out[13]:
In [14]:
fromSkyline(m.kv, m.kdiag)
Out[14]:
In [15]:
dis_df = DataFrame(
x_deflections = m.displacements[:, 1],
)
fm_df = DataFrame(
normal_force_1 = m.actions[:, 1],
normal_force_2 = m.actions[:, 2],
);
In [16]:
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] -= round(t[2][i], 2)
end
end
end
In [17]:
dis_df
Out[17]:
In [18]:
fm_df
Out[18]:
In [19]:
gr(size=(400,600))
Out[19]:
In [20]:
p = Vector{Plots.Plot{Plots.GRBackend}}(2)
titles = ["EEM fig 1.1 u(x)", "EEM fig 1.1 N(x)"]
fors = vcat(fm_df[:, :normal_force_1], q*l);
In [21]:
p[1] = plot(data[:x_coords], -fors,
xlim=(0,l), ylim=(-6.0, 0.0),
xlabel="x [m]", ylabel="Normal force [N]", color=:blue,
line=(:dash,1), marker=(:dot,1,0.8,stroke(1,:black)),
title=titles[2], leg=false
)
Out[21]:
In [22]:
p[2] = plot(data[:x_coords], m.displacements[:, 1],
xlim=(0, l), ylim=(0.0, 0.00003),
xlabel="x [m]", ylabel="Deflection [m]", color=:red,
line=(:dash,1), marker=(:circle,1,0.8,stroke(1,:black)),
title=titles[1], leg=false
)
Out[22]:
In [23]:
plot(p..., layout=(2, 1))
Out[23]:
In [24]:
if :struc_el in keys(data)
struc_el = data[:struc_el]
end
Out[24]:
In [25]:
ndim = 1
nst = struc_el.np_types;
In [26]:
fin_el = struc_el.fin_el
Out[26]:
In [27]:
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[27]:
In [28]:
nodof = fin_el.nodof # Degrees of freedom per node
Out[28]:
In [29]:
ndof = fin_el.nod * nodof # Degrees of freedom per fin_el, in this case each finite element has 2 nodes
Out[29]:
In [30]:
penalty = 1e20 # used to fix fixed_freedoms
if :penalty in keys(data)
penalty = data[:penalty]
end
In [31]:
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 [32]:
prop
Out[32]:
In [33]:
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[33]:
In [34]:
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[34]:
In [35]:
collect(x_coords)
Out[35]:
In [36]:
#
# 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 [37]:
?CSoM.formnf!
Out[37]:
In [38]:
CSoM.formnf!(nodof, nn, nf)
nf
Out[38]:
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 [39]:
neq = maximum(nf)
Out[39]:
In [40]:
kdiag = zeros(Int64, neq);
In [41]:
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[41]:
In [42]:
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 [43]:
g_g
Out[43]:
In [44]:
for i in 2:neq
kdiag[i] = kdiag[i] + kdiag[i-1]
end
kdiag
Out[44]:
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 [45]:
m.kv
Out[45]:
The function fromSkyline(m.kv, m.kdiag) uses the indices in kdiag to reconstruct the (symmetrical) stiffness matrix.
In [46]:
sm = fromSkyline(m.kv, m.kdiag)
Out[46]:
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 [67]:
?AbstractSparseArray
Out[67]:
In [47]:
sparse(sm)
Out[47]:
In [48]:
kv = zeros(kdiag[neq])
gv = zeros(kdiag[neq])
print("There are $(neq) equations,")
println(" and the skyline storage is $(kdiag[neq]).\n")
In [49]:
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[49]:
In [50]:
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 [51]:
km
Out[51]:
In [52]:
kv
Out[52]:
In [53]:
# Add radial stress if fin_el is 3d and axisymmetric
if ndim == 3 && struc_el.axisymmetric
nst = 4
end
In [54]:
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[54]:
In [55]:
CSoM.sparin!(kv, kdiag)
In [68]:
kv
Out[68]:
In [56]:
loads[2:end] = CSoM.spabac!(kv, loads[2:end], kdiag)
loads
Out[56]:
In [57]:
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[57]:
In [58]:
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[58]:
In [59]:
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 [60]:
m.displacements
Out[60]:
In [61]:
m.actions
Out[61]:
In [ ]: