<commit_before>
In this Chapter we implement a finite element method (FEM) for an unstructured triangular mesh on a general polygonal domain.  Our code solves certain linear and nonlinear Poisson equation problems with reasonably-general boundary conditions.  The approach is to
\begin{center}
\emph{first do element-wise assembly of the residual equations,}
\end{center}
and thereby get a functioning \pSNES-based code.  We will not even think about matrices in writing the initial code.  Only after making sure it works do we
\begin{center}
\emph{write additional code to assemble an approximate Jacobian matrix.}
\end{center}

We do the following new tasks in this Chapter:
\begin{itemize}
\item read a triangular mesh from a file into \PETSc \pVecs and \pISs,
\item index that mesh unstructured way,
\item pre-allocate a \PETSc \pMat, and
\item implement Neumann boundary conditions.
\end{itemize}

Similarly to the example in the last Chapter, the residual equations $\bF(\bu)=0$, which is to say the nonlinear system seen by the \pSNES solver, is the discretized weak form of the PDE.  Our direct construction of these equations roughly follows the residual implementations of Chapters \ref{chap:nl}--\ref{chap:of}.

The FEM ``stiffness'' matrix $A$, a major piece in typical FEM introductions \citep{Braess2007,Elmanetal2005}, arises here as the Jacobian when solving $\bF(\bu)=0$.  In our first runs, done before writing any matrix-assembly code, the \pSNES solver approximates this linear system through finite-differencing evaluations of $\bF$.  Eventually we do assemble a Jacobian, but one that is only correct for the linear case.  For the smooth nonlinear case tested here, however, this inaccurate linearization suffices.

The example here contrasts with the rectangular-element ($Q^1$) FEM of the last Chapter by not using a \pDMDA and a structured grid.  Instead we implement a minimal, though naive, mesh-topology infrastructure.  (Doing so substantially increases our workload!)  The relatively-simple \Triangle program generates ASCII files describing the mesh.  An ``index set'' \pIS  type holds the global indices of nodes and triangles (elements).  \pVec and \pMat objects, used for the solution and Jacobian respectively, are indexed through the \pIS values.

In contrast with both earlier and later Chapters, we generate a code that only works on one MPI process (i.e.~a serial code).  However, at the end of the Chapter we briefly consider, though we do not implement, the additional constructs needed to distribute a mesh across many processes.  These considerations motivate the \pDMPlex example in Chapter \ref{chap:dp}.  In fact, the benefits of \PETSc's mesh-topology \pDM tools, including \pDMDA in previous Chapters and \pDMPlex in Chapter \ref{chap:dp}, become clearer by forgoing them here.


\section{A 2D nonlinear Poisson problem}

\begin{marginfigure}
\input{tikz/generalpoissondomain.tex}
\caption{Problem \eqref{eq:un:poissonstrong} on a domain.}
\label{fig:un:generalpoissondomain}
\end{marginfigure}

Let $\Omega \subset \RR^2$ be a bounded (open) region as in Figure \ref{fig:un:generalpoissondomain}.  We suppose its boundary $\partial\Omega$ is well-behaved, for example polygonal or Lipschitz-continuous \citep[section 1.2]{Ciarlet2002}.  Also we decompose the boundary into (measurable) disjoint subsets $\partial_D \Omega$ and $\partial_N \Omega$.

Chapter \ref{chap:st} solves the Poisson equation $- \grad^2 u = f$.  Here we allow a more general nonlinear form.  Let $a(u,x,y)$ and $f(u,x,y)$ be continuous functions, and assume there is $\eps$ so that
\begin{equation}
a(u,x,y) \ge \eps > 0, \label{eq:un:uniformelliptic}
\end{equation}
that is, \emph{uniform ellipticity} \citep{Evans2010}.  Boundary conditions are also needed, and here we consider nonhomogeneous Dirichlet and Neumann type.  Thus the  nonlinear Poisson \emph{strong form} we solve is to find $u(x,y)$ so that
\begin{align}
- \Div \left(a(u,x,y) \grad u\right) &= f(u,x,y) \quad \text{ on } \Omega, \label{eq:un:poissonstrong} \\
u &= g_D \quad \text{ on } \partial_D \Omega, \notag \\
\frac{\partial u}{\partial n} &= g_N \quad \text{ on } \partial_N \Omega. \notag
\end{align}
By definition, $\partial u/\partial n = \bn \cdot \grad u$ where $\bn$ is the outward unit normal on $\partial \Omega$ (Figure \ref{fig:un:generalpoissondomain}).

