I am learning to use bempp to solve acoustic problems. I test it with a simple monopole source example:

I use the Fredholm integral equation of the first kind to calculate **Dirichlet data** according to **Neumann data**:

where **Neumann data** g_{\mathrm{N}}:=\gamma^{1, \mathrm{ext}} u and **Dirichlet data** g_{\mathrm{D}}:=\gamma^{0, \mathrm{ext}} u. K_{\kappa} and V_{\kappa} represents **Double Layer Potential Operator** and **Single Layer Potential Operator**.

I expect that the Dirichlet data calculated by the integral is equal to the monopole Dirichlet data, but the result is different. I can’t understand why this happens. Can anyone point out the problem? I have been troubled for a long time.

Thank you!

I installed the latest bempp from github according to the documentation. Here is my test code:

```
import bempp.api
import numpy as np
grid = bempp.api.shapes.sphere(h = 0.1)
dp0_space = bempp.api.function_space(grid, "DP", 0)
p1_space = bempp.api.function_space(grid, "P", 1)
k = 10
@bempp.api.complex_callable()
def monopole_neumann_data(x, n, domain_index, result):
r = (x[0]**2+x[1]**2+x[2]**2)**0.5
result[0] = (-1j*k - 1/r)*(np.exp(-1j*k*r))/(4*np.pi*r)
@bempp.api.complex_callable()
def monopole_dirichlet_data(x, n, domain_index, result):
r = (x[0]**2+x[1]**2+x[2]**2)**0.5
result[0] = (np.exp(-1j*k*r))/(4*np.pi*r)
monopole_neumann_fun = bempp.api.GridFunction(dp0_space, fun=monopole_neumann_data)
monopole_dirichlet_fun = bempp.api.GridFunction(p1_space, fun=monopole_dirichlet_data)
identity = bempp.api.operators.boundary.sparse.identity(
p1_space, p1_space, p1_space)
double_layer = bempp.api.operators.boundary.helmholtz.double_layer(
p1_space, p1_space, p1_space, k)
single_layer = bempp.api.operators.boundary.helmholtz.single_layer(
dp0_space, p1_space, p1_space, k)
left_side = -.5*identity+double_layer
right_side = single_layer*monopole_neumann_fun
output_dirichlet_fun, info = bempp.api.linalg.gmres(left_side, right_side,tol=1e-5)
print(monopole_dirichlet_fun.coefficients[0])
print(output_dirichlet_fun.coefficients[0])
```

I get the result:

```
(-0.06738962350437123+0.042509072890037296j)
(0.05712675186005465-0.05505303182725231j)
```

PS: If I change the formulation to p(x)=\frac{e^{i k r}}{4 \pi r}, I get a result agree with the original data. But this is not reasonable in theory.