# Compute arbitrary terms - low-level generic assembly procedures (deprecated)¶

This section present the first version of generic assembly procedure which has been implemented in *GetFEM* and is now considered as deprecated. It allows to make the assembly of arbitrary matrices in the linear case. In the nonlinear case, some special “non_linear_term” object have to be implemented, which could be a bit tricky and obliges to use very low-level internal tools of *GetFEM*. The generic weak form language (GWFL) has been developed to circumvent these difficulties (see Compute arbitrary terms - high-level generic assembly procedures - Generic Weak-Form Language (GWFL)).

As it can be seen in the file `getfem/getfem_assembling.h`

, all the
previous assembly procedures use a `getfem::generic_assembly`

object and provide it an adequate
description of what must be done. For example, the assembly of a volumic source
term for a scalar FEM is done with the following excerpt of code:

```
getfem::generic_assembly assem;
assem.push_im(mim);
assem.push_mf(mf);
assem.push_mf(mfdata);
assem.push_data(F);
assem.push_vec(B);
assem.set("Z=data(#2);"
"V(#1)+=comp(Base(#1).Base(#2))(:,j).Z(j);");
assem.assembly();
```

The first instructions declare the object, and set the data that it will use: a
mesh_im object which holds the integration methods, two mesh_fem objects, the input data
`F`

, and the destination vector `B`

.

The input data is the vector \(F\), defined on `mfdata`

. One wants to
evaluate \(\sum_{j} f_j (\int_\Omega \phi^i \psi^j)\). The instruction must be
seen as something that will be executed for each convex `cv`

of the mesh. The
terms `#1`

and `#2`

refer to the first mesh_fem and the second one (i.e. `mf`

and `mfdata`

). The instruction `Z=data(#2);`

means that for each convex, the
“tensor” `Z`

will receive the values of the first data argument provided with
`push_data`

, at indexes corresponding to the degrees of freedom attached to the
convex of the second (`#2`

) mesh_fem (here, ```
Z =
F[mfdata.ind_dof_of_element(cv)]
```

).

The part `V(#1)+=...`

means that the result of the next expression will be
accumulated into the output vector (provided with `push_vec`

). Here again,
`#1`

means that we will write the result at indexes corresponding to the degrees
of freedom of the current convex with respect to the first (`#1`

) mesh_fem.

The right hand side `comp(Base(#1).Base(#2))(:,j).Z(j)`

contains two operations.
The first one is a computation of a tensor on the convex:
`comp(Base(#1).Base(#2))`

is evaluated as a 2-dimensions tensor,
\(\int\phi^i \psi^j\), for all degrees of freedom \(i\) of `mf`

and
\(j\) of `mfdata`

attached to the current convex. The next part is a
reduction operation, `C(:,j).Z(j)`

: each named index (here \(j\)) is summed,
i.e. the result is \(\sum_j c_{i,j} z_j\).

The integration method used inside `comp(Base(#1).Base(#2))`

is taken from
`mim`

. If you need to use integration methods from another mesh_im object, you can
specify it as the first argument of `comp`

, for example ```
comp(\%2,
Base(#1).Grad(#2))
```

will use the second mesh_im object (New in getfem++-2.0).

An other example is the assembly of the stiffness matrix for a vector Laplacian:

```
getfem::generic_assembly assem;
assem.push_im(mim);
assem.push_mf(mf);
assem.push_mf(mfdata);
assem.push_data(A);
assem.push_mat(SM);
assem.set("a=data$1(#2);"
"M$1(#1,#1)+=sym(comp(vGrad(#1).vGrad(#1).Base(#2))(:,j,k,:,j,k,p).a(p))");
assem.assembly();
```

Now the output is written in a sparse matrix, inserted with
`assem.push_mat(SM)`

. The `$1`

in `M$1(#1,#1)`

just indicates that we refer
to the first matrix “pushed” (it is optional, but if the assembly builds two
matrices, the second one must be referred this way). The `sym`

function ensure
that the result is symmetric (if this is not done, some round-off errors may
cancel the symmetricity, and the assembly will be a little bit slower). Next, the
`comp`

part evaluates a 7D tensor,

where \(\varphi^i_j\) is a \(jth\) component of the \(ith\) base
function of `mf`

and \(\psi^p\) is a (scalar) base function of the second
mesh_fem. Since we want to assemble

the reduction is:

In the `comp`

function, `vGrad`

was used instead of `Grad`

since we said
that we were assembling a *vector* Laplacian: that is why each `vGrad`