The data of problem \eqref{eq:un:poissonstrong}, besides the region $\Omega$ and its boundary, includes a \emph{diffusion coefficient} $a$, a \emph{source term} $f$, \emph{Dirichlet data} $g_D$, and \emph{Neumann data} $g_N$.  For the purpose of numerical solutions we will assume these data are continuous.  Poisson problem \eqref{poissonsquare} in Chapter \ref{chap:st} is the homogeneous Dirichlet case where $\Omega$ is a square, $a\equiv 1$, $f=f(x,y)$, $\partial_D \Omega = \partial \Omega$, $\partial_N \Omega = \emptyset$, and $g_D=0$.

One could allow more general boundary conditions in \eqref{eq:un:poissonstrong}, in particular \emph{Robin} conditions $\alpha u + \beta \frac{\partial u}{\partial n} = \gamma$ where $\alpha,\beta,\gamma$ could vary along the boundary \citep{Elmanetal2005}; see Exercise \ref{chap:un}.\ref{exer:un:robin}.  One could also allow $a$ and/or $f$ to depend on the gradient of $u$; recall $a=|\grad u|^{p-2}$ in the $p$-Laplacian equation of Chapter \ref{chap:of}, for example.

As \eqref{eq:un:poissonstrong} is stated there may be no solution where ``$\Div(a\grad u)$'' makes sense as a continuous function, even for polygonal regions and continuous data.  That is, there may be no $u\in C^2(\Omega) \cap C(\overline \Omega)$ so that \eqref{eq:un:poissonstrong} applies at all points.  In the linear case, at least, a solution exists if we convert \eqref{eq:un:poissonstrong} to a \emph{weak form}.  Such a weak formulation already arose in Chapter \ref{chap:of} as the gradient of an objective function, but the weak form here generally does not have a minimization origin; see Exercise \ref{chap:un}.\ref{exer:un:notminimization}.  In this Chapter we will derive the weak form from the strong form in the traditional way, by multiplying by a test function and integrating.


\section{Weak form with general boundary values}

Recall $L^2(\Omega)$ is the space of square-integrable real functions on $\Omega$.  Define the following Hilbert space:
    $$H^1(\Omega) = \left\{u \,:\, u \in L^2(\Omega) \,\&\, \grad u \in L^2(\Omega)\right\}.$$
Definition \eqref{eq:of:sobolevdefn} of Sobolev spaces $W^{1,p}(\Omega)$ includes this space: $H^1(\Omega) = W^{1,2}(\Omega)$.

We use two subsets of $H^1(\Omega)$: \emph{trial functions} come from $H_{g}^1(\Omega)$, the functions with value $g_D$ along $\partial_D \Omega$, and \emph{test functions} come from $H_{0}^1(\Omega)$, with value $0$ along $\partial_D \Omega$.  Note that $H_{0}^1(\Omega)$ is a linear subspace of $H^1(\Omega)$ while generally $H_{g}^1(\Omega)$ is not.

Now choose any $v\in H_{0}^1(\Omega)$, multiply the first equation in \eqref{eq:un:poissonstrong} by $v$, and integrate by parts:
\begin{equation*}
\int_\Omega a(u) \grad u \cdot \grad v - \int_{\partial\Omega} \frac{\partial u}{\partial n} v = \int_\Omega f(u) v.
\end{equation*}
In writing such integrals we will generally show dependence on the solution $u$ but suppress it for $x,y$ so that $a(u,x,y)=a(u)$, and similarly for $f$.

Next we use the boundary information, that $v=0$ on $\partial_D\Omega$ and $\partial u/\partial n=g_N$ on $\partial_N\Omega$:
\begin{equation}
\int_\Omega a(u) \grad u \cdot \grad v = \int_\Omega f(u) v + \int_{\partial_N\Omega} g_N v\quad \text{ for any } v\in H_{0}^1(\Omega). \label{eq:un:poissonweak}
\end{equation}
Equation \eqref{eq:un:poissonweak} is the \emph{weak formulation} of \eqref{eq:un:poissonstrong}.  Any $u \in H_{g}^1(\Omega)$ satisfying \eqref{eq:un:poissonweak} is a \emph{weak solution}.  Observe that the Dirichlet boundary conditions $g_D$ are incorporated into defining $H_{g}^1(\Omega)$ while the Neumann boundary data $g_N$ appears in the weak form \eqref{eq:un:poissonweak}.

Here are the ``rules'' for passing between the strong \eqref{eq:un:poissonstrong} and weak \eqref{eq:un:poissonweak} formulations:\begin{itemize}
\item A well-behaved function $u \in C^2(\Omega) \cap C(\overline \Omega)$ which satisfies \eqref{eq:un:poissonstrong} also solves \eqref{eq:un:poissonweak}.  The derivation above shows this.
\item If $u \in H_{g}^1(\Omega)$ solves \eqref{eq:un:poissonweak} then we accept it, by definition, as a solution.   If it is also well-behaved, $u \in C^2(\Omega) \cap C(\overline \Omega)$, then we may reverse the derivation to show it solves \eqref{eq:un:poissonstrong}.
\end{itemize}

For the linear case, where functions $a$ and $f$ are independent of $u$, the theory is relatively complete.  If $\partial_D \Omega$ has positive measure then it is a theorem that a solution to \eqref{eq:un:poissonweak} exists and is unique (\emph{well-posedness}; \citep{Ciarlet2002,Evans2010}).  Furthermore there exist conditions on the domain (e.g.~it is a convex polygon) and the boundary data, from which one can show that $u$ solving \eqref{eq:un:poissonweak} is in $C^2(\Omega) \cap C(\overline \Omega)$ (\emph{regularity}; \citep{Evans2010}).  Such theoretical matters are, however, beyond our scope.

In nonlinear cases each particular problem must be analyzed for well-posedness.  In terms of practical computation, however, some nonlinear cases are covered by our method and are solvable using our code.  These include the 2D versions of
\begin{itemize}
\item the \emph{Bratu} equation wherein $a\equiv 1$ and $f=\lambda e^u$, and
\item ``uniformized'' \emph{porous medium} equations \citep{Ockendonetal2003}, in which, for example, $a=u^{m-1}+\eps$ for some $m\ge 1$ and $\eps>0$.
\end{itemize}
For the former see Exercise \ref{chap:un}.\ref{exer:un:bratu}.\sidenote{Exercise \ref{chap:nl}.\ref{exer:nl:bratu} solves the 1D version.}  We will test our code on an easy case of the latter, namely with $a(u)=1+u^2$; see the exact solution \eqref{eq:un:exactsolution} in ``case 1,'' later in this Chapter.


\section{A $P^1$ finite element method}

Our FEM for the Poisson problem requires weak formulation \eqref{eq:un:poissonweak} to be true for $u$ from a finite-dimensional subspace of $H_{g}^1(\Omega)$, and for all test functions $v$ in a finite-dimensional subspace of $H_{0}^1(\Omega)$.  In the \emph{Galerkin} method here, these subspaces, both built from a triangulation of $\Omega$, are nearly the same.

We require $\Omega$ to be polygonal, so that $\partial\Omega$ is a closed polygon (Figure \ref{fig:un:polygon}), and we require that the segments forming $\partial\Omega$ are either entirely in $\partial_D\Omega$ or entirely in $\partial_N\Omega$.  Furthermore we assume $\partial_D\Omega$ is closed and contains at least one segment of positive length, leading to well-posedness of the continuum problem, at least in the linear case.

\begin{marginfigure}
\input{tmp/blob.tikz}
\caption{A polygonal domain $\Omega$.  The Dirichlet portion of the boundary $\partial_D\Omega$ is shown in bold.}
\label{fig:un:polygon}
\end{marginfigure}

By definition, a \emph{triangulation} is a finite set of non-overlapping, non-empty open triangles $\triangle_k$ which tile $\Omega$:
\begin{equation}
\Th = \left\{\triangle_k \,\Big|\, \cup_k \overline{\triangle}_k = \overline{\Omega} \, \text{ and} \,\, \triangle_k \cap \triangle_l = \emptyset \, \text{ if } k\ne l\right\}. \label{eq:un:triangulation}
\end{equation}
Let $K$ be the number of triangles, which are indexed $k=0,\dots,K-1$.  The $N$ nodes are located at $\bx_i=(x_i,y_i)$ for $i=0,\dots,N-1$.  (In contrast to many references, such as \citet{Elmanetal2005} which we follow in many ways, numbering here is zero-based, suitable for a C implementation.)  An example triangulation $\Th$ is shown in Figure \ref{fig:un:number-mesh}.  The subscript ``$h$'' in ``$\Th$'' denotes the typical or maximum size $h$ (e.g.~diameter) of the triangles, for now just a reminder of our desired limit $h\to 0$.  

We will use $P^1$ finite elements \citep{Elmanetal2005}, and thus the functions in use are continuous on the whole of $\overline\Omega$ while being linear on each triangle $\triangle_k$.  Such functions have three degrees of freedom on each $\triangle_k$, which one may think of as the coefficients in the linear formula $a + b x + c y$.  In fact there are more convenient local bases than $\{1,x,y\}$, as follows.  For each node $j$ there is a basis (``hat'') function  $\psi_j(x,y)$ which is linear on each triangle, continuous on all of $\overline{\Omega}$, and equal to one at only one node $j$,
\begin{equation*}
\psi_j(\bx_i) = \delta_{ij}.
\end{equation*}
See Figure \ref{fig:un:hatfunction} and compare Figure \ref{fig:of:q1hat}.  Functions $\psi_j$ are in $H^1(\Omega)$ \citep{Braess2007}, with piecewise-constant partial derivatives $\partial\psi_j/\partial x$ and $\partial\psi_j/\partial y$.  The set $\{\psi_j\}_{j=0,\dots,N-1}$ is linearly-independent.

\begin{marginfigure}
\input{tmp/blob.1.elenum.tikz}

\medskip

\input{tmp/blob.1.nodenum.tikz}
\caption{A triangulation of the polygon in Figure \ref{fig:un:polygon}, with element (top) and node (bottom) numbering.  There are $K=15$ elements, $N=13$ nodes, and $n_D=4$ nodes in $\partial_D\Omega$.}
\label{fig:un:number-mesh}
\end{marginfigure}

Hat functions allow us to interpolate and extend the Dirichlet data $g_D$ from its nodal values on $\partial_D \Omega$ to the whole of $\overline\Omega$.  We number the $n_D$ nodes which are in the Dirichlet boundary by $\bx_{j_l} \in \partial_D\Omega$ for $l=0,\dots,n_D-1$.  Let $\hat g_D \in C(\overline\Omega)$ be the piecewise-linear interpolant of $g_D$ which has value zero at all the nodes $\bx_j$ which are not in the Dirichlet boundary $\partial_D \Omega$:
\begin{equation}
\hat g_D(x,y) = \sum_{l=0}^{n_D-1} g_D(\bx_{j_l}) \psi_{j_l}(x,y). \label{eq:un:hatgdefine}
\end{equation}

Using \eqref{eq:un:hatgdefine} we can specify the finite-dimensional subspaces used in the method in terms of $\hat g_D$ and the span of certain basis functions $\psi_j$.  First, the \emph{test functions} are zero along $\partial_D \Omega$,
\begin{align*}
S_{0}^h &= \Span\left<\psi_j \,:\, \bx_j \notin \partial_D \Omega\right> = \Span\left<\psi_j \,:\, j \neq j_l\right>.
\end{align*}
Second, the \emph{trial functions} have value $g_D$ along $\partial_D \Omega$,
\begin{equation}
S_{g}^h = \left\{\hat g_D + w \,:\, w \in S_{0}^h\right\}.
\end{equation}
Note that $S_{0}^h$ is a linear subspace of $H_{0}^1(\Omega)$ while $S_{g}^h$ is merely an affine subspace of $H^1(\Omega)$.  Also
\begin{equation}
\dim(S_{0}^h)=\dim(S_{g}^h)=N-n_D.
\end{equation}

\begin{marginfigure}
\input{tikz/hatfunction.tex}
\caption{Hat functions $\psi_j$.}
\label{fig:un:hatfunction}
\end{marginfigure}

The FEM itself can now be stated.  It requires \eqref{eq:un:poissonweak} to be true for some $u^h\in S_{g}^h$ and for all $v^h\in S_{0}^h$.  By linearity it suffices to consider only a basis of test functions $S_{0}^h$, so we require
\begin{equation}<commit_msg>more text improvements in Chapter 7<commit_after>
In this Chapter we implement a finite element method (FEM) for an unstructured triangular mesh on a general polygonal domain.  Our code solves certain linear and nonlinear Poisson equation problems with reasonably-general boundary conditions.  The approach is to
\begin{center}
\emph{first do element-wise assembly of the residual equations,}
\end{center}
and thereby get a functioning \pSNES-based code.  We will not even think about matrices in writing the initial code.  Only after making sure it works do we
\begin{center}
\emph{write additional code to assemble an approximate Jacobian matrix.}
\end{center}

We do the following new tasks in this Chapter:
\begin{itemize}
\item read a triangular mesh from a file into \PETSc \pVecs and \pISs,
\item index that mesh unstructured way,
\item pre-allocate a \PETSc \pMat, and
\item implement Neumann boundary conditions.
\end{itemize}

Similarly to the example in the last Chapter, the residual equations $\bF(\bu)=0$, which is to say the nonlinear system seen by the \pSNES solver, is the discretized weak form of the PDE.  Our direct construction of these equations roughly follows the residual implementations of Chapters \ref{chap:nl}--\ref{chap:of}.

The FEM ``stiffness'' matrix $A$, a major piece in typical FEM introductions \citep{Braess2007,Elmanetal2005}, arises here as the Jacobian when solving $\bF(\bu)=0$.  In our first runs, done before writing any matrix-assembly code, the \pSNES solver approximates this linear system through finite-differencing evaluations of $\bF$.  Eventually we do assemble a Jacobian, but one that is only correct for the linear case.  For the smooth nonlinear case tested here, however, this inaccurate linearization suffices.

The example here contrasts with the rectangular-element ($Q^1$) FEM of the last Chapter by not using a \pDMDA and a structured grid.  Instead we implement a minimal, though naive, mesh-topology infrastructure.  (Doing so substantially increases our workload!)  The relatively-simple \Triangle program generates ASCII files describing the mesh.  An ``index set'' \pIS  type holds the global indices of nodes and triangles (elements).  \pVec and \pMat objects, used for the solution and Jacobian respectively, are indexed through the \pIS values.

In contrast with both earlier and later Chapters, we generate a code that only works on one MPI process (i.e.~a serial code).  However, at the end of the Chapter we briefly consider, though we do not implement, the additional constructs needed to distribute a mesh across many processes.  These considerations motivate the \pDMPlex example in Chapter \ref{chap:dp}.  In fact, the benefits of \PETSc's mesh-topology \pDM tools, including \pDMDA in previous Chapters and \pDMPlex in Chapter \ref{chap:dp}, become clearer by forgoing them here.


\section{A 2D nonlinear Poisson problem}

\begin{marginfigure}
\input{tikz/generalpoissondomain.tex}
\caption{Problem \eqref{eq:un:poissonstrong} on a domain.}
\label{fig:un:generalpoissondomain}
\end{marginfigure}

Let $\Omega \subset \RR^2$ be a bounded (open) region as in Figure \ref{fig:un:generalpoissondomain}.  We suppose its boundary $\partial\Omega$ is well-behaved, for example polygonal or Lipschitz-continuous \citep[section 1.2]{Ciarlet2002}.  Also we decompose the boundary into (measurable) disjoint subsets $\partial_D \Omega$ and $\partial_N \Omega$.

Chapter \ref{chap:st} solves the Poisson equation $- \grad^2 u = f$.  Here we allow a more general nonlinear form.  Let $a(u,x,y)$ and $f(u,x,y)$ be continuous functions, and assume there is $\eps$ so that
\begin{equation}
a(u,x,y) \ge \eps > 0, \label{eq:un:uniformelliptic}
\end{equation}
that is, \emph{uniform ellipticity} \citep{Evans2010}.  Boundary conditions are also needed, and here we consider nonhomogeneous Dirichlet and Neumann type.  Thus the  nonlinear Poisson \emph{strong form} we solve is to find $u(x,y)$ so that
\begin{align}
- \Div \left(a(u,x,y) \grad u\right) &= f(u,x,y) \quad \text{ on } \Omega, \label{eq:un:poissonstrong} \\
u &= g_D \quad \text{ on } \partial_D \Omega, \notag \\
\frac{\partial u}{\partial n} &= g_N \quad \text{ on } \partial_N \Omega. \notag
\end{align}
By definition, $\partial u/\partial n = \bn \cdot \grad u$ where $\bn$ is the outward unit normal on $\partial \Omega$ (Figure \ref{fig:un:generalpoissondomain}).

