PtFEM: 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 [49]:

    
#Pkg.update()



In [50]:

    
versioninfo()
using PtFEM









    



Julia Version 0.5.1
Commit 6445c82 (2017-03-05 13:25 UTC)
Platform Info:
  OS: macOS (x86_64-apple-darwin13.4.0)
  CPU: Intel(R) Core(TM) i5-6267U CPU @ 2.90GHz
  WORD_SIZE: 64
  BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Haswell)
  LAPACK: libopenblas64_
  LIBM: libopenlibm
  LLVM: libLLVM-3.7.1 (ORCJIT, broadwell)

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. Ex4.1.1.jl
A program skeleton as described in PtFEM, e.g. Program 4.1 is p4_1.jl
The skeleton programs call the appropriate Julia functions in subdirectory PtFEM

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

Discussion of the initial implementation of Example 4.1.

This version (p4_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 [51]:

    
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 [52]:

    
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[52]:





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 [53]:

    
?p4_1









    



search: p4_1 p4_7 p4_6 p4_5 p4_4 p4_3 p4_2







    Out[53]:




p4_1
Method for static equilibrium analysis of a rod.
Constructors
p4_1(data::Dict)
p4_1(m::jFEM, data::Dict) # Used to re-use factiored global stiffness matrix

Arguments
* `m`    : Previously created jFEM model
* `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 PtFEM

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 = p4_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 [54]:

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



In [55]:

    
m = p4_1(data);









    



There are 4 equations.



In [56]:

    
m.displacements









    Out[56]:





5×1 Array{Float64,2}:
  0.0       
 -1.09375e-5
 -1.875e-5  
 -2.34375e-5
 -2.5e-5



In [57]:

    
m.actions









    Out[57]:





4×2 Array{Float64,2}:
 4.375  -4.375
 3.125  -3.125
 1.875  -1.875
 0.625  -0.625

Step by step through Ex4_1

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



In [58]:

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









    Out[58]:





PtFEM.Rod(4,1,1,PtFEM.Line(2,1))



In [59]:

    
?Rod









    



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







    Out[59]:




Rod
Concrete 1D structural element with only axial stresses.
Constructor
Rod(nels, np_types, nip, fin_el)

Arguments
* nels::Int64             : Number of fin_els (stored in field nxe)
* np_types::Int64         : Number of different property types
* nip::Int64              : Number of integration points
* fin_el::FiniteElement   : Line(nod, nodof)

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



In [60]:

    
ndim = 1
nst = struc_el.np_types;



In [61]:

    
fin_el = struc_el.fin_el









    Out[61]:





PtFEM.Line(2,1)



In [62]:

    
?PtFEM.mesh_size









    Out[62]:




mesh_size
Function mesh_size returns the number of fin_els (nels) and the number of nodes (nn) in a 1, 2 or 3-d geometry-created mesh.
Method
(nels, nn) = mesh_size(fin_el, nxe, [nye[, nze]])

Arguments
* fe::Element           : Shape of finite element
                          1D: Line
                          2D: Trangle or Quadrilateral
                          3D: Hexahedron
* nxe::Int64            : Number of fin_els in x direction
* nye::Int64            : Number of fin_els in y direction (for 2D and 3D)
* nze::Int64            : Number of fin_els in z direction (3D only)



In [63]:

    
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) = PtFEM.mesh_size(fin_el, struc_el.nxe, struc_el.nye)
elseif typeof(fin_el) == Hexahedron
    (nels, nn) = PtFEM.mesh_size(fin_el, struc_el.nxe, struc_el.nye, struc_el.nze)
end









    Out[63]:





(4,5)



In [64]:

    
nodof = fin_el.nodof         # Degrees of freedom per node









    Out[64]:





1



In [65]:

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









    Out[65]:





2



In [66]:

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



In [67]:

    
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 [68]:

    
prop









    Out[68]:





1×1 Array{Float64,2}:
 100000.0



In [69]:

    
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[69]:





((1,5),
[0 1 … 1 1])



In [70]:

    
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[70]:





5-element LinSpace{Float64}:
 0.0,0.25,0.5,0.75,1.0



In [71]:

    
collect(x_coords)









    Out[71]:





5-element Array{Float64,1}:
 0.0 
 0.25
 0.5 
 0.75
 1.0



In [72]:

    
#
# 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 [73]:

    
?PtFEM.formnf!









    Out[73]:




formnf!
Returns nodal freedom numbering array nf
Function
formnf!(nodof::Int64, nn::Int64, nf::Matrix{Int64})

Arguments
* nodof::Int64       : Number of degrees of freedom for each node
* nn::Int64          : Number of nodes in mesh
* nf::Array{Int64,2} : Nodal freedom matrix (updated)



In [74]:

    
PtFEM.formnf!(nodof, nn, nf)
nf









    Out[74]:





1×5 Array{Int64,2}:
 0  1  2  3  4

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



In [75]:

    
neq = maximum(nf)









    Out[75]:





4



In [76]:

    
kdiag = zeros(Int64, neq);



In [77]:

    
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[77]:





4-element Array{Float64,1}:
 0.25
 0.25
 0.25
 0.25



In [78]:

    
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: num_to_g! not defined

 in macro expansion; at ./In[78]:3 [inlined]
 in anonymous at ./<missing>:?



In [79]:

    
kdiag









    Out[79]:





4-element Array{Int64,1}:
 0
 0
 0
 0



In [80]:

    
g_g









    Out[80]:





2×4 Array{Int64,2}:
 0  0  0  0
 0  0  0  0



In [81]:

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









    Out[81]:





4-element Array{Int64,1}:
 0
 0
 0
 0



In [82]:

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









    



There are 4 equations, and the skyline storage is 0.



In [83]:

    
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[83]:





1×5 Array{Int64,2}:
 0  1  2  3  4



In [84]:

    
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 [85]:

    
km









    Out[85]:





2×2 Array{Float64,2}:
  400000.0  -400000.0
 -400000.0   400000.0



In [86]:

    
kv









    Out[86]:





0-element Array{Float64,1}



In [87]:

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



In [88]:

    
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[88]:





5-element Array{Float64,1}:
 -0.625
 -1.25 
 -1.25 
 -1.25 
 -0.625



In [89]:

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









    



UndefVarError: sparin! not defined



In [90]:

    
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[90]:





5×1 Array{Float64,2}:
  0.0  
 -1.25 
 -1.25 
 -1.25 
 -0.625

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



In [96]:

    
eld









    Out[96]:





2-element Array{Float64,1}:
 0.0
 0.0



In [92]:

    
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 [93]:

    
m.displacements









    Out[93]:





5×1 Array{Float64,2}:
  0.0  
 -1.25 
 -1.25 
 -1.25 
 -0.625



In [94]:

    
m.actions









    Out[94]:





4×2 Array{Float64,2}:
 0.0  0.0
 0.0  0.0
 0.0  0.0
 0.0  0.0



In [ ]: