BEM++ overview


  • Overview presentation here.
  • Applicable only for 3D problems
  • Support Laplace, Helmholtz and Maxwell equations with Dirichlet and Neumann boundary conditions
  • Support H-matrices and itertive solvers for linear systems

Conda environment preparation

  • Create conda environment
    conda create -n bempp python=3.6
  • Activate this environment
    source activate bempp
  • Install in this environment NumPy stack + MPI packages
    • Numpy
    • Scipy
    • Matplotlib
    • mpi4py
    • Cython
    • Jupyter Notebook kernel
      conda install jupyter notebook

Threading building block

BEM++ uses TBB package, so you should install it before bulding BEM++.

For OS X:

brew install tbb

For Ubuntu:

sudo apt-get install libtbb-dev


Probably you will need install boost library if it had not already installed.

For OS X:

brew install boost

For Ubuntu:

sudo apt-get install libboost-all-dev

Cloning source code and build from source

  • Clone BEM++ source code in your folder with lib source codes

    cd your-folder
    git clone
  • Got to the cloned folder and build C/C++ core of BEM++

    cd ./bempp/
    mkdir build
    cd ./build
    cmake ..
    make -j4

    The last command runs build process and parallel it in 4 threads.

Create and install Python package

After successfully finished building in the previous step go to the home folder

cd ..

and run

python install

This command initiates converting C/C++ core of BEM++ to Python package and installs it to the target place inside your environment.

Check correctness of installation

To check that BEM++ is installed correctly, run python inside environment, where you have installed BEM++, and import it

import bempp.api

Next steps...

Now we are ready to try BEM++ for solving different integral equations, but how to give integral equation for this package?

Main ingredients

Remember first-kind integral equation $$ \int_{\partial \Omega} \frac{q(y)}{\Vert x - y \Vert} dy = f(x), \quad x \in \partial \Omega. $$ and list steps required to solve it numerically.

From IE to LS

  • Discretization
  • Local basis functions
  • Test functions
  • Operator
  • Linear system solver


In [1]:
import bempp.api
import numpy as np
grid = bempp.api.shapes.regular_sphere(3)

Q: what means argument of the function?

Plot the following objects

  • cube
  • ellipsoid
  • rectangle with hole


  • You can import your mesh in Gmsh file
  • You can export mesh from Gmsh file
    bempp.api.export(grid=grid, file_name=filename)

Function space

After setting grid, we can define functions which are used for discretization in this mesh.

In [2]:
space = bempp.api.function_space(grid, "DP", 0)

What mean these arguments?

  • The first argument is always grid object
  • The second argument can be discontinious polynomial ("DP"), polynomial ("P") or some special function space ("DUAL")
  • The third argument is the order of polynomial

Study degrees of freedom for different types and order of function spaces.

Grid function

  • After introducing space we can obtain grid function, which is representation of data on given grid.
  • This object consists of a set of basis function coefficients and a corresponding space object

In [3]:
import numpy as np
def fun(x, normal, domain_index, result):    
    result[0] = np.exp(1j * x[0])
grid_fun = bempp.api.GridFunction(space, fun=fun)


  • Boundary operators $$ A: D \to R $$ is defined by domain space ($D$), range space ($R$) and dual-to-range ($V$). The last space is neccessary for weak reformulation and it will be discussed later.
  • Potential operators map from given space to set of external points

Operators in BEM++

  • All available operators are described here
  • Main groups of operators are operators for Laplace equation, for Maxwell equations and for Helmholtz equation
  • Algebra in operators is implemented, so you can perform the following operatios with operators: sum, multiply by scalar, squared and others...
  • Lazy evaluation: discretization is performed not after defining operators, but only at the moment of solving linear system

In [4]:
slp = bempp.api.operators.boundary.laplace.single_layer(space, space,
scaled_operator = 1.5 * slp
sum_operator = slp + slp
squared_operator = slp * slp
  • General operator $A: D \to R$ mapping has the form $$ Au= f, $$ where $u \in D$, $f \in R$
  • We can inner multiply both side on some function from dual space (here we get dual-to-rage space!) and get $$ \langle Au, v\rangle = \langle f, v\rangle $$
  • In case of integral equation $$ \langle Au, v\rangle = \int_{\Gamma} [Au](y)\bar{v}(y)dy. $$

In [5]:
slp_discrete = slp.weak_form()
print("Shape of the matrix: {0}".format(slp_discrete.shape))
print("Type of the matrix: {0}".format(slp_discrete.dtype))
x = np.random.rand(slp_discrete.shape[1])
y = slp_discrete * x

Shape of the matrix: (512, 512)
Type of the matrix: float64

Linear system solvers

  • This is the final step for solving integral equation
  • But to get it we need to do the following steps:

    • Find weak form of operator:

      discrete_op = A.weak_form()
    • Compute projection of the right-hand side onto the dual-to-range:

      p = f.projections(dual_to_range)
    • Now we have a matrix and a right-hand side, we can call linear system solver from SciPy (with appropriate preconditioner) and get vector of coefficients
    • After that we come back to the domain space and compute solution function
      sol_fun = bempp.api.GridFunction(A.domain, coefficients=x)

Can we do more convenient?

  • Yes! BEM++ has its own wrapper for all above steps
  • Module
    has solvers that require operator and grid function instead of matrix and vector.
  • Automatic LU preconditioner is available


Now let's see how it works for particular problem...