The data of problem \eqref{eq:un:poissonstrong}, besides the region $\Omega$ and its boundary, includes a \emph{diffusion coefficient} $a$, a \emph{source term} $f$, \emph{Dirichlet data} $g_D$, and \emph{Neumann data} $g_N$.  For the purpose of numerical solutions we will assume these data are continuous.  Poisson problem \eqref{poissonsquare} in Chapter \ref{chap:st} is the homogeneous Dirichlet case where $\Omega$ is a square, $a\equiv 1$, $f=f(x,y)$, $\partial_D \Omega = \partial \Omega$, $\partial_N \Omega = \emptyset$, and $g_D=0$.

One could allow more general boundary conditions in \eqref{eq:un:poissonstrong}, in particular \emph{Robin} conditions $\alpha u + \beta \frac{\partial u}{\partial n} = \gamma$ where $\alpha,\beta,\gamma$ could vary along the boundary \citep{Elmanetal2005}; see Exercise \ref{chap:un}.\ref{exer:un:robin}.  One could also allow $a$ and/or $f$ to depend on the gradient of $u$; recall $a=|\grad u|^{p-2}$ in the $p$-Laplacian equation of Chapter \ref{chap:of}, for example.

As \eqref{eq:un:poissonstrong} is stated there may be no solution where ``$\Div(a\grad u)$'' makes sense as a continuous function, even for polygonal regions and continuous data.  That is, there may be no $u\in C^2(\Omega) \cap C(\overline \Omega)$ so that \eqref{eq:un:poissonstrong} applies at all points.  In the linear case, at least, a solution exists if we convert \eqref{eq:un:poissonstrong} to a \emph{weak form}.  Such a weak formulation already arose in Chapter \ref{chap:of} as the gradient of an objective function, but the weak form here generally does not have a minimization origin; see Exercise \ref{chap:un}.\ref{exer:un:notminimization}.  In this Chapter we will derive the weak form from the strong form in the traditional way, by multiplying by a test function and integrating.


\section{Weak form with general boundary values}

Recall $L^2(\Omega)$ is the space of square-integrable real functions on $\Omega$ and that $W^{1,2}(\Omega)$ is the Sobolev space in which the gradient is also square integrable.\sidenote{See definition \eqref{eq:of:sobolevdefn}.}  We use two subsets of $W^{1,2}(\Omega)$: \emph{trial functions} from $W^{1,2}_g(\Omega)$, the functions with value $g_D$ along $\partial_D \Omega$, and \emph{test functions} from $W^{1,2}_0(\Omega)$, with value $0$ along $\partial_D \Omega$.

Now choose any $v\in W^{1,2}_0(\Omega)$, multiply the first equation in \eqref{eq:un:poissonstrong} by $v$, and integrate by parts:
\begin{equation*}
\int_\Omega a(u) \grad u \cdot \grad v - \int_{\partial\Omega} \frac{\partial u}{\partial n} v = \int_\Omega f(u) v.
\end{equation*}
(In writing such integrals we will generally show dependence on the solution $u$ but suppress it for $x,y$ so that $a(u,x,y)=a(u)$, and similarly for $f$.)  Next we use the boundary information, that $v=0$ on $\partial_D\Omega$ and $\partial u/\partial n=g_N$ on $\partial_N\Omega$:
\begin{equation}
\int_\Omega a(u) \grad u \cdot \grad v = \int_\Omega f(u) v + \int_{\partial_N\Omega} g_N v\quad \text{ for any } v\in W^{1,2}_0(\Omega). \label{eq:un:poissonweak}
\end{equation}

Equation \eqref{eq:un:poissonweak} is the \emph{weak formulation} of \eqref{eq:un:poissonstrong}.  Any $u \in W^{1,2}_g(\Omega)$ satisfying \eqref{eq:un:poissonweak} is a \emph{weak solution}.  Observe that the Dirichlet boundary conditions $g_D$ are incorporated into defining $W^{1,2}_g(\Omega)$ while the Neumann boundary data $g_N$ appears in the weak form \eqref{eq:un:poissonweak}.

