In this tutorial, we will solve a simple Laplace problem inside the unit sphere with Dirichlet boundary conditions. Let be the unit sphere with boundary . Let be the outward pointing normal on . The PDE and boundary conditions are given by
The boundary data is a source located at the point .
For this example, we will use a direct integral equation of the first kind. Let
be the Green’s function for Laplace in three dimensions. From Green’s representation theorem, it follows that every harmonic function in satisfies
where and are the single and double layer potential operators.
Taking the limit we obtain the boundary integral equation
Here, and are the single and double layer boundary operators.
We now demonstrate how to solve this problem with Bempp. We begin by importing Bempp and NumPy.
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.
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.
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.
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.
result = 1./(4 * np.pi * ((x - .9)**2 + x**2 + x**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
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 plane for . We first define points at which to plot the solution.
plot_grid = np.mgrid[-1:1:n_grid_points*1j, -1:1:n_grid_points*1j]
points = np.vstack((plot_grid.ravel(),
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.
dlp_pot = bempp.api.operators.potential.laplace.double_layer(
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.
u_evaluated = u_evaluated.reshape((n_grid_points,n_grid_points))
radius = np.sqrt(plot_grid**2 + plot_grid**2)
u_evaluated[radius>1] = np.nan
# Plot the image
from matplotlib import pyplot as plt