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
For this method, we need to use an integral representation of the equation, that, in this case, is
Then, we can form a system of equations
that we obtain by discretization of the boundary. If we take constant variables over the discretization, the integral can be computed analytically by
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
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.
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 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] - coords[elem1]) 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] - coords[elem] dcos = dcos / norm(dcos) rotmat = np.array([[dcos, -dcos], [dcos, dcos]]) r_A = rotmat.dot(coords[elem] - pt_col) r_B = rotmat.dot(coords[elem] - pt_col) theta_A = arctan2(r_A, r_A) theta_B = arctan2(r_B, r_B) if norm(r_A) <= 1e-6: G_coeff = r_B*(log(norm(r_B)) - 1) + theta_B*r_B elif norm(r_B) <= 1e-6: G_coeff = -(r_A*(log(norm(r_A)) - 1) + theta_A*r_A) else: G_coeff = r_B*(log(norm(r_B)) - 1) + theta_B*r_B -\ (r_A*(log(norm(r_A)) - 1) + theta_A*r_A) 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 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 mesh = meshio.read("disk.msh") 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 = fig.add_subplot(111, projection='3d') 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.