# Problems

The role of problems in JuliaFEM is to work as a container for a set of elements. They contain elements and an information how the elements are assembled to the global assembly. The key point is that thanks to multiple dispatch, each problem defines also the physics behind discretization and all problems are assembled using a command `assemble_elements!`

.

As an example, a heat equation in two dimensions is discretized. Mathematically known also as Poisson problem. Strong form of the problem is

and corresponding variational problem is to find $u\in\mathcal{U}$ such that for all $v\in\mathcal{V}$

Let's call $k$ thermal conductivity, $f$ heat source and $g$ heat flux.

Each new problem must be a subtype of `FieldProblem`

or `BoundaryProblem`

. The main difference between these two are that `FieldProblem`

is assembling local matrices for domain $\Omega$ whereas `BoundaryProblem`

is creating (in general) constraint matrices for boundary $\Gamma_{\mathrm{D}}$. The general structure to solve in JuliaFEM is currently described by four different matrices and two vectors, i.e.

We believe that this construction is general enough to describe all possible situations in future. Quite often $\boldsymbol{C}_1 = \boldsymbol{C}_2^\mathrm{T}$ and $\boldsymbol{D} = \boldsymbol{0}$ so that we have a typical saddle point problem

which is equivalent to minimization problem

## Discretizing field problem

So, first we must define a new type, e.g. `Heat`

, which is a subclass of `FieldProblem`

. Problem-wide parameters can be defined into struct if needed.

```
type Heat <: FieldProblem
end
```

In principle it's possible to assemble one element at a time, but way more memory efficient is to assemble all elements with same kind (basis) at same time and preallocate memory only one time before looping through element list. Implementation for assembling local stiffness matrix is

```
function assemble_elements!{B}(problem::Problem{Heat}, assembly::Assembly,
elements::Vector{Element{B}}, time::Float64)
println("Assembling elements of kind $B")
bi = BasisInfo(B)
ndofs = length(B)
Ke = zeros(ndofs, ndofs)
for element in elements
fill!(Ke, 0.0)
for ip in get_integration_points(element)
J, detJ, N, dN = element_info!(bi, element, ip, time)
k = element("thermal conductivity", ip, time)
Ke += ip.weight * k*dN'*dN * detJ
end
gdofs = get_gdofs(problem, element)
add!(assembly.K, gdofs, gdofs, Ke)
end
end
nothing # hide
```

Here, first some memory is allocated to calculate Jacobian, gradients etc. to `BasisInfo`

. `Ke`

is used to store local stiffness matrix. Then iterate over all elements and integration points, evaluate basis and add contribution to local stiffness matrix. Finally, get global degrees of freedom for element by using command `get_gdofs`

and finally add contribution to global stiffness matrix `K`

.

From performance point of view, it's important that memory allocations inside at least the innermost loop is minimized, although assembling global stiffness matrix can be parallelized (at least almost) perfectly and is not considered as a bottleneck when models get big enough. It's anyway a good idea to pay attention to the memory allocations.

### Setting and getting global degrees of freedom for element

`get_gdofs`

is returning the global degrees of freedom for element. They can be set manually using `set_gdofs(problem, element, dofs)`

. Otherwise they are calculated automatically based on the problem dimension using formula `g(i,j,d) = d*(i-1)+j`

, where `i`

is local node number, `j`

is the number of degree of freedom and `d`

is the maximum number of degrees of freedom for node. With this formula dofs are ordered so that first comes all dofs for node 1, then for node 2 and so on. For 3 dofs/node we get $(u_{11}, u_{12}, u_{13}, u_{21}, u_{22}, u_{23}, \ldots, u_{N1}, u_{N2}, u_{N3})$, where the first index is node number and second index is dof number.

Let's create some test problem and test our implementation:

```
el1 = Element(Quad4, [1, 2, 3, 4])
X = Dict(1 => [0.0,0.0], 2 => [1.0,0.0], 3 => [1.0,1.0], 4 => [0.0,1.0])
update!(el1, "geometry", X)
update!(el1, "thermal conductivity", 6.0)
elements = [el1]
assembly = Assembly()
problem = Problem(Heat, "test problem", 1)
nothing # hide
```

Now the struct `Heat`

is empty. If we need to store some problem-wide settings to that struct, they can be found from `problem.properties`

. When creating a new `Problem`

, the syntax is `Problem(type, name, field_dimension)`

, where two first arguments are self descriptive. The third one is the information, how many degrees of freedom is in this problem. As Poisson problem is scalar equation, there is only 1 degrees of freedom per node. For example in continuum mechanics, where the unknown field is displacement, there is usually 2-6 degrees of freedom per node, depending on problem type. Next we do the actual assembling of elements into global stiffness matrix:

```
time = 0.0
assemble_elements!(problem, assembly, elements, time)
```

