The block assembly
The block matrix assembly in FEAT2 is an alternative assembly strategy for block matrices. Instead of setting up each submatrix in a block matrix separately, this technique directly works on the complete block matrix, thus saving a lot of computational time for the precomputation of all the basis functions. However, for applying this technique, one has to know a bit more about the assembly of matrices and vectors.
In contrast to the "scalar" assembly, the "block" assembly does not have separate structures for a linear form, bilinear form or trilinear form. Instead, appropriate callback routines directly calculate local matrix/vector entries which are imposed by the block assembly routines into the global matrices/vectors. This involves some additional loops in the callback routines, and the callback routines directly have to deal with finite element basis functions.
Involved modules and functionality
The following modules are involved in the block assembly:
Module  Description 

blockmatassemblybase.f90  Definition of structures 
blockmatassembly.f90  The actual assembly routines 
blockmatassemblystdop.f90  Implementation of standard operators with block assembly routines 
feevaluation2.f90  Evaluation of finite element functions for nonlinearities 
Currently, the block assembly routines in blockmatassembly.f90
realise the following functionalities:
Assembly of block matrices: Subroutine
bma_buildMatrix
This subroutine realises the assembly of block matrices. A special callback routine has to be provided that calculates local matrix entries.
Assembly of block vectors: Subroutine
bma_buildVector
This subroutine realises the assembly of block vectors. A special callback routine has to be provided that calculates local vector entries.
Assembly of integrals: Subroutines
bma_buildIntegral
andbma_buildIntegrals
These subroutines can be used to calculate one or multiple integrals or other distributed operators. They are usually used to compute $L_2$ or $H^1$ norms but can also be used to compute the $sup$ norm of a vector or vector field.
The assembly routines support scalarvalues as well as vectorvalues finite element functions and can also deal with interleaved matrices and vectors. However, each functionality needs a special implementation in the corresponding callback routines, see below.
Basics
When a block assembly routine is invoked, the following steps are done:
The routine loops over sets of cells. On every cell, there exist a set of cubature points.
On every cell and in each cubature point on every cell, the assembly routine calculates the values of the basis functions of all finite element spaces that are involved in the assembly. This includes the basis functions that come from the submatrices in the block matrix as well as from the basis functions of all involved nonlinearities.
If there are finite element nonlinearities specified (e.g., a convection given by a finite element function), their values in all the cubature points are calculated.
Next, a callback routine is invoked. This callback routine has to compute the local matrices/vectors/integrals in all the cubature points on all the elements given.
Finally, the assembly routines collect the calculated local matrix/vector/integral contributions and form a global matrix/vector/integral from them.
Simple block matrix assembly
In the simplest form, calling the matrix assembly reads as follows:
type(t_matrixBlock) :: rmatrix
...
call bma_buildMatrix (rmatrix,BMA_CALC_STANDARD,fcalcLocalMatrices)
The routine fcalcLocalMatrices
is a userdefined callback routine that calculates the local matrix contributions.
Simple block vector assembly
In the simplest form, calling the vector assembly reads as follows:
type(t_vectorBlock) :: rrhs
...
call bma_buildVector (rrhs,BMA_CALC_STANDARD,fcalcLocalVectors)
The routine fcalcLocalVectors
is a userdefined callback routine that calculates the local vector contributions.
Simple integral assembly: General integral
In the simplest form, calling the integral assembly reads as follows:
type(t_scalarCubatureInfo), target :: rmyCubature
real(DP) :: dintvalue
...
call bma_buildIntegral (dintvalue,BMA_CALC_STANDARD,fcalcLocalIntegral,&
rtriangulation,rcubatureInfo=rmyCubature)
The routine fcalcLocalIntegral
is a callback routine that calculates the local integral contribution on a set of cells. rtriangulation
specifies the underlying mesh and rmyCubature
the cubature formula to be used. The result of the integration is saved to dintvalue
.
The alternative routine bma_buildIntegrals
cab be used to calculate multiple integrals at once. For example,
type(t_scalarCubatureInfo), target :: rmyCubature
real(DP), dimension(3) :: Dintvalues
...
call bma_buildIntegrals (Dintvalue,BMA_CALC_STANDARD,fcalcLocalIntegral,&
rtriangulation,rcubatureInfo=rmyCubature)
calculates three domain integrals at once into Dintvalues(:). The corresponding local contributions are calculated by fcalcLocalIntegral
. rmyCubature
defines the cubature rule to be used on the mesh rtriangulation
.
Remark: This form of the integration does not necessitate any finite element function and is rather seldom used. The original idea of the integral assembly is to calculate norms which involves finite element functions. For such an application, the routine is used in a slightly different manner.
Providing nonlinearities / FE functions / temporary memory
All assembly routines can be called with an optional parameter revalVectors
. In combination with the module feevaluation2.f90
, this parameter offers the possibility to specify finite element functions that live on the same mesh as on which the assembly is applied. All finite element functions given in this structure are automatically evaluated in all cubature points, and the calculated values are provided to the callback routines for further processing.
Example: A discrete stationary Oseen equation has the form $$Delta u_h + v_h nabla u_h + nabla p_h = f_h$$ with $u_h$, $v_h$ and $f_h$ vector fields. The operator $n(v_h,cdot):=v_h nabla cdot$ needs a finite element function $v_h$ for the assembly, which specifies the direction of a flow. Such a function can be provided to the callback routines via
revalVectors
.
Let a finite element function $v_h$ be given as a block vector structure rvectorVh
. This vector can be provided to the assembly routines as follows. We demonstrate the approach for the matrix assembly, however, for vector or integral assembly, the approach is exactly the same.
Providing a scalar function: Let $v_h=(v_1)$ be a block vector containing only one subvector. Such a subvector is provided to the assembly routines as follows:
type(t_fev2Vectors) :: rmyVectors
type(t_vectorBlock) :: rvectorVh
type(t_matrixBlock) :: rmatrix
...
call fev2_addVectorToEvalList(rmyVectors,rvectorVh%RvectorBlock(1),0)
call bma_buildMatrix (rmatrix,BMA_CALC_STANDARD,&
fcalcLocalMatrices,revalVectors=rmyVectors)
call fev2_releaseVectorList (rmyVectors)
There is no special initialisation routine for the structure rmyVectors
. One simply adds all scalar components of the vector $v_h$ to rmyVectors
and provides this structure to the bma_buildXXX
routine. The assembly routine will evaluate the provided functions automatically.
Derivatives: The last parameter
nmaxDerivative
in the routinefev2_addVectorToEvalList
is a derivative quantifier. It advises the evaluation routines what to evaluate.nmaxDerivative=0
(as specified above) calculates the function values in the cubature points, =1 the 1st derivatives, =2 the 2nd derivatives and so on. This can be used for more complicated operators which also involve derivatives of nonlinearities.
The following example assumes that $v_h=(v_1,v_2,v_3)$ has multiple components. It adds all components, one by one, and advises the assembly routine to calculate also their derivatives in the cubature points. So in the callback function the values of $v_h$ as well as the values of $nabla v_h$ are available.
type(t_fev2Vectors) :: rmyVectors
type(t_vectorBlock) :: rvectorVh
type(t_matrixBlock) :: rmatrix
...
call fev2_addVectorToEvalList(rmyVectors,rvectorVh%RvectorBlock(1),1)
call fev2_addVectorToEvalList(rmyVectors,rvectorVh%RvectorBlock(2),1)
call fev2_addVectorToEvalList(rmyVectors,rvectorVh%RvectorBlock(3),1)
call bma_buildMatrix (rmatrix,BMA_CALC_STANDARD,&
fcalcLocalMatrices,revalVectors=rmyVectors)
call fev2_releaseVectorList (rmyVectors)
Providing vector fields: Let $v_h=(v_1,v_2)$ be a block vector describing a vector field (for example, a velocity vector field of a flow in 2D). The following example provides the vector field to the evaluation routine via the function fev2_addVectorFieldToEvalList
. This subroutine expects all components of the vector field in the parameters and allows 1D, 2D and 3D vector fields:
type(t_fev2Vectors) :: rmyVectors
type(t_vectorBlock) :: rvectorVh
type(t_matrixBlock) :: rmatrix
...
call fev2_addVectorFieldToEvalList(rmyVectors,0,&
rvectorVh%RvectorBlock(1),rvectorVh%RvectorBlock(2))
call bma_buildMatrix (rmatrix,BMA_CALC_STANDARD,&
fcalcLocalMatrices,revalVectors=rmyVectors)
call fev2_releaseVectorList (rmyVectors)
Again, the parameter nmaxDerivative
(which is set to nmaxDerivative=0
in this example) specifies in how far derivatives from the vector field should be calculated. A call in the form
...
call fev2_addVectorFieldToEvalList(rmyVectors,1,&
rvectorVh%RvectorBlock(1),rvectorVh%RvectorBlock(2))
...
would calculate $v_h$ as well as $nabla v_h$ and provide it to the callback routine.
Providing temporary memory: Additionally to vectors and vector fields, the user can also specify "dummy vectors" and "dummy vector fields" in the t_fev2Vectors
structure. With this technique, it is possible to provide temporary memory in all cubature points on all elements to callback routines. Such temporary memory is needed in complicated operators for intermediate calculations which cannot be done in one simple step.
Example: A nonNewtonian fluid following the Power law model is given by the equation $$nu(x)Delta u + u nabla u + nabla p = f$$ with a nonlinear viscosity $$nu(x) = nu_0 z^{e/2  1},$$ $e,,nu_0inmathbb{R}$ a constant and $z$ be given by $$z=D(u)^2+varepsilon,qquad D(u)=frac{1}{2} (nabla u + nabla u^T),$$ with $varepsiloninmathbb{R}$ a constant. Computation routines for this model at first evaluate $nabla u$ and $nabla u^T$, then $D(u)$ and at the end $nu(x)$. For the intermediate steps, some temporary memory is necessary which can be provided by "dummy vectors".
The following routines provide temporary memory:
fev2_addDummyVectorToEvalList

provides one or multiple temporary arrays with memory in every cubature point on every element. The following example adds one temporary array to the structure
rmyVectors
.type(t_fev2Vectors) :: rmyVectors type(t_matrixBlock) :: rmatrix ... call fev2_addDummyVectorToEvalList(rmyVectors) call bma_buildMatrix (rmatrix,BMA_CALC_STANDARD,& fcalcLocalMatrices,revalVectors=rmyVectors) ...
fev2_addDummyVecFieldToEvalList

provides temporary memory for one or multiple vector fields with memory in every cubature point on every element. The following example adds one temporary vector field with 2 vector components to the structure
rmyVectors
.type(t_fev2Vectors) :: rmyVectors type(t_matrixBlock) :: rmatrix ... call fev2_addDummyVecFieldToEvalList(rmyVectors,2) call bma_buildMatrix (rmatrix,BMA_CALC_STANDARD,& fcalcLocalMatrices,revalVectors=rmyVectors) ...
Using the parameter
nsubarrays
allows to reserve also some subarrays associated to this temporary array. This is used, e.g., to allocate memory for a vector field as well as its derivatives. The following example allocates temporary memory for a temporary vector field $v_h=(v_1,v_2)$ and its derivative $nabla v_h$ in all cubature points on all elements. We neednsubarrays=3
, corresponding to the three terms $v_h$, $partial_x v_h$ and $partial_y v_h$.type(t_fev2Vectors) :: rmyVectors type(t_matrixBlock) :: rmatrix ... call fev2_addDummyVecFieldToEvalList(rmyVectors,2,3) call bma_buildMatrix (rmatrix,BMA_CALC_STANDARD,& fcalcLocalMatrices,revalVectors=rmyVectors) ...
fev2_addCellDummyToEvalList

provides one or multiple temporary arrays with memory for one information on every element (cell based information). The following example adds one temporary cell based array to the structure
rmyVectors
.type(t_fev2Vectors) :: rmyVectors type(t_matrixBlock) :: rmatrix ... call fev2_addCellDummyToEvalList(rmyVectors) call bma_buildMatrix (rmatrix,BMA_CALC_STANDARD,& fcalcLocalMatrices,revalVectors=rmyVectors) ...
fev2_addDummyVecFieldToEvalList

provides temporary memory for one or multiple vector fields with memory for one information on every element (cell based information). The following example adds one temporary vector field with 2 cell based vector components to the structure
rmyVectors
.type(t_fev2Vectors) :: rmyVectors type(t_matrixBlock) :: rmatrix ... call fev2_addCellVecFieldDummy(rmyVectors,2) call bma_buildMatrix (rmatrix,BMA_CALC_STANDARD,& fcalcLocalMatrices,revalVectors=rmyVectors) ...
Standard operators
The module blockmatassemblystdop.f90
provides a lot of standard implementations for the most important standard operators  matrices, vectors and integrals. These implemetations are given as callback routines which can either be provided to the bma_buildXXXX
routines directly, or which can be called inside of a callback function. The following naming conventions are used here:
bma_fcalc_XXXX
 Routines with this name can be used as callback function in a call to
bma_buildXXXX
. Parameters are passed via at_collection
structure and the vector evaluation structurerevalVectors
. There is no need for writing a userdefined callback function. It is a wrapper for the correspondingbma_docalc_XXXX
subroutine. bma_docalc_XXXX
 This is the immediate version of the
bma_fcalc_XXXX
routine from above. Routines of typebma_docalc_XXXX
can be called inside of a callback function ofbma_buildXXXX
to calculate a specific operator. This allows to easily compute multiple operators in one step.
Using bma_fcalc_XXXX
for simple calculations
With the routines bma_fcalc_XXXX
from blockmatassemblystdop.f90
one can easily assemble all types of matrices with very few commands. Parameters for these routines are passed via the "QuickAccess" arrays of the collection structure in a specified manner. In the following, we give some examples how these routines are used.
 Assembly of a Laplace matrix

