BCs.tex

Unlike Lagrange FEs \emph{boundary conditions} (BCs) for $C^1$ FEs can be
challenging to implement. We will demonstrate the difficulty of implementing BCs for
$C^1$ FEs through an example. The example we will use is the \emph{Poisson
problem}, which would normally have a $C^0$ FE discretization. The motivation
for using this problem is that it should be a very familiar implementation in
$C^0$ and thus a useful demonstration of the difficulties of implementing
$C^1$ FEs. Therefore, consider the following \emph{Poisson problem}
\begin{equation}
  \begin{split}
    -\Delta &u = f \quad \text{on } \Omega \\
    &u = 0 \quad \text{on } \partial \Omega.
  \end{split}
  \label{eqn:Poisson}
\end{equation}
The weak form of \eqref{eqn:Poisson} is given by
\begin{equation}
  \begin{split}
    \text{Fi}&\text{nd }u \in X := H^1_0(\Omega) \text{ such that} \\
    &(\nabla u, \nabla v) = (f, v) \quad \forall v \in X.
  \end{split}
  \label{eqn:PoissonWeak}
\end{equation}
Thus, for conforming FEs the FE formulation is given by
\begin{equation}
  \begin{split}
    &\text{Find }u^h \in X^h \subset X \text{ such that} \\
    (\nabla &u^h, \nabla v^h) = (f, v^h) \quad \forall v^h \in X^h.
  \end{split}
  \label{eqn:PoissonFE}
\end{equation}
Since $u \in H^1_0(\Omega)$ and $u^h \in X^h \subset H^1_0(\Omega)$ only $C^0$
Lagrange FEs need be used. However, for demonstration purposes we assume $u^h
\in X^h \subset X \subset C^1(\Omega)$. Therefore, Lagrange FEs are no longer
an appropriate choice for our FE discretization.

With Lagrange FEs the most common way to deal with BCs is to enforce them by
explicitly setting the degrees of freedom (DoFs) on the boundary to zero and
reduce the linear system to contain only the internal nodes. This can be seen
graphically in \autoref{fig:ReducedBCs}. One significant benefit of treating BCs in
this way is that one reduces the size of the linear system that must be solved, and
therefore one reduce computational time.

\input{ReducedBC}

This method, however, can be problematic when using $C^1$ FEs. Consider the
boundary $\partial \Omega = \Gamma_1 \bigcup \Gamma_2$, where $\Gamma_1$ and
$\Gamma_2$ correspond to the vertical and horizontal boundaries, respectively
(see \autoref{fig:SplitBoundary}).  On $\Gamma_1$ we know $u=0$ which implies
$u_y = u_{yy} = 0$ on $\Gamma_1$.  On $\Gamma_2$ we know $u=0$ which implies
$u_x = u_{xx} = 0$ on $\Gamma_2$. Thus, simply setting $u^h = 0$ on the boundary
results in the DoFs corresponding to $u_x,\,u_y,\,u_{xx}$, and $u_{yy}$ not
necessarily being set to zero. Therefore, one must force these DoFs to be zero
too or another scheme must be implemented to deal with this case. However,
setting the appropriate DoFs to zero in this manner requires some analysis of
the problem domain before implementation, which is problematic when one wants to
implement a generalized code. To read more about implementation of boundary
conditions for Hermite type FEs see Remark 2.3.8 in Ciarlet \cite{Ciarlet}. To
avoid analyzing boundary conditions for each domain one can implement a Lagrange
Multiplier scheme
\cite{Babuska1973,Barbosa1991,Barbosa1992,Bramble1981,Pitkaranta1980,Pitkaranta1981}.

\input{SplitBoundary}

We begin by rewriting the FE formulation, \eqref{eqn:PoissonFE}, in the
following way
\begin{equation}
  \begin{split}
    A u &= \ell \\
    \Lambda u &= b,
  \end{split}
  \label{eqn:ConstrainedPoisson}
\end{equation}
where is $\Lambda u = b$ is the constraint equation descibing the boundary
conditions and $A u = \ell$ is the FE formulation of the Poisson equation.
Then the \emph{variational form} of \eqref{eqn:ConstrainedPoisson} is
\begin{equation}
  L(u,\lambda) = \frac{1}{2}u^T A u - u^T \ell + \lambda^T (\Lambda u - b),
  \label{eqn:Variational}
\end{equation}
where $\lambda$ is the Lagrange multiplier. Using the first-order optimality
condition results in
\begin{equation}
  \begin{split}
    \frac{d L}{du} &= Au - \ell + \Lambda^t \lambda = 0 \\
    \frac{d L}{d\lambda} &= \Lambda u - b = 0.
  \end{split}
  \label{eqn:Condition}
\end{equation}
Thus, given $b = \mathbf{0}$ the new system of equations using Lagrange
multipliers for boundary constraints is given by
\begin{equation}
  \begin{bmatrix}
    A & \Lambda^T \\
    \Lambda &  \mathbf{0} \\
  \end{bmatrix} \begin{bmatrix}
    u \\ \lambda
  \end{bmatrix} = \begin{bmatrix}
    \ell \\ \mathbf{0}
  \end{bmatrix}.
  \label{eqn:Lagrange}
\end{equation}
Since the homogeneous boundary conditions were not enforced in $A$, we note that
$A$ is singular.  To avoid poor conditioning of the system \eqref{eqn:Lagrange},
we add a small perturbation matrix, $\varepsilon \mathbf{I}$, in the lower right
corner of the matrix of the linear system \eqref{eqn:Lagrange}. We note that a
similar approach (the penalty method) is used in the Navier-Stokes equations
context (see, e.g., \cite{Gunzburger89}).