In [1]:
#Pkg.update()
PtFEM: The toolkit as described in the book "Programming the Finite Element Method" by I. M. Smith & D. V. Griffiths (5th edition).
In [2]:
using PtFEM, DataFrames, Plots
As in the book, the input are subdivided by chapter name (e.g. "4 Static Equilibrium") in the "examples" sub-directory. In this inital example components 1 and 2 are combined.
In [3]:
path = joinpath("/Users/rob/.julia/v0.5/PtFEM", "src", "Chap04")
include(path*"/FE4_1.jl")
Out[3]:
In [4]:
struct_el = :Rod;
fin_el = :Line;
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;
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_translation = 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
(nels, nn) = PtFEM.mesh_size(fin_el, struc_el.nxe)
elseif typeof(fin_el) == Triangle || typeof(fin_el) == Quadrilateral
(nels, nn) = mesh_size(fin_el, struc_el.nxe, struc_el.nye)
elseif typeof(fin_el) == Hexahedron
(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
Out[29]:
In [30]:
penalty = 1e20
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 all dynamic arrays storen 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, all arrays have been initialized. Time to start the real work. First determine the global numbering:
In [37]:
?PtFEM.formnf!
Out[37]:
In [38]:
PtFEM.formnf!(nodof, nn, nf)
nf
Out[38]:
Node N (=10) is fixed, nodes 1 to 9 represents 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)
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]
PtFEM.num_to_g!(fin_el.nod, nodof, nn, ndof, num, nf, g)
println(g)
g_g[:, i] = g
PtFEM.fkdiag!(ndof, neq, g, kdiag)
println(kdiag)
end
In [43]:
kdiag
Out[43]:
In [44]:
g_g
Out[44]:
In [45]:
for i in 2:neq
kdiag[i] = kdiag[i] + kdiag[i-1]
end
kdiag
Out[45]:
In [46]:
kv = zeros(kdiag[neq])
gv = zeros(kdiag[neq])
print("There are $(neq) equations,")
println(" and the skyline storage is $(kdiag[neq]).\n")
In [47]:
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[47]:
In [48]:
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 [49]:
km
Out[49]:
In [50]:
kv
Out[50]:
In [51]:
# Add radial stress if fin_el is 3d and axisymmetric
if ndim == 3 && struc_el.axisymmetric
nst = 4
end
In [52]:
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[52]:
In [53]:
PtFEM.sparin!(kv, kdiag)
In [54]:
loads[2:end] = PtFEM.spabac!(kv, loads[2:end], kdiag)
loads
Out[54]:
In [55]:
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[55]:
In [56]:
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[56]:
In [57]:
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 [58]:
m.displacements
Out[58]:
In [59]:
m.actions
Out[59]:
In [ ]: