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=0results in the original cubature rule.nlevels=1results in a sum about 2 intervals.nlevels=2results in a sum about 4 intervals.nlevels=4results in a sum about 8 intervals.
etc. Similarly in 2D with quadrilateral elements,
nlevels=0results in the original cubature rule.nlevels=1results in a sum about 4 subelements in each element.nlevels=2results in a sum about 16 subelements in each element.nlevels=3results 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.
