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, "B-P", 1)
dual0 = bempp.api.function_space(grid, "DUAL", 0)

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

\begin{pmatrix}u\\\lambda\end{pmatrix}

= \begin{pmatrix}g_\text{D}\\\beta g_\text{D}\end{pmatrix}, $$ 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)
multi[0,1] = bempp.api.operators.boundary.laplace.single_layer(dual0, p1, dual0)
multi[1,0] = bempp.api.operators.boundary.laplace.hypersingular(p1, dual0, p1)
multi[1,1] = bempp.api.operators.boundary.laplace.adjoint_double_layer(dual0, dual0, p1)

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]

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("Solution took",it_count,"iterations")
Solution took 8 iterations

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

In [6]:
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)
hyp = bempp.api.operators.boundary.laplace.hypersingular(p1, dual0, p1)
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("Error:",error)
Error: 4.247437124527201