In [1]:

    

using PtFEM


    

In [2]:

    

l = 1.0       # Total length [m]
N = 5         # 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
np_types = 1  # Number of proerty types
EA = 1.0e5    # Strain stiffness
nip = 1       # Number of integration points


    

    
Out[2]:





1

In [3]:

    

struct_el = :Rod
fin_el = :Line;


    

In [4]:

    

data = Dict(
  # Rod(nxe, np_types, nip, fin_el(nod, nodof))
  :struc_el => getfield(Main, Symbol(struct_el))(els, np_types, nip,
    getfield(Main, Symbol(fin_el))(nod, nodof)),
  :properties => [EA;],
  :x_coords => 0.0:l/els:l,
  :support => [(N, [0])],
  :loaded_nodes => [
      (1, [-0.625]),
      (2, [-1.25]),
      (3, [-1.25]),
      (4, [-1.25]),
      (5, [-0.625])
    ]
)


    

    
Out[4]:





Dict{Symbol,Any} with 5 entries:
  :support      => Tuple{Int64,Array{Int64,1}}[(5,[0])]
  :struc_el     => PtFEM.Rod(4,1,1,PtFEM.Line(2,1))
  :properties   => [100000.0]
  :x_coords     => 0.0:0.25:1.0
  :loaded_nodes => Tuple{Int64,Array{Float64,1}}[(1,[-0.625]),(2,[-1.25]),(3,[-…

In [5]:

    

@time m = FE4_1(data);


    

    




There are 4 equations and the skyline storage is 7.

 

In [6]:

    

using DataTables


    

In [7]:

    

dis_dt = DataTable(
    x_translation = m.displacements[:, 1],
)


    

    
Out[7]:




x_translation
1
-2.5e-5
2
-2.34375e-5
3
-1.875e-5
4
-1.09375e-5
5
0.0

In [8]:

    

fm_dt = DataTable(
    normal_force_1 = m.actions[:, 1],
    normal_force_2 = m.actions[:, 2]
)


    

    
Out[8]:




normal_force_1
normal_force_2
1
-0.625
0.625
2
-1.875
1.875
3
-3.125
3.125
4
-4.375
4.375

In [9]:

    

using Plots
gr(size=(400,500))


    

    
Out[9]:





Plots.GRBackend()

In [10]:

    

x = 0.0:l/els:l
u = convert(Array, dis_dt[:x_translation])
Nf = vcat(
    convert(Array, fm_dt[:normal_force_1])[1],
    convert(Array, fm_dt[:normal_force_2])
)


    

    
Out[10]:





5-element Array{Float64,1}:
 -0.625
  0.625
  1.875
  3.125
  4.375

In [11]:

    

p = Vector{Plots.Plot{Plots.GRBackend}}(2)
titles = ["p4.1.1 u(x)", "p4.1.1 N(x)"];


    

In [12]:

    

p[1]=plot(ylim=(0.0, 1.0), xlim=(0.0, 5.0),
    yflip=true, xflip=true, xlab="Normal force [N]",
    ylab="x [m]", title=titles[2]
);


    

In [13]:

    

vals = convert(Array, fm_dt[:normal_force_2]);


    

In [14]:

    

for i in 1:els
    plot!(p[1], [vals[i], vals[i]], [(i-1)*l/els, i*l/els], color=:blue,
        color=:blue, fill=true, fillalpha=0.1, leg=false)
    
    delta = abs(((i-1)*l/els) - (i*l/els)) / 20.0
    y1 = collect(((i-1)*l/els):delta:(i*l/els))
    for j in 1:length(y1)
        plot!(p[1], [vals[i], 0.0], [y1[j], y1[j]], color=:blue, alpha=0.5)
  end
end


    

In [15]:

    

p[2] = plot(u, x, xlim=(-0.00003, 0.0), yflip=true,
    xlab="Displacement [m]", ylab="x [m]",
    fillto=0.0, fillalpha=0.1, leg=false, title=titles[1]);


    

In [16]:

    

plot(p..., layout=(1, 2))


    

    
Out[16]:

In [ ]: