Hello

First of all, thanks for making this software.

I’m a masters student currently working with BEMPP to solve acoustic scattering problems off of a human head model to numerically calculate HRTF’s.

As I have understood it, i can formulate my BIE as \frac{1}{2}p - M[p] = p_{inc}, with a sound-hard boundary condition \frac{\partial}{\partial n} p = 0.

I get results out using:

```
k = 5 # wave number o
space = bempp.api.function_space(grid, "DP", 0)
I = bempp.api.operators.boundary.sparse.identity(space, space, space)
M = bempp.api.operators.boundary.helmholtz.double_layer(space, space, space,k,assembler = "fmm") # double layer potential
lhs = 0.5 * I - M
@bempp.api.complex_callable
def f(x, n, domain_index, result):
result[0] = np.exp(1j * k * x[0])
rhs = bempp.api.GridFunction(space, fun=f)
sol, info = bempp.api.linalg.gmres(lhs, rhs,tol=1E-4)
```

in a couple of minutes.

However when i try to add the Burton-Miller term:

```
I = bempp.api.operators.boundary.sparse.identity(space, space, space)
M = bempp.api.operators.boundary.helmholtz.double_layer(space, space, space,k,assembler = "fmm") # double layer potential
E = bempp.api.operators.boundary.helmholtz.hypersingular(space, space, space, k, assembler= "fmm") # derivative of dlp
lhs = 0.5 * I - M - 1j/k * E # burton miller
```

The gmres solver hangs forever (or at least becomes really slow), and i am forced to use `space = bempp.api.function_space(grid, "DP", 1)`

which will give me three values per element, even though I’m only interested in a single value for each element.

I’m pretty new to BEM so I’m certainly just doing something wrong, but i’m hoping that the smart people here may have a solution.

Best regards

Rasmus