BEM++ overview

https://bempp.com/

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

Installation

https://bitbucket.org/bemppsolutions/bempp

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

Boost

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 git@bitbucket.org:bemppsolutions/bempp.git

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 setup.py 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

Grids



In [1]:

    
import bempp.api
import numpy as np
grid = bempp.api.shapes.regular_sphere(3)
grid.plot()

Q: what means argument of the function?

Plot the following objects

cube
ellipsoid
rectangle with hole

Import/export

You can import your mesh in Gmsh file
```
bempp.api.import_grid(filename)
```

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)
grid_fun.plot()

Operators

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,
                                                        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
```
bempp.api.linalg
```
has solvers that require operator and grid function instead of matrix and vector.
Automatic LU preconditioner is available

Example

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