```
julia> full(assembly.K)
ERROR: UndefVarError: full not defined
```

There is actually one `Assembly`

inside `Problem`

and elements are defined to problem using `add_elements!`

, so a more realistic use case to create global assembly would be to use `assemble!(problem, time)`

as shown below:

```
el2 = Element(Tri3, [3, 2, 5])
X[5] = [2.0, 1.0]
elements = [el1, el2]
update!(elements, "geometry", X)
update!(elements, "thermal conductivity", 6.0)
problem = Problem(Heat, "test problem", 1)
add_elements!(problem, elements)
assemble!(problem, time)
```

```
julia> full(problem.assembly.K)
ERROR: UndefVarError: problem not defined
```

Now, function defined above is actually executed two times, first for elements using `Tri3`

basis and after that for `Quad4`

. That is, `assemble!(problem, time)`

is grouping elements by their type and calling function for each element type separately. It also does some initializations and gives possibility to mangle matrices before and after assembly. These hacks may be useful if one needs to add some discrete values to the matrices after assembly or e.g. save matrices to disk for later diagnoses.

We also need to deal with the integrals on the right hand side. The first integral is done over the domain and can be included to the same `assemble_elements!`

-function than stiffness matrix. Boundary term can be handled in different ways. One option is to define it yet in same function and search for fields like `surface traction S1`

, where `S1`

is one side of the element. This is how it is done in ABAQUS. Another way is to use lower dimensional "boundary element" in assembly and add surface term to that element. This is how it is done in Code Aster.

`assembly_elements!`

-function defined above can be overridden by restricting the type of elements list, `elements::Vector{Element{B}}`

to a some spesific elements. This allows, for example, to optimize assembly for some certain element type what is commonly used. Another use case is to define different assembly function for boundary elements. In 2d setting, voluminal elements like Tet3, Tet6, Quad4, Quad8, Quad9 are integrated over volume and they one dimensional counterparts Seg2, Seg3 can be used to assign boundary fluxes.

The updated `assemble_elements!`

-function, which can also handle volume load from right hand side of the equation, i.e.

looks like following:

```
function assemble_elements!{B}(problem::Problem{Heat}, assembly::Assembly,
elements::Vector{Element{B}}, time::Float64)
println("Assembling volume elements of kind $B")
bi = BasisInfo(B)
ndofs = length(B)
Ke = zeros(ndofs, ndofs)
fe = zeros(ndofs)
for element in elements
fill!(Ke, 0.0)
fill!(fe, 0.0)
for ip in get_integration_points(element)
J, detJ, N, dN = element_info!(bi, element, ip, time)
k = element("thermal conductivity", ip, time)
Ke += ip.weight * k*dN'*dN * detJ
if haskey(element, "heat source")
f = element("heat source", ip, time)
fe += ip.weight * N'*f * detJ
end
end
gdofs = get_gdofs(problem, element)
add!(assembly.K, gdofs, gdofs, Ke)
add!(assembly.f, gdofs, fe)
end
end
nothing # hide
```

At last, implement boundary elements to handle heat flux. To choose what elements should use this assembly function, `elements::Vector{Element{B}}`

must be restricted only to group where `B<:Union{Seg2, Seg}`

.

```
function assemble_elements!{B<:Union{Seg2,Seg3}}(problem::Problem{Heat},
assembly::Assembly,
elements::Vector{Element{B}},
time::Float64)
println("Assembling boundary elements of kind $B")
bi = BasisInfo(B)
ndofs = length(B)
fe = zeros(ndofs)
for element in elements
haskey(element, "heat flux") || continue
fill!(fe, 0.0)
for ip in get_integration_points(element)
J, detJ, N, dN = element_info!(bi, element, ip, time)
g = element("heat flux", ip, time)
fe += ip.weight * N'*g * detJ
end
gdofs = get_gdofs(problem, element)
add!(assembly.f, gdofs, fe)
end
end
nothing # hide
```

To test everything implemented, create some small test problem:

```
el1 = Element(Quad4, [1, 2, 3, 4])
el2 = Element(Tri3, [3, 2, 5])
el3 = Element(Seg2, [3, 5])
X = Dict(1 => [0.0, 0.0],
2 => [1.0, 0.0],
3 => [1.0, 1.0],
4 => [0.0, 1.0],
5 => [2.0, 1.0])
elements = [el1, el2, el3]
update!(elements, "geometry", X)
update!(elements, "thermal conductivity", 6.0)
update!(el3, "heat flux", 264.0)
update!(el1, "heat source", 132.0)
problem = Problem(Heat, "test problem", 1)
add_elements!(problem, elements)
assemble!(problem, time)
```

Global stiffness matrix $\boldsymbol{K}$ and force vector $\boldsymbol{f}$ are

