CSoM: Julia experiments with the toolkit as described in the book "Programming the Finite Element Method" by I. M. Smith & D. V. Griffiths (5th edition).



In [1]:

    
#Pkg.update()



In [2]:

    
versioninfo()
using CSoM









    



Julia Version 0.5.0

Derived from example 4.1 of "Programming the FEM" (PtFEM).

This is the first example program in the book.

Currently each example in PtFEM consists of 3 Julia component3:

The input data, e.g. p4.1.1.jl
A program skeleton as described in PtFEM, e.g. Program 4.1 is FE4_1.jl
The skeleton programs call the appropriate Julia functions in subdirectory CSoM

As in the book, the input are subdivided by chapter name (e.g. "4 Static Equilibrium") in the "examples" sub-directory.



In [3]:

    
path = joinpath("/Users/rob/.julia/v0.5/CSoM", "src", "Chap04")
include(path*"/FE4_1.jl")









    Out[3]:





FE4_1

Discussion of the initial implementation of Example 4.1.

This version (FE4_1.jl) is an almost straight translation from Fortran to Julia (maybe with the exception of the inputs). This is true for all FEx_x.jl files.

These experiments are aimed at a next version (P4.1.jl) which contains a more Julia-based implementation, e.g. using Julia sparse matrices vs. the skyline format and using Julia (and underlying libraries) for solving the global matrix equations, e.g. y = A \ x. All Px.x.jl will be similarly be Julia based scripts, and Px.x.x.jl contain the corresponding input files.

Made the number of fin_els and distributed load parameters

N: Number of fin_els (nxe) in :x direction (left to right) F: Total of distributed load



In [4]:

    
N = 4;    
F = 5.0;

Distribute the load

Distribute the total fin_el loads over N+1 nodes. Assign half-load to both end points. Last line just ensures the type of dist_loads is concrete.



In [5]:

    
dist_loads = [[(i, [-F/N]) for i in 1:(N+1)];]
dist_loads[1] = (1, [-F/(2*N)])
dist_loads[size(dist_loads,1)] = (N+1, [-F/(2*N)])
dist_loads = convert(Vector{Tuple{Int64, Vector{Float64}}}, dist_loads)









    Out[5]:





5-element Array{Tuple{Int64,Array{Float64,1}},1}:
 (1,[-0.625])
 (2,[-1.25]) 
 (3,[-1.25]) 
 (4,[-1.25]) 
 (5,[-0.625])

Notice that in this example the x axis goes from left to right (not as in the book) and the distributed load is compressive. Clamped end is node 1.



In [6]:

    
?FE4_1









    



search:







    Out[6]:




FE4_1
Backbone method for static equilibrium analysis of a rod.
Constructors
FE4_1(data::Dict)

Arguments
* `data` : Dictionary containing all input data

Dictionary keys
* struc_el::StructuralElement                            : Type of  structural fin_el
* support::Array{Tuple{Int64,Array{Int64,1}},1}        : Fixed-displacements vector
* loaded_nodes::Array{Tuple{Int64,Array{Float64,1}},1} : Node load vector
* properties::Vector{Float64}                          : Material properties
* x_coords::LinSpace{Float64}                          : Xcoordinate vector

Optional dictionary keys
* etype::Vector{Int64}                                 : Element material vector

Examples
using CSoM

data = Dict(
  # Rod(nels, np_types, nip, finite_element(nod, nodof))
  :struc_el => Rod(4, 1, 1, Line(2, 1)),
  :properties => [1.0e5;],
  :x_coords => linspace(0, 1, 5),
  :support => [(1, [0])],
  :loaded_nodes => [(1,[-0.625]),(2,[-1.25]),(3,[-1.25]),(4,[-1.25]),(5,[-0.625])]
)

m = FE4_1(data)

println("Displacements:")
m.displacements |> display
println()

println("Actions:")
m.actions |> display
println()

Related help
?StructuralElement  : Help on structural elements
?Rod                : Help on a Rod structural fin_el
?FiniteElement      : Help on finite element types



In [7]:

    
data = Dict(
# Rod(nxe, np_types, nip, fin_el(nod, nodof))
:struc_el => Rod(N, 1, 1, Line(2, 1)),
  :properties => [1.0e5;],
  :x_coords => linspace(0, 1, (N+1)),
  :support => [(1, [0])],
  :loaded_nodes => dist_loads
);









    



UndefVarError: Rod not defined



In [8]:

    
m = FE4_1(data);









    



UndefVarError: data not defined



In [9]:

    
m.displacements









    



UndefVarError: m not defined



In [10]:

    
plot(x=data[:x_coords], y=m.displacements[:,1],Geom.line())









    



UndefVarError: data not defined



In [11]:

    
m.actions









    



UndefVarError: m not defined



In [12]:

    
plot(x=data[:x_coords], y=[F; m.actions[:,1]],Geom.line())









    



UndefVarError: data not defined

Step by step through FE4_1

Above the results using FE4_1.jl. Below goes through the steps in FE4_1.jl. Initial part just checks the input Dict.