The following code calculates a Laplace matrix at position (ix,iy) of a block matrix. We specify "1.0" as multiplier in front of the Laplace matrix:
type(t_matrixBlock) :: rmatrix type(t_collection) :: rcoll ... rcoll%DquickAccess(1) = 1.0_DP ! multiplier rcoll%IquickAccess(1) = ix rcoll%IquickAccess(2) = iy call bma_buildMatrix (& rmatrix,BMA_CALC_STANDARD,bma_fcalc_laplace, rcoll) ...
 Assembly of a tensor Laplace matrix

The following code calculates a tensor Laplace matrix at position (1,1) of a block matrix  i.e., a $dtimes d$ block matrix with the Laplace operator on the diagonal. The dimension $d$ of the tensor is taken from the underlying triangulation. We specify "1.0" as multiplier in front of the tensor Laplace matrix and use a tensor size of 2. This is typically used for the Stokes equations in 2D where the Stokes operator is exactly a $2times 2$ block matrix with Laplace on the diagonal:
type(t_matrixBlock) :: rmatrix type(t_collection) :: rcoll ... rcoll%DquickAccess(1) = 1.0_DP ! multiplier rcoll%IquickAccess(1) = 1 rcoll%IquickAccess(2) = 1 rcoll%IquickAccess(3) = 2 call bma_buildMatrix (& rmatrix,BMA_CALC_STANDARD,bma_fcalc_laplaceTensor, rcoll) ...
 Assembly of a Mass matrix

The following code calculates a mass matrix at position (ix,iy) of a block matrix. We specify "1.0" as multiplier in front of the Laplace matrix:
type(t_matrixBlock) :: rmatrix type(t_collection) :: rcoll ... rcoll%DquickAccess(1) = 1.0_DP ! multiplier rcoll%IquickAccess(1) = ix rcoll%IquickAccess(2) = iy call bma_buildMatrix (& rmatrix,BMA_CALC_STANDARD,bma_fcalc_mass, rcoll) ...
A tensor mass matrix can be assembled via
bma_fcalc_massTensor
, similar to the tensor Laplace matrix above.  Assembly of a "scalar" convection operator

The following code assembles the convection operator $n(u_h,v_h,w_h)=( (u_h nabla) v_h, w_h )$ at position (1,1) of a block matrix with convection (2.0,3.0). We apply the operator only to the scalar matrix at position (1,1) which realises the convection in a convectiondiffusionreaction equation of type $$Delta u + u nabla u + u = f.$$
type(t_matrixBlock) :: rmatrix type(t_collection) :: rcoll ... rcoll%IquickAccess(1) = 1 ! xposition rcoll%IquickAccess(2) = 1 ! yposition rcoll%IquickAccess(3) = 0 ! scalar destination rcoll%IquickAccess(4) = 0 ! constant direction rcoll%DquickAccess(1) = 1.0_DP ! multiplier rcoll%DquickAccess(2) = 2.0_DP ! xdirection rcoll%DquickAccess(3) = 3.0_DP ! ydirection call bma_buildMatrix (& rmatrix,BMA_CALC_STANDARD,bma_fcalc_convec_ugradvw, rcoll) ...
 Assembly of a convection operator

The following code assembles the convection operator $n(u_h,v_h,w_h)=( (u_h nabla) v_h, w_h )$ at position (1,1) of a block matrix. The nonlinearity $u_h$ is given as a finite element function, realised by the block vector
rvelocity
. We apply the nonlinearity to the full velocity tensor; the dimension is derived from the underlying triangulation, which is assumed to be 2D here.type(t_matrixBlock) :: rmatrix type(t_vectorBlock) :: rvelocity type(t_fev2Vectors) :: rmyVectors type(t_collection) :: rcoll ... ! Parameters rcoll%IquickAccess(1) = 1 ! xposition rcoll%IquickAccess(2) = 1 ! yposition rcoll%IquickAccess(3) = 1 ! tensor destination rcoll%IquickAccess(4) = 1 ! nonconstant direction rcoll%DquickAccess(1) = 1.0_DP ! multiplier ! Provide the vector field call fev2_addVectorFieldToEvalList(rmyVectors,0,& rvelocity%RvectorBlock(1),rvelocity%RvectorBlock(2)) call bma_buildMatrix (& rmatrix,BMA_CALC_STANDARD,bma_fcalc_convec_ugradvw, rcoll,& revalVectors = rmyVectors) call fev2_releaseVectorList (rmyVectors) ...
 Assembly of a constant RHS vector

The following code assembles a RHS vector to the constant RHS function $f=1$. The RHS is assumed to have only one component:
type(t_matrixBlock) :: rrhs type(t_collection) :: rcoll ... ! Parameters rcoll%IquickAccess(1) = 1 ! Number of components rcoll%DquickAccess(1) = 1.0_DP ! The value of the RHS function call bma_buildVector (& rmatrix,BMA_CALC_STANDARD,bma_fcalc_rhsConst, rcoll) ...
 Assembly of a $u_h_{L_2}$

The following code assembles the $L_2$ norm of a finite element function $u_h=(u_1)$ given by
rx
.type(t_vectorBlock) :: rx type(t_fev2Vectors) :: rmyVectors real(DP) :: dintegral ... ! Provide the vector call fev2_addVectorToEvalList(rmyVectors,rx%RvectorBlock(1),0) call bma_buildIntegral (& dintvalue,BMA_CALC_STANDARD,bma_fcalc_L2norm,revalVectors=rmyVectors) call fev2_releaseVectorList (rmyVectors) ...
Remark: If a finite element function is specified via
revalVectors
, a triangulation, boundary definition and cubature rule do not have to be specified. Although it may be advisable top specify a cubature rule, triangulation and boundary description is taken from the first vector in the vector list. This simplifies the call tobma_buildIntegral
.
 Assembly of a $u_h_{H^1}$

The following code assembles the $H^1$ norm of a finite element function $u_h=(u_1,u_2)$ given by
rx
. The callback routine needs the 1st derivatives of the function $u_h$, and thus,nmaxDerivative=1
is specified infev2_addVectorToEvalList
:type(t_vectorBlock) :: rx type(t_fev2Vectors) :: rmyVectors real(DP) :: dintegral ... ! Provide the vector and its derivative call fev2_addVectorToEvalList(rmyVectors,rx%RvectorBlock(1),1) call fev2_addVectorToEvalList(rmyVectors,rx%RvectorBlock(2),1) call bma_buildIntegral (& dintvalue,BMA_CALC_STANDARD,bma_fcalc_H1norm,revalVectors=rmyVectors) call fev2_releaseVectorList (rmyVectors) ...
 Assembly of $u_h_{infty}$

