Dimensionality of input and output arrays

quadpy is fully vectorized and it gives the user the opportunity to evaluate the function at many points at once. This dramatically increases the speed at which integrals are computed.

The input x to the user function f is typically of shape (d, n, p) where the leading dimension d is the geometrical dimensionality of the domain (e.g. 2 for domains in 2D space), n is the number of domains over which you want to integrate (e.g., 250 triangles at once), and p is the number of integration points of the scheme. The output of f must be of shape (..., n, p), e.g., (n, p) for scalar-valued functions. You can realize that by simply looping, e.g.,

def f(x):
    out = numpy.empty(x.shape[1:])
    for i in range(x.shape[1]):
        for j in range(x.shape[2]):
            out[i, j] = numpy.sin(x[0, i, j]) * numpy.sin(x[1, i, j])
    return out

but the vectorized version

def f(x):
    return numpy.sin(x[0]) * numpy.sin(x[1])

is easier to read and much faster.

If you return an output of shape (a, n, p) ("vector"), (a, b, n, p) ("matrix") or any length before (n, p), quadpy returns an output of shape (a,), (a, b), etc. as one would expect.