← Back to top

Solving the Laplacian on an unstructered grid using the finite element method in Python

One of the classical problems in numerical computing is solving the Poission equation: \nabla^{{2}}{u}={1} Here I will show how the finite element method together with open source Python packages can be used to solve this equation on an unstructured grid. All the code is available on Github.

Weak formulation

The idea behind the finite element method is to take a set of basis functions and expressing the solution {u} as a linear combination of these basis functions. The coefficients determining {u} can be cast in a linear system by testing (taking the inner product) with each individual basis function. \int_{\Omega}\phi_{{i}}\nabla^{{2}}{u}{d}\Omega=\int_{\Omega}\phi_{{i}}{d}\Omega \int_{\Omega}\phi_{{i}}\nabla^{{2}}{u}{d}\Omega=\int_{\Omega}\phi_{{i}}\nabla\cdot\nabla{u}{d}\Omega=\int_{\Gamma}\phi_{{i}}\nabla{u}\cdot\vec{{{n}}}{d}\Gamma-\int_{\Gamma}\nabla\phi_{{i}}\cdot\nabla{u}{d}\Gamma=\int_{\Omega}\phi_{{i}}{d}\Omega where the last equality follows from Green's theorem.

To keep things simple for now let's demand {u}={0} on the boundaries. To fulfill this we take \phi_{{i}}={0} on the boundary nodes and remove the boundary integral. This is actually equivalent to removing these basis functions from the linear system and taking {c}_{{i}}={0} .

To find the linear system we express {u} itself as a linear system of the basis functions {u}={\sum_{{{j}={1}}}^{{n}}}{c}_{{j}}\phi_{{j}} (Galerkin method). Substitution in the previous integral gives {\sum_{{{j}={1}}}^{{{n}}}}{c}_{{j}}\int_{\Gamma}\nabla\phi_{{i}}\cdot\nabla\phi_{{j}}{d}\Gamma=-\int_{\Omega}\phi_{{i}}{d}\Omega\qquad{i}\in{\left\lbrace{1},{2},{3},\ldots,{n}\right\rbrace} Since we choose the basis functions ourselves the integrals can be computed. Thus we can we can read these equations as a linear system, where we have {n} unknowns represented by {c}_{{j}} and {n} equations to find them.

Geometry and basis functions

Let's choose a simple but not trivial geometry, represented by two dimensional disk with a rectangular shaped interior border. We create this geometry in gmsh and save the resulting mesh as a '.msh' file. Gmsh has the concept of 'physical groups' to tag vertices in the mesh so we can find them easily later in our program. As we wan't to exclude the border vertices from our matrix we give them the physical names 'inner-border' and 'outer-border'. Each of these physical names is assigned a number, and every vertex in the mesh file gets assigned such a number, also referenced as a 'tag'. We can read this mesh file in our Python program using meshio.
Figure 1. the tool gmsh is used to create a geometry consiting of a disk with a rectangular shaped interior border.
mesh = meshio.read('assets/domain.msh')

vertices = mesh.points

(inner_border_tag, _) = mesh.field_data['inner-border']
(outer_border_tag, _) = mesh.field_data['outer-border']
(domain_tag, _) = mesh.field_data['domain']

line_tags = mesh.cell_data_dict['gmsh:physical']['line']
triangle_tags = mesh.cell_data_dict['gmsh:physical']['triangle']

lines = mesh.cells_dict['line']
triangles = mesh.cells_dict['triangle']

# We do not want to include the sea vertices
inactive = np.full(len(vertices), False)
#inactive[ np.ndarray.flatten(lines[line_tags == border_tag]) ] = True
inactive[ np.ndarray.flatten(lines[line_tags == inner_border_tag]) ] = True
inactive[ np.ndarray.flatten(lines[line_tags == outer_border_tag]) ] = True

# To map the index in de 'vertices' array to an index
# into our matrix we need to subtract all inactive vertices
# that have come before the vertex we're currently 'looking' at.
subtract = np.cumsum(inactive)

map_index = np.arange(len(vertices)) - subtract
Let's assign a basis function \phi_{{i}} to every node {i} and require that \phi_{{i}}={1} at the node and falls linearly to zero at the neighboring nodes. At all other nodes the function is zero. These basis functions can be written as \phi_{{i}}={1}-{C}_{{1}}{x}-{C}_{{2}}{y} on the neighoubring triangles for suitable constants {C}_{{1}},{C}_{{2}} which depend only on the triangle dimensions. See pages 116-119 of [1] for a good description of these basis elements and information on how to compute the relevant integrals.

We loop over every triangle, and for every node in the triangle that is not part of the border we add a contribution to the linear system. We then solve the linear system using a sparse matrix solver from the scipy package, and visualize the results using matplotlib.

Python solution




[1] J. van Kan, A. Segal, F. Vermolen. Numerical Methods in Scientific Computing. 2014.