# Some Standard assembly procedures (low-level generic assembly)¶

Procedures defined in the file `getfem/getfem_assembling.h`

allow the
assembly of stiffness matrices, mass matrices and boundary conditions for a few
amount of classical partial differential equation problems. All the procedures
have vectors and matrices template parameters in order to be used with any matrix library.

CAUTION: The assembly procedures do not clean the matrix/vector at the begining of the assembly in order to keep the possibility to perform several assembly operations on the same matrix/vector. Consequently, one has to clean the matrix/vector before the first assembly operation.

## Laplacian (Poisson) problem¶

An assembling procedure is defined to solve the problem:

where \(\Omega\) is an open domain of arbitrary dimension, \(\Gamma_{D}\) and \(\Gamma_{N}\) are parts of the boundary of \(\Omega\), \(u(x)\) is the unknown, \(a(x)\) is a given coefficient, \(f(x)\) is a given source term, \(U(x)\) the prescribed value of \(u(x)\) on \(\Gamma_{D}\) and \(F(x)\) is the prescribed normal derivative of \(u(x)\) on \(\Gamma_{N}\). The function to be called to assemble the stiffness matrix is:

```
getfem::asm_stiffness_matrix_for_laplacian(SM, mim, mfu, mfd, A);
```

where

`SM`

is a matrix of any type having the right dimension (i.e.`mfu.nb_dof()`

),`mim`

is a variable of type`getfem::mesh_im`

defining the integration method used,`mfu`

is a variable of type`getfem::mesh_fem`

and should define the finite element method for the solution,`mfd`

is a variable of type`getfem::mesh_fem`

(possibly equal to`mfu`

) describing the finite element method on which the coefficient \(a(x)\) is defined,`A`

is the (real or complex) vector of the values of this coefficient on each degree of freedom of`mfd`

.

Both mesh_fem should use the same mesh (i.e. ```
&mfu.linked_mesh() ==
&mfd.linked_mesh()
```

).

It is important to pay attention to the fact that the integration methods stored
in `mim`

, used to compute the elementary matrices, have to be chosen of
sufficient order. The order has to be determined considering the polynomial
degrees of element in `mfu`

, in `mfd`

and the geometric transformations for
non-linear cases. For example, with linear geometric transformations, if `mfu`

is a \(P_{K}\) FEM, and `mfd`

is a \(P_{L}\) FEM, the integration will
have to be chosen of order \(\geq 2(K-1) + L\), since the elementary integrals
computed during the assembly of `SM`

are
\(\int\nabla\varphi_i\nabla\varphi_j\psi_k\) (with \(\varphi_i\) the basis
functions for `mfu`

and \(\psi_i\) the basis functions for `mfd`

).

To assemble the source term, the function to be called is:

```
getfem::asm_source_term(B, mim, mfu, mfd, V);
```

where `B`

is a vector of any type having the correct dimension (still
`mfu.nb_dof()`

), `mim`

is a variable of type `getfem::mesh_im`

defining the integration
method used, `mfd`

is a variable of type `getfem::mesh_fem`

(possibly equal to `mfu`

)
describing the finite element method on which \(f(x)\) is defined, and `V`

is the vector of the values of \(f(x)\) on each degree of freedom of `mfd`

.

The function `asm_source_term`

also has an optional argument, which is a
reference to a `getfem::mesh_region`

(or just an integer `i`

, in which case
`mim.linked_mesh().region(i)`

will be considered). Hence for the Neumann
condition on \(\Gamma_{N}\), the same function:

```
getfem::asm_source_term(B, mim, mfu, mfd, V, nbound);
```

is used again, with `nbound`

is the index of the boundary \(\Gamma_{N}\) in
the linked mesh of `mim`

, `mfu`

and `mfd`

.

There is two manner (well not really, since it is also possible to use Lagrange multipliers, or to use penalization) to take into account the Dirichlet condition on \(\Gamma_{D}\), changing the linear system or explicitly reduce to the kernel of the Dirichlet condition. For the first manner, the following function is defined:

```
getfem::assembling_Dirichlet_condition(SM, B, mfu, nbound, R);
```

where `nbound`

