Volumetric Integration Scheme for PDEs

Updated 26 May 2026

Volumetric integration scheme is a numerical technique for computing integrals over complex geometries in PDEs and remapping applications.
It employs flexible domain decomposition with mapped reference cells and high-order singular quadrature to handle challenging geometrical features.
Integration with fast multipole methods and local orthogonal polynomial expansions ensures O(N) performance and strict conservation of physical quantities.

A volumetric integration scheme refers to a numerical technique for computing integrals of functions over domains with potentially complex geometry, typically in the context of partial differential equations (PDEs), boundary value problems, or mesh-based discretizations. In modern applied mathematics and computational physics, volumetric integration is essential for solving inhomogeneous PDEs, evaluating fundamental solutions (volume potentials), and ensuring accurate conservation of mass, momentum, or energy in mesh remapping or mesh-to-voxel projection.

1. Volume Potential Formulation and Applications

The core objective is to compute integrals of the form

$\mathcal V[f](\mathbf r_0) = \int_\Omega G(\mathbf r, \mathbf r_0)\, f(\mathbf r) \, dA(\mathbf r), \qquad \mathbf r_0\in\Omega,$

where $\Omega \subset \mathbb{R}^2$ is a (possibly multiply-connected) domain, $f$ is a source density (solution to $L u = f$ ), and $G$ is the free-space Green’s function associated with the constant-coefficient elliptic operator $L$ (Anderson et al., 2022). This framework underpins the solution of inhomogeneous elliptic PDEs via decomposition into a particular solution (volume potential) plus a homogeneous (boundary-driven) term.

These volume potentials also appear in boundary integral methods for the Stokes and Navier–Stokes equations, as well as in the construction of exact conservative remeshings for mesh-to-voxel projection in computational physics (Gray et al., 2018, Powell et al., 2014).

2. Domain Decomposition and Mappings

Efficient volumetric integration in complex domains requires flexible domain decomposition. The standard approach is to cover $\Omega$ by a union of (i) regular “box-code” grid cells (squares), and (ii) a narrow strip of boundary-fitted triangular cells that may be straight-edged or curved. Each cell $\mathcal{C}_k$ is mapped from a reference shape ( $\widehat{T}_0$ for triangles, $\widehat{B}_0$ for squares) using a $\Omega \subset \mathbb{R}^2$ 0-diffeomorphism $\Omega \subset \mathbb{R}^2$ 1.

Affine mappings are used for straight triangles.
Gordon–Hall blending is used for curved triangles, with parametric representation of curved edges.
Regular boxes employ a simple linear scaling and shift map.

This decomposition aligns with standard meshing tools such as Gmsh (Anderson et al., 2022).

3. Area-to-Boundary Conversion via the Poincaré Lemma

A central innovation for high-order volumetric integration is the reduction of area integrals to boundary integrals through the Poincaré lemma. For any star-shaped region $\Omega \subset \mathbb{R}^2$ 2 and kernel $\Omega \subset \mathbb{R}^2$ 3, the formula

$\Omega \subset \mathbb{R}^2$ 4

replaces a 2D singular volume potential by a contour integral around the cell boundary, each coupled to a 1D quadrature in the interpolation parameter $\Omega \subset \mathbb{R}^2$ 5 (Anderson et al., 2022). The "star-point" $\Omega \subset \mathbb{R}^2$ 6 is chosen as the cell point closest to the target $\Omega \subset \mathbb{R}^2$ 7.

4. Singular and Near-Singular Quadrature

Despite the conversion to boundary integrals, the $\Omega \subset \mathbb{R}^2$ 8-integral can still be weakly singular. High-order accuracy is achieved by employing:

Alpert's endpoint-corrected trapezoid rules for logarithmic or weak algebraic singularities when the target is inside the cell.
Modified Gaussian quadrature (Kolm–Rokhlin) to resolve near-singular behavior for targets close to, but outside, the cell.

Composite Gauss–Legendre quadrature with interval refinement, especially near star-points on the boundary, supports uniform accuracy even when targets approach cell edges or vertices (Anderson et al., 2022).

5. Local Orthogonal-Polynomial Expansions and Fast Multipole Coupling

Within each cell, the source term $\Omega \subset \mathbb{R}^2$ 9 is approximated by a high-order orthogonal polynomial expansion:

On triangles, Koornwinder polynomials $f$ 0 are employed and coefficients are computed using well-conditioned Vioreanu–Rokhlin nodes.
On boxes, tensor-product Chebyshev polynomials and nodes are used.

This expansion enables separation of smooth, far-field interactions (handled by global quadrature and fast multipole methods—FMM) from singular or near-singular local interactions (handled by precomputed local corrections). For each target,

$f$ 1

where “far” cell contributions are rapidly evaluated via FMM and smooth quadrature, while “near”/“self” cells invoke local high-order singular corrections. The composite approach ensures $f$ 2 scaling with the number of degrees of freedom, with nearly all computational time concentrated in the FMM (Anderson et al., 2022).

6. Implementation Details and Numerical Performance

The scheme is compatible with adaptive Cartesian box codes and common mesh generators. All local quadrature corrections (on reference triangles or boxes) can be precomputed and reused across the mesh. Empirical results demonstrate:

8th-order spatial convergence for the maximum error as mesh is refined.
Achievable $f$ 3– $f$ 4 absolute errors even in nearly touching complex geometries.
For static geometries, over $f$ 5 of computational effort is in the FMM stage; singular and near-singular corrections are negligible in cost.

The algorithm handles both box and triangle boundary cells, achieves high-order accuracy, and maintains strict moment conservation in remapping contexts (Anderson et al., 2022, Powell et al., 2014).

The high-order volumetric integration method based on area-to-line conversion and orthogonal-polynomial cell expansions is particularly distinguished by:

Compatibility with fast multipole methods for $f$ 6 complexity.
Uniform applicability on union-of-triangles + boxes domains, with exact adherence to adaptive, possibly curved, domain boundaries.
Integration with standard boundary integral equation solvers for both homogeneous and inhomogeneous elliptic PDEs.

Previous approaches, including boundary plus remainder formulations for the Stokes equations (Gray et al., 2018) and physically conservative remeshing via analytic polyhedron intersection (Powell et al., 2014), emphasize exact moment preservation and accurate flux transfer. As such, these schemes form a toolkit for high-order, physically robust volumetric integration across disparate computational physics applications. Future research directions may extend these methods to 3D, dynamically adaptive meshes, or more general kernels beyond classical Green’s functions.