Essential boundary conditions

$\newcommand{\bm}[1]{\boldsymbol{\mathbf{#1}}} \newcommand{\bvec}[1]{\mathbf{#1}} \newcommand{\bx}{\mathbf{x}} \newcommand{\bu}{\mathbf{u}} \newcommand{\bmv}{\mathbf{v}} \newcommand{\bmw}{\mathbf{w}} \newcommand{\Nh}[1][]{N_{h,#1}} \newcommand{\bU}{\mathbf{U}} \newcommand{\bp}{\bvec{p}} \newcommand{\VV}{\mathbf{V}} \newcommand{\WW}{\mathbf{W}} \newcommand{\QQ}{\mathbf{Q}} \renewcommand{\AA}{\mathbf{A}} \newcommand{\BB}{\mathbf{B}} \newcommand{\FF}{\mathbf{F}} \newcommand{\GG}{\mathbf{G}} \newcommand{\CC}{\mathbf{C}} \newcommand{\MM}{\mathbf{M}} \newcommand{\XX}{\mathbf{X}} \newcommand{\KK}{\mathbf{K}} \newcommand{\UU}{\mathbf{U}} \newcommand{\II}{\mathbf{I}} \renewcommand{\SS}{\mathbf{S}} \newcommand{\HH}{\mathbf{H}} \newcommand{\DD}{\mathbf{D}} \newcommand{\NN}{\mathbf{N}} \newcommand{\bmg}{\mathbf{g}} \newcommand{\R}{\mathbb{R}} \newcommand{\JJ}{\mathbf{J}} \newcommand{\RR}{\mathbf{R}}$

We describe here the way essential (also called Dirichlet) boundary conditions are imposed in the redbKIT library.

Poisson problem

For the sake of simplicity, we restrict the discussion to the following Poisson problem:

$\begin{equation} \left\{ \begin{aligned} - \Delta u &= f & \qquad \quad \text{in } & \Omega \\ u &= g & \qquad \quad \text{on } &\Gamma = \partial \Omega. \end{aligned} \right. \label{eq:LaplStrong} \end{equation}$

Its weak formulation reads: find $u \in V = \{ v \in H^1(\Omega) : v|_{\Gamma} = g\}$ such that

$\begin{equation} \underbrace{\int_{\Omega} \nabla u \cdot \nabla v \, d\Omega}_{ a(u,v)} = \underbrace{\int_{\Omega} fv \, d\Omega}_{F(v)} \qquad \forall v \in V_0 = H^1_0(\Omega). \label{eq:LaplWeak} \end{equation}$

Finite element approximation

We seek an approximate solution to \eqref{eq:LaplWeak} in a discrete space $X_h \subset H^1(\Omega)$ of dimension $N$ made of Lagrange finite elements. We denote by ${\phi_1, \dots, \phi_N}$ the nodal basis of $X_h$ and by $\{ \mathbf{n}_1, \dots, \mathbf{n}_N \}$ the associated nodes. We assume without loss of generality that the nodes located at the boundary are $\{\mathbf{n}_1, \dots, \mathbf{n}_{N_\Gamma}\},$ so that

$V_{0h} = \text{span} \{ \phi_{N_\Gamma+1}, \dots, \phi_N \} = \{ v_h \in X_h \, : \, v_h|_{\Gamma} = 0 \}.$

Finally, we introduce the Lagrange interpolant of the essential data $g$, i.e.

$g_h = \sum_{k=1}^{N_\Gamma} g(\mathbf{n}_k) \, \phi_k|_{\Gamma}$

The finite element approximation reads:

find $ u_h \in V_h = \{ v_h \in X_h : v_h|_{\Gamma} = g_h\} $ such that $$ \begin{equation} a(u_h,v_h) = F(v_h) \qquad \forall v \in V_{0h}. \label{eq:LaplWeakFE} \end{equation} $$

Algebraic formulation

At the algebraic level, we obtain the following linear system

$\begin{equation} \AA \bu = \mathbf{f}, \label{eq:LaplAlg} \end{equation}$

where $\AA \in \R^{N \times N}$ is the stiffness matrix, $\mathbf{f} \in \R^{N}$ is a vector encoding the action of the forcing term, $\bu \in \R^{N}$ is the vector of FE degrees of freedom (dofs), i.e.

$u = \sum_{i=1}^N \bu_i \phi_i.$

To distinguish the interior degrees of freedom from those associated with Dirichlet data, we denote by

$\bu = \begin{pmatrix} \bu_{\Gamma} \\[2pt] \bu_{I} \end{pmatrix}$

the vector of unknowns split in two subvectors, the one ($\bu_\Gamma$) related with the nodes lying on the boundary $\Gamma$, and that ($\bu_I$) related with the internal nodes. At the algebraic level, problem \eqref{eq:LaplWeakFE} reads

$\begin{pmatrix} \II_{\Gamma \Gamma} & 0 \\ \AA_{II} & \AA_{I\Gamma} \\[2pt] \end{pmatrix} \begin{pmatrix} \bu_{\Gamma} \\[2pt] \bu_{I} \end{pmatrix} = \begin{pmatrix} \mathbf{g} \\[2pt] \mathbf{f}_{I} \end{pmatrix}$

where $\mathbf{g}_i = g(\mathbf{n}_i)$ for $1 \leq i \leq N_{\Gamma}$.

As the subsytem associated with the Dirichlet dofs is trivial ($\bu_{\Gamma} = \mathbf{g}$), it is natural to eliminate the essential dofs $\bu_{\Gamma}$, yielding

$\AA_{II} \bu_I = \mathbf{f}_I - \AA_{I\Gamma} \mathbf{g}.$

The resulting linear system enjoys the following remarkable properties:

if $\AA$ is symmetric, then also $\AA_{II}$ is symmetric;
it has minimal size $N-N_{\Gamma}$, since only internal dofs are unknowns;
$\AA$ and $\AA_{II}$ have different sparsity pattern.

Extension to nonlinear problems

We now consider the case of nonlinear PDEs, whose finite element approximation results in a nonlinear system of equations

$\RR(\bu) = \mathbf{0}.$

If such a system is solved by means of Newton’s method, at each iteration we are required to solve a linear system of the form

$\JJ(\bu^{k}) \delta \bu = - \RR(\bu^k),$

and then set $\bu^{k+1} = \bu^k + \delta \bu$. In this case we impose essential boundary conditions as follows:

we first set $$ \bu^{0} = \begin{pmatrix} \mathbf{g} \\[2pt] \bu^0_{I} \end{pmatrix}; $$ then, for $k=0,1,\dots$ until convergence, we solve $$ \JJ_{II}(\bu^{k}) \, \delta \bu_I = - \RR_I(\bu^k), $$ and update $\bu^{k+1}$ as follows $$ \bu^{k+1} = \bu^k + \begin{pmatrix} \mathbf{0} \\[2pt] \delta \bu_I \end{pmatrix}. $$

References

Ern, Alexandre, and Jean-Luc Guermond. Theory and practice of finite elements. Vol. 159. Springer Science & Business Media, 2013.