is the index of the boundary \(\Gamma_D\) where the Dirichlet
condition is applied, `R`

is the vector of the values of \(R(x)\) on each
degree of freedom of `mfu`

. This operation should be the last one because it
transforms the stiffness matrix `SM`

. It works only for Lagrange elements. At
the end, one obtains the discrete system:

where \(U\) is the discrete unknown.

For the second manner, one should use the more general:

```
getfem::asm_dirichlet_constraints(H, R, mim, mf_u, mf_mult,
mf_r, r, nbound).
```

See the Dirichlet condition as a general linear constraint that must satisfy the
solution \(u\). This function does the assembly of Dirichlet conditions of
type \(\int_{\Gamma} u(x)v(x) = \int_{\Gamma}r(x)v(x)\) for all \(v\) in
the space of multiplier defined by `mf_mult`

. The fem `mf_mult`

could be often
chosen equal to `mf_u`

except when `mf_u`

is too “complex”.

This function just assemble these constraints into a new linear system \(H u=R\), doing some additional simplification in order to obtain a “simple” constraints matrix.

Then, one should call:

```
ncols = getfem::Dirichlet_nullspace(H, N, R, Ud);
```

which will return a vector \(U_d\) which satisfies the Dirichlet condition, and an orthogonal basis \(N\) of the kernel of \(H\). Hence, the discrete system that must be solved is:

and the solution is $U=N U_{int}+U_d$. The output matrix \(N\) should be a
\(nbdof \times nbdof\) (sparse) matrix but should be resized to `ncols`

columns. The output vector \(U_d\) should be a \(nbdof\) vector. A big
advantage of this approach is to be generic, and do not prescribed for the finite
element method `mf_u`

to be of Lagrange type. If `mf_u`

and `mf_d`

are
different, there is implicitly a projection (with respect to the \(L^2\) norm)
of the data on the finite element `mf_u`

.

If you want to treat the more general scalar elliptic equation \(\mbox{div}(A(x)\nabla u)\), where \(A(x)\) is square matrix, you should use:

```
getfem::asm_stiffness_matrix_for_scalar_elliptic(M, mim, mfu,
mfdata, A);
```

The matrix data `A`

should be defined on `mfdata`

. It is expected as a vector
representing a \(n \times n \times nbdof\) tensor (in Fortran order), where
\(n\) is the mesh dimension of `mfu`

, and \(nbdof\) is the number of dof
of `mfdata`

.

## Linear Elasticity problem¶

The following function assembles the stiffness matrix for linear elasticity:

```
getfem::asm_stiffness_matrix_for_linear_elasticity(SM, mim, mfu,
mfd, LAMBDA, MU);
```

where `SM`

is a matrix of any type having the right dimension (i.e. here
`mfu.nb_dof()`

), `mim`

is a variable of type `getfem::mesh_im`

defining the integration
method used, `mfu`

is a variable of type `getfem::mesh_fem`

and should define the finite
element method for the solution, `mfd`

is a variable of type `getfem::mesh_fem`

(possibly
equal to `mfu`

) describing the finite element method on which the Lamé
coefficient are defined, `LAMBDA`

and `MU`

are vectors of the values of Lamé
coefficients on each degree of freedom of `mfd`

.

Caution

Linear elasticity problem is a vectorial problem, so the target dimension of
`mfu`

(see `mf.set_qdim(Q)`

) should be the same as the dimension of the
mesh.

In order to assemble source term, Neumann and Dirichlet conditions, same functions as in previous section can be used.

## Stokes Problem with mixed finite element method¶

The assembly of the mixed term \(B = - \int p\nabla.v\) is done with:

```
getfem::asm_stokes_B(MATRIX &B, const mesh_im &mim,
const mesh_fem &mf_u, const mesh_fem &mf_p);
```

## Assembling a mass matrix¶

Assembly of a mass matrix between two finite elements:

```
getfem::asm_mass_matrix(M, mim, mf1, mf2);
```

It is also possible to obtain mass matrix on a boundary with the same function:

getfem::asm_mass_matrix(M, mim, mf1, mf2, nbound);

where `nbound`

is the region index in `mim.linked_mesh()`

, or a
`mesh_region`

object.