cubatures Module

Numerical integration rules (cubatures) on the reference elements of the finite element method.

A cubature holds the abscissae $\xi_q$ and weights $w_q$ of a rule such that $$\int_{\hat\Omega} f(\xi)\,d\xi \approx \sum_{q=1}^{n} w_q\, f(\xi_q)$$ on the reference element $\hat\Omega$.

Element types and reference domains

Constant	Element	Reference domain	Measure
CUB_LIN	line	$[-1,1]$	2
CUB_TRI	triangle	$x,y \ge 0,\ x+y \le 1$	1/2
CUB_QUA	quadrilateral	$[-1,1]^2$	4
CUB_TET	tetrahedron	$x,y,z \ge 0,\ x+y+z \le 1$	1/6
CUB_HEX	hexahedron	$[-1,1]^3$	8
CUB_WED	wedge / prism	triangle $\times\ [-1,1]$	1

Triangle and tetrahedron abscissae are the Cartesian coordinates on the unit simplex, which equal the leading barycentric (area / volume) coordinates; the last barycentric coordinate is $1 - \sum_i \xi_i$.

Degree of exactness

Rules are selected by their degree of exactness $p$:

• CUB_TRI, CUB_TET: all polynomials of total degree $\le p$.
• CUB_LIN, CUB_QUA, CUB_HEX: all polynomials of degree $\le p_i$ in each coordinate $i$ (tensor-product Gauss–Legendre with $p_i/2 + 1$ points per direction). One degree or one per direction.
• CUB_WED: total degree $\le p_1$ in the triangle coordinates and degree $\le p_2$ in the axial coordinate. One or two degrees.

All rules have positive weights and interior abscissae. Triangle rules of degree 2, 4 and 5 and the tetrahedron rule of degree 2 are the tabulated symmetric rules; all other simplex rules are collapsed (conical product) Gauss–Jacobi rules, available for any degree.

Point ordering

Tensor-product rules vary the first coordinate fastest, and the Gauss points of each direction run from +1 down to -1. Wedge rules vary the triangle point fastest and the axial point slowest.

Conventions

Abscissae and weights are real64; kinds are not exported, so take real64 from iso_fortran_env. Invalid arguments to the constructor or set stop with error stop and a message that starts with cubature%set:.

References

1. Abramowitz, M., Stegun, I. A. (1964). Handbook of Mathematical Functions, §22 and §25.4. National Bureau of Standards.
2. Stroud, A. H. (1971). Approximate Calculation of Multiple Integrals. Prentice-Hall.
3. Radon, J. (1948). Zur mechanischen Kubatur. Monatshefte für Mathematik, 52, 286–300.
4. Dunavant, D. A. (1985). High degree efficient symmetrical Gaussian quadrature rules for the triangle. International Journal for Numerical Methods in Engineering, 21(6), 1129–1148.
5. Karniadakis, G. E., Sherwin, S. J. (2005). Spectral/hp Element Methods for Computational Fluid Dynamics, 2nd ed. Oxford University Press.
6. Zienkiewicz, O. C., Taylor, R. L., Zhu, J. Z. (2005). The Finite Element Method: Its Basis and Fundamentals, 6th ed. Butterworth-Heinemann.

Variables

Type	Visibility	Attributes		Name		Initial
integer,	public,	parameter	::	CUB_LIN	=	1	Line element
integer,	public,	parameter	::	CUB_TRI	=	2	Triangle element
integer,	public,	parameter	::	CUB_QUA	=	3	Quadrilateral element
integer,	public,	parameter	::	CUB_TET	=	4	Tetrahedron element
integer,	public,	parameter	::	CUB_HEX	=	5	Hexahedron element
integer,	public,	parameter	::	CUB_WED	=	6	Wedge (prism) element

Interfaces

public interface cubature

Construct a rule, e.g. cubature(CUB_TRI, 4) or cubature(CUB_HEX, [3, 3, 1])

• private pure function new_iso(elm, degree) result(self)

  Construct a rule of degree degree in every direction

  Arguments

  Type	Intent	Optional	Attributes		Name
  integer,	intent(in)			::	elm	Element type, CUB_*
  integer,	intent(in)			::	degree	Degree of exactness

  Return Value type(cubature)

• private pure function new_aniso(elm, degree) result(self)

  Construct a rule with per-direction degrees

  Arguments

  Type	Intent	Optional	Attributes		Name
  integer,	intent(in)			::	elm	Element type, CUB_*
  integer,	intent(in)			::	degree(:)	Degrees of exactness, size 1 or number of directions

  Return Value type(cubature)

Derived Types

type, public ::  cubature

Abscissae and weights of an integration rule on a reference element.

Read more…

Components

Type	Visibility	Attributes		Name		Initial
integer,	public		::	elm	=	0	Element type, CUB_*
integer,	public		::	dim	=	0	Spatial dimension
integer,	public		::	degree(3)	=	0	Degree of exactness, see module documentation
integer,	public		::	npoints	=	0	Number of points
real(kind=rk),	public,	allocatable	::	abscissae(:,:)			Abscissae, shape [dim, npoints]
real(kind=rk),	public,	allocatable	::	weights(:)			Weights, shape [npoints]

Constructor

Construct a rule, e.g. cubature(CUB_TRI, 4) or cubature(CUB_HEX, [3, 3, 1])

private pure function new_iso (elm, degree)	Construct a rule of degree degree in every direction
private pure function new_aniso (elm, degree)	Construct a rule with per-direction degrees

Type-Bound Procedures

generic, public :: set => set_iso, set_aniso	Build the rule in place
procedure, public :: is_valid	True if the rule is set and consistent
procedure, public :: summary	Write element, dimension, degrees, points
procedure, public :: show	Write the summary, points and weights
procedure, public :: destroy	Deallocate and reset