# Coming back to the Boundary element method

During October (2017) I wrote a program per day for some well-known numerical methods in both Python and Julia. It was intended to be an exercise to learn some Julia. You can see a summary here. I succeeded in 30 of the challenges, except for the BEM (Boundary Element Method), where I could not figure out what was wrong that day. The original post is here.

Thomas Klimpel found the mistake and wrote an email where he described my mistakes. So, I am creating a new post with a correct implementation of the BEM.

## The Boundary Element Method

We want so solve the equation

\begin{equation*} \nabla^2 u = -f(x, y)\quad \forall (x, y) \in \Omega\, , \end{equation*}

with

\begin{equation*} u(x) = g(x, y), \quad \forall (x, y)\in \partial \Omega \, . \end{equation*}

For this method, we need to use an integral representation of the equation, that, in this case, is

\begin{equation*} u(\boldsymbol{\xi}) = \int_{S} [u(\mathbf{x}) F(\mathbf{x}, \boldsymbol{\xi}) - q(\mathbf{x})G(\mathbf{x}, \boldsymbol{\xi})]\mathrm{d}S(\mathbf{x}) + \int_{V} f(\mathbf{x}) G(\mathbf{x}, \boldsymbol{\xi}) \mathrm{d}V(\mathbf{x}) \end{equation*}

with

\begin{equation*} G(\mathbf{x}, \boldsymbol{\xi})= -\frac{1}{2\pi}\ln|\mathbf{x} - \xi| \end{equation*}

and

\begin{equation*} F(\mathbf{x}, \boldsymbol{\xi}) = -\frac{1}{2\pi |\mathbf{x} - \boldsymbol{\xi}|^2} (\mathbf{x} - \boldsymbol{\xi})\cdot\hat{\mathbf{n}} \end{equation*}

Then, we can form a system of equations

\begin{equation*} [G]\{q\} = [F]\{u\}\, , \end{equation*}

that we obtain by discretization of the boundary. If we take constant variables over the discretization, the integral can be computed analytically by

\begin{equation*} G_{nm} = -\frac{1}{2\pi}\left[r \sin\theta\left(\ln|r| - 1\right) + \theta r\cos\theta\right]^{\theta_B, r_B}_{\theta_A, r_A} \end{equation*}

and

\begin{equation*} F_{nm} = \left[\frac{\theta}{2\pi}\right]^{\theta_B}_{\theta_A} \end{equation*}

for points $n$ and $m$ in different elements, and the subindices $A,B$ refer to the endpoints of the evaluation element. We should be careful evaluating this expression since both $r_A$ and $r_B$ can be (close to) zero and make it explode. Also, here it was the main problem since I forgot to compute the angles with respect to elements that are, in general, not aligned with horizontal or vertical axes.

For diagonal terms the integrals evaluate to

\begin{equation*} G_{nn} = -\frac{L}{2\pi}\left(\ln\left\vert\frac{L}{2}\right\vert - 1\right) \end{equation*}

and

\begin{equation*} F_{nn} = - \frac{1}{2} \end{equation*}

with $L$ the size of the element.

Following is the code. Keep in mind that this code works for purely Dirichlet problems. For mixed Dirichlet-Neumann the influence matrices would need rearranging to separate known and unknowns in opposite sides of the equation.

You can download the files for this project here. It includes a YML file to create a conda environment with the dependencies listed. For example, it uses version 3.0 of Meshio.

import numpy as np
from numpy import log, arctan2, pi, mean, arctan
from numpy.linalg import norm, solve
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import meshio

