The Boundary Element Method

Model problem

The model problem we consider here is the scattering of an incoming acoustic wave \(u_\text{inc}\) by an obstacle \(\Omega\). Thus, we want to solve the following homogeneous Helmholtz equation written in terms of the scattered field \(u\):

with the Sommerfeld radiation condition at infinity, which states that there can be only outgoing waves at infinity:

and where the total field \(u_\text{tot} = u_\text{inc} + u\).

If the wavenumber \(k\) is constant in \(\mathbb{R}^3 \backslash \Omega\), the boundary element method can be applied. It consists in reformulating the problem in terms of unknowns on the boundary \(\Gamma\) of \(\Omega\).

First, let us introduce the Green kernel \(\mathcal{G}\), which for the helmholtz equation in 3D is

Let us also introduce the Single Layer Potential \(\operatorname{SL}\), which for any \(q \in H^{-1/2}(\Gamma)\) is defined as

An interesting property of \(\text{SL}\) is that it produces solutions of the PDE at hand in \(\mathbb{R}^3 \backslash \Gamma\) which satisfy the necessary conditions at infinity (here the Helmholtz equation and the Sommerfeld radiation condition).

Thus, we now need to find a so-called ansatz \(p \in H^{-1/2}(\Gamma)\) such that \(\forall \boldsymbol{x} \in \mathbb{R}^3 \backslash \Omega\)

where \(u\) also verifies the Dirichlet boundary condition \(u = - u_\text{inc}\) on \(\Gamma\).

In order to find \(p\), we define a variational problem by multiplying (42) by a test function q and integrating over \(\Gamma\):

Using the Dirichlet boundary condition \(u = - u_\text{inc}\) on \(\Gamma\), we end up with the following variational problem to solve: find \(p : \Gamma \rightarrow \mathbb{C}\) such that

Note that knowing \(p\) on \(\Gamma\), we can indeed compute \(u\) anywhere using the potential formulation (42). Thus, we essentially gained one space dimension, as we only have to solve for \(p : \Gamma \rightarrow \mathbb{C}\) in (43). Another advantage of the boundary element method is that for a given mesh size, it is usually more accurate than the finite element method.

Of course, these benefits of the boundary element method come with a drawback: after discretization of (43), for example with piecewise linear continuous (P1) functions on \(\Gamma\), we end up with a linear system whose matrix is full: because \(\mathcal{G}(\boldsymbol{x}-\boldsymbol{y})\) never vanishes, every interaction coefficient is nonzero. Thus, the matrix \(A\) of the linear system can be very costly to store (\(N^2\) coefficients) and invert (factorization in \(\mathcal{O}(N^3)\)) (\(N\) is the size of the linear system). Moreover, compared to the finite element method, the matrix coefficients are much more expensive to compute because of the double integral and the evaluation of the Green function \(\mathcal{G}\). Plus, the choice of the quadrature formulas has to be made with extra care because of the singularity of \(\mathcal{G}\).

Boundary Integral Operators

In order to formulate our model Dirichlet problem, we have used the Single Layer Potential \(\operatorname{SL}\):

Depending on the choice of the boundary integral formulation or boundary condition, the Double Layer Potential \(\operatorname{DL}\) can also be of use:

Similarly, we have used the Single Layer Operator \(\mathcal{SL}\) in our variational problem

There are three other building blocks that can be of use in the boundary element method, and depending on the problem and the choice of the formulation a boundary integral method makes use of one or a combination of these building blocks:

the Double Layer Operator \(\mathcal{DL}\):

the Transpose Double Layer Operator \(\mathcal{TDL}\):

the Hypersingular Operator \(\mathcal{HS}\):

the BEMTool library

In order to compute the coefficients of the BEM matrix, FreeFEM is interfaced with the boundary element library BEMTool. BEMTool is a general purpose header-only C++ library written by Xavier Claeys, which handles

• BEM Potentials and Operators for Laplace, Yukawa, Helmholtz and Maxwell equations

• both in 2D and in 3D

• 1D, 2D and 3D triangulations

• \(\mathbb{P}_k\)-Lagrange for \(k = 0,1,2\) and surface \(\mathbb{RT}_0\)

Hierarchical matrices

Although BEMTool can compute the BEM matrix coefficients by accurately and efficiently evaluating the boundary integral operator, it is very costly and often prohibitive to compute and store all \(N^2\) coefficients of the matrix. Thus, we have to rely on a matrix compression technique. To do so, FreeFEM relies on the Hierarchical Matrix, or H-Matrix format.

Hierarchical block structure