Here are the ``rules'' for passing between the strong \eqref{eq:un:poissonstrong} and weak \eqref{eq:un:poissonweak} formulations:\begin{itemize}
\item A well-behaved function $u \in C^2(\Omega) \cap C(\overline \Omega)$ which satisfies \eqref{eq:un:poissonstrong} also solves \eqref{eq:un:poissonweak}.  The derivation above shows this.
\item If $u \in W^{1,2}_g(\Omega)$ solves \eqref{eq:un:poissonweak} then we accept it, by definition, as a solution.   If it is also well-behaved, $u \in C^2(\Omega) \cap C(\overline \Omega)$, then we may reverse the derivation to show it solves \eqref{eq:un:poissonstrong}.
\end{itemize}

For the linear case, where functions $a$ and $f$ are independent of $u$, the theory of \eqref{eq:un:poissonweak} is relatively complete.  If $\partial_D \Omega$ has positive measure then a solution to \eqref{eq:un:poissonweak} exists and is unique (\emph{well-posedness}; \citep{Ciarlet2002,Evans2010}).  Furthermore there exist conditions on the domain (e.g.~it is a convex polygon) and the boundary data, from which one can show that $u$ solving \eqref{eq:un:poissonweak} is in $C^2(\Omega) \cap C(\overline \Omega)$ (\emph{regularity}; \citep{Evans2010}).  Such theoretical matters are, however, beyond our scope.

In nonlinear cases each particular problem must be analyzed for well-posedness.  In terms of practical computation, however, some important nonlinear cases are covered by our method and are solvable using the code we construct in this Chapter.  These include the 2D versions of
\begin{itemize}
\item the \emph{Bratu} equation wherein $a\equiv 1$ and $f=\lambda e^u$, and
\item ``uniformized'' \emph{porous medium} equations \citep{Ockendonetal2003}, in which, for example, $a=u^{m-1}+\eps$ for some $m\ge 1$ and $\eps>0$.
\end{itemize}
For the former see Exercise \ref{chap:un}.\ref{exer:un:bratu}.\sidenote{Exercise \ref{chap:nl}.\ref{exer:nl:bratu} solves the 1D version.}  Later in the Chapter we test our code on an easy case of the porous medium equation, namely with $a(u)=1+u^2$.


\section{A $P^1$ finite element method}

Our FEM for the Poisson problem seeks $u$, from a finite-dimensional subspace of $W^{1,2}_g(\Omega)$, for which \eqref{eq:un:poissonweak} holds for all test functions $v$ in a finite-dimensional subspace of $W^{1,2}_0(\Omega)$.  In the \emph{Galerkin} method here, these subspaces, both built from a triangulation of $\Omega$, are nearly the same.

We require $\Omega$ to be polygonal, so that $\partial\Omega$ is a closed polygon (Figure \ref{fig:un:polygon}), and we require that the segments forming $\partial\Omega$ are either entirely in $\partial_D\Omega$ or entirely in $\partial_N\Omega$.  Furthermore we assume $\partial_D\Omega$ is closed and contains at least one segment of positive length, leading to well-posedness of the continuum problem, at least in the linear case.

\begin{marginfigure}
\input{tmp/blob.tikz}
\caption{A polygonal domain $\Omega$.  The Dirichlet portion of the boundary $\partial_D\Omega$ is shown in bold.}
\label{fig:un:polygon}
\end{marginfigure}

By definition, a \emph{triangulation} is a finite set of non-overlapping, non-empty open triangles $\triangle_k$ which tile $\Omega$:
\begin{equation}
\Th = \left\{\triangle_k \,\Big|\, \cup_k \overline{\triangle}_k = \overline{\Omega} \, \text{ and} \,\, \triangle_k \cap \triangle_l = \emptyset \, \text{ if } k\ne l\right\}. \label{eq:un:triangulation}
\end{equation}
Let $K$ be the number of triangles, which are indexed $k=0,\dots,K-1$.  The $N$ nodes are located at $\bx_i=(x_i,y_i)$ for $i=0,\dots,N-1$.  (In contrast to many references, such as \citet{Elmanetal2005} which we follow in many ways, numbering here is zero-based, suitable for a C implementation.)  An example triangulation $\Th$ is shown in Figure \ref{fig:un:number-mesh}.  The subscript ``$h$'' in ``$\Th$'' denotes the typical or maximum size $h$ (e.g.~diameter) of the triangles, for now just a reminder of our desired limit $h\to 0$.  

