- 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

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.

```
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

```
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.

```
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()

- 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

```
```

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

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