How to properly apply non-homogeneous Dirichlet boundary conditions with FEM?

Question

In general, Dirichlet boundary conditions won't be satisfied exactly for FEM for non-homogeneous boundary conditions. The FEM codes I've seen set the degrees of freedom to interpolate the Dirichlet boundary condition but I haven't found any mathematical justification for this. It seems to me that setting essential boundary conditions should probably minimize some functional of the error (e.g. minimize $ ||u -u_h|| $ over the portion of the boundary that the Dirichlet BC is applied ) even though this would be more computationally expensive.

Is there any justification for setting the BC like this and if so, what would the proper norm be?

Answer 1

There is mathematical justification for setting Dirichlet boundary degrees of freedom to a value. However, you should adjust your variational form accordingly. If you are looking at a general problem, say:

Find $u\in\mathcal{U}$ such that

$a(u,w)=l(w) \ \ \forall w\in\mathcal{V}$

where

$\mathcal{U}=\{u:\int \nabla u^2 < \infty, u=g\text{ on }\Gamma_D\}$

$\mathcal{V}=\{u:\int \nabla u^2 < \infty, u=0\text{ on }\Gamma_D\}$

Instead we can write $u = v + g$ where $v\in\mathcal{V}$ and $g$ is the Dirichlet condition. Then the variational form becomes

$a(v+g,w)=l(w)$

or by using the linearity of $a(.,.)$

$a(v,w)=l(w)-a(g,w)$

In a finite element code, you can form your element stiffness matrix as if there were no boundary conditions. Then you take the column of the local matrix which corresponds to the Dirichlet boundary condition, scale it by the coefficient you want to enforce, and subtract it from the right-hand-side. This is the discrete form of what I wrote above, $-a(g,w)$. Then you zero out that column and the corresponding Dirichlet row, placing a 1 in the diagonal and the coefficient you wish to enforce. This decouples the equation from the system and yet sets the value you wish to enforce.

I recommend The Finite Element Method: Linear Static and Dynamic Finite Element Analysis, by Tom Hughes. He has an expanded discussion of this issue starting on page 8.

Answer 2

Here's my justification for setting the nodal values to $g(\mathbf{x})$. There might be some details missing because this is an online answer and not a math paper.

Let's start with the homogeneous variational problem of finding $u \in V$ such that

$$ a(u,v) = \mathcal{B}(v) \quad \forall \quad v \in V $$ where $V$ is the space of $L_2$ functions that vanish at $\Gamma_D$. Let's assume for simplicity the problem is 2D:

Then the Galerkin approximation involves writing both $u$ and $v$ as linear combinations of 32 shape functions $h_j(\mathbf{x})$ times the nodal values $u_j$ and $v_j$ $$ u(\mathbf{x}) = \sum_{j=1}^{32} h_j(\mathbf{x}) \cdot u_j = H \cdot \mathbf{u} $$

$$ v(\mathbf{x}) = \sum_{j=1}^{32} h_j(\mathbf{x}) \cdot v_j = H \cdot \mathbf{v} $$

for a matrix $H \in \mathbb{R}^{1 \times 32}$ such that the problem can be cast in matrix form as

$$ A \cdot \mathbf{u} = \mathbf{b} $$ where the elements of the matrix $A$ are

$$ a_{i,j} = a(h_i,h_j) $$ and the elements of vector $\mathbf{b}$ are

$$ b_i = \mathcal{B}(h_i) $$

Note that since the $h_j(\mathbf{x})$ form a basis of $V$ they all vanish on $\Gamma_D$ so the solution $u(\mathbf{x})$ also vanishes on $\Gamma_D$. This procedure is right but it is inconvenient because it is not easy to find the right shape functions that vanish on $\Gamma_D$. We can simplify it by allowing nodes to exist on $\Gamma_D$

and "extending" the expansion of $u$ and $v$ with zeros for $j=33,34,35$:

$$ u_h(\mathbf{x}) = \sum_{j=1}^{32} h_j(\mathbf{x}) \cdot u_j + \sum_{j=33}^{35} h_j(\mathbf{x}) \cdot 0 = \tilde{H} \cdot \begin{bmatrix}\mathbf{u} \\ \mathbf{0}\end{bmatrix} $$

$$ v(\mathbf{x}) = \sum_{j=1}^{32} h_j(\mathbf{x}) \cdot v_j + \sum_{j=33}^{35} h_j(\mathbf{x}) \cdot 0 = \tilde{H} \cdot \begin{bmatrix}\mathbf{v} \\ \mathbf{0}\end{bmatrix} $$ for an extended matrix $\tilde{H} \in \mathbb{R}^{1 \times 35}$.

First, we note that these "extended" objects