The following call uses an extended syntax of
bma_buildIntegral
to calculate the $sup$ norm of a finite element function. Normally,bma_buildIntegral
calculates integrals by summing up the contributions of all subdomains. Modifying the optional parametercoperation
,bma_buildIntegral
can be told to compute the maximum of the values of the subdomains. In combination withbma_fcalc_MAXnorm
, this allows to compute the maximum norm:type(t_vectorBlock) :: rx type(t_fev2Vectors) :: rmyVectors real(DP) :: dmaxnorm ... ! Provide the vector call fev2_addVectorToEvalList(rmyVectors,rx%RvectorBlock(1),0) call bma_buildIntegral (& dintvalue,BMA_CALC_STANDARD,bma_fcalc_MAXnorm,& revalVectors=rmyVectors,coperation=BMA_INT_MAX) call fev2_releaseVectorList (rmyVectors) ...
Using bma_docalc_XXXX
in callback routines
For more complex matrices and vectors, the user has to provide a userdefined callback routine that does the assembly. The actual calculation routines bma_docalc_XXXX
are designed to be called in such userdefined callback routines for the assembly of standard operators. Using these routines allows to assemble a lot of standard operators with high efficiency.
The following example demonstrates how to set up the operator $text{DR}:=Delta + I$ at position (1,1) of a block matrix. This operator realises the lefthand side of the diffusionreaction equation $Delta u + u = f$:
subroutine fcalc_DR(RmatrixData,rassemblyData,rmatrixAssembly,&
npointsPerElement,nelements,revalVectors,rcollection)
...
call bma_docalc_laplace(RmatrixData,rassemblyData,rmatrixAssembly,&
npointsPerElement,nelements,1.0_DP,1,1)
call bma_docalc_mass(RmatrixData,rassemblyData,rmatrixAssembly,&
npointsPerElement,nelements,1.0_DP,1,1)
end subroutine
...
type(t_matrixBlock) :: rmatrix
...
call bma_buildMatrix (rmatrix,BMA_CALC_STANDARD,fcalc_DR)
...
The callback routine just calls the bma_docalc_XXXX
routines for all operators to be assembled. For more complicated operators which probably involve a nonlinearity, parameters must be passed to the userdefined callback routine via the collection structure and the vector evaluation structure revalVectors
. The callback routine can then use this information and call the predefined calculation routines.
The following example demonstrates this approach. The callback routine assembles the operator for the diffusionconvectionequation $$nuDelta u_h + v_h nabla u_h = f_h.$$ The convection is given as a finite element function $v_h$ which is available inside of the callback function via revalVectors%p_RvectorData(1)
. The constant $nu$ is passed via the DquickAccess
array of the collection.
subroutine fcalc_DC(RmatrixData,rassemblyData,rmatrixAssembly,&
npointsPerElement,nelements,revalVectors,rcollection)
...
real(DP) :: dnu
! Get the parameter
dnu = rcollection%DquickAccess(1)
! Calculate nu*Laplace
call bma_docalc_laplace(RmatrixData,rassemblyData,rmatrixAssembly,&
npointsPerElement,nelements,dnu,1,1)
! Calculate the convection
call bma_docalc_convec_ugradvw(RmatrixData,rassemblyData,rmatrixAssembly,&
npointsPerElement,nelements,1.0_DP,1,1,&
rvectorField=revalVectors%p_RvectorData(1))
end subroutine
type(t_matrixBlock) :: rmatrix
type(t_vectorBlock) :: rvectorVh
type(t_fev2Vectors) :: rmyVectors
type(t_collection) :: rcoll
...
! Parameters
rcoll%DquickAccess(1) = nu ! multiplier
! Provide the nonlinearity
call fev2_addVectorToEvalList(rmyVectors,rvectorVh,0)
call bma_buildMatrix (&
rmatrix,BMA_CALC_STANDARD,fcalc_DC, rcoll,&
revalVectors = rmyVectors)
call fev2_releaseVectorList (rmyVectors)
...
Block matrices: Advanced operator assembly
For complex operators, it is necessary to use an appropriately designed callback routine for the matrix assembly. While standard operators can be assembled with the techniques above, the implementation of complex operators involves loops over all elements, cubature points as well as trial and test basis functions  and probably, the dimension of the FE space. However, there is no such a "general" case as it is for the "scalar" assembly. For optimal efficiency and a maximum of generality, one has to implement an operator four times, once for each of the following cases:
 Scalar valued trial/test FE space
 Vector valued trial/test FE space
 Scalar valued, interleaved trial/test FE space
 Vector valued, interleaved trial/test FE space
The four implementations are rather similar, however, differ in some small parts of the loops. One should note firthermore that some operators cannot be implemented in all cases. The user has to decide on which implementation is to choose for the desired case.
Remark: Also the default implementations in
blockmatassemblystdop.f90
do not cover all the above four cases for all types of finite element spaces. Most implementations cover only case 1 while only the most important standard operators (like Laplace and Mass) are designed for all cases. For the implementation of a new operator, the user is advised to do copy&paste the code for case 1 and later generalise this if necessary.
Input/output of the callback routine
The callback routine used for the calculation of matrices has the following interface:
subroutine fcalcLocalMatrices(RmatrixData,rassemblyData,rmatrixAssembly,&
npointsPerElement,nelements,revalVectors,rcollection)
The variables have the following meaning:
Variable  Description 

RmatrixData  An $mtimes n$ array of t_bmaMatrixData structures, one for each submatrix in the global matrix 
rmatrixData  Data necessary for the assembly. Contains determinants and cubature weights for the cubature,... 
rmatrixAssembly  Structure with all data about the assembly 
npointsPerElement  Number of points per element 
nelements  Number of elements 
revalVectors  Provided FEM routines, encapsules nonlinearities and temporary memory 
rcollection  Collection structure for userdefined parameters 
The callback routine has to use the values of the finite element trial and test spaces, cubature points/weights and the values of all nonlinearities/coefficients to compute local matrices for a set of nelements
elements. The main assembly routine will later on impose these local matrices into the global matrix. Depending on which of the above four cases the callback routine is used four, different variables of the above parameters can be used and have to be written to. Here a small overview about input/output variables in the four cases.
 Case 1: Scalar valued trial/test FE space

The following parameters hold:
Variable Input/Output Description RmatrixData(:,:)%p_Dentry Output The local matrix entries, to be computed rassemblyData%p_DcubWeight Input Cubature weights in all cubature points RmatrixData(:,:)%p_DbasTrial For all matrices: values of the corresponding trial basis functions in all cubature points RmatrixData(:,:)%p_DbasTest For all matrices: values of the corresponding test basis functions in all cubature points RmatrixData(:,:)%ndofTrial Number of degrees of freedom in the trial space RmatrixData(:,:)%ndofTest Number of degrees of freedom in the test space RmatrixData(:,:)%ndimfeTrial Input =1, indicates scalarvalued trial FE space RmatrixData(:,:)%ndimfeTest =1, indicates scalarvalued test FE space RmatrixData(:,:)%bisInterleaved =.false. here, the matrix is not interleaved p_Dentry (1:ndofTrial, 1:ndofTest, 1:nelements)
has to be filled by the callback routine with the entries of the local matrices.p_DcubWeight (1:ncubp, 1:nelements)
specify the cubature weights (including the Jacobian determinant of the mapping, etc.)p_DbasTrial (1:ndofTrial, :, 1:ncubp, 1:nelements)
andp_DbasTest (1:ndofTest, :, 1:ncubp, 1:nelements)
specify the values of the local basis functions. The 2nd dimension specifies a derivative quantifier DER_xxxx for the desired derivative.
 Case 2: Vector valued trial/test FE space