def assem(coords, elems):
"""Assembly matrices for the BEM problem

Parameters
----------
coords : ndarray, float
Coordinates for the nodes.
elems : ndarray, int
Connectivity for the elements.

Returns
-------
Gmat : ndarray, float
Influence matrix for the flow.
Fmat : ndarray, float
Influence matrix for primary variable.
"""
nelems = elems.shape[0]
Gmat = np.zeros((nelems, nelems))
Fmat = np.zeros((nelems, nelems))
for ev_cont, elem1 in enumerate(elems):
for col_cont, elem2 in enumerate(elems):
pt_col = mean(coords[elem2], axis=0)
if ev_cont == col_cont:
L = norm(coords[elem1[1]] - coords[elem1[0]])
Gmat[ev_cont, ev_cont] = - L/(2*pi)*(log(L/2) - 1)
Fmat[ev_cont, ev_cont] = - 0.5
else:
Gij, Fij = influence_coeff(elem1, coords, pt_col)
Gmat[ev_cont, col_cont] = Gij
Fmat[ev_cont, col_cont] = Fij
return Gmat, Fmat

def influence_coeff(elem, coords, pt_col):
"""Compute influence coefficients

Parameters
----------
elems : ndarray, int
Connectivity for the elements.
coords : ndarray, float
Coordinates for the nodes.
pt_col : ndarray
Coordinates of the colocation point.

Returns
-------
G_coeff : float
Influence coefficient for flows.
F_coeff : float
Influence coefficient for primary variable.
"""
dcos = coords[elem[1]] - coords[elem[0]]
dcos = dcos / norm(dcos)
rotmat = np.array([[dcos[1], -dcos[0]],
[dcos[0], dcos[1]]])
r_A = rotmat.dot(coords[elem[0]] - pt_col)
r_B = rotmat.dot(coords[elem[1]] - pt_col)
theta_A = arctan2(r_A[1], r_A[0])
theta_B = arctan2(r_B[1], r_B[0])
if norm(r_A) <= 1e-6:
G_coeff = r_B[1]*(log(norm(r_B)) - 1) + theta_B*r_B[0]
elif norm(r_B) <= 1e-6:
G_coeff = -(r_A[1]*(log(norm(r_A)) - 1) + theta_A*r_A[0])
else:
G_coeff = r_B[1]*(log(norm(r_B)) - 1) + theta_B*r_B[0] -\
(r_A[1]*(log(norm(r_A)) - 1) + theta_A*r_A[0])
F_coeff = theta_B - theta_A
return -G_coeff/(2*pi), F_coeff/(2*pi)

def eval_sol(ev_coords, coords, elems, u_boundary, q_boundary):
"""Evaluate the solution in a set of points

Parameters
----------
ev_coords : ndarray, float
Coordinates of the evaluation points.
coords : ndarray, float
Coordinates for the nodes.
elems : ndarray, int
Connectivity for the elements.
u_boundary : ndarray, float
Primary variable in the nodes.
q_boundary : [type]
Flows in the nodes.

Returns
-------
solution : ndarray, float
Solution evaluated in the given points.
"""
npts = ev_coords.shape[0]
solution = np.zeros(npts)
for k in range(npts):
for ev_cont, elem in enumerate(elems):
pt_col = ev_coords[k]
G, F = influence_coeff(elem, coords, pt_col)
solution[k] += u_boundary[ev_cont]*F - q_boundary[ev_cont]*G
return solution

#%% Simulation
elems = mesh.cells["line"]
bound_nodes = list(set(elems.flatten()))
coords = mesh.points[bound_nodes, :2]
x, y = coords.T
x_m, y_m = 0.5*(coords[elems[:, 0]] + coords[elems[:, 1]]).T
theta = np.arctan2(y_m, x_m)
u_boundary = 3*np.cos(6*theta)

#%% Assembly
Gmat, Fmat = assem(coords, elems)

#%% Solution
q_boundary = solve(Gmat, Fmat.dot(u_boundary))

#%% Evaluation
ev_coords =  mesh.points[:, :2]
ev_x, ev_y = ev_coords.T
solution = eval_sol(ev_coords, coords, elems, u_boundary, q_boundary)

#%% Visualization
tris = mesh.cells["triangle"]
fig = plt.figure()
ax.plot_trisurf(ev_x, ev_y, solution, cmap="RdYlBu", lw=0.3,
edgecolor="#3c3c3c")
plt.xticks([])
plt.yticks([])
ax.set_zticks([])
plt.savefig("bem_solution.png", bbox_inches="tight", transparent=True,
dpi=300)


The result in this case is the following.