Solving a Laplace problem with Dirichlet boundary conditions

Download this tutorial as a Jupyter notebook

Background

In this tutorial, we will solve a simple Laplace problem inside the unit sphere with Dirichlet boundary conditions. Let \displaystyle \Omega be the unit sphere with boundary \displaystyle \Gamma. Let \displaystyle \nu be the outward pointing normal on \displaystyle \Gamma. The PDE and boundary conditions are given by

\Delta u = 0\quad\text{in }\Omega,\\  u = f\quad\text{on }\Gamma.

The boundary data is a source \displaystyle \hat{u} located at the point \displaystyle (0.9,0,0).

\displaystyle   \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

\displaystyle   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 \displaystyle u in \displaystyle \Omega satisfies

\displaystyle   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

\displaystyle   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 \displaystyle \mathcal{V} and \displaystyle \mathcal{K} are the single and double layer potential operators.

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

\displaystyle   \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, \displaystyle \mathsf{V} and \displaystyle \mathsf{K} are the single and double layer boundary operators.

Implementation

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

import bempp.api
import numpy as np

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.

grid = bempp.api.shapes.sphere(h=0.1)

We now define the spaces. For this example we will use two spaces: the space of continuous, piecewise linear functions; and 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.

dp0_space = bempp.api.function_space(grid, "DP", 0)
p1_space = bempp.api.function_space(grid, "P", 1)

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

identity = bempp.api.operators.boundary.sparse.identity(
    p1_space, p1_space, dp0_space)
dlp = bempp.api.operators.boundary.laplace.double_layer(
    p1_space, p1_space, dp0_space)
slp = bempp.api.operators.boundary.laplace.single_layer(
    dp0_space, p1_space, dp0_space)

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

def dirichlet_data(x, n, domain_index, result):
    result[0] = 1./(4 * np.pi * ((x[0] - .9)**+ x[1]**+ x[2]**2)**(0.5))
    
dirichlet_fun = bempp.api.GridFunction(p1_space, fun=dirichlet_data)

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

rhs = (.5 * identity + dlp) * dirichlet_fun

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

neumann_fun, info = bempp.api.linalg.cg(slp, rhs, tol=1E-3)

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

n_grid_points = 150
plot_grid = np.mgrid[-1:1:n_grid_points*1j, -1:1:n_grid_points*1j]
points = np.vstack((plot_grid[0].ravel(),
                    plot_grid[1].ravel(),
                    np.zeros(plot_grid[0].size)))

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.

slp_pot = bempp.api.operators.potential.laplace.single_layer(
    dp0_space, points)
dlp_pot = bempp.api.operators.potential.laplace.double_layer(
    p1_space, points)
u_evaluated = slp_pot * neumann_fun - dlp_pot * dirichlet_fun

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.

# Filter out solution values associated with points outside the unit circle.
u_evaluated = u_evaluated.reshape((n_grid_points,n_grid_points))
radius = np.sqrt(plot_grid[0]**+ plot_grid[1]**2)
u_evaluated[radius>1] = np.nan

# Plot the image
from matplotlib import pyplot as plt

plt.imshow(np.log(np.abs(u_evaluated.T)), extent=(-1,1,-1,1))
plt.title('Computed solution')
plt.colorbar()