The following parameters hold:
Variable Input/Output Description RmatrixData(:,:)%p_Dentry Output The local matrix entries, to be computed rassemblyData%p_DcubWeight Input Cubature weights in all cubature points RmatrixData(:,:)%p_DbasTrial For all matrices: values of the corresponding trial basis functions in all cubature points RmatrixData(:,:)%p_DbasTest For all matrices: values of the corresponding test basis functions in all cubature points RmatrixData(:,:)%ndofTrial Number of degrees of freedom in the trial space RmatrixData(:,:)%ndofTest Number of degrees of freedom in the test space RmatrixData(:,:)%ndimfeTrial >1, Dimension of the trial FE space RmatrixData(:,:)%ndimfeTest >1, Dimension of the test FE space RmatrixData(:,:)%bisInterleaved Input =.false. here, the matrix is not interleaved p_Dentry (1:ndofTrial, 1:ndofTest, 1:nelements)
has to be filled by the callback routine with the entries of the local matrices.p_DcubWeight (1:ncubp, 1:nelements)
specify the cubature weights (including the Jacobian determinant of the mapping, etc.)p_DbasTrial (1:ndofTrial*ndimfeTrial, :, 1:ncubp, 1:nelements)
andp_DbasTest (1:ndofTest*ndimfeTest, :, 1:ncubp, 1:nelements)
specify the values of the local basis functions. The 2nd dimension specifies a derivative quantifier DER_xxxx for the desired derivative. More preciselyp_DbasTrial (1:ndofTrial + 0*ndofTrial , :,:,:)
= 1st dimension, trial spacep_DbasTrial (1:ndofTrial + 1*ndofTrial , :,:,:)
= 2nd dimension, trial spacep_DbasTrial (1:ndofTrial + 2*ndofTrial , :,:,:)
= 3rd dimension, trial space, ...p_DbasTest (1:ndofTest + 0*ndofTest , :,:,:)
= 1st dimension, test spacep_DbasTest (1:ndofTest + 1*ndofTest , :,:,:)
= 2nd dimension, test spacep_DbasTest (1:ndofTest + 2*ndofTest , :,:,:)
= 3rd dimension, test space, ...
 Case 3: Scalar valued, interleaved trial/test FE space

The following parameters hold:
Variable Input/Output Description RmatrixData(:,:)%p_DentryIntl Output The local matrix entries, to be computed rassemblyData%p_DcubWeight Input Cubature weights in all cubature points RmatrixData(:,:)%p_DbasTrial For all matrices: values of the corresponding trial basis functions in all cubature points RmatrixData(:,:)%p_DbasTest For all matrices: values of the corresponding test basis functions in all cubature points RmatrixData(:,:)%ndofTrial Number of degrees of freedom in the trial space RmatrixData(:,:)%ndofTest Number of degrees of freedom in the test space RmatrixData(:,:)%nvar Number of variables per matrix entry RmatrixData(:,:)%ndimfeTrial Input =1, indicates scalarvalued trial FE space RmatrixData(:,:)%ndimfeTest =1, indicates scalarvalued test FE space RmatrixData(:,:)%bisInterleaved =.true. here, the matrix is interleaved p_DentryIntl (1:nvar, 1:ndofTrial, 1:ndofTest, 1:nelements)
has to be filled by the callback routine with the entries of the local matrices.p_DcubWeight (1:ncubp, 1:nelements)
specify the cubature weights (including the Jacobian determinant of the mapping, etc.)p_DbasTrial (1:ndofTrial, :, 1:ncubp, 1:nelements)
andp_DbasTest (1:ndofTest, :, 1:ncubp, 1:nelements)
specify the values of the local basis functions. The 2nd dimension specifies a derivative quantifier DER_xxxx for the desired derivative.
 Case 4: Vector valued, interleaved trial/test FE space

The following parameters hold:
Variable Input/Output Description RmatrixData(:,:)%p_DentryIntl Output The local matrix entries, to be computed rassemblyData%p_DcubWeight Input Cubature weights in all cubature points RmatrixData(:,:)%p_DbasTrial For all matrices: values of the corresponding trial basis functions in all cubature points RmatrixData(:,:)%p_DbasTest For all matrices: values of the corresponding test basis functions in all cubature points RmatrixData(:,:)%ndofTrial Number of degrees of freedom in the trial space RmatrixData(:,:)%ndofTest Number of degrees of freedom in the test space RmatrixData(:,:)%ndimfeTrial >1, Dimension of the trial FE space RmatrixData(:,:)%ndimfeTest >1, Dimension of the test FE space RmatrixData(:,:)%nvar Number of variables per matrix entry RmatrixData(:,:)%bisInterleaved Input =.true. here, the matrix is interleaved p_DentryIntl (1:nvar, 1:ndofTrial, 1:ndofTest, 1:nelements)
has to be filled by the callback routine with the entries of the local matrices.p_DcubWeight (1:ncubp, 1:nelements)
specify the cubature weights (including the Jacobian determinant of the mapping, etc.)p_DbasTrial (1:ndofTrial*ndimfeTrial, :, 1:ncubp, 1:nelements)
andp_DbasTest (1:ndofTest*ndimfeTest, :, 1:ncubp, 1:nelements)
specify the values of the local basis functions. The 2nd dimension specifies a derivative quantifier DER_xxxx for the desired derivative. More preciselyp_DbasTrial (1:ndofTrial + 0*ndofTrial , :,:,:)
= 1st dimension, trial spacep_DbasTrial (1:ndofTrial + 1*ndofTrial , :,:,:)
= 2nd dimension, trial spacep_DbasTrial (1:ndofTrial + 2*ndofTrial , :,:,:)
= 3rd dimension, trial space, ...p_DbasTest (1:ndofTest + 0*ndofTest , :,:,:)
= 1st dimension, test spacep_DbasTest (1:ndofTest + 1*ndofTest , :,:,:)
= 2nd dimension, test spacep_DbasTest (1:ndofTest + 2*ndofTest , :,:,:)
= 3rd dimension, test space, ...
Basic structure of the callback routine
The actual structure of the callback routine is now demonstrated on an example:
Example: The following example demonstrates how to set up a Laplace operator in 2D at position (1,1) of a block matrix, here for scalarvalued FEM spaces like $Q_1$, $Q_2$, etc:
subroutine fcalc_laplace(RmatrixData,rassemblyData,rmatrixAssembly,&
npointsPerElement,nelements,revalVectors,rcollection)
...
real(DP) :: dbasIx, dbasJx, dbasIy, dbasJy
integer :: iel, icubp, idofe, jdofe, ndofTrial, ndofTest
real(DP), dimension(:,:,:), pointer :: p_DlocalMat
real(DP), dimension(:,:,:,:), pointer :: p_DbasTrial,p_DbasTest
real(DP), dimension(:,:), pointer :: p_DcubWeight
! Get cubature weights data
p_DcubWeight => rassemblyData%p_DcubWeight
! Get local data
ndofTrial = RmatrixData(1,1)%ndofTrial
ndofTest = RmatrixData(1,1)%ndofTrial
p_DbasTrial => RmatrixData(1,1)%p_DbasTrial
p_DbasTest => RmatrixData(1,1)%p_DbasTest
p_DlocalMat => RmatrixData(1,1)%p_Dentry
! Loop over the elements in the current set.
do iel = 1,nelements
! Loop over all cubature points on the current element
do icubp = 1,npointsPerElement
! Loop over the test functions
do idofe=1,ndofTest
! Get the values of the test basis functions
dbasIx = p_DbasTest(idofe,DER_DERIV2D_X,icubp,iel)
dbasIy = p_DbasTest(idofe,DER_DERIV2D_Y,icubp,iel)
! Loop over the trial basis functions
do jdofe=1,ndofTrial
! Get the values of the trial basis functions
dbasJx = p_DbasTrial(jdofe,DER_DERIV2D_X,icubp,iel)
dbasJy = p_DbasTrial(jdofe,DER_DERIV2D_Y,icubp,iel)
! Calculate the local matrix entries
p_DlocalMat(jdofe,idofe,iel) = p_DlocalMat(jdofe,idofe,iel) + &
p_DcubWeight(icubp,iel) * ( dbasJx*dbasIx + dbasJy*dbasIy )
end do ! jdofe
end do ! idofe
end do ! icubp
end do ! iel
end subroutine
In the above example it can be seen that setting up a Laplace matrix needs a callback routine with altogether four nested do loops:
 An outer DO loop over elements
 An inner DO loop over cubature points
 For every cubature point, an inner loop over the DOFs in the test space
 For every DOF in the test space, an inner loop over the DOFs in the trial space
