Ways to approximate singular surface integral in BEMPP

Consider the following integral

I(\vec z) = \int_{\partial\Omega} \frac{e^{0.1(x_1+2x_2+3x_3)}}{|\vec z-\vec x|} \,dL(x), \hspace{3mm}\\
\vec x = (x_1,x_2,x_3) \in \Omega, \hspace{3mm} \vec z = (z_1,z_2,z_3) \in \partial\Omega


\Omega = \left\{(x,y,z) \mid x^2+\left(\frac{y}{2}\right)^2+\left(\frac{z}{3}\right)^2\le 1 \right\}

L basically means that integration is over boundary (\,dS)
I am trying to approximate the singular integral I for various \vec z. Is there any way to compute approximate this singular surface integral in BEMPP? Is there any source code available for similar problem to the integral above?
What are some approximation schemes implemented in BEMPP for such integrals?


The integral you want can be computed using the Laplace single layer potential operator. This operator computes:


where \vec{y} is a point not on \partial\Omega, and G(\vec{x},\vec{y})=1/(4\pi|\vec{x}-\vec{y}|). So in your case, you should get the result you want if you take:

\displaystyle f(\vec{x})=4\pi\mathrm{e}^{0.1(x_1+2x_2+3x_3)}.

The following code (using Bempp-cl) does this for the points (1,1,0), (2, 2, 0), (5, 6, 0), and (0, 0, 0):

import bempp.api
import numpy as np

grid = bempp.api.shapes.ellipsoid(1, 2, 3)
space = bempp.api.function_space(grid, "DP", 0)

def f(x, n, domain_index, result):
    result[0] = 4 * np.pi * np.exp(0.1 * (x[0] + 2 * x[1] + 3 * x[2]))

function = bempp.api.GridFunction(space, fun=f)

points = np.array([
    [1., 2., 5., 0.], [1., 2., 6., 0.], [0., 0., 0., 0.],

single = bempp.api.operators.potential.laplace.single_layer(space, points)


@timo will have to let you know which approximation schemes are used to compute these integrals, as I’m not sure.

1 Like

Wow! Thanks for your reply.
I have run the code in BEMPP and I also get this warning
SparseEfficiencyWarning: splu requires CSC matrix format warn('splu requires CSC matrix format', SparseEfficiencyWarning)
Is there any paper out there explaining how BEMPP computes the above integral using Laplace single layer potential operator?
In addition, print(single.evaluate(function)) always returns a numpy array with three numbers. It seems it is evaluating grid function on a single element (Triangle). How can I get overall integral value?
For instance I changed the points to

points = np.array([
    [1., 2., 5.],

I get as output

[[12.33608467 29.76894033 29.76894033]]

Why do I get a numpy list? what is this list? I want the value of surface integral

I have read this
It seems the points are considered as the column of the array(Matrix)?
Consider this

points = np.array([
    [1., 2., 5.], 
    [1., 2., 6.], 
    [0., 0., 0.],

Points are (1, 1, 0), (2, 2, 0), (5, 6, 0) and each column of output is the value of surface integral for each point(column)

[[12.33608467 29.76894033 29.76894033]]

Am I wrong?

oops, yes looks like I mixed up the definition of the points array. I’ll edit my original reply so that it’s correct

1 Like

However, when you only give it something like this

points = np.array([
    [1., 2., 5.],

Still outputs something like this

[[12.33608467 29.76894033 29.76894033]]

I wonder what sort of point it calculates according to the input given?!
Points (1,?,?), (2,?,?) and (5,?,?)
What does the library fill question marks (?) with?