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 set-up

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

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

+ O(N)

kernel evaluations.

The

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.