In [13]:

    
if :struc_el in keys(data)
    struc_el = data[:struc_el]
end









    



UndefVarError: data not defined



In [14]:

    
?Rod









    



search: prod prod! produce cumprod cumprod! broadcast broadcast!

Couldn't find Rod
Perhaps you meant mod, prod, rol, ror, Ref, cor, for, Col, Cmd, cd, cld or cond






    Out[14]:




No documentation found.
Binding Rod does not exist.



In [15]:

    
ndim = 1
nst = struc_el.np_types;









    



UndefVarError: struc_el not defined



In [16]:

    
fin_el = struc_el.fin_el









    



UndefVarError: struc_el not defined



In [17]:

    
?mesh_size









    



search:

Couldn't find mesh_size
Perhaps you meant filesize






    Out[17]:




No documentation found.
Binding mesh_size does not exist.



In [18]:

    
if typeof(fin_el) == Line
    (nels, nn) = 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









    



UndefVarError: fin_el not defined



In [19]:

    
nodof = fin_el.nodof         # Degrees of freedom per node









    



UndefVarError: fin_el not defined



In [20]:

    
ndof = fin_el.nod * nodof    # Degrees of freedom per fin_el









    



UndefVarError: fin_el not defined



In [21]:

    
penalty = 1e20
if :penalty in keys(data)
    penalty = data[:penalty]
end









    



UndefVarError: data not defined



In [22]:

    
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









    



UndefVarError: data not defined

 in anonymous at ./<missing>:?



In [23]:

    
prop









    



UndefVarError: prop not defined



In [24]:

    
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)









    



UndefVarError: nodof not defined



In [25]:

    
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









    



UndefVarError: nn not defined



In [26]:

    
collect(x_coords)









    



UndefVarError: x_coords not defined



In [27]:

    
#
# 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);









    



UndefVarError: struc_el not defined

Ok, all arrays have been initialized. Time to start the real work. First determine the global numbering:



In [28]:

    
?formnf!









    



search: UniformScaling

Couldn't find formnf!
Perhaps you meant foreach






    Out[28]:




No documentation found.
Binding formnf! does not exist.



In [29]:

    
formnf!(nodof, nn, nf)
nf









    



UndefVarError: formnf! not defined

Node 1 is fixed, nodes 2 to 5 represents degrees of freedom. We need 4 equations to solve for 4 displacements.



In [30]:

    
neq = maximum(nf)









    



UndefVarError: nf not defined



In [31]:

    
kdiag = zeros(Int64, neq);









    



UndefVarError: neq not defined



In [32]:

    
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









    



UndefVarError: nels not defined



In [33]:

    
for i in 1:nels
    num = [i; i+1]
    num_to_g!(fin_el.nod, nodof, nn, ndof, num, nf, g)
    println(g)
    g_g[:, i] = g
    fkdiag!(ndof, neq, g, kdiag)
    println(kdiag)
end









    



UndefVarError: nels not defined

 in anonymous at ./<missing>:?



In [34]:

    
kdiag









    



UndefVarError: kdiag not defined



In [35]:

    
g_g









    



UndefVarError: g_g not defined



In [36]:

    
for i in 2:neq
    kdiag[i] = kdiag[i] + kdiag[i-1]
end
kdiag









    



UndefVarError: neq not defined

 in anonymous at ./<missing>:?



In [37]:

    
kv = zeros(kdiag[neq])
gv = zeros(kdiag[neq])
print("There are $(neq) equations,")
println(" and the skyline storage is $(kdiag[neq]).\n")









    



UndefVarError: kdiag not defined



In [38]:

    
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









    



UndefVarError: neq not defined



In [39]:

    
for i in 1:nels
    km = rod_km!(km, prop[etype[i], 1], ell[i])
    g = g_g[:, i]
    fsparv!(kv, km, g, kdiag)
end









    



UndefVarError: nels not defined

 in anonymous at ./<missing>:?



In [40]:

    
km









    



UndefVarError: km not defined



In [41]:

    
kv









    



UndefVarError: kv not defined



In [42]:

    
# Add radial stress if fin_el is 3d and axisymmetric
if ndim == 3 && struc_el.axisymmetric
    nst = 4
end



In [43]:

    
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









    



UndefVarError: data not defined



In [44]:

    
sparin!(kv, kdiag)
loads[2:end] = spabac!(kv, loads[2:end], kdiag)
loads









    



UndefVarError: sparin! not defined



In [45]:

    
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'









    



UndefVarError: nf not defined



In [46]:

    
loads[1] = 0.0
for i in 1:nels
    km = rod_km!(km, prop[etype[i], 1], ell[i])
    g = g_g[:, i]
    eld = loads[g+1]
    actions[i, :] = km * eld
end
actions









    



UndefVarError: loads not defined



In [47]:

    
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);









    



UndefVarError: FEM not defined



In [48]:

    
m.displacements









    



UndefVarError: m not defined



In [49]:

    
m.actions









    



UndefVarError: m not defined