# The scalar assembly

The "scalar" assembly is the "original" way for assembling finite element matrices in FEAT2. It is fast, comfortable and designed to assemble "scalar" matrices, i.e., matrices for operators $A_h : V_h to W_h$ between two finite element spaces $V_h$ and $W_h$. "Block" matrices are assembled by assembling all scalar submatrices, one after the other. There exist assembly routines for domain integrals and boundary integrals, so this type of assembly can be used to assemble weak forms as well as, e.g., inhomogeneous Neumann boundary conditions.

Contents | |
---|---|

1 | Involved modules |

2 | Nomenclature |

3 | Derivative quantifiers |

4 | Assembling linear forms |

5 | Assembling trilinear forms |

6 | Advanced assembly |

## Involved modules

The design of the "scalar" assembly focuses on the weak form. There exist special structures that express a linear, bilinear or even trilinear form, the three most common types of operators. There are five main modules involved in this approach:

Module | Content |
---|---|

scalarpde.f90 | Contains structures that describe linear, bilinear and trilinear forms |

linearformevaluation.f90 | Contains subroutines to assembles linear forms. Domain integration and boundary integration |

bilinearformevaluation.f90 | Contains subroutines to assemble bilinear forms |

trilinearformevaluation.f90 | Contains subroutines to assemble trilinear forms |

stdoperators.f90 | Support module. Uses the above modules to realise standard operators, e.g., Mass and Laplace operators. |

derivatives.f90 | Definition of derivatives of basis functions |

collection.f90 | Collection structure, used to pass parameters to callback functions |

feevaluation.f90 | Evaluation of finite element functions in cubature points. |

## Nomenclature

The nomenclature is as follows:

**Linear form:**A "linear form" is the integral on the right-hand side of an equation. Linear forms assemble*right-hand side vectors*. Example:

$$F(v_h) := int_{Omega} f v_h $$

**Bilinear form:**A "bilinear form" describes a set of domain integrals on the left-hand side of an equation. Bilinear forms assemble*matrices*. A bilinear form may have constant coefficients (e.g., Laplace, Mass) or nonconstant coefficients (e.g., anisotropic difusion).- Example for a bilinear form: $$a(v_h,w_h) = int_{Omega} v_h w_h $$
- Example for a bilinear form, with $g:Omegatomathbb{R}$ an additional function: $$b(v_h,w_h) = int_{Omega} g v_h w_h $$

**Trilinear form:**A "trilinear form" is a special type of bilinear form, where the coefficient in front of the bilinear form is expessed as a nonconstant finite element function. Example:

$$n(u_h,v_h,w_h) = int_{Omega} u_h v_h w_h $$

## Derivative quantifiers (derivative.f90)

The module `derivatives.f90`

declares a set of quantifiers to configure derivatives of basis functions. The following table gives an overview about these quantifiers. The operators $partial_{x_1}$, $partial_{x_2}$ and $partial_{x_3}$ refer to a derivative in x, y or z direction, respectively.

Id | Dimension | Type of derivative |
---|---|---|

DER_FUNC | 1D, 2D, 3D | Function value, no derivative |

DER_DERIVX_1D | 1D | Identifies $partial_{x}$ |

DER_DERIVXX_1D | Identifies $partial^2_{xx}$ | |

DER_DERIVX_2D | 2D | Identifies $partial_{x_1}$ |

DER_DERIVY_2D | Identifies $partial_{x_2}$ | |

DER_DERIVXX_2D | Identifies $partial^2_{x_1x_1}$ | |

DER_DERIVYY_2D | Identifies $partial^2_{x_2x_2}$ | |

DER_DERIVXY_2D | Identifies $partial^2_{x_1x_2}$ | |

DER_DERIVX_3D | 3D | Identifies $partial_{x_1}$ |

DER_DERIVY_3D | Identifies $partial_{x_2}$ | |

DER_DERIVZ_3D | Identifies $partial_{x_3}$ | |

DER_DERIVXX_2D | Identifies $partial^2_{x_1x_1}$ | |

DER_DERIVYY_2D | Identifies $partial^2_{x_2x_2}$ | |

DER_DERIVZZ_2D | Identifies $partial^2_{x_3}$ | |

DER_DERIVXY_2D | Identifies $partial^2_{x_1x_2}$ | |

DER_DERIVXZ_2D | Identifies $partial^2_{x_1x_3}$ | |

DER_DERIVYZ_2D | Identifies $partial^2_{x_2x_3}$ |

## Assembling linear forms (right-hand side vectors)

### Domain integration

A "linear form for domain integration" has the following rather general shape:

$$F(v_h) := int_{Omega} f v_h + g_1 partial_{x_1} v_h + g_2 partial_{x_2} v_h + ... $$

The sum below the integral contains combinations of theoretically arbitrary functions $f$ and $g_i$, i=1,2,..., multiplied with a test function $v_h$ and/or its derivatives. There can be up to `SCPDE_NNAB=21`

many elements in this sum. This type of linear form is the most commonly used type of integration used for right-hand side assembly.

For the scalar assembly of this domain integral, this linear form must be mapped by the application to the structure `t_linform`

from `scalarpde.f90`

as follows:

Variable | Content |
---|---|

(t_linearForm)%itermCount | Number of terms in the linear form |

(t_linearForm)%Idescriptors(1) | Derivative quantifier for the test function, 1st term |

(t_linearForm)%Idescriptors(2) | Derivative quantifier for the test function, 2nd term |

(t_linearForm)%Idescriptors(3) | Derivative quantifier for the test function, 3rd term |

... | ... |

The descriptors `Idescriptors(:)`

take one of the values `DER_xxxx`

and identify which derivative to apply to the test function in the linear form.

Example:The following code identifies the linear form $F(v_h) := int_{Omega} f v_h $:

```
type(t_linform) :: rlinform
...
rlinform%itermCount = 1
rlinform%Idescriptors(1) = DER_FUNC
```

Example:The following code identifies the linear form $$F(v_h) := int_{Omega} f partial_{x_1} v_h $$ in 2D:

```
type(t_linform) :: rlinform
...
rlinform%itermCount = 1
rlinform%Idescriptors(1) = DER_DERIVX_2D
```

Example:The following code identifies the linear form $$F(v_h) := int_{Omega} g_1 partial_{x_1} v_h + g_2 partial_{x_2} v_h + g_3 partial_{x_3} v_h $$ in 3D:

```
type(t_linform) :: rlinform
...
rlinform%itermCount = 3
rlinform%Idescriptors(1) = DER_DERIVX_3D
rlinform%Idescriptors(2) = DER_DERIVY_3D
rlinform%Idescriptors(3) = DER_DERIVZ_3D
```

To invoke the actual assembly, the subroutine `linf_buildVectorScalar`

from `linearformevaluation.f90`

has to be called. This subroutine expects on the one hand the `t_linform`

structure configuring the linear form, on the other hand a special callback routine that calculates the values of all the functions $f$, $g_1$, $g_2$,... in the cubature points. The vector where the data should be calculated to has to be created in advance.

Example:The following piece of code demonstrates how to assemble the linear form $$F(v_h) := int_{Omega} f v_h $$ for a function $$f(x)=32x_1x_2$$ in 2D:

```
subroutine frhs (&
..., nelements, npointsPerElement, Dpoints, ..., Dcoefficients...)
...
do iel = 1,nelements
do ipt = 1,npointsPerElement
dx = Dpoints(1,ipt,iel)
dy = Dpoints(2,ipt,iel)
Dcoefficients(1,ipt,iel) = 32.0 * dx * dy
end do
end do
end subroutine
...
type(t_linform) :: rlinform
type(t_vectorScalar) :: rrhs
...
rlinform%itermCount = 1
rlinform%Idescriptors(1) = DER_FUNC
call linf_buildVectorScalar (&
rlinform, .false., rrhs, fcoeff_buildVectorSc_sim = frhs)
...
```

This example demonstrates a special design approach of the complete FEAT2 library. The callback function `frhs`

introduced above has to calculate the values of the function $f$ *simultaneously in all cubature points on all given elements*. This is a difference to many other finite element codes which call such callback functions for each cubature point or for each element separately. The kernel internally divides the set of all elements into smaller chunks of elements and invokes the callback routine for the corresponding element/point-sets. The loops inside of the callback routine are usually perfectly parallisable which is the heart of efficiency in the FEAT2 package.

Example:The following example demonstrates how to assemble the linear form $$F(v_h) := int_{Omega} g_1 partial_{x_1} v_h + g_2 partial_{x_2} v_h + g_3 partial_{x_3} v_h $$ in 3D for functions $$ g_1(x) = x_1^2,quad g_2(x) = x_2^2,quad g_3(x) = x_3^2 $$

```
subroutine frhs (&
..., nelements, npointsPerElement, Dpoints, ..., Dcoefficients...)
...
do iel = 1,nelements
do ipt = 1,npointsPerElement
dx = Dpoints(1,ipt,iel)
dy = Dpoints(2,ipt,iel)
dz = Dpoints(3,ipt,iel)
Dcoefficients(1,ipt,iel) = dx**2
Dcoefficients(2,ipt,iel) = dy**2
Dcoefficients(3,ipt,iel) = dz**2
end do
end do
end subroutine
...
type(t_linform) :: rlinform
type(t_vectorScalar) :: rrhs
...
rlinform%itermCount = 3
rlinform%Idescriptors(1) = DER_DERIVX_3D
rlinform%Idescriptors(2) = DER_DERIVY_3D
rlinform%Idescriptors(3) = DER_DERIVZ_3D
call linf_buildVectorScalar (&
rlinform, .false., rrhs, fcoeff_buildVectorSc_sim = frhs)
...
```

As can be seen, the array `Dcoefficients`

is filled with the values for $g_1$, $g_2$ and $g_3$ in all cubature points on all elements. `Dcoefficients(1,:,:)`

corresponds to `Idescriptors(1)`

(so $g_1$), `Dcoefficients(2,:,:)`

to `Idescriptors(2)`

(so $g_2$) and `Dcoefficients(3,:,:)`

to `Idescriptors(3)`

(so $g_3$).

## Assembling bilinear forms

Bilinear forms are designed for the assembly of linear operators. They can be assembled with constant and nonconstant coefficients. COnstant coefficients are used, e.g., for operators like Laplace and Mass while nonconstant coefficients allow to impose nonlinearities. Note that the situation where a nonconstant coefficient is given by a finite element function is a special case which can more easily be addressed with trilinear forms, see below.

### Domain integration

A "bilinear form for domain integration" has the following rather general shape:

$$a(v_h,w_h) := int_{Omega} g_1 v_h w_h + g_2 partial_{x_1} v_h w_h + g_3 partial_{x_1} v_h w_h + g_4 partial_{x_1} v_h partial_{x_1} w_h + ... $$

The sum below the integral contains combinations of trial functions $v_h$ and test functions $w_h$, multiplied by coefficients $g_i$, i=1,2,... .
If all $g_i$ are constants, the bilinear form is termed "bilinear form with constant coefficients", otherwise "with nonconstant coefficients". There can be up to `SCPDE_NNAB=21`

many elements in this sum.

For the scalar assembly of this domain integral, this bilinear form must be mapped by the application to the structure `t_bilinearForm`

from `scalarpde.f90`

. The structure and setup is similar to a linear form:

Variable | Content |
---|---|

(t_bilinearForm)%itermCount | Number of terms in the linear form |

(t_bilinearForm)%ballCoeffConstant | .true. if the bilinear form has constant coefficients |

(t_bilinearForm)%Dcoefficients(:) | Constant coefficients |

(t_bilinearForm)%Idescriptors(1,1) | Derivative quantifier for the trial function, 1st term |

(t_bilinearForm)%Idescriptors(2,1) | Derivative quantifier for the test function, 1st term |

(t_bilinearForm)%Idescriptors(1,2) | Derivative quantifier for the trial function, 2nd term |

(t_bilinearForm)%Idescriptors(2,2) | Derivative quantifier for the test function, 2nd term |

(t_bilinearForm)%Idescriptors(1,3) | Derivative quantifier for the trial function, 3rd term |

(t_bilinearForm)%Idescriptors(2,3) | Derivative quantifier for the test function, 3rd term |

... | ... |

The descriptors `Idescriptors(:,:)`

take one of the values `DER_xxxx`

and identify which derivative to apply to the trial and test functions in the bilinear form. `Idescriptors(1,:)`

refers to the trial function, `Idescriptors(2,:)`

to the test function.

Example:The following code sets up a structure for the bilinear form $$a(v_h,w_h) := int_{Omega} v_h w_h$$ which is used to assemble mass matrices. This bilinear form has one term with constant coefficient 1.0:

```
type(t_bilinearForm) :: rform
...
rform%itermCount = 1
rform%ballCoeffConstant = .true.
rform%Dcoefficients(1) = 1.0
rform%Idescriptors(1,1) = DER_FUNC
rform%Idescriptors(2,1) = DER_FUNC
```

Example:The following code sets up a structure for the bilinear form $$a(v_h,w_h) := int_{Omega} partial_{x_1}v_h partial_{x_1}w_h + partial_{x_2}v_h partial_{x_2}w_h$$ which is used to assemble a Laplace matrix in 2D. This bilinear form has two term with constant coefficient 1.0:

```
type(t_bilinearForm) :: rform
...
rform%itermCount = 2
rform%ballCoeffConstant = .true.
rform%Dcoefficients(1) = 1.0
rform%Idescriptors(1,1) = DER_DERIVX_2D
rform%Idescriptors(2,1) = DER_DERIVX_2D
rform%Dcoefficients(2) = 1.0
rform%Idescriptors(1,2) = DER_DERIVY_2D
rform%Idescriptors(2,2) = DER_DERIVY_2D
```

Example:The following code sets up a structure for the bilinear form $$a(v_h,w_h) := int_{Omega} 2 partial_{x_1}v_h w_h + 1 partial_{x_2}v_h w_h$$ which corresponds to the operator $begin{pmatrix} 2 \ 1 end{pmatrix} nabla u$ in 2D.

```
type(t_bilinearForm) :: rform
...
rform%itermCount = 2
rform%ballCoeffConstant = .true.
rform%Dcoefficients(1) = 2.0
rform%Idescriptors(1,1) = DER_DERIVX_2D
rform%Idescriptors(2,1) = DER_FUNC
rform%Dcoefficients(2) = 1.0
rform%Idescriptors(1,2) = DER_DERIVY_2D
rform%Idescriptors(2,2) = DER_FUNC
```

**Assembly with constant coefficients:** For the actual assembly, the subroutine `bilf_buildMatrixScalar`

from the module `bilinearformevaluation.f90`

has to be invoked. The routine accepts the configured `t_bilinearForm`

structure and the matrix to be build. **Remark:** For sparse matrices (e.g. in CSR format), the sparsity pattern / matrix structure has to be calculated in advance before calculating the content. Furthermore, memory for the content has to be allocated in advance as well.

Example:The following code calculates a mass matrix in CSR format, based on a discretisation structure`rspatialDiscr`

that describes the underlying finite element space $V_h$:

```
type(t_spatialDiscretisation) :: rspatialDiscr
type(t_matrixScalar) :: rmatrix
type(t_bilinearForm) :: rform
...
! Create the matrix structure, allocate memory for the entries
call bilf_createMatrixStructure (rspatialDiscr,LSYSSC_MATRIX9,rmatrix)
call lsyssc_allocEmptyMatrix (rmatrix,.true.)
...
! Set up the bilinear form
rform%itermCount = 1
rform%ballCoeffConstant = .true.
rform%Dcoefficients(1) = 1.0
rform%Idescriptors(1,1) = DER_FUNC
rform%Idescriptors(2,1) = DER_FUNC
! Assemble the matrix entries
call bilf_buildMatrixScalar (rform, .false., rmatrix)
```

**Assembly with nonconstant coefficients:** For nonlinear equations, it is often necessary to assemble an operator with a nonconstant coefficient $g_i$. Similar to a linear form, nonconstant coefficients are realised by a callback function. The callback function returns for every cubature point on a set of elements the coefficients in front of all terms in the bilinear form. Nonconstant coefficients are activated by

```
rform%ballCoeffConstant = .false.
```

in the setup of the bilinear form. The coefficients `Dcoefficients(:,:,:)`

are not used in this case.

Example:The following example demonstrates how to set up the matrix for the bilinear form of the operator $begin{pmatrix} 2 \ 1 end{pmatrix} nabla u$ in 2D. The coefficients "2" and "1" only hold for all coordinates $x=(x_1,x_2)$ with $x_1 > 0$ and $x_2 > 0$:

```
subroutine fcoeff (...,nelements,npointsPerElement,Dpoints,&
..., Dcoefficients, ...)
...
integer :: ipt,iel
real(DP) :: dx,dy
! Loop over the points and elements, set up the coefficients
do iel = 1,nelements
do ipt = 1,npointsPerElement
dx = Dpoints(1,ipt,iel)
dy = Dpoints(2,ipt,iel)
if ((dx .gt. 0.0_DP) .and. (dy .gt. 0.0_DP)) then
Dcoefficients(1,ipt,iel) = 2.0_DP
Dcoefficients(2,ipt,iel) = 1.0_DP
else
Dcoefficients(1,ipt,iel) = 0.0_DP
Dcoefficients(2,ipt,iel) = 0.0_DP
end if
end do
end do
end subroutine
...
type(t_matrixScalar) :: rmatrix
type(t_bilinearForm) :: rform
...
! Set up a bilinear form. Nonconstant coefficients.
rform%itermCount = 2
rform%ballCoeffConstant = .false.
rform%Idescriptors(1,1) = DER_DERIVX_2D
rform%Idescriptors(2,1) = DER_FUNC
rform%Idescriptors(1,1) = DER_DERIVY_2D
rform%Idescriptors(2,1) = DER_FUNC
! Assemble the matrix
call bilf_buildMatrixScalar (rform, .false., rmatrix, &
fcoeff_buildMatrixSc_sim = fcoeff)
```

As can be seen from the example, the additional coefficient function `fcoeff`

returns the coefficients for all the terms in the bilinear form. It has to be specified in `bilf_buildMatrixScalar`

in accordance to `ballCoeffConstant = .false.`

, which activates nonconstant coefficients. The coefficient array `Dcoefficients`

is not used in this case.

## Assembling trilinear forms

Trilinear forms are rather similar to bilinear forms. A trilinear form can be seen as a bilinear form with a nonconstant coefficient, given by a finite element function. FEAT2 automatically evaluates the the additional finite element function in all cubature points and multiplies the values with the trial and test functions. Additional constant or nonconstant coefficients can be specified as well. The assembly routines for trilinear forms can be found in the module `trilinearform.f90`

.

### Domain integration

A "trililinear form for domain integration" has the following shape:

$$b(u_h,v_h,w_h) := int_{Omega} g_1 u_h v_h w_h + g_2 partial_{x_1} u_h v_h w_h + g_3 u_h partial_{x_1} v_h w_h + g_4 partial_{x_1} u_h partial_{x_1} v_h w_h + ... $$

The sum below the integral contains combinations of the finite element function $u_h$, trial function $v_h$ and test function $w_h$, multiplied by coefficients $g_i$, i=1,2,... . Similar to the bilinear form, if all $g_i$ are constants, the trilinear form is termed "trilinear form with constant coefficients", otherwise "with nonconstant coefficients". There can be up to `SCPDE_NNAB=21`

many elements in this sum.

For the scalar assembly of this domain integral, this trilinear form must be mapped by the application to the structure `t_trilinearForm`

from `scalarpde.f90`

. The structure and setup is similar to a linear form:

Variable | Content |
---|---|

(t_trilinearForm)%itermCount | Number of terms in the linear form |

(t_trilinearForm)%ballCoeffConstant | .true. if the bilinear form has constant coefficients |

(t_trilinearForm)%Dcoefficients | Constant coefficients |

(t_trilinearForm)%Idescriptors(1,1) | Derivative quantifier for the FE function $u_h$, 1st term |

(t_trilinearForm)%Idescriptors(2,1) | Derivative quantifier for the trial function, 1st term |

(t_trilinearForm)%Idescriptors(3,1) | Derivative quantifier for the test function, 1st term |

(t_trilinearForm)%Idescriptors(1,2) | Derivative quantifier for the FE function $u_h$, 2nd term |

(t_trilinearForm)%Idescriptors(2,2) | Derivative quantifier for the trial function, 2nd term |

(t_trilinearForm)%Idescriptors(3,2) | Derivative quantifier for the test function, 2nd term |

(t_trilinearForm)%Idescriptors(1,3) | Derivative quantifier for the FE function $u_h$, 3rd term |

(t_trilinearForm)%Idescriptors(2,3) | Derivative quantifier for the trial function, 3rd term |

(t_trilinearForm)%Idescriptors(3,3) | Derivative quantifier for the test function, 3rd term |

... | ... |

The descriptors `Idescriptors(:,:)`

take one of the values `DER_xxxx`

and identify which derivative to apply to the function $u_h$, the trial functions and the test functions in the trilinear form. `Idescriptors(1,:)`

