Meshless quadrature formulas arising from numerical differentiation

We suggest a method for simultaneously generating high order quadrature weights for integrals over Lipschitz domains and their boundaries that requires neither meshing nor moment computation. The weights are determined on pre-defined scattered nodes as a minimum norm solution of a sparse underdetermined linear system arising from a discretization of a suitable boundary value problem by either collocation or meshless finite differences. The method is easy to implement independently of the domain's representation, since it only requires as inputs the position of all quadrature nodes and the direction of outward-pointing normals at each node belonging to the boundary. Numerical experiments demonstrate the robustness and high accuracy of the method on a number of smooth and piecewise smooth domains in 2D and 3D, including some with reentrant corners and edges.