---

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