$$ \tilde{\mathbf{v}} = \begin{bmatrix} \mathbf{v} \\ \mathbf{0} \end{bmatrix} \quad \tilde{{A}} = \begin{bmatrix} {A} & {C} \\ {D} & {E} \\ \end{bmatrix} \quad \tilde{\mathbf{u}} = \begin{bmatrix} \mathbf{u} \\ \mathbf{0} \end{bmatrix} \quad \tilde{\mathbf{b}} = \begin{bmatrix} \mathbf{b} \\ \mathbf{e} \end{bmatrix} $$ still represent the original Galerkin problem. In effect,

$$ \begin{aligned} \tilde{\mathbf{v}}^T \cdot \tilde{{A}} \cdot \tilde{\mathbf{u}} &= \tilde{\mathbf{v}}^T \cdot \tilde{\mathbf{b}} \\ \begin{bmatrix} \mathbf{v}^T & \mathbf{0}^T \end{bmatrix} \cdot \begin{bmatrix} {A} & {C} \\ {D} & {E} \\ \end{bmatrix} \cdot \begin{bmatrix} \mathbf{u} \\ \mathbf{0} \end{bmatrix} & = \begin{bmatrix} \mathbf{v}^T & \mathbf{0}^T \end{bmatrix} \begin{bmatrix} \mathbf{b} \\ \mathbf{e} \end{bmatrix} \\ \begin{bmatrix} \mathbf{v}^T & \mathbf{0}^T \end{bmatrix} \cdot \begin{bmatrix} {A} \cdot \mathbf{u} + {C} \cdot \mathbf{0} \\ {D} \cdot \mathbf{u} + {E} \cdot \mathbf{0} \\ \end{bmatrix} &= \begin{bmatrix} \mathbf{v}^T & \mathbf{0}^T \end{bmatrix} \begin{bmatrix} \mathbf{b} \\ \mathbf{e} \end{bmatrix} \\ \mathbf{v}^T \cdot {A} \cdot \mathbf{u} + \mathbf{0}^T \cdot {D} \cdot \mathbf{u} &= \mathbf{v}^T \cdot \mathbf{b} + \mathbf{0}^T \cdot \mathbf{e}\\ \mathbf{v}^T \cdot {A} \cdot \mathbf{u} &= \mathbf{v}^T \cdot \mathbf{b}\\ \end{aligned} $$

Since this identity must hold $\forall \mathbf{v}$, then ${A} \cdot \mathbf{u} - \mathbf{b} = 0$, which is the original homogeneous problem. We now need to prove that if we have

$$ {K} = \begin{bmatrix} {A} & {C} \\ {0} & {I} \\ \end{bmatrix} \quad \mathbf{f} = \begin{bmatrix} \mathbf{b} \\ \mathbf{0} \\ \end{bmatrix} $$ such that ${A} \cdot \mathbf{u} = \mathbf{b}$, where ${I}$ is the identity matrix of size $3 \times 3$, then the vector $\mathbf{\varphi}$ such that ${K} \cdot \mathbf{\varphi} = \mathbf{f}$ is equal to

$$ \mathbf{\varphi} = \begin{bmatrix} \mathbf{u} \\ \mathbf{0} \\ \end{bmatrix} $$

Indeed, let $\mathbf{\varphi} = \begin{bmatrix} \mathbf{\varphi}_1 & \mathbf{\varphi}_2 \end{bmatrix}^T$. Then ${K} \cdot \mathbf{\varphi}$ is

$$ \begin{bmatrix} {A} & {C} \\ {0} & {I} \end{bmatrix} \cdot \begin{bmatrix} \mathbf{\varphi}_1 \\ \mathbf{\varphi}_2 \end{bmatrix} = \begin{bmatrix} {A} \cdot \mathbf{\varphi}_1 + {C} \cdot \mathbf{\varphi}_2\\ {0} \cdot \mathbf{\varphi}_1 + {I} \cdot \mathbf{\varphi}_2 \end{bmatrix} = \begin{bmatrix} \mathbf{b}\\ \mathbf{0} \end{bmatrix} $$

From the second row, $\mathbf{\varphi}_2 = \mathbf{0}$. Replacing this result in the first row, ${A} \cdot \mathbf{\varphi}_1 = \mathbf{b}$. Therefore $\mathbf{\varphi}_1 = {A}^{-1} \cdot \mathbf{b} = \mathbf{u}$.

This part proves that putting a one in the diagonal of the stiffness matrix and a zero in the RHS vector in the rows corresponding to the nodes at $\Gamma_D$ is correct.

Now let's investigate non-homogeneous Dirichlet BCs. The "lifting" procedure from textbooks translates into finding $u_h \in V$ such that