In the innerst loop, the actual matrix entries are computed using cubature. The variables dbasIx
, dbasIy
, dbasJx
and dbasJy
fetch the values of the basis functions in the cubature points. The sum in the innerst loop
...
p_DlocalMat(...) = p_DlocalMat(...) + &
p_DcubWeight(...) * ( dbasJx*dbasIx + dbasJy*dbasIy )
...
realises the summation $$sum_{x_k}sum_{i,j} omega(x_k) (partial_{x_1}psi_j(x_k) partial_{x_1}varphi_i(x_k) + partial_{x_2}psi_j(x_k) partial_{x_2}varphi_i(x_k) )$$ with $x_k$ the cubature points, $psi_j$ the local test functions and $varphi_i$ the local trial functions on the cell $T$. This is the discrete counterpart to the integral $$(nabla psi_j, nabla varphi_i)T = int_T partial{x_1}psi_j partial_{x_1}varphi_i + partial_{x_2}psi_j partial_{x_2}varphi_i$$ which realises the local Laplace operator on a cell $T$ in 2D. Summing up the contributions of all elements gives the global Laplace operator.
Remark: For the actual assembly, the order of the summation is changed. The above loops do not calculate $(nabla psi_j, nabla varphi_i)_T$ but only its contribution in one cubature point for one combination of trial/test functions, summing up later. In this order, the assembly is faster.
Block vectors: Advanced vector assembly
THe assembly of complex righthand side vectors is rather similar to the assembly of complex vectors. One needs an appropriately designed callback routine which applies loops over all elements, cubature points as well as trial and test basis functions. For optimal efficiency and a maximum of generality, also here one has to implement an operator four times, once for each of the following cases:
 Scalar valued trial/test FE space
 Vector valued trial/test FE space
 Scalar valued, interleaved trial/test FE space
 Vector valued, interleaved trial/test FE space
Input/output of the callback routine
The callback routine used for the calculation of matrices has the following interface:
subroutine fcalcLocalVectors(rvectorData,rassemblyData,rvectorAssembly,&
npointsPerElement,nelements,revalVectors,rcollection)
The variables have the following meaning:
Variable  Description 

RvectorData  An array of t_bmaVectorData structures, one for each subvector in the global vector 
rassemblyData  Data necessary for the assembly. Contains determinants and cubature weights for the cubature,... 
rvectorAssembly  Structure with all data about the assembly 
npointsPerElement  Number of points per element 
nelements  Number of elements 
revalVectors  Provided FEM routines, encapsules nonlinearities and temporary memory 
rcollection  Collection structure for userdefined parameters 
The callback routine has to use the values of the finite element test spaces, cubature points/weights and the values of all nonlinearities/coefficients to compute local vectors for a set of nelements
elements. The main assembly routine will later on impose these local vectors into the global vector. Depending on which of the above four cases the callback routine is used four, different variables of the above parameters can be used and have to be written to. Here a small overview about input/output variables in the four cases.
 Case 1: Scalar valued test FE space

The following parameters hold:
Variable Input/Output Description RvectorData(:)%p_Dentry Output The local vector entries, to be computed rassemblyData%p_DcubWeight Input Cubature weights in all cubature points RvectorData(:)%p_DbasTest For all vectors: values of the corresponding test basis functions in all cubature points RvectorData(:)%ndofTest Number of degrees of freedom in the test space RvectorData(:)%ndimfeTest Input =1, indicates scalarvalued test FE space RvectorData(:)%bisInterleaved =.false. here, the vector is not interleaved p_Dentry (1:ndofTest, 1:nelements)
has to be filled by the callback routine with the entries of the local vectors.p_DcubWeight (1:ncubp, 1:nelements)
specify the cubature weights (including the Jacobian determinant of the mapping, etc.)p_DbasTest (1:ndofTest, :, 1:ncubp, 1:nelements)
specify the values of the local basis functions. The 2nd dimension specifies a derivative quantifier DER_xxxx for the desired derivative.
 Case 2: Vector valued test FE space

The following parameters hold:
Variable Input/Output Description RvectorData(:)%p_Dentry Output The local vector entries, to be computed rassemblyData%p_DcubWeight Input Cubature weights in all cubature points RvectorData(:)%p_DbasTest For all vectors: values of the corresponding test basis functions in all cubature points RvectorData(:)%ndofTest Number of degrees of freedom in the test space RvectorData(:)%ndimfeTest >1, Dimension of the test FE space RvectorData(:)%bisInterleaved Input =.false. here, the vector is not interleaved p_Dentry (1:ndofTest, 1:nelements)
has to be filled by the callback routine with the entries of the local vectors.p_DcubWeight (1:ncubp, 1:nelements)
specify the cubature weights (including the Jacobian determinant of the mapping, etc.)p_DbasTest (1:ndofTest*ndimfeTest, :, 1:ncubp, 1:nelements)
specify the values of the local basis functions. The 2nd dimension specifies a derivative quantifier DER_xxxx for the desired derivative. More preciselyp_DbasTest (1:ndofTest + 0*ndofTest , :,:,:)
= 1st dimension, test spacep_DbasTest (1:ndofTest + 1*ndofTest , :,:,:)
= 2nd dimension, test spacep_DbasTest (1:ndofTest + 2*ndofTest , :,:,:)
= 3rd dimension, test space, ...
 Case 3: Scalar valued, interleaved test FE space

