# Cubature

All assembly routines in FEAT2 have a build-in logic to choose a proper cubature formula, which is usually a strong Gauss Formula. In many cases however, the Gauss formula is much to strong, so a lot of computational time can be saved by choosing an appropriate cubature formula for the specific situation and finite element space.

## Basic cubature formulas

The basic cubature support is realised in the module `cubature.f90` which provides a list of IDs for all implemented cubature formulas. The following tables contain a brief overview about IDs which are commonly used:

1D cubature formulas:

ID cell type degree #cub.pts Formula
CUB_G1_1D interval 2 1 Gauss 1-pt / Midpoint rule
CUB_TRZ_1D 2 2 Trapezoidal rule
CUB_G2_1D 4 2 Gauss 2-pt
CUB_G3_1D 6 3 Gauss 3-pt
CUB_G4_1D 8 4 Gauss 4-pt
CUB_G5_1D 10 5 Gauss 5-pt

2D cubature formulas:

ID cell type degree #cub.pts Formula
CUB_G1_T triangle 2 1 Gauss 1-pt
CUB_TRZ_T 2 3 Trapezoidal rule
CUB_G3_T 3 3 Gauss 3-pt
CUB_G1_2D quad 2 1 Gauss 1-pt
CUB_TRZ_2D 2 4 Trapezoidal rule
CUB_MID_2D 2 4 (Edge) Midpoint rule
CUB_G2_2D 4 4 Gauss 2x2
CUB_G3_2D 6 9 Gauss 3x3
CUB_G4_2D 8 16 Gauss 4x4
CUB_G5_2D 10 25 Gauss 5x5

3D cubature formulas:

ID cell type degree #cub.pts Formula
CUB_G1_3D_T tetra 2 1 Gauss 1-pt
CUB_TRZ_3D_T 2 4 Trapezoidal rule
CUB_S2_3D_T 2 4 Stroud 4-pt
CUB_S3_3D_T 3 10 Stroud 10-pt
CUB_S5_3D_T 5 15 Stroud 15-pt
CUB_G1_3D hex 2 1 Gauss 1-pt
CUB_MIDAREA_3D 2 6 (Face) Midpoint rule
CUB_TRZ_3D 2 8 Trapezoidal rule
CUB_G2_3D 4 8 Gauss 2x2x2
CUB_G3_3D 6 27 Gauss 3x2x3
CUB_G4_3D 8 64 Gauss 4x4x4
CUB_G5_3D 10 125 Gauss 5x5x5

## Automatic cubature formulas

Additionally to the above specific cubature formulas, there exist a number of "automatic" cubature rules. These cubature rules are chosen automatically depening on the underlying finite element space and the object which is to be assembled. The set of automatic cubature rules is restricted to those cubature rules which are element and dimension independent. The following list gives an overview about the corresponding IDs and formulas:

ID Formula
CUB_GEN_AUTO_G1 1-point Gauss formula
CUB_GEN_AUTO_G2 2-point Gauss formula
CUB_GEN_AUTO_G3 3-point Gauss formula
CUB_GEN_AUTO_G4 4-point Gauss formula
CUB_GEN_AUTO_G5 5-point Gauss formula
CUB_GEN_AUTO_G6 6-point Gauss formula
CUB_GEN_AUTO_TRZ Trapezoidal rule
CUB_GEN_AUTO_MID Midpoint rule
CUB_GEN_AUTO_SIMPSON Simpson rule
CUB_GEN_AUTO_LUMPMASS Cubature rule that lumps a mass matrix (if available)
CUB_GEN_AUTO General automatic cubature formula, chosen depending on the dimension and type of operator to be assembled. This is the default.

For most cases, the user will specify one of the above automatic cubature formulas for the assembly. This gives a large simplification in particular for conformal meshes, where different cubature rules have to be chosen depending on the type of the underlying cell. Automatic cubature formulas choose the corresponding cubature rule automatically depending on the cell type.

## Using cubature formulas in the assembly

To enforce a specific cubature rule during the assembly of a right-hand side or a matrix, a structure of type `t_scalarCubatureInfo` from the module `spatialdiscretsation.f90` can be specified to any assembly routine. This structure defines for every utilised finite element space a corresponding cubature formula. For creating the structure, the routine `spdiscr_createDefCubStructure` is used, which sets up a cubature structure depending on a finite element space, a mesh and (possibly) the type of the operator to be assembled. Any of the above cubature IDs can be specified here. The structure can be destroyed using `spdiscr_releaseCubStructure`.

The following example demonstrates this approach. It creates a cubature structure for the Gauss 3-pt rule and passes this to the assembly routine for right-hand side vectors. At the end, the structure is released:

``````type(t_spatialDiscr) :: rdiscr
type(t_vectorScalar) :: rrhs
type(t_scalarCubatureInfo) :: rcubature

...
! Define the cubature rule
call spdiscr_createDefCubStructure (rdiscr,rcubature,CUB_GEN_AUTO_G3)

! Assemble the vector
call linf_buildVectorScalar (...,rrhs,rcubatureInfo = rcubature,...)

! Release the cubature structure
call spdiscr_releaseCubStructure (rcubatureInfo)
``````

As can be seen from the example, the structure `rcubature` is specified as additional (optional) parameter `rcubatureInfo` in the call to the assembly routine `linf_buildVectorScalar`. The cubature formula `CUB_GEN_AUTO_G3` tells `linf_buildVectorScalar` to automatically choose an appropriate 3-point Gauss formula (3x3 for quads in 2D, 3x3x3 for hexas in 3D) for the cubature.

## Summed cubature rules

The routine `spdiscr_createDefCubStructure` accepts an optional parameter `nlevels` which allows to specify "summed cubature formulas". For a summed cubature formula, the element on which the cubature is applied is regularly refined into subelements, and the actual cubature formula is aplied piecewise. For example, in 1D:

• `nlevels=0` results in the original cubature rule.
• `nlevels=1` results in a sum about 2 intervals.
• `nlevels=2` results in a sum about 4 intervals.
• `nlevels=4` results in a sum about 8 intervals.

etc. Similarly in 2D with quadrilateral elements,

• `nlevels=0` results in the original cubature rule.
• `nlevels=1` results in a sum about 4 subelements in each element.
• `nlevels=2` results in a sum about 16 subelements in each element.
• `nlevels=3` results in a sum about 64 subelements in each element.

etc. This technique can be used to improve integrals in cells with large gradients, e.g., in the region of shocks.

The level parameter is directly specified in the call to `spdiscr_createDefCubStructure`. The following example demonstrates how to set p a right-hand side vector with a summed Gauss 3-point rule where each element is refined 4 times:

``````type(t_spatialDiscr) :: rdiscr
type(t_vectorScalar) :: rrhs
type(t_scalarCubatureInfo) :: rcubature

...
! Define the summed G3 cubature rule, 4 levels of refinement
call spdiscr_createDefCubStructure (rdiscr,rcubature,CUB_GEN_AUTO_G3,4)

! Assemble the vector
call linf_buildVectorScalar (...,rrhs,rcubatureInfo = rcubature,...)

! Release the cubature structure
call spdiscr_releaseCubStructure (rcubatureInfo)
``````

The following picture demonstrates a summed 2x2 Gauss formula for level `nlevels=2`. Every cell is locally refined two times, thus containing 16 subelements. The Gauss formula is applied on each of the subelements, resulting in 64 cubature points in total on each element. 