Hi

I’m using BEMPP to solve the exterior Helmholtz problem a single point on a surface but for many incident directions.

I would like to, instead of solving the linear system N times for all source position and keeping only the pressure at x_1, i would rather like to place the source at x_1 and evaluate N field points by applying the reciprocity principle.

In the direct case I am computing the BIE \frac{1}{2}p - M[p] + \frac{j}{k} E[p] = p_{inc}

where M[\cdot] is the double layer potential and E[\cdot] is the hypersingular operator, which should correspond to the Burton-Miller formulation for scattering off a surface S with neumann BC \frac{\partial}{\partial \boldsymbol{n}_r}p(\boldsymbol{r}) = 0 , \ \boldsymbol{r} \in S

The code is as follows:

```
k = 1
grid = bempp.api.import_grid('ball5k.stl') # a unit sphere with 5000 triangular elements
space = bempp.api.function_space(grid, "P", 1)
# burton miller sound hard scattering
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")
lhs = 0.5 * I - M - 1j/k * E
# point source
r = 2.5
theta= 0 * np.pi/180
phi = 0 * np.pi/180
r_s = r * np.array([np.cos(theta) * np.cos(phi) , np.sin(theta) * np.cos(phi) , np.sin(phi) ])
@bempp.api.complex_callable
def f(x, n, domain_index,
dist = np.sqrt((r_s[0] - x[0])**2 + (r_s[1] - x[1])**2 + (r_s[2] - x[2])**2)
result[0] = np.exp(1j * k * dist)/(4 * np.pi * dist)
rhs = bempp.api.GridFunction(space, fun=f)
sol, info = bempp.api.linalg.gmres(lhs, rhs,tol=1E-5)
```

I can the evaluate the solution at a point near the boundary

```
r2 = 1 + 1e-3
theta2 = 45 * np.pi/180
phi2 = 20 * np.pi/180
r_r = r2 * np.array([[np.cos(theta2) * np.cos(phi2)],[ np.sin(theta2) * np.cos(phi2)] , [ np.sin(phi2) ]])
dlp_pot = bempp.api.operators.potential.helmholtz.double_layer(space,r_r , k)
dist2 = np.sqrt((r_s[0] - r_r[0])**2 + (r_s[1] - r_r[1])**2 + (r_s[2] - r_r[2])**2)
p = np.exp(1j *k * dist2)/(4 * np.pi * dist2) + dlp_pot.evaluate(sol)
```

However exchanging the points r_s and r_r leads to different results unless r_r is moved far away from the boundary.

Is this due to singularities near the surface? Is there a solution to this problem? Do i have to reformulate my BIE for the reciprical case?

Any help would be really appreciated.

Best Regards

Rasmus

EDIT: A potential solution could be to make an element into the source by making it vibrate with constant velocity. Can anyone help explain, how i would go about adding a separate boundary condition \frac{\partial}{\partial \boldsymbol{n}} p = v_s on a single element of the grid?