Unfortunately, BEM matrices generally do not have fast decreasing singular values. However, they can exhibit sub-blocks with rapidly decreasing singular values, thanks to the asymptotically smooth nature of the BEM kernel. Let us look for example at the absolute value of the matrix coefficients in the 2D (circle) case below:

• blocks near the diagonal contain information about the near-field interactions, which are not low-rank in nature

• blocks away from the diagonal corresponding to the interaction between two clusters of geometric points \(X\) and \(Y\) satisfying the so-called admissibility condition

(44)\[\max(\text{diam}(X),\text{diam}(Y)) \leq \eta \text{ dist}(X,Y)\]

are far-field interactions and have exponentially decreasing singular values. Thus, they can be well-approximated by low-rank matrices.

The idea is then to build a hierarchical representation of the blocks of the matrix, then identify and compress admissible blocks using low-rank approximation.

We can then build the H-Matrix by taking the following steps:

1. build a hierarchical partition of the geometry, leading to a cluster tree of the unknowns. It can for example be defined using bisection and principal component analysis.

2. from this hierarchical clustering, define and traverse the block cluster tree representation of the matrix structure, identifying the compressible blocks using admissibility condition (44)

3. compute the low-rank approximation of the identified compressible blocks using e.g. partial ACA ; the remaining leaves corresponding to near-field interactions are computed as dense blocks.

The Htool library

the H-Matrix format is implemented in the C++ library Htool. Htool is a parallel header-only library written by Pierre Marchand and Pierre-Henri Tournier. It is interfaced with FreeFEM and provides routines to build hierarchical matrix structures (cluster trees, block trees, low-rank matrices, block matrices) as well as efficient parallel matrix-vector and matrix-matrix product using MPI and OpenMP. Htool is interfaced with BemTool to allow the compression of BEM matrices using the H-Matrix format in FreeFEM.

Build the geometry

The geometry of the problem (i.e. the boundary \(\Gamma\)) can be discretized by a line (2D) or surface (3D) mesh:

1border b(t = 0, 2*pi){x=cos(t); y=sin(t);}
2meshL ThL = buildmeshL(b(100));

1load "msh3"
2mesh Th = square(10,10);
3meshL ThL = extract(Th);

1load "msh3"
2meshS ThS = square3(10,10);
3int[int] labs = [1,2];
4meshL ThL = extract(ThS, label=labs);

3D

In 3D, the geometry of the boundary surface can be discretized with a surface mesh of type meshS, which can be built by several ways, for example using the square3 constructor:

1load "msh3"
2real R = 3, r=1, h=0.2;
3int nx = R*2*pi/h, ny = r*2*pi/h;
4func torex = (R+r*cos(y*pi*2))*cos(x*pi*2);
5func tore y= (R+r*cos(y*pi*2))*sin(x*pi*2);
6func torez = r*sin(y*pi*2);
7meshS ThS = square3(nx,ny,[torex,torey,torez],removeduplicate=true);

or from a 2D mesh using the movemesh23 keyword:

1load "msh3"
2mesh Th = square(10,10);
3meshS ThS = movemesh23(Th, transfo=[x,y,cos(x)^2+sin(y)^2]);

We can also extract the boundary of a mesh3:

1load "msh3"
2mesh3 Th3 = cube(10,10,10);
3int[int] labs = [1,2,3,4];
4meshS ThS = extract(Th3, label=labs);

You can find much more information about surface mesh generation here.

Orientation of normal vector

Depending on whether your problem is posed on a bounded or unbounded domain, you may have to set the orientation of the outward normal vector \(\boldsymbol{n}\) to the boundary. You can use the OrientNormal function with the parameter unbounded set to 0 or 1 (the normal vector \(\boldsymbol{n}\) will then point to the exterior of the domain you are interested in):

1border b(t = 0, 2*pi){x=cos(t); y=sin(t);}
2meshL ThL = buildmeshL(b(100));
3ThL = OrientNormal(ThL,unbounded=1);
4plot(ThL,dim=2);

You can use shift + t on a plot of a boundary mesh to display the outward normal vector \(\boldsymbol{n}\):

Define the type of operator

For now, FreeFEM allows to solve the following PDE with the boundary element method:

with

• \(k = 0\) (Laplace)

• \(k \in \mathbb{R}^*_+\) (Helmholtz)

• \(k \in \imath \mathbb{R}^*_+\) (Yukawa)

First, the BEM plugin needs to be loaded:

1load "bem"

The information about the type of operator and the PDE can be specified by defining a variable of type BemKernel:

1BemKernel Ker("SL",k=2*pi);

You can choose the type of operator depending on your formulation (see Boundary Integral Operators):

• "SL": Single Layer Operator \(\mathcal{SL}\)

• "DL": Double Layer Operator \(\mathcal{DL}\)

• "TDL": Transpose Double Layer Operator \(\mathcal{TDL}\)

• "HS": Hyper Singular Operator \(\mathcal{HS}\)

Define the variational problem

We can then define the variational form of the BEM problem. The double BEM integral is represented by the int1dx1d keyword in the 2D case, and by int2dx2d for a 3D problem. The BEM keyword inside the integral takes the BEM kernel operator as argument:

1BemKernel Ker("SL", k=2*pi);
2varf vbem(u,v) = int2dx2d(ThS)(ThS)(BEM(Ker,u,v));

You can also specify the BEM kernel directly inside the integral:

1varf vbem(u,v) = int2dx2d(ThS)(ThS)(BEM(BemKernel("SL",k=2*pi),u,v));

Depending on the choice of the BEM formulation, there can be additional terms in the variational form. For example, Second kind formulations have an additional mass term:

1BemKernel Ker("HS", k=2*pi);
2varf vbem(u,v) = int2dx2d(ThS)(ThS)(BEM(Ker,u,v)) - int2d(ThS)(0.5*u*v);

We can also define a linear combination of two BEM kernels, which is useful for Combined formulations:

1complex k=2*pi;
2BemKernel Ker1("HS", k=k);
3BemKernel Ker2("DL", k=k);
4BemKernel Ker = 1./(1i*k) * Ker1 + Ker2;
5varf vbem(u,v) = int2dx2d(ThS)(ThS)(BEM(Ker,u,v)) - int2d(ThS)(0.5*u*v);

As a starting point, you can find how to solve a 2D scattering problem by a disk using a First kind, Second kind and Combined formulation, for a Dirichlet (here) and Neumann (here) boundary condition.

Assemble the H-Matrix

Assembling the matrix corresponding to the discretization of the variational form on an fespace Uh is similar to the finite element case, except that we end up with an HMatrix instead of a sparse matrix:

1fespace Uh(ThS,P1);
2HMatrix<complex> H = vbem(Uh,Uh);

Behind the scenes, FreeFEM is using Htool and BEMTool to assemble the H-Matrix.

1ff-mpirun -np 4 script.edp -wg

You can specify the different Htool parameters as below. These are the default values:

1HMatrix<complex> H = vbem(Uh,Uh,
2  compressor = "partialACA", // or "fullACA", "SVD"
3  eta = 10.,                 // parameter for the admissibility condition
4  eps = 1e-3,                // target compression error for each block
5  minclustersize = 10,       // minimum block side size min(n,m)
6  maxblocksize = 1000000,    // maximum n*m block size
7  commworld = mpiCommWorld); // MPI communicator

You can also set the default parameters globally in the script by changing the value of the global variables htoolEta, htoolEpsilon, htoolMinclustersize and htoolMaxblocksize.

Once assembled, the H-Matrix can also be plotted with

1display(H, wait=true);

FreeFEM can also output some information and statistics about the assembly of H:

1if (mpirank == 0) cout << H.infos << endl;

Solve the linear system

Generally, the right-hand-side of the linear system is built as the discretization of a standard linear form:

1Uh<complex> p, b;
2varf vrhs(u,v) = -int2d(ThS)(uinc*v);
3b[] = vrhs(0,Uh);

We can then solve the linear system to obtain \(p\), with the standard syntax:

1p[] = H^-1*b[];

Under the hood, FreeFEM solves the linear system with GMRES with a Jacobi (diagonal) preconditioner.

Compute the solution

Finally, knowing \(p\), we can compute the solution \(u\) of our initial problem (39) using the Potential as in (42). As for the BemKernel, the information about the type of potential can be specified by defining a variable of type BemPotential:

1BemPotential Pot("SL", k=2*pi);

In order to benefit from low-rank compression, instead of using (42) to sequentially compute the value \(u(\boldsymbol{x})\) at each point of interest \(\boldsymbol{x}\), we can compute the discretization of the Potential on a target finite element space UhOut defined on an output mesh ThOut with an H-Matrix.

First, let us define the variational form corresponding to the potential that we want to use to reconstruct our solution. Similarly to the kernel case, the POT keyword takes the potential as argument. Note that we have a single integral, and that v plays the role of \(\boldsymbol{x}\).

1varf vpot(u,v) = int2d(ThS)(POT(Pot,u,v));

We can then assemble the rectangular H-Matrix from the potential variational form:

1fespace UhOut(ThOut,P1);
2HMatrix<complex> HP = vpot(Uh,UhOut);

Computing \(u\) on UhOut is then just a matter of performing the matrix-vector product of HP with p:

1UhOut<complex> u;
2u[] = HP*p[];
3plot(u);