Frugal numerical integration scheme for polytopal domains

This paper introduces a novel numerical integration scheme tailored for polytopic domains, circumventing the need for sub-tessellation or sub-tetrahedralization. Our method involves defining integration points on a Cartesian bounding box surrounding the polytopic domain and computing integration weights through moment matching with analytically computed integrals of monomials using Euler’s homogeneous function theorem. The fact that points are defined across the bounding box renders the scheme particularly suited for methods where the variable of interest is defined on the bounding box, i.e. the polytopal version of Interior Penalty Discontinuous Galerkin Method. We demonstrate the method’s typical accuracy, achieving an error of \(\sim 10^{-15}\) $\sim {10}^{-15}$ for simple integrands and a maximum error of \(5\times 10^{-14}\) $5\times {10}^{-14}$ for the computation of mass-like matrices over complex non-convex polyhedra. Additionally, we highlight the method’s efficiency, as it requires only a small number of integration points, precisely matching the dimension of the polynomial space with a total degree equal to that of the integrand. The key advantage of the method is that since the distribution of points is defined on the bounding box it is shape independent. As a result the most expensive operation (QR decomposition) can be done once for all. Then for each polytopic shape one simply needs to compute the analytical integral of monomials and solve a small matrix equation. Alongside this paper comes a python implementation of the proposed method freely available at github.com/LMSD-KULeuven/polyquad.