The following parameters hold:
Variable Input/Output Description RvectorData(::)%p_DentryIntl Output The local vector entries, to be computed rassemblyData%p_DcubWeight Input Cubature weights in all cubature points RvectorData(:)%p_DbasTest For all vectors: values of the corresponding test basis functions in all cubature points RvectorData(:)%ndofTest Number of degrees of freedom in the test space RvectorData(:)%nvar Number of variables per vector entry RvectorData(:)%ndimfeTest Input =1, indicates scalarvalued test FE space RvectorData(:)%bisInterleaved =.true. here, the vector is interleaved p_DentryIntl (1:nvar, 1:ndofTest, 1:nelements)
has to be filled by the callback routine with the entries of the local vectors.p_DcubWeight (1:ncubp, 1:nelements)
specify the cubature weights (including the Jacobian determinant of the mapping, etc.)p_DbasTest (1:ndofTest, :, 1:ncubp, 1:nelements)
specify the values of the local basis functions. The 2nd dimension specifies a derivative quantifier DER_xxxx for the desired derivative.
 Case 4: Vector valued, interleaved test FE space

The following parameters hold:
Variable Input/Output Description RvectorData(:)%p_DentryIntl Output The local vector entries, to be computed rassemblyData%p_DcubWeight Input Cubature weights in all cubature points RvectorData(:)%p_DbasTest For all vectors: values of the corresponding test basis functions in all cubature points RvectorData(:)%ndofTest Number of degrees of freedom in the test space RvectorData(:)%ndimfeTest >1, Dimension of the test FE space RvectorData(:)%nvar Number of variables per vector entry RvectorData(:)%bisInterleaved Input =.true. here, the vector is interleaved p_DentryIntl (1:nvar, 1:ndofTest, 1:nelements)
has to be filled by the callback routine with the entries of the local vectors.p_DcubWeight (1:ncubp, 1:nelements)
specify the cubature weights (including the Jacobian determinant of the mapping, etc.)p_DbasTest (1:ndofTest*ndimfeTest, :, 1:ncubp, 1:nelements)
specify the values of the local basis functions. The 2nd dimension specifies a derivative quantifier DER_xxxx for the desired derivative. More preciselyp_DbasTest (1:ndofTest + 0*ndofTest , :,:,:)
= 1st dimension, test spacep_DbasTest (1:ndofTest + 1*ndofTest , :,:,:)
= 2nd dimension, test spacep_DbasTest (1:ndofTest + 2*ndofTest , :,:,:)
= 3rd dimension, test space, ...
Basic structure of the callback routine
The actual structure of the callback routine is now demonstrated on an example:
Example: The following example demonstrates how to set up the righthand side $f=x_1 x_2$ in 2D for the first subvector of a block vector, here for scalarvalued FEM spaces like $Q_1$, $Q_2$, etc:
subroutine fcalc_rhsxy(rvectorData,rassemblyData,rvectorAssembly,&
npointsPerElement,nelements,revalVectors,rcollection)
...
! Local variables
integer :: iel, icubp, idofe, dbasI, dval
real(DP), dimension(:,:), pointer :: p_DlocalVector
real(DP), dimension(:,:,:,:), pointer :: p_DbasTest
real(DP), dimension(:,:), pointer :: p_DcubWeight
type(t_bmaVectorData), pointer :: p_rvectorData
real(DP), dimension(:,:,:), pointer :: p_Dpoints
! Get cubature weights data
p_DcubWeight => rassemblyData%p_DcubWeight
! Get the coordinates of the cubature points
p_Dpoints => rassemblyData%revalElementSet%p_DpointsReal
! Get the data arrays of the subvector
p_rvectorData => RvectorData(icomp)
p_DlocalVector => p_rvectorData%p_Dentry
p_DbasTest => p_rvectorData%p_DbasTest
! Loop over the elements in the current set.
do iel = 1,nelements
! Loop over the cubature points
do icubp = 1,npointsPerElement
! Get the coordinates of the cubature point.
dx = p_Dpoints(1,icubp,iel)
dy = p_Dpoints(2,icubp,iel)
! Loop over the test space
do idofe=1,p_rvectorData%ndofTest
! Get the value of the basis function
dbasI = p_DbasTest(idofe,DER_FUNC,icubp,iel)
! Value of the function f=x_1 x_2
dval = dx * dy
! Multiply the values of the basis functions
! (1st derivatives) by the cubature weight and sum up
! into the local vectors.
p_DlocalVector(idofe,iel) = p_DlocalVector(idofe,iel) + &
p_DcubWeight(icubp,iel) * dval * dbasI
end do ! jdofe
end do ! icubp
end do ! iel
end subroutine
In the above example it can be seen that setting up a righthand side needs a callback routine with altogether three nested do loops:
 An outer DO loop over elements
 An inner DO loop over cubature points
 For every cubature point, an inner loop over the DOFs in the test space
In the innerst loop, the actual vector entries are computed using cubature. The variables dbasI
fetches the values of the basis functions in the cubature points. The sum in the innerst loop
...
p_DlocalVector(...) = p_DlocalVector(...) + &
p_DcubWeight(...) * dval * dbasI
...
realises the summation $$sum_{x_k}sum_{i} omega(x_k) f(x_k) psi_i(x_k) $$ with $x_k$ the cubature points and $psi_i$ the test functions on the cell $T$. This is the discrete counterpart to the integral $$(f, psi_i)_T = int_T f psi_i $$ which realises the local righthand side on a cell $T$ in 2D. Summing up the contributions of all elements gives the global righthand side.
Remark: For the actual assembly, the order of the summation is changed. The above loops do not calculate $(f, psi_i)_T$ but only its contribution in one cubature point for one test functions, summing up later. In this order, the assembly is faster.
Accessing FE functions and temporary arrays
As introduced above, it is possible to pass finite element functions an temporary arrays to the callback routine in order to do complex calculations there. To demonstrate the exact passing strategy, we start with an example based on a matrix assembly.
Example: Assume that there are two finite element functions to be passed to a callback routine as well as a set of temporary arrays.
The first finite element function, $v_h=(v_1)$, is a scalar finite element function, denoted by
rvectorVh
for which we need $v_h$ and $nabla v_h$.The second one, $w_h=(w_1,w_2,w_3)$, is a vector field for which we also need $w_h$ and $nabla w_h$.
Furthermore, we need a temporary arrays where we store intermediate information for processing $v_h$, $nabla v_h$, $w_h$ and $nabla w_h$.
Finally, we need a temporary arrays that store intermediate information for processing $v_h$, $nabla v_h$, $w_h$ and $nabla w_h$.
The main program sets up a vector evaluation structure as follows. Herein, the dummy vector fields use "4" temporary vectors for each dimension, one for the function values and three for the x/y/z derivatives of every component  resulting in 12 temporary arrays in total.
...
type(t_vectorBlock) :: rvectorVh
type(t_vectorBlock) :: rvectorWh
type(t_fev2Vectors) :: rmyVectors
...
! V_h, nabla V_h
call fev2_addVectorToEvalList(&
rmyVectors,rvectorVh%RvectorBlock(1),1)
! W_h, nabla W_h
call fev2_addVectorFieldToEvalList(rmyVectors,1,&
rvectorWh%RvectorBlock(1),&
rvectorWh%RvectorBlock(2),&
rvectorWh%RvectorBlock(3))
! Temporary memory, corresponding to V_h, nabla V_h.
call fev2_addDummyVectorToEvalList (rmyVectors,4)
! Temporary memory, corresponding to W_h, nabla W_h.
call fev2_addDummyVecFieldToEvalList (rmyVectors,3,4)
! Cell based temp memory, corresponding to V_h, nabla V_h.
call fev2_addDummyVectorToEvalList (rmyVectors,4)
! Cell based temp memory, corresponding to W_h, nabla W_h.
call fev2_addDummyVecFieldToEvalList (rmyVectors,3,4)
! Assembly
call bma_buildMatrix (rmatrix,BMA_CALC_STANDARD,&
fcalcMatrix,revalVectors=rmyVectors)
! Cleanup
call fev2_releaseVectorList (rmyVectors)
...
The callback routine fcalcMatrix
which is called in bma_buildMatrix
receives a parameter revalVectors=rmyVectors
which contains the vectors and temporary arrays in exactly the order they were added to rmyVectors
. The data is accessible via revalVectors%p_RvectorData
which is an array of all data added to rmyVectors
. revalVectors%p_RvectorData
contains data arrays for all cubature points and all elements that are currently in process of the assembly.
The content of the structure revalVectors%p_RvectorData(i)
depends on the type of the data which is represented. The callback routine has to "know" which data is added to revalVectors
and use it in the appropriate order. In the following tables, there will be an overview which data is accessible in which situation:
 Structural data