refers to the finite element function $u_h$, `Idescriptors(2,:)`

refers to the trial function, `Idescriptors(3,:)`

to the test function.

Example:The following code assembles the trilinear form $$b(u_h,v_h,w_h) := int_{Omega} u_h v_h w_h$$ which is used to assemble mass matrices. This trilinear form has one term with constant coefficient 1.0. The finite element function $u_h$ is specified by the structure`rvectorUh`

:

```
type(t_matrixScalar) :: rmatrix
type(t_vectorScalar) :: rvectorUh
type(t_trilinearForm) :: rform
...
rform%itermCount = 1
rform%ballCoeffConstant = .true.
rform%Dcoefficients(1) = 1.0
rform%Idescriptors(1,1) = DER_FUNC
rform%Idescriptors(2,1) = DER_FUNC
rform%Idescriptors(3,1) = DER_FUNC
call trilf_buildMatrixScalar (rform,.false.,rmatrix,rvectorUh)
```

Example:The following code assembles the trilinear form $$b(u_h,v_h,w_h) := int_{Omega} u_h partial_{x_1} v_h w_h.$$ This trilinear form has one term with constant coefficient 1.0. The finite element function $u_h$ is specified by the structure`rvectorUh`

which is passed as additional parameter in the call to`trilf_buildMatrixScalar`

:

```
type(t_matrixScalar) :: rmatrix
type(t_vectorScalar) :: rvectorUh
type(t_trilinearForm) :: rform
...
rform%itermCount = 1
rform%ballCoeffConstant = .true.
rform%Dcoefficients(1) = 1.0
rform%Idescriptors(1,1) = DER_FUNC
rform%Idescriptors(2,1) = DER_DERIVX_2D
rform%Idescriptors(3,1) = DER_FUNC
call trilf_buildMatrixScalar (rform,.false.,rmatrix,rvectorUh)
```

**Nonconstant coefficients:** Nonconstant coefficients are specified by setting

```
rform%ballCoeffConstant = .false.
```

and specifying an appropriate callback function in the call to `trilf_buildMatrixScalar`

, similar to the bilinear form. The coefficient array `Dcoefficients(:)`

is not used in this case.

## Advanced assembly

### Passing parameters

For the assembly of the right-hand side as well as for the assembly of bilinear/triliear forms with nonconstant coefficients, callback functions are used. Usually, these callback functions are parametrised, i.e., parameters have to be passed to these functions. For this purpose, the usual "collection strategy" of the FEAT2 package is used. The assembly routine accepts an optional parameter `rcollection`

of type `t_collection`

which is passed to the callback function:

Type of data | passed via |
---|---|

Double precision data can be passed via | `(t_collection)%DquickAccess(:)` |

Integers can be passed via | `(t_collection)%IquickAccess(:)` |

Strings can be passed via | `(t_collection)%SquickAccess(:)` |

Block vectors can be passed via | `(t_collection)%p_rvectorQuickAccessX` with N=1,2,... |

Parsers can be passed via | `(t_collection)%p_rfparserQuickAccessN` with N=1,2,... |

Parameter lists can be passed via | `(t_collection)%p_rparlistQuickAccessN` with N=1,2,... |

Complex parameters can be specified via | named items in the collection. |

etc. Using only the quick access arrays, the collection does not even have to be initialised.

Example:The following example shows how to pass variables`dradius`

and`itype`

to the callback routine:

```
subroutine frhs (..., Dcoefficients, rcollection)
...
real(DP) :: dradius
integer :: itype
! Get parameters from rcollection
dradius = rcollection%DquickAccess(1)
itype = rcollection%IquickAccess(1)
...
end subroutine
...
real(DP) :: dradius
integer :: itype
type(t_collection) :: rcoll
...
! Save parameters to rcoll, pass to linf_buildVectorScalar.
rcoll%DquickAccess(1) = dradius
rcoll%IquickAccess(1) = itype
call linf_buildVectorScalar (..., frhs, rcoll)
...
```

The parameters `dradius`

and `itype`

are written to `rcoll%DquickAccess(1)`

and `rcoll%IquickAccess(1)`

. `rcoll`

is passed to `linf_buildVectorScalar`

which passes it to `frhs`

. In `frhs`

`dradius`

and `itype`

can be recovered from the quick access arrays.

