This is the first part of series of posts about programming new features to JuliaFEM. Simple examples demonstrating main concepts is given. As an example problem, we aim to program new truss element to JuliaFEM.

In principle, all one needs to do is to define a new type of problem, define what kind of field is expected to return and then implement assemble!-function which takes the assembly, problem, and a list of elements to assemble. All development can be done using Jupyter Notebooks if wanted.

In [1]:

using JuliaFEM
import JuliaFEM: assemble!, get_unknown_field_name
Pkg.status("JuliaFEM")

WARNING: Method definition info(Any...) in module Base at util.jl:532 overwritten in module Logging at C:\Users\jahx06\.julia\v0.6\Logging\src\Logging.jl:115.
WARNING: Method definition warn(Any...) in module Base at util.jl:585 overwritten in module Logging at C:\Users\jahx06\.julia\v0.6\Logging\src\Logging.jl:115.
01-Sep 11:36:33:WARNING:root:replacing docs for 'JuliaFEM.field :: Tuple{Float64}' in module 'JuliaFEM'.
01-Sep 11:36:33:WARNING:root:replacing docs for 'JuliaFEM.field :: Tuple{Float64}' in module 'JuliaFEM'.
01-Sep 11:36:37:WARNING:root:replacing docs for 'JuliaFEM.solver :: Tuple{}' in module 'JuliaFEM'.

 - JuliaFEM                      0.3.2+             master

The first thing needs to be done is to define new type which is a subtype of FieldProblem or BoundaryProblem. Boundary problems are affecting to the boundaries of domain while field problems is the actual field equation. A new type can have some internal properties defined if needed, as long as good defaults are provided. Our problem, in this case is simply:

In [2]:

type Truss <: FieldProblem
end

It's mandatory to define what kind of result can be expected from assembly procedure. This is more likely to be changed in future, but for now the return type is defined by a function get_unknown_field_name. In this particular case, we want to solve equations of elasticity and the unknown field is thus displacement:

In [3]:

function get_unknown_field_name(problem::Problem{Truss})
    return "displacement"
end

Out[3]:

get_unknown_field_name (generic function with 6 methods)

The last thing is to provide a function which takes care of assembling element for problem global matrices defined inside problem as problem.assembly. To keep things simple, a very basic assembly procedure is implemented and it is improved in later posts. For now it's enough to know that local stiffness matrix for truss element is

\begin{align} K_{\mathrm{e}} & =\frac{EA}{L}\left[\begin{array}{rr} 1 & -1\\ -1 & 1 \end{array}\right] & f_{e} & =\frac{qL}{2}\left[\begin{array}{c} 1\\ 1 \end{array}\right] \end{align}

In [4]:

function assemble!(assembly::Assembly, problem::Problem{Truss},
                   element::Element{Seg2}, time)
    X = element("geometry", time)  # get geometry
    L = norm(X[2] - X[1])  # calculate length of rod
    E = 1.0
    A = 1.0
    q = 1.0
    Ke = E*A/L*[1.0 -1.0; -1.0 1.0]
    fe = q*L/2*[1.0, 1.0]
    gdofs = get_gdofs(problem, element)  # find global dofs of element
    add!(assembly.K, gdofs, gdofs, Ke)
    add!(assembly.f, gdofs, fe)
end

Out[4]:

assemble! (generic function with 36 methods)

Now, to test our implementation we create a simple 1d truss problem with two elements, assemble that and examine the global stiffness matrix and force vector:

In [5]:

X = Dict(1 => [0.0], 2 => [1.0], 3 => [2.0])
element1 = Element(Seg2, [1, 2])
element2 = Element(Seg2, [2, 3])
elements = [element1, element2]
update!(elements, "geometry", X)
problem = Problem(Truss, "test problem", 1)
add_elements!(problem, elements)
assemble!(problem)

01-Sep 11:37:08:WARNING:root:Assemble problem test problem: seems that problem is uninitialized.
01-Sep 11:37:08:INFO:root:Initializing problem test problem at time 0.0 automatically.
01-Sep 11:37:08:WARNING:root:assemble!() this is default assemble operation, decreased performance can be expected without preallocation of memory!

Out[5]:

true

We can then find the results from problem.assembly, i.e.

In [6]:

full(problem.assembly.K)

Out[6]:

3×3 Array{Float64,2}:
  1.0  -1.0   0.0
 -1.0   2.0  -1.0
  0.0  -1.0   1.0

In [7]:

full(problem.assembly.f)

Out[7]:

3×1 Array{Float64,2}:
 0.5
 1.0
 0.5