Weakly imposing a Dirichlet boundary condition

This tutorial shows how to implement the weak imposition of a Dirichlet boundary condition, as proposed in the paper Boundary Element Methods with Weakly Imposed Boundary Conditions (2019).

First, we import Bempp and NumPy.

In [1]:
import bempp.api
import numpy as np

Next, we define the grid for our problem, and the function spaces that we will use. In this example, we use a sphere with P1 and DUAL0 function spaces.

In [2]:
h = 0.3
grid = bempp.api.shapes.sphere(h=h)
p1 = bempp.api.function_space(grid, "P", 1)
dual0 = bempp.api.function_space(grid, "DUAL", 0)

Next, we define the blocked operators proposed in the paper: $$\left(\left(\begin{array}{cc}-\mathsf{K}&\mathsf{V}\\\mathsf{W}&\mathsf{K}'\end{array}\right)+ \left(\begin{array}{cc}\tfrac12\mathsf{Id}&0\\\beta\mathsf{Id}&-\tfrac12\mathsf{Id}\end{array}\right)\right)


\left(\begin{array}{c}g_\text{D}\\\beta g_\text{D}\end{array}\right),$$ where $\beta>0$ is a parameter of our choice. In this example, we use $\beta=0.1$.

In [3]:
beta = 0.1
multi = bempp.api.BlockedOperator(2,2)
multi[0,0] = -bempp.api.operators.boundary.laplace.double_layer(p1, p1, dual0, assembler="fmm")
multi[0,1] = bempp.api.operators.boundary.laplace.single_layer(dual0, p1, dual0, assembler="fmm")
multi[1,0] = bempp.api.operators.boundary.laplace.hypersingular(p1, dual0, p1, assembler="fmm")
multi[1,1] = bempp.api.operators.boundary.laplace.adjoint_double_layer(dual0, dual0, p1, assembler="fmm")

diri = bempp.api.BlockedOperator(2,2)
diri[0,0] = 0.5 * bempp.api.operators.boundary.sparse.identity(p1, p1, dual0)
diri[1,0] = beta * bempp.api.operators.boundary.sparse.identity(p1, dual0, p1)
diri[1,1] = -0.5 * bempp.api.operators.boundary.sparse.identity(dual0, dual0, p1)

Next, we define the function $g_\text{D}$, and define the right hand side.

Here, we use $$g_\text{D}=\sin(\pi x)\sin(\pi y)\sinh(\sqrt2\pi z),$$ as in section 5 of the paper.

In [4]:
def f(x, n, d, res):
    res[0] = np.sin(np.pi*x[0]) * np.sin(np.pi*x[1]) * np.sinh(np.sqrt(2)*np.pi*x[2])

f_fun = bempp.api.GridFunction(p1, fun=f)

rhs = [2*diri[0,0]*f_fun, diri[1,0]*f_fun]
Numba kernel time: 3.8878002166748047
Numba kernel time: 2.7627720832824707
Numba kernel time: 0.00045013427734375

Now we solve the system. We set use_strong_form=True to activate mass matrix preconditioning.

In [5]:
sol, info, it_count = bempp.api.linalg.gmres(multi+diri, rhs, return_iteration_count=True, use_strong_form=True)
print(f"Solution took {it_count} iterations")
Numba kernel time: 0.002165555953979492
Numba kernel time: 1.853149175643921
Numba kernel time: 0.0011739730834960938
Numba kernel time: 0.0011909008026123047
Numba kernel time: 0.001819610595703125
Solution took 18 iterations

For this problem, we know the analytic solution. We compute the error in the $\mathcal{B}_\text{D}$ norm:

In [8]:
def g(x, n, d, res):
    grad = np.array([
            np.cos(np.pi*x[0]) * np.sin(np.pi*x[1]) * np.sinh(np.sqrt(2)*np.pi*x[2]) * np.pi,
            np.sin(np.pi*x[0]) * np.cos(np.pi*x[1]) * np.sinh(np.sqrt(2)*np.pi*x[2]) * np.pi,
            np.sin(np.pi*x[0]) * np.sin(np.pi*x[1]) * np.cosh(np.sqrt(2)*np.pi*x[2]) * np.pi * np.sqrt(2)
    res[0] = np.dot(grad, n)

g_fun = bempp.api.GridFunction(dual0, fun=g)

e_fun = [sol[0]-f_fun,sol[1]-g_fun]

error = 0
# V norm
slp = bempp.api.operators.boundary.laplace.single_layer(dual0, p1, dual0, assembler="fmm")
hyp = bempp.api.operators.boundary.laplace.hypersingular(p1, dual0, p1, assembler="fmm")
error += np.sqrt(np.dot(e_fun[1].coefficients.conjugate(),(slp * e_fun[1]).projections(dual0)))
error += np.sqrt(np.dot(e_fun[0].coefficients.conjugate(),(hyp * e_fun[0]).projections(p1)))
# D part
error += beta**.5 * e_fun[0].l2_norm()

print(f"Error: {error}")
Numba kernel time: 0.0009734630584716797
Numba kernel time: 0.0009589195251464844
Error: 7.761257842012365