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")
```

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

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

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

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

Out[5]:

We can then find the results from `problem.assembly`

, i.e.

In [6]:

```
full(problem.assembly.K)
```

Out[6]:

In [7]:

```
full(problem.assembly.f)
```

Out[7]: