FreeFEM Documentation on GitHub

stars - forks

The Boundary Element Method

Introduction to the Boundary Element Method (BEM)

This is a short overview of the Boundary Element Method, introduced through a simple model example. For a thorough mathematical description of the Boundary Element Method you can refer to the PhD thesis of Pierre Marchand. Some figures used in this documentation are taken from the manuscript.

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\):

(39)\[\begin{split}\left \{ \begin{aligned} - \Delta u - k^2 u &= 0 \;\; &\text{in} \;\; &\mathbb{R}^3 \backslash \Omega \\ u &= - u_\text{inc} \;\; &\text{on} \;\; &\Gamma\\ &\text{+ radiation condition}\hspace{-2.8cm} \end{aligned} \right .\end{split}\]

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

\[\lim_{|\boldsymbol{x}| \rightarrow \infty} |\boldsymbol{x}| \left( \frac{\partial}{\partial |\boldsymbol{x}|} - \imath k \right) u(\boldsymbol{x}) = 0\]

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

(40)\[\mathcal{G}(\boldsymbol{x}) = \exp(\imath k |\boldsymbol{x}|) / (4 \pi |\boldsymbol{x}|).\]

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

(41)\[\operatorname{SL}(q)(\boldsymbol{x}) = \int_{\Gamma} \mathcal{G}(\boldsymbol{x}-\boldsymbol{y}) q(\boldsymbol{y}) d\sigma(\boldsymbol{y}), \quad \forall \boldsymbol{x} \in \mathbb{R}^3 \backslash \Gamma.\]

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\)

(42)\[u(\boldsymbol{x}) = \operatorname{SL}(p)(\boldsymbol{x}) = \int_{\Gamma} \mathcal{G}(\boldsymbol{x}-\boldsymbol{y}) p(\boldsymbol{y}) d\sigma(\boldsymbol{y}),\]

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\):

\[\int_{\Gamma} u(\boldsymbol{x}) q(\boldsymbol{x}) d\sigma(\boldsymbol{x}) = \int_{\Gamma \times \Gamma} \frac{\exp(\imath k |\boldsymbol{x}-\boldsymbol{y}|)}{4 \pi |\boldsymbol{x}-\boldsymbol{y}|} p(\boldsymbol{y}) q(\boldsymbol{x}) d\sigma(\boldsymbol{x,y}) \quad \forall q : \Gamma \rightarrow \mathbb{C}.\]

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

(43)\[\int_{\Gamma \times \Gamma} \frac{\exp(\imath k |\boldsymbol{x}-\boldsymbol{y}|)}{4 \pi |\boldsymbol{x}-\boldsymbol{y}|} p(\boldsymbol{y}) q(\boldsymbol{x}) d\sigma(\boldsymbol{x,y}) = - \int_{\Gamma} u_\text{inc}(\boldsymbol{x}) q(\boldsymbol{x}) d\sigma(\boldsymbol{x}) \quad \forall q : \Gamma \rightarrow \mathbb{C}.\]

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}\):

\[q \mapsto \operatorname{SL}(q)(\boldsymbol{x}) = \int_{\Gamma} \mathcal{G}(\boldsymbol{x}-\boldsymbol{y}) q(\boldsymbol{y}) d\sigma(\boldsymbol{y}).\]

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

\[q \mapsto \operatorname{DL}(q)(\boldsymbol{x}) = \int_{\Gamma} \frac{\partial}{\partial \boldsymbol{n} (\boldsymbol{y})} \mathcal{G}(\boldsymbol{x}-\boldsymbol{y}) q(\boldsymbol{y}) d\sigma(\boldsymbol{y}).\]

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

\[p, q \mapsto \mathcal{SL}(p,q) = \int_{\Gamma \times \Gamma} p(\boldsymbol{x}) q(\boldsymbol{y}) \mathcal{G}(\boldsymbol{x - y}) d \sigma(\boldsymbol{x,y}).\]

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}\):

\[p, q \mapsto \mathcal{DL}(p,q) = \int_{\Gamma \times \Gamma} p(\boldsymbol{x}) q(\boldsymbol{y}) \frac{\partial}{\partial \boldsymbol{n} (\boldsymbol{y})} \mathcal{G}(\boldsymbol{x - y}) d \sigma(\boldsymbol{x,y})\]

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

\[p, q \mapsto \mathcal{TDL}(p,q) = \int_{\Gamma \times \Gamma} p(\boldsymbol{x}) q(\boldsymbol{y}) \frac{\partial}{\partial \boldsymbol{n} (\boldsymbol{x})} \mathcal{G}(\boldsymbol{x - y}) d \sigma(\boldsymbol{x,y})\]

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

\[p, q \mapsto \mathcal{HS}(p,q) = \int_{\Gamma \times \Gamma} p(\boldsymbol{x}) q(\boldsymbol{y}) \frac{\partial}{\partial \boldsymbol{n} (\boldsymbol{x})} \frac{\partial}{\partial \boldsymbol{n} (\boldsymbol{y})} \mathcal{G}(\boldsymbol{x - y}) d \sigma(\boldsymbol{x,y})\]

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.

Low-rank approximation

Let \(\textbf{B} \in \mathbb{C}^{N \times N}\) be a dense matrix. Assume that \(\textbf{B}\) can be written as follows:

\[\textbf{B} = \sum_{j=1}^r \textbf{u}_j \textbf{v}_j^T\]

where \(r \leq N, \textbf{u}_j \in \mathbb{C}^{N}, \textbf{v}_j \in \mathbb{C}^{N}.\)

If \(r < \frac{N^2}{2 N}\), the computing and storage cost is reduced to \(\mathcal{O}(r N) < \mathcal{O}(N^2)\). We say that \(\textbf{B}\) is low rank.

Usually, the matrices we are interested in are not low-rank, but they may be well-approximated by low-rank matrices. We may start by writing their Singular Value Decomposition (SVD):

\[\textbf{B} = \sum_{j=1}^N \sigma_j \textbf{u}_j \textbf{v}_j^T\]

where \((\sigma_j)_{j=1}^N\) are the singular values of \(\textbf{B}\) in decreasing order, and \((\textbf{u}_j)_{j=1}^N\) and \((\textbf{v}_j)_{j=1}^N\) its left and right singular vectors respectively.

Indeed, if \(\textbf{B}\) has fast decreasing singular values \(\sigma_j\), we can obtain a good approximation of \(\textbf{B}\) by truncating the SVD sum, keeping only the first \(r\) terms. Although the truncated SVD is actually the best low-rank approximation possible (Eckart-Young-Mirsky theorem), computing the SVD is costly (\(\mathcal{O}(N^3)\)) and requires computing all \(N^2\) coefficients of the matrix, which we want to avoid.

Thankfully, there exist several techniques to approximate a truncated SVD by computing only some coefficients of the initial matrix, such as randomized SVD, or Partially pivoted Adaptive Cross Approximation (partial ACA), which requires only \(2 r N\) coefficients.

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.


Solve a BEM problem with 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:


In 2D, the geometry of the boundary can be defined with the border keyword and discretized by constructing a line or curve mesh of type meshL using buildmeshL:

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

With the extract keyword, we can also extract the boundary of a 2D mesh (need to load "msh3"):

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

or of a meshS ; we can also specify the boundary labels we want to extract:

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

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


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:

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

or from a 2D mesh using the movemesh23 keyword:

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

We can also extract the boundary of a mesh3:

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

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

Define the type of operator

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

\[-\Delta u - k^2 u = 0, \quad k \in \mathbb{C},\]


  • \(k = 0\) (Laplace)

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

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

First, the BEM plugin needs to be loaded:

1 load "bem"

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

1 BemKernel 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:

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

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

1 varf 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:

1 BemKernel Ker("HS", k=2*pi);
2 varf 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:

1 complex k=2*pi;
2 BemKernel Ker1("HS", k=k);
3 BemKernel Ker2("DL", k=k);
4 BemKernel Ker = 1./(1i*k) * Ker1 + Ker2;
5 varf 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:

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

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


Since Htool is a parallel library, you need to use FreeFem++-mpi or ff-mpirun to be able to run your BEM script. The MPI parallelism is transparent to the user. You can speed up the computation by using multiple cores:

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

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

1 HMatrix<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

1 display(H, wait=true);

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

1 if (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:

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

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

1 p[] = 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:

1 BemPotential 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}\).

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

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

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

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

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

2D example script

Let us summarize what we have learned with a 2D version of our model problem where we study the scattering of a plane wave by a disc:

 1 load "bem"
 3 real k = 10;
 5 int n = 100;
 7 border circle(t = 0, 2*pi){x=cos(t); y=sin(t);}
 8 meshL ThL = buildmeshL(circle(n));
10 varf vbem(u,v) = int1dx1d(ThL)(ThL)(BEM(BemKernel("SL",k=k),u,v));
12 fespace Uh(ThL,P1);
13 HMatrix<complex> H = vbem(Uh,Uh);
15 func uinc = exp(1i*k*x);
16 Uh<complex> p, b;
17 varf vrhs(u,v) = -int1d(ThL)(uinc*v);
18 b[] = vrhs(0,Uh);
20 p[] = H^-1*b[];
22 varf vpot(u,v) = int1d(ThL)(POT(BemPotential("SL",k=k),u,v));
24 int np = 200;
25 int R = 4;
26 border b1(t=-R, R){x=t; y=-R;}
27 border b2(t=-R, R){x=R; y=t;}
28 border b3(t=-R, R){x=-t; y=R;}
29 border b4(t=-R, R){x=-R; y=-t;}
30 mesh ThOut = buildmesh(b1(np)+b2(np)+b3(np)+b4(np)+circle(-n));
32 fespace UhOut(ThOut,P1);
33 HMatrix<complex> HP = vpot(Uh,UhOut);
35 UhOut<complex> u, utot;
36 u[] = HP*p[];
38 utot = u + uinc;
39 plot(utot,fill=1,value=1,cmm="u_total");
Table of content