The same technique can be used for bilinear forms and trilinear forms. In all cases, the assembly routine passes the specified `t_collection`

structure through to the callback routine, which can then extract parameters saved to this structure.

### Complex coefficients with multiple FE functions

Simple nonlinearities involving only one finite element function as nonconstant coefficient can be assembled using a trilinear form. However, for equations with complex nonlinearities, it is often necessary to pass multiple finite element functions to callback routines and calculate a complex nonlinear coefficient from this information. For this purpose, the "scalar assembly" can be advised to allocate some temporary memory for values in cubature points, while callback routines can use functions from `feevaluation.f90`

to evaluate values of finite element functions in the cubature points. In the following section, this technique is addressed for the assembly of a nonlinear coefficient in a bilinear form. The technique can, however, be used in exactly the same way for right-hand side vectors and trilinear forms.

To calculate such complex coefficients, the following rules apply:

- Finite element vectors should be passed via the collection to the callback routine. This can be done either by using the "Quick access" arrays or named variables in a collection.
- The user has to specify
`ntempArrays>0`

in the call to the assembly routine, which defines the number of temporary arrays to be allocated. - Inside of the callback routine,
`fevl_evaluate_sim`

from`feevaluation.f90`

can be used to evaluate the finite element function(s) in the cubature points.

The following example demonstrates this approach. We want to assemble a modified mass matrix with nonconstant coefficients, given by the bilinear form
$$a_{u_h}(v_h,w_h) = int_Omega 3(u_1^2 + u_2^2) v_h w_h.$$ This bilinear form has only one term in the bilinear form which, however, depends on a 2D vector field $u_h=(u_1,u_2)$. The following code realises the vector field by a block vector `rvectorUh`

with two components. The vector is passed to the callback routine via a collection and evaluated there. In the call to the assembly method, we specify `ntemparrays=2`

as we need two temporary arrays for $u_1$ and $u_2$.

```
subroutine fcoeff (...,nelements,npointsPerElement,&
..., rdomainIntSubset, Dcoefficients, rcollection)
...
integer :: ipt,iel
real(DP) :: du1, du2
type(t_vectorBlock), pointer :: p_rvectorUh
real(dp), dimension(:,:), pointer :: p_Du1,p_Du2
! Get the vector and pointers to the two temp memory arrays
rvectorUh => rcollection%p_rvectorQuickAccess1
p_Du1 => rdomainIntSubset%p_DtempArrays(:,:,1)
p_Du2 => rdomainIntSubset%p_DtempArrays(:,:,2)
! Calculate u_1 and u_2
call fevl_evaluate_sim (p_rvectorUh%RvectorBlock(1), &
rdomainIntSubset, DER_FUNC, p_Du1)
call fevl_evaluate_sim (p_rvectorUh%RvectorBlock(2), &
rdomainIntSubset, DER_FUNC, p_Du2)
! Loop over the points and elements, set up the coefficients
do iel = 1,nelements
do ipt = 1,npointsPerElement
! Get the calculated values, calculate the coefficient
du1 = p_Du1(ipt,iel)
du2 = p_Du2(ipt,iel)
Dcoefficients(1,ipt,iel) = 3.0_DP*(du1**2 + du2**2)
end do
end do
end subroutine
...
type(t_vectorBlock), target :: rvectorUh
type(t_matrixScalar) :: rmatrix
type(t_bilinearForm) :: rform
type(t_collection) :: rcoll
...
! Set up the bilinear form, nonconstant coefficients
rform%itermCount = 1
rform%ballCoeffConstant = .false.
rform%Dcoefficients(1) = 1.0
rform%Idescriptors(1,1) = DER_FUNC
rform%Idescriptors(2,1) = DER_FUNC
! Pass the FE function via the collection
rcoll%p_rvectorQuickAccess1 => rvectorUh
! Assemble the matrix entries
call bilf_buildMatrixScalar (rform, .false., rmatrix, &
fcoeff_buildMatrixSc_sim = fcoeff,&
rcollection = rcoll,&
ntempArrays = 2)
...
```

The routine `bilf_buildMatrixScalar`

provides `ntempArrays=2`

temporary arrays in `rdomainIntSubset%p_DtempArrays`

for the callback routine. This array provides memory in all cubature points on all elements of the current element set. The callback routine can fill this memory with arbitrary information to calculate the nonlinear coefficient(s).