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.
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.
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$.
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.
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.
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:
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