Solving a Laplace problem with Dirichlet boundary conditions

Background

• Laplace problem inside the unit sphere with Dirichlet boundary conditions. Let $\Omega$ be the unit sphere with boundary $\Gamma$. Let $\nu$ be the outward pointing normal on $\Gamma$. The PDE and boundary conditions are given by

\begin{align} \Delta u &= 0&&\text{in }\Omega,\\ u &= f&&\text{on }\Gamma. \end{align}

The boundary data is a source $\hat{u}$ located at the point $(0.9,0,0)$. $$ \hat{u}(\mathbf x)=\frac{1}{4\pi\sqrt{(x-0.9)^2+y^2+z^2}}. $$

For this example, we will use a direct integral equation of the first kind. Let $$ g(\mathbf x,\mathbf y) = \frac{1}{4\pi |\mathbf x-\mathbf y|} $$ be the Green's function for Laplace in three dimensions. From Green's representation theorem, it follows that every harmonic function $u$ in $\Omega$ satisfies

$$ u(\mathbf x) = \int_{\Gamma} g(\mathbf x,\mathbf y)\frac{\partial u(\mathbf y)}{\partial \nu(\mathbf{y})}\mathrm{d}\mathbf y-\int_{\Gamma}\frac{\partial g(\mathbf x,\mathbf y)}{\partial \nu(\mathbf{y})}u(\mathbf y)\mathrm{d}\mathbf y,~\mathbf x\in\Omega\setminus\Gamma $$

or equivalantly

$$ u(\mathbf x) = \left[\mathcal{V}\frac{\partial u(\mathbf y)}{\partial \nu(\mathbf{y})}\right] (\mathbf{x}) - \left[\mathcal{K}u\right] (\mathbf{x}),~\mathbf x\in\Omega\setminus\Gamma, $$

where $\mathcal{V}$ and $\mathcal{K}$ are the single and double layer potential operators.

Taking the limit $\mathbf x\rightarrow \Gamma$ we obtain the boundary integral equation

$$ \left[\mathsf{V}\frac{\partial u}{\partial n}\right] (\mathbf x)=\left[(\tfrac12\mathsf{Id}+\mathsf{K})u\right] (\mathbf x),~\mathbf x\in\Gamma. $$

Here, $\mathsf{V}$ and $\mathsf{K}$ are the single and double layer boundary operators.

This equation is the target equation that we will solve with respect to unknown norrmal derivative on the boundary $\frac{\partial u}{\partial n}$.

When we get normal derivative on the boundary $$ \frac{\partial u}{\partial n} $$ we can use Green representation theorem to compute $u(x)$ in any points inside the domain $\Omega$.

Implementation

We now demonstrate how to solve this problem with Bempp. We begin by importing Bempp and NumPy.

We next define a mesh or grid. For this problem, we will use the built-in function sphere that defines a simple spherical grid. Details of other available grids can be found in the grids tutorial.

We now define the spaces. We will use two spaces:

• the space of continuous, piecewise linear functions
• the space of piecewise constant functions.

The space of piecewise constant functions has the right smoothness for the unknown Neumann data. We will use continuous, piecewise linear functions to represent the known Dirichlet data.

We can now define the operators. We need the identity, single layer, and double layer boundary operator.

We now define the GridFunction object on the sphere grid that represents the Dirichlet data.

We next assemble the right-hand side of the boundary integral equation, given by $$ (\tfrac12\mathsf{Id}+\mathsf{K})u. $$

We use implemented operators algebra to define right-hand side.

We now solve the linear system using a conjugate gradient (CG) method.

We now want to provide a simple plot of the solution in the $(x,y)$ plane for $z=0$. We first define points at which to plot the solution.

The variable points now contains in its columns the coordinates of the evaluation points. We can now use Green's representation theorem to evaluate the solution on these points. Note in particular the last line of the following code. It is a direct implementation of Green's representation theorem.

We now plot the 2D slice of the solution. For a full three dimensional visualization, Bempp can export the data to Gmsh. Since the solution decays quickly, we use a logarithmic plot.