part has
three dimensions (dof number, component number, and derivative number). For a
scalar Laplacian, we could have used
`comp(Grad(#1).Grad(#1).Base(#2))(:,k,:,k,p).a(p)`

. But the vector form has the
advantage to work in both vector and scalar case.

The last instruction, `assem.assembly()`

, does evaluate the expression on each
convex. For an assembly over a boundary just call `assem.assembly(rg)`

, where
`rg`

is a `getfem::mesh_region`

object. `rg`

might also be a number, in that case the mesh
region taken into account is `mim.linked_mesh().region(rg)`

.

The third example shows how to compute the \(L^2\) norm of a scalar or vector field on a mesh boundary:

```
assem.push_im(mim);
assem.push_mf(mf);
assem.push_data(U);
std::vector<scalar_type> v(1);
assem.push_vec(v);
assem.set("u=data(#1);"
"V()+=u(i).u(j).comp(vBase(#1).vBase(#1))(i,k,j,k)");
assem.assembly(boundary_number);
```

This one is easy to read. When `assembly`

returns, `v[0]`

will contain

The fourth and last example shows an (sub-optimal) assembly of the linear elasticity problem with a complete Hooke tensor:

```
assem.set("h=data$1(qdim(#1),qdim(#1),qdim(#1),qdim(#1),#2);"
"t=comp(vGrad(#1).vGrad(#1).Base(#2));"
"e=(t{:,2,3,:,5,6,:}+t{:,3,2,:,5,6,:}+t{:,2,3,:,6,5,:}+t{:,3,2,:,6,5,:})/4;"
"M(#1,#1)+= sym(e(:,j,k,:,m,n,p).h(j,k,m,n,p))");
```

The original equations are:

where \(h\) is the Hooke tensor, and \(:\) means the scalar product
between matrices. Since we assume it is not constant, \(h\) is given on the
second mesh_fem: \(h_{ijkl}(x)=\sum_p h_{ijkl}^p \psi^p\). Hence the first line
declares that the first data “pushed” is indeed a five-dimensions tensor, the
first fourth ones being all equal to the target dimension of the first mesh_fem, and
the last one being equal to the number of degrees of freedom of the second mesh_fem.
The `comp`

part still computes the same 7D tensor than for the vector Laplacian
case. From this tensor, one evaluates
\(\varepsilon(\varphi^i)_{jk}\varepsilon(\phi^l)_{mn}\psi^p\) via
permutations, and finally the expression is reduced against the hook tensor.

## available operations inside the `comp`

command¶

`Base(#i)`

: evaluate the value of the base functions of the*ith*mesh_fem`Grad(#i)`

: evaluate the value of the gradient of the base functions of the*ith*mesh_fem`Hess(#i)`

: evaluate the value of the Hessian of the base functions of the*ith*mesh_fem`Normal()`

: evaluate the unit normal (should not be used for volumic integrations !)`NonLin$x(#mf1,... #mfn)`

: evaluate the*xth*non-linear term (inserted with`push_nonlinear_term(pnonlinear_elem_term)`

) using the listed mesh_fem objects.`GradGT()`

,`GradGTInv()`

: evaluate the gradient (and its inverse) of the geometric transformation of the current convex.

Note

you may reference any data object inside the `comp`

command, and perform
reductions inside the `comp()`

. This feature is mostly interesting for
speeding up assembly of nonlinear terms (see the file
`getfem/getfem_nonlinear_elasticity.h`

for an example of use).

## others operations¶

Slices may be mixed with reduction operations `t(:,4,i,i)`

takes a slice at
index 4 of the second dimension, and reduces the diagonal of dimension 3 and 4.
*Please note that index numbers for slices start at 1 and not 0 !!*

`mdim(#2)`

is evaluated as the mesh dimension associated to the second mesh_fem,
while `qdim(#2)`

is the target dimension of the mesh_fem.

The diagonal of a tensor can be obtained with `t{:,:,3,3}`

(which is strictly
equivalent to `t{1,2,3,3}`

: the colon is just here to improve the readability).
This is the same operator than for permutation operations. Note that
`t{:,:,1,1}`

or `t{:,:,4,4}`

are not valid operations.

The `print`

command can be used to see the tensor: `"print comp(Base(#1));"`

will print the integrals of the base functions for each convex.

If there is more than one data array, output array or output sparse
matrix, one can use `data$2`

, `data$3`

, `V$2`

, `M$2`

,…