Fast integration in 3D boundary element methods

The boundary integral method is a common component in solution strategies for elliptic (and sometimes parabolic) PDEs, especially when classical PDEs have to be solved on domains with complicated and/or infinite boundaries. For more general nonlinear problems the boundary integral method can often be coupled with standard finite elements as part of an efficient composite solution strategy. The main advantage of the boundary integral equation approach is that one has only to perform discretisation over the 2D boundary of a 3D domain, whereas the main disadvantage is that the resulting stiffness matrix is full and contains (singular) surface integrals. The latter fact poses problems both for matrix setup and system solution. In recent years much progress has been made on the derivation of fast iterative methods for solving the system but less has been done on the computation of the matrix itself. In many modern codes the formation of the stiffness matrix can form a dominant part of the computation time. For N degrees of freedom the cost of assembling the complete matrix is CN² with C moderately large and depending on the number of kernel evaluations used by the integration routine.

In this work we describe a new fast quadrature method for approximating the stiffness matrix which reduces this cost. For example, in the case of continuous piecewise linear basis functions on a mesh with N nodes, the new method computes the stiffness matrix using N² + O(N) kernel evaluations. The N² term comes from the fact that most of the matrix is computed using novel quadrature rules which just use kernel evaluations at pairs of nodes, with the lower order term coming from conventional rules applied to the small remaining part of the matrix. The computed (approximate) stiffness matrix is guaranteed to be sufficiently accurate so as to preserve the known theoretical convergence properties of the numerical method.

In this talk we will describe the essential theoretical ideas behind this method allowing us to prove a number of theorems governing its stability, convergence and complexity. The analytical tools required include novel inverse estimates for finite element functions on general mesh sequences and the theory of polynomial interpolation at arbitrary points in the plane.

We also describe the practical implementation of the method and illustrate this with examples involving several boundary integral equations on a selection of smooth and non-smooth surfaces. The numerical experiments show that the theoretical predictions are realised even in the case of the strongly anisotropic meshes commonly used for handling edge singularities in 3D.