We will use $P^1$ finite elements \citep{Elmanetal2005}, and thus the functions in use are continuous on the whole of $\overline\Omega$ while being linear on each triangle $\triangle_k$.  Such functions have three degrees of freedom on each $\triangle_k$, which one may think of as the coefficients in the linear formula $a + b x + c y$.  In fact there are more convenient local bases than $\{1,x,y\}$, as follows.  For each node $j$ there is a basis (``hat'') function  $\psi_j(x,y)$ which is linear on each triangle, continuous on all of $\overline{\Omega}$, and equal to one at only one node $j$,
\begin{equation*}
\psi_j(\bx_i) = \delta_{ij}.
\end{equation*}
See Figure \ref{fig:un:hatfunction} and compare Figure \ref{fig:of:q1hat}.  Functions $\psi_j$ are in $W^{1,2}(\Omega)$ \citep{Braess2007}, with piecewise-constant partial derivatives $\partial\psi_j/\partial x$ and $\partial\psi_j/\partial y$.  The set $\{\psi_j\}_{j=0,\dots,N-1}$ is linearly-independent.

\begin{marginfigure}
\input{tmp/blob.1.elenum.tikz}

\medskip

\input{tmp/blob.1.nodenum.tikz}
\caption{A triangulation of the polygon in Figure \ref{fig:un:polygon}, with element (top) and node (bottom) numbering.  There are $K=15$ elements, $N=13$ nodes, and $n_D=4$ nodes in $\partial_D\Omega$.}
\label{fig:un:number-mesh}
\end{marginfigure}

Hat functions allow us to interpolate and extend the Dirichlet data $g_D$ from its nodal values on $\partial_D \Omega$ to the whole of $\overline\Omega$.  We number the $n_D$ nodes which are in the Dirichlet boundary by $\bx_{j_l} \in \partial_D\Omega$ for $l=0,\dots,n_D-1$.  Let $\hat g_D \in C(\overline\Omega)$ be the piecewise-linear interpolant of $g_D$ which has value zero at all the nodes $\bx_j$ which are not in the Dirichlet boundary $\partial_D \Omega$:
\begin{equation}
\hat g_D(x,y) = \sum_{l=0}^{n_D-1} g_D(\bx_{j_l}) \psi_{j_l}(x,y). \label{eq:un:hatgdefine}
\end{equation}

Using \eqref{eq:un:hatgdefine} we can specify the finite-dimensional subspaces used in the method in terms of $\hat g_D$ and the span of certain basis functions $\psi_j$.  First, the \emph{test functions} are zero along $\partial_D \Omega$,
\begin{align*}
S_{0}^h &= \Span\left<\psi_j \,:\, \bx_j \notin \partial_D \Omega\right> = \Span\left<\psi_j \,:\, j \neq j_l\right>.
\end{align*}
Second, the \emph{trial functions} have value $g_D$ along $\partial_D \Omega$,
\begin{equation}
S_{g}^h = \left\{\hat g_D + w \,:\, w \in S_{0}^h\right\}.
\end{equation}
Note that $S_{0}^h$ is a linear subspace of $W^{1,2}_0(\Omega)$ while $S_{g}^h$ is merely an affine subspace of $W^{1,2}(\Omega)$.  Also
\begin{equation}
\dim(S_{0}^h)=\dim(S_{g}^h)=N-n_D.
\end{equation}

\begin{marginfigure}
\input{tikz/hatfunction.tex}
\caption{Hat functions $\psi_j$.}
\label{fig:un:hatfunction}
\end{marginfigure}

The FEM itself can now be stated.  It requires \eqref{eq:un:poissonweak} to be true for some $u^h\in S_{g}^h$ and for all $v^h\in S_{0}^h$.  By linearity it suffices to consider only a basis of test functions $S_{0}^h$, so we require
\begin{equation}
{'fraction_non_alphanumeric': 0.12605237762470467, 'fraction_numerical': 0.008984725965858042, 'mean_word_length': 4.579650946899369, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '2821533', 'n_tokens_mistral': 9721, 'n_tokens_neox': 9301, 'n_words': 3832}