I will in this post attempt to summarise some different approaches to handling Dirichlet boundary conditions in implementations of FEM codes. We shall be concerned with the (weak) Poisson equation with inhomogeneous boundary condition,

\begin{equation} -\nabla\cdot\nabla u = f \text{ in } \Omega,\ u = g \text{ on } \partial\Omega. \end{equation}

Here \(f\in L^2(\Omega),g\in L^{\infty}(\partial\Omega)\) represent data which we use to solve for \(u\in H^1(\Omega)\). Since then the trace of \(u\) is equal to \(g\), we can say that \(g\in H^{1/2}(\partial\Omega)\). Recast in a variational form, one typically wants to find \(u\in H^1_g(\Omega)= \\\{ v\in H^1(\Omega) : v\lvert_{\partial\Omega} = g \\\}\) solving

\begin{equation} \int_{\Omega} \nabla u \cdot \nabla v = \int_{\Omega} fv \ \ \forall v\in H^1_0(\Omega). \label{eq1}\end{equation}

To put the condition \(u\lvert_{\partial\Omega} = g\) into the computer is the tricky part once we know how to discretise. This formulation is also difficult to analyse with respect to well-posedness¹. Fortunately we may reduce to the homogeneous case in the following way. Set \(w=u-\hat{g}\), where \(\hat{g}\in H^1(\Omega)\) is the extension of \(g\) to the entire domain (this exists if our domain is Lipschitz). We shall call this extension simply \(g\) in the following. Then \(w\in H^1_0(\Omega)\) and the equation for \(w\) becomes homogeneous; its variational form nets us a unique solution. We may also find a unique extension of \(g\) by taking the extension with the smallest (infimum) \(H^1\) norm, thus getting a unique solution \(u\). The modified variational form for \(w\) becomes

\begin{equation} \int_{\Omega} \nabla w \cdot \nabla v = \int_{\Omega} fv - \int_{\Omega} \nabla g \cdot \nabla v \ \ \forall v\in H^1_0(\Omega). \end{equation}

From here it is still not entirely clear how we should discretise. The only unknown data lies in our function \(w\) so we should respect that when solving the discrete system. But what do we do about \(\nabla g\), in particular \(g\) could be sampled from some physical system so its gradient might not be something we can uniquely construct. We could try projecting it, but we can actually avoid the problem entirely.

The natural choice of finite elements for the Poisson equation are the Lagrange \(\mathbf{P}^k\) elements whose dual space consists in part of nodal evaluations. A subspace will be the piecewise linears \(\mathbf{P}^1\) so consider the following. Set \(V_h\subset H^1(\Omega),\ V_{h,g}\subset H^1_g(\Omega)\) as finite² dimensional \(\mathbf{P}^k\) subspaces. Decompose the discrete solution \(u_h\in V_{h,g}\) as \(u_h=w_h+u_{h,g}\) where \(w_h \in V_{h,0}\) and \(u_{h,g}\in V_{h,g}\). Then we can pick \(u_{h,g}\) to be the function which is \(0\) on the interior of \(\Omega\), meaning that we can solve the discrete version of \eqref{eq1} for \(w_h\) as if our Dirichlet boundary condition was \(g=0\). The nodal values along the boundary \(x_i\) will be \(g(x_i)\) regarded as known quantities.

We will now investigate different philosophies regarding the implementation of this in code. Let \(n=\text{dim}(V_{h,g})\). In all sections below we will consider the linear \(n\times n\) symmetric system arising from discretising \eqref{eq1}, \(Ax=b\) with our unknowns in \(x\).

Note 20/11/23: For reasons I will omit the non-recommended approaches for now.

Nice way: restriction and extension

Back to the drawing board then. As is the case in many situations of a mathematical nature, it helps to consider what is actually happening on a higher level of mathematical abstraction. When the solution is known on the boundary, what we are solving for is the restriction of the solution on the interior. When we then have the interior solution we want to extend back to the entire domain. Restriction and extension are in this context linear operators which can be represented by matrices (our spaces are finite-dimensional).

Let \(n_I\) be the number of dof’s without the boundary nodes. Let \(R:\mathbb{R}^{n}\to \mathbb{R}^{n_I}\) be the restriction to the interior and let \(R_{\partial}:\mathbb{R}^{n}\to \mathbb{R}^{n-n_I}\) be the restriction to the boundary. Then if \(x_{\partial}\in\mathbb{R}^{n-n_I}\) is the dof vector of boundary evaluations,

\begin{equation} RAR^Tx_I = R(b - A R_{\partial}^Tx_{\partial}) \label{eq:2}\end{equation}

is solvable for the interior dof’s \(x_I\), after which we get the entire solution from \(R^T x_I + R_{\partial}^T x_{\partial}\).

This procedure changes nothing about the system, it retains its scaling and is still symmetric. Basically \eqref{eq:2} is exactly the system we want to solve. Pretty nice I would say.

¹ Note that the space \(H^1_g(\Omega)\) is not a vector space!

² If \(g\) is outside the space, this space is defined with respect to subspace projection using the Hilbert inner product in \(H^1\).

³ We have ordered the dof’s so that the node evaluations occur first.