# Sound-hard Burton miller acoustic scattering

**URL:** <https://bempp.discourse.group/t/sound-hard-burton-miller-acoustic-scattering/45>\
**Category:** Applications\
**Created:** [23 November 2020 10:30 UTC](https://bempp.discourse.group/t/sound-hard-burton-miller-acoustic-scattering/45 "2020-11-23T10:30:08Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![Rasmus](https://avatars.discourse-cdn.com/v4/letter/r/7993a0/32.png) [@Rasmus](https://bempp.discourse.group/u/Rasmus)\
**Post date:** [23 November 2020 10:30 UTC](https://bempp.discourse.group/t/sound-hard-burton-miller-acoustic-scattering/45/1 "2020-11-23T10:30:08Z")

</div>

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:

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

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

---

<div class="post-metadata">

**Author:** ![Elwin](https://yyz2.discourse-cdn.com/free1/user_avatar/bempp.discourse.group/elwin/32/38_2.png) [@Elwin](https://bempp.discourse.group/u/Elwin)\
**Post date:** [23 November 2020 12:16 UTC](https://bempp.discourse.group/t/sound-hard-burton-miller-acoustic-scattering/45/2 "2020-11-23T12:16:18Z")

</div>

Rasmus,

Please try: space = bempp.api.function\_space(grid, “P”, 1).  
The hypersingular operator requires piecewise linear functions. These are the continuous versions instead of the discontinuous versions you used.

Best,  
Elwin

---

<div class="post-metadata">

**Author:** ![Rasmus](https://avatars.discourse-cdn.com/v4/letter/r/7993a0/32.png) [@Rasmus](https://bempp.discourse.group/u/Rasmus)\
**Post date:** [23 November 2020 13:29 UTC](https://bempp.discourse.group/t/sound-hard-burton-miller-acoustic-scattering/45/3 "2020-11-23T13:29:25Z")

</div>

Thanks a ton, Elwin!

Now at least it seems to give me results in a decent amount of time.

In this way I seem to get a result for every vertex. If i want to approximate the surface solution on the midpoint of each element, can i simply take the arithmetic mean of the pressure at the 3 vertices of the element?

Best regards,  
Rasmus

EDIT: Also to the moderators, if i run into other problems using BEMPP in this project, would you prefer if i ask for help in this thread or make a new one for separate issues?

---

<div class="post-metadata">

**Author:** ![Elwin](https://yyz2.discourse-cdn.com/free1/user_avatar/bempp.discourse.group/elwin/32/38_2.png) [@Elwin](https://bempp.discourse.group/u/Elwin)\
**Post date:** [23 November 2020 17:28 UTC](https://bempp.discourse.group/t/sound-hard-burton-miller-acoustic-scattering/45/4 "2020-11-23T17:28:22Z")

</div>

Rasmus,

The continuous P1 elements are node based: value 1 at a given vertex in the mesh, value 0 at all other vertices and piecewise linear at each triangle element. Hence, the solution you get are values of the surface pressure at each vertex. You can evaluate the surface pressure inside the triangles with linear interpolation of the values at the three vertices.

Best,  
Elwin

---

<div class="post-metadata">

**Author:** ![LCheng](https://avatars.discourse-cdn.com/v4/letter/l/9fc348/32.png) [@LCheng](https://bempp.discourse.group/u/LCheng)\
**Post date:** [1 July 2021 06:25 UTC](https://bempp.discourse.group/t/sound-hard-burton-miller-acoustic-scattering/45/5 "2021-07-01T06:25:37Z")

</div>

Hi, Rasmus, have you compared this solution with the analytical one? It seems that large error is caused by utilizing the BIE `0.5p - M[p] = p_{inc}` for this wavenumber. The following figure is the calculation for ka = 5 (a - radius of the sphere).

![image](https://global.discourse-cdn.com/free1/uploads/bempp/original/1X/e5ac9f56c6ac0226b20a8b84bc262e629e3c5bc3.png)

(DoF - 8192, dp0\_space = bempp.api.function\_space(grid, “DP”, 0) )

---

<div class="post-metadata">

**Author:** ![Gustavo](https://yyz2.discourse-cdn.com/free1/user_avatar/bempp.discourse.group/gustavo/32/36_2.png) [@Gustavo](https://bempp.discourse.group/u/Gustavo)\
**Post date:** [10 November 2021 11:59 UTC](https://bempp.discourse.group/t/sound-hard-burton-miller-acoustic-scattering/45/6 "2021-11-10T11:59:13Z")

</div>

Sorry, I am pretty new to BEMPP. How did you generated this graph?

Thanks in advance,  
Gustavo

---

<div class="post-metadata">

**Author:** ![LCheng](https://avatars.discourse-cdn.com/v4/letter/l/9fc348/32.png) [@LCheng](https://bempp.discourse.group/u/LCheng)\
**Post date:** [10 November 2021 12:18 UTC](https://bempp.discourse.group/t/sound-hard-burton-miller-acoustic-scattering/45/7 "2021-11-10T12:18:44Z")

</div>

Hi, [Jupyter Notebook Viewer (nbviewer.org)](https://nbviewer.org/github/bempp/bempp-cl/blob/master/notebooks/helmholtz/helmholtz_combined_exterior.ipynb) is a good reference. The sound pressure directivity in the upper figure is ploted using the azimuthal angle and pressure magnitude on the probe points. Actually, here I exported the values and then using the ‘polar’ function in Octave. Definitely, you can also do this in Python.

Best wishes,  
Long