$$ a(u_h,v) = \mathcal{B}(v) - a(u_g,v) \quad \forall \quad v \in V $$

where functions in $V$ vanish at $\Gamma_D$, and $u_g$ is a known function (i.e. a known lifting) that satisfies $u_g(\mathbf{x}) = g(\mathbf{x})$ for $\mathbf{x} \in \Gamma_D$. Note that the right-hand side of the formulation contains known functions only.

Now, let's take back the $u_g$ into the left-hand side (we have already assumed $a$ was bi-linear) then we have to solve

$$ a(u_h+u_g,v) = \mathcal{B}(v) $$

Let's first write the functions $u_h$ and $v$ over the whole set of nodes

$$ u_h(\mathbf{x}) = \sum_{j=1}^{32} h_j(\mathbf{x}) \cdot u_j + \sum_{j=33}^{35} h_j(\mathbf{x}) \cdot 0 = \tilde{H} \cdot \begin{bmatrix}\mathbf{u} \\ \mathbf{0}\end{bmatrix} $$

$$ v(\mathbf{x}) = \sum_{j=1}^{32} h_j(\mathbf{x}) \cdot v_j + \sum_{j=33}^{35} h_j(\mathbf{x}) \cdot 0 = \tilde{H} \cdot \begin{bmatrix}\mathbf{v} \\ \mathbf{0}\end{bmatrix} $$ for the "extended" matrix $\tilde{H}$ of size $1 \times 35$.

We then write the non-homogeneous part $u_g$ as

$$ u_h(\mathbf{x}) = \sum_{j=1}^{32} h_j(\mathbf{x}) \cdot 0 + \sum_{j=33}^{35} h_j(\mathbf{x}) \cdot g(\mathbf{x}_j) = \tilde{H} \cdot \begin{bmatrix}\mathbf{0} \\ \mathbf{g}\end{bmatrix} $$ so the sum $u_h + u_g$ is

$$ u_h(\mathbf{x}) + u_g(\mathbf{x}) = \tilde{H} \cdot \begin{bmatrix}\mathbf{u} \\ \mathbf{g}\end{bmatrix} $$

Now, the stiffness matrix of the first homogeneous problem was $A \in \mathbb{R}^{32 \times 32}$. We then extend it to a new matrix of size $35 \times 35$

$$ \begin{bmatrix} A & C \\ D & E \end{bmatrix} $$ such that the discretized problem is now

$$ \begin{aligned} \begin{bmatrix} \mathbf{v}^T & \mathbf{0}^T \end{bmatrix} \cdot \begin{bmatrix} A & C \\ D & E \end{bmatrix} \cdot \begin{bmatrix} \mathbf{u} \\ \mathbf{g} \end{bmatrix} &= \begin{bmatrix} \mathbf{v}^T & \mathbf{0}^T \end{bmatrix} \cdot \begin{bmatrix} \mathbf{b} \\ \mathbf{e} \end{bmatrix} \\ \begin{bmatrix} \mathbf{v}^T & \mathbf{0}^T \end{bmatrix} \cdot \begin{bmatrix} A \cdot \mathbf{u} + C \cdot \mathbf{g} \\ D \cdot \mathbf{u} + E \cdot \mathbf{g} \end{bmatrix} &= \begin{bmatrix} \mathbf{v}^T & \mathbf{0}^T \end{bmatrix} \cdot \begin{bmatrix} \mathbf{b} \\ \mathbf{e} \end{bmatrix} \\ \mathbf{v}^T \cdot A \cdot \mathbf{u} + \mathbf{v}^T \cdot C \cdot \mathbf{g} &= \mathbf{v}^T \cdot \mathbf{b} \end{aligned} $$ for all $\mathbf{u} \in \mathbb{R}^J$. That is to say,

$$ A \cdot \mathbf{u} + C \cdot \mathbf{g} = \mathbf{b} $$

Now you can prove as an exercise that if we have

$$ K = \begin{bmatrix} A & C \\ 0 & I \\ \end{bmatrix} \quad \mathbf{f} = \begin{bmatrix} \mathbf{b} \\ \mathbf{g} \\ \end{bmatrix} $$ such that $A \cdot \mathbf{u} + C \cdot \mathbf{g} = \mathbf{b}$, then the vector $\mathbf{\varphi}$ such that $K \cdot \mathbf{\varphi} = \mathbf{f}$ is equal to

$$ \mathbf{\varphi} = \begin{bmatrix} \mathbf{u} \\ \mathbf{g} \\ \end{bmatrix} $$

Answer 3

To add up to Nathan's great answer with the variational reasoning, one often needs algorithmic details when implementing finite elements. For example,

I also have a more detailed explanation on the subject in my personal notes. Please see the chapter "Constrained Linear Systems".