Gaussian Quadrature for an n-dimensional simplex

• GitHub: Source code repository
• 

This code was originally published by Greg von Winckel (Contact: gregvw(at)math(dot)unm(dot)edu, http://math.unm.edu/~gregvw) on the MathWorks File Exchange https://www.mathworks.com/matlabcentral/fileexchange/9435-n-dimensional-simplex-quadrature. (The given email and web addresses seems now defunct; however, a web search for Greg von Winckel easily finds up-to-date pointers.) The code in this package is a fairly literal translation from Matlab to Julia.

Construct Gauss points and weights for an n-dimensional simplex domain with vertices specified by the n*(n-1) matrix vert, where each row contains the coordinates (x1,...,xn) for a vertex. The order of the quadrature scheme in each dimension must be the same in this implementation.

X, W = simplexquad(n, vert)     # Specify the vertices
X, W = simplexquad(T, n, dim)   # Specify the dimension and use unit simplex

X will be a matrix for which the jth column are the grid points in each coordinate xj.

To integrate a function f, use e.g.

sum(W[i] * f(X[i,:]) for i in 1:length(W))

I tested the package for up to D=5 dimensions and order N=10, and found the integration error for polynomials of order P≤N (which should have only floating-point round-off error) to be less than 10eps. This is tested by the test suite.

Integration points and weights for various numbers of points:

The area of the red disks corresponds to their integration weights. It is evident that the location of the integration points is not a tensor product and does not respect the symmetry of the domain.

The first four simplexes are

Examine JuAFEM.jl and its generate_quadrature.jl file https://kristofferc.github.io/JuAFEM.jl/dev/reference/quadrature/.