The following table depicts general information that describe the type of data represented by
revalVectors%p_RvectorData(i)
:Variable Description bisInterleaved Specifies if the data is interleaved bisVectorField Specifies if the data is a vector field bisCellBased Specifies if the data is cell based ndimVectorField Dimension of the vector field. nvar Number of entries per vector if the vector is interleaved ndimfe Dimension of the underlying FEM space (for vector valued FEM spaces. =1 for standard FEM spaces) nmaxDerivativeIdx Number of allocated subarrays. This corresponds to the maximum derivative quantifier DER_xxxx
that can be used  Scalar data in every cubature point

If the structure represents a scalar FE function in every cubature point:
Variable Description p_Ddata Array with values in every cubature point on every element  Vector data in every cubature point

If the structure represents a fector field FE function in every cubature point:
Variable Description p_DdataVec Array with values in every dimension, cubature point and on every element  Interleaved scalar data in every cubature point

If the structure represents an interleaved scalar FE function in every cubature point:
Variable Description p_DdataIntl Array with values in every dimension, cubature point and on every element  Scalar data on every element

If the structure represents one data entry per element:
Variable Description p_DcellData Array with one value per element  Vector data on every element

If the structure represents one vector field data entry per element:
Variable Description p_DcellDataVec Array with one value in every dimension on every element  Interleaved scalar data on every element

If the structure represents an one interleaved data entry per element:
Variable Description p_DcellDataIntl Array with one value per element
For vector data added by fev2_addVector_XXXX
, the data is initialised with the values in the cubature points on all elements. "Dummy" arrays are treated the same way as data arrays but are not initialised; they represent uninitialised temporary memory the callback routine can use arbitrarily.
Example: A callback routine
fcalcMatrix
corresponding to the above code may look like this:
subroutine fcalcMatrix (RmatrixData,rassemblyData,rmatrixAssembly,&
npointsPerElement,nelements,revalVectors,rcollection)
...
real(DP), dimension(:,:,:), pointer :: p_DdataVh
real(DP), dimension(:,:,:,:), pointer :: p_DdataWh
real(DP), dimension(:,:,:), pointer :: p_DdataTempVh
real(DP), dimension(:,:,:,:), pointer :: p_DdataTempWh
real(DP), dimension(:,:), pointer :: p_DdataCellVh
real(DP), dimension(:,:,:), pointer :: p_DdataCellWh
! Get the data arrays
p_DdataVh => revalVectors%p_RvectorData(1)%p_Ddata
p_DdataWh => revalVectors%p_RvectorData(2)%p_DdataVec
p_DdataTempVh => revalVectors%p_RvectorData(1)%p_Ddata
p_DdataTempWh => revalVectors%p_RvectorData(2)%p_DdataVec
p_DdataCellVh => revalVectors%p_RvectorData(1)%p_DcellData
p_DdataCellWh => revalVectors%p_RvectorData(2)%p_DcellDataVec
...
do iel = 1,nelements
do icubp = 1,npointsPerElement
! Get the value of V_h in the cubature point
dvh = p_DdataVh(icubp,iel,DER_FUNC)
! Save for later use
p_DdataTempVh(icubp,iel,DER_FUNC) = dvh**2
! Sum up to norm over the cell
p_DdataCellVh(iel,DER_FUNC) = p_DdataCellVh(iel,DER_FUNC) + dvh**2
do idofe=1,ndofTest
dbasI = p_DbasText(idofe,DER_FUNC,icubp,iel)
do jdofe=1,ndofTrial
dbasJ = p_DbasTrial(jdofe,DER_FUNC,icubp,iel)
! Calculate the local matrix entries
p_DlocalMat(jdofe,idofe,iel) = p_DlocalMat(jdofe,idofe,iel) + &
p_DcubWeight(icubp,iel) * ( dvh*dbasJ*dbasI )
end do ! jdofe
end do
end do
end do
...
end subroutine
This code calculates the bilinear form $$m(v_h)(varphi_h,psi_h) = ( psi_h, v_h varphi_h )$$ which corresponds to the operator $v_h I$ with $I$ the identity operator (resulting in a mass matrix with nonconstant coefficient).
Using nonconstant data: The data arrays p_DdataXXXX
in this subroutine are extracted from the revalVectors%p_RvectorData(:)
structures in exactly the order they were added to rmyVectors
above. The shape of the data arrays correspond to the type of data that was added and is designed for fast data access inside of the loops. In the above example, p_DdataVh
contains the values of $v_1$ in all cubature points and all elements. As soon as iel
and icubp
identify one cubature point on one element, the value dvh
of $v_1$ can be obtained via p_DdataVh
and be used in all the loops over the test and trial functions that calculate the local matrices.
Using temporary arrays: The temp arrays can be arbitrarily used in the callback routine. In the above example, the value of $v_h^2$ is saved to p_DdataTempVh
which has the same layout as p_DdataVh
. The data can then be accessed in a later step of the assembly.
Cell based data is used in a similar way. Such temporary data provides one data entry per cell. The above example uses, e.g., p_DdataCellVh
to compute the norm of $v_h$ on a cell $T$, $v_h_{L_2(T)}$, which may be used later, for example in a stabilisation technique.