using PtFEM

ProjDir = @__DIR__

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

struct_el = :Rod
fin_el = :Line

data = Dict(
  # Rod(nels, np_types, nip, finite_element(nod, nodof))
  :struc_el => Rod(4, 1, 1, Line(2, 1)),
  :properties => [2.0e3; 1.0e3],
  :etype => [2, 2, 1, 1],
  :x_coords => collect(range(0, stop=1, length=5)),
  :support => [
    (1, [0])
    ],
  :fixed_freedoms => [
    (5, 1, 0.05)
    ]
)

data |> display
println()

@time fem, dis_df, fm_df = p41(data)
println()

display(fem.actions)
println()
display(fem.displacements)
println()

display(fm_df)
println()
display(dis_df)
  
using Plots
gr(size=(400,500))

x = 0.0:l/els:l
u = convert(Array, dis_df[:, :x_translation])
  
p = Vector{Plots.Plot{Plots.GRBackend}}(undef, 2)
titles = ["PtFEM Ex42.1 u(x)", "PtFEM Ex42.1 N(x)"]
 
p[1]=plot(ylim=(0.0, 1.0), xlim=(0.0, 70.0),
  yflip=true, xflip=false, xlab="Normal force [N]",
  ylab="x [m]", title=titles[2]
)
vals = fm_df[:, :normal_force_2]
for i in 1:els
    plot!(p[1], 
      [vals[i], vals[i]],
      [(i-1)*l/els, i*l/els], 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

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

plot(p..., layout=(1, 2))
savefig(ProjDir*"/Ex41.2.png")