```
julia> K = full(problem.assembly.K)
ERROR: UndefVarError: problem not defined
julia> f = full(problem.assembly.f)
ERROR: UndefVarError: problem not defined
```

To get unique solution, some essential boundary condition must be given, e.g. set dofs 1 and 4 fixed, $u_1 = u_4 = 0$.

```
julia> u = zeros(5)
5-element Array{Float64,1}:
0.0
0.0
0.0
0.0
0.0
julia> all_dofs = collect(1:5)
5-element Array{Int64,1}:
1
2
3
4
5
julia> fixed_dofs = [1, 4]
2-element Array{Int64,1}:
1
4
julia> free_dofs = setdiff(all_dofs, fixed_dofs)
3-element Array{Int64,1}:
2
3
5
julia> u[free_dofs] = K[free_dofs,free_dofs] \ f[free_dofs]
ERROR: UndefVarError: K not defined
julia> u
5-element Array{Float64,1}:
0.0
0.0
0.0
0.0
0.0
```

## Discretizing boundary problem

Can be e.g. Dirichlet boundary, contact / mesh tie problem between two domains, BEM formulation, kinematic coupling (MPC) and so on.

A good question is to determine how to map element local degrees of freedom with global ones. As the plan is to develop a multiphysical FEM platform, it's a hard question how to create this kind of mapping in a dynamic way. Yet another question is how to set boundary conditions for different type of physics. For now, these questions are answered by defining the name of unknown field using function `get_unknown_field_name`

and giving the dimension of unknown field when creating a problem, so that change of information between two problem is possible. This solution has already identified shortcomings and can be expected to change in future.

```
function get_unknown_field_name(::Problem{Heat})
return "temperature"
end
nothing # hide
```

```
type DirichletBC <: BoundaryProblem
end
```

Our implementation to handle boundary condition $u = u_0$ looks following:

```
function assemble_elements!{E}(problem::Problem{DirichletBC},
assembly::Assembly,
elements::Vector{Element{E}},
time::Float64)
name = get_parent_field_name(problem)
dim = get_unknown_field_dimension(problem)
println("Assembling Dirichlet boundary condition")
println("Field name = $name, dofs/node = $dim")
data = Dict{Int64,Float64}()
for element in elements
for i=1:dim
haskey(element, "$name $dim") || continue
gdofs = get_gdofs(problem, element)
ldofs = gdofs[i:dim:end]
xis = get_reference_element_coordinates(E)
for (ldof, xi) in zip(ldofs, xis)
data[ldof] = interpolate(element, "$name $dim", xi, time)
end
end
end
for (dof, val) in data
add!(assembly.C1, dof, dof, 1.0)
add!(assembly.C2, dof, dof, 1.0)
add!(assembly.g, dof, val)
end
end
```

To fix dofs 1 and 4 like before:

```
bel1 = Element(Seg2, [1, 4])
update!(bel1, "geometry", X)
update!(bel1, "temperature 1", 0.0)
bc = Problem(DirichletBC, "fixed", 1, "temperature")
add_elements!(bc, [bel1])
assemble!(bc, time)
```

Now we have all matrices needed:

```
julia> C1 = full(bc.assembly.C1, 5, 5)
ERROR: UndefVarError: bc not defined
julia> C2 = full(bc.assembly.C2, 5, 5)
ERROR: UndefVarError: bc not defined
julia> D = full(bc.assembly.D, 5, 5)
ERROR: UndefVarError: bc not defined
julia> g = full(bc.assembly.g, 5, 1)
ERROR: UndefVarError: bc not defined
```

Together with already calculated matrices, we can now form saddle point problem $\boldsymbol{A}\boldsymbol{x} = \boldsymbol{b}$:

```
julia> A = [K C1; C2 D]
ERROR: UndefVarError: K not defined
julia> b = [f; g]
ERROR: UndefVarError: f not defined
```

Solution:

```
julia> nz = [1, 2, 3, 4, 5, 6, 9]
7-element Array{Int64,1}:
1
2
3
4
5
6
9
julia> x = zeros(10)
10-element Array{Float64,1}:
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
julia> x[nz] = A[nz,nz] \ b[nz]
ERROR: UndefVarError: A not defined
```

As a result we have found $\boldsymbol{u}$ and $\boldsymbol{\lambda}$:

```
julia> u = x[1:5]
5-element Array{Float64,1}:
0.0
0.0
0.0
0.0
0.0
julia> la = x[6:10]
5-element Array{Float64,1}:
0.0
0.0
0.0
0.0
0.0
julia> u' * la
0.0
```

Both field problem and boundary problems can of course have all four matrices and two vectors. For example, in finite sliding contact algorithms all four matrices are used as algorithms are contributing to stiffness matrix also when linearized properly.