Boundary elements method — calculation of solid angle

Question

I am developing a BEM code based on a deal.ii tutorial, see https://www.dealii.org/8.3.0/doxygen/deal.II/step_34.html . Consider the Poisson equation

$$ \Delta u=-f\,, $$ and its Green's function $G\left(\mathbf{x},\mathbf{x}'\right)$ with the property $$ \Delta G=-\delta\left(\mathbf{x}-\mathbf{x}'\right)\,. $$ Applying Green's 2nd identity leads to the well known boundary integral equation (BIE) $$ \int\limits_{\Omega}u\left(\mathbf{x}\right)\delta\left(\mathbf{x}-\mathbf{x}'\right)\mathrm{d}V =\int\limits_{\Omega}f\left(\mathbf{x}\right)G\left(\mathbf{x},\mathbf{x}'\right)\mathrm{d}V +\oint\limits_{\partial\Omega}\mathbf{n}\cdot\left[\cdots\right]\mathrm{d}\Gamma $$ I am considering the integral on left hand side of the equation. It is well known that $$ \int\limits_{\Omega}u\left(\mathbf{x}\right)\delta\left(\mathbf{x}-\mathbf{x}'\right)\mathrm{d}V=\alpha\left(\mathbf{x}'\right)u\left(\mathbf{x}'\right)\,, $$ where as the factor $\alpha$, i.e. the fraction of the solid angle, depends on the location of the point $\mathbf{x}'$. Typically it is given by $$ \alpha\left(\mathbf{x}\right)=\begin{cases} 0 & \mathbf{x}\notin\Omega\\ 1 & \mathbf{x}\in\Omega\\ \frac{1}{2} & \mathbf{x}\in\Gamma & \text{smooth boundary}\\ \frac{\theta}{4\pi} & \mathbf{x}\in\Gamma & \text{3D corner} \end{cases} $$ You get these expressions by simple analytic evaluations of the integral. If the BIE is discretized, these analytic expressions may be not be appropiate any more. Therefore, the author of the tutorial states that the fraction of solid angle may be calculate by $$ \alpha\left(\mathbf{x}'\right)=1+\oint\mathbf{n}\cdot\nabla G\left(\mathbf{x},\mathbf{x}'\right)\mathrm{d}\Gamma\,, $$ and in fact uses this term in the discretization.

Finally my question: How can be the equation above can be derived?

Answer 1

If a homogeneous Neumann problem is considered, i.e. $f=0$ in $\Omega$ and $\mathbf{n}\cdot\nabla\phi$ on $\Gamma$, one solution is given by the constant function

$$ \phi\left(\mathbf{x}\right)=\phi_0\,. $$

Then the boundary integral equation reads

$$\int\limits_{\Omega}\phi\left(\mathbf{x}' \right)\delta\left(\mathbf{x}'-\mathbf{x}\right)\mathrm{d}V' =\phi_0\int\limits_{\Omega}\delta\left(\mathbf{x}'-\mathbf{x}\right)\mathrm{d}V' =-\phi_0\oint\limits_{\Gamma}\tilde{\mathbf{n}}\cdot \left[\nabla' G\left(\mathbf{x}'-\mathbf{x}\right)\right]\mathrm{d}\Gamma'\,,$$

$\tilde{\mathbf{n}}$ is the outward normal of the considered domain. Let the considered domain be the whole space substracted with an interior domain $\Gamma_0$. Then the boundary integral reads $$ -\oint\limits_{\Gamma}\tilde{\mathbf{n}}\cdot \left[\nabla' G\left(\mathbf{x}'-\mathbf{x}\right)\right]\mathrm{d}\Gamma' =-\oint\limits_{\Gamma_\infty}\tilde{\mathbf{n}}\cdot \left[\nabla' G\left(\mathbf{x}'-\mathbf{x}\right)\right]\mathrm{d}\Gamma'+\\ +\oint\limits_{\Gamma_0}\mathbf{n}\cdot \left[\nabla' G\left(\mathbf{x}'-\mathbf{x}\right)\right]\mathrm{d}\Gamma'\,,$$ where $\mathbf{n}$ is the inward normal. It can be shown that the integral on $\Gamma_\infty$ is unity (use spherical coordinates). Using this result and the 2nd equation, we have $$ \alpha\left(\mathbf{x}\right) =1+\oint\limits_{\Gamma_0}\mathbf{n}\cdot \left[\nabla' G\left(\mathbf{x}'-\mathbf{x}\right)\right]\mathrm{d}\Gamma'\,. $$