Integrating function of surface normals

Applications
1

I am trying to compute various quantities from the scattering amplitude, in the setting of acoustic scattering from an obstacle \Omega with Dirichlet boundary conditions.

One way to compute the scattering amplitude (cf. Taylor 2023, PDE II, Eq. (3.22)) is via the formula

a(\omega, \theta, \lambda)= \frac{i \lambda}{4 \pi} \int_{\partial \Omega}(\nu(y) \cdot \theta) e^{i \lambda(\omega-\theta) \cdot y} d S(y) +\frac{1}{4 \pi} \int_{\partial \Omega} e^{-i \lambda \theta \cdot y} \operatorname{DtN}_\lambda(e^{i\lambda \omega \cdot y}) d S(y) .

Here, \nu(y) is the unit surface normal at the point y \in \partial \Omega. The second term (involving the Dirichlet-to-Neumann operator) I can compute easily.

I am having some trouble with the first term. I can’t work out how to define a grid function that takes the normal vector of a grid as input. I can probably use grid.normals to extract the normal vectors, and then integrate by hand. But I would prefer to use the convenience of grid_function.integrate() for the boundary integral.

Any help would be much appreciated.

2

I think I figured it out.

def normal_fun(omega, space):
    @bempp.api.real_callable
    def normal_fun_callable(x, normal, domain_index, result):
        result[0] = np.dot(normal, omega)
    
    return bempp.api.GridFunction(space, fun=normal_fun_callable)

sphere = bempp.api.shapes.sphere(r=1., origin=(0, 0, 0), h=.2)
sphspace = bempp.api.function_space(sphere, "DP", 0)
omega = np.array([1., 0, 0])
integrand = normal_fun(omega, sphspace)
integrand.integrate()

This produces a value around E-15, which I can easily accept as 0. Seems to work fine.