Placing monopole source at or near boundary

Applications
1

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?

2

Dear Rasmus, the operators stored under bempp.api.potential are not suited for evaluation on the boundary. They are only highly accurate away from the boundary. The closer you are to the boundary the less accurate they become due to the singularity.

To evaluate potentials at the boundary you need to choose the operators in bempp.api.boundary. For the single layer potential operator as you move to the boundary you can just use the corresponding boundary operator. For the double layer potential operator you need to take the jump relationship into account and add/subtract a half identity term depending from what side of the boundary you are coming.

1 Like
3

Thank you so much, Timo.

I will adjust my code accordingly. So am I correct in assuming that the boundary solution given by sol is actually correct, and the errors are only added when evaluating the field points?

Should it be positive or negative for the outer side of the boundary?