In [2]:
using CSoM, DataFrames, Plots
The Julia DataFrames package is used to display the final results. The Julia 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]: