\section*{Finite Element Methods}

The most general numerical techniques for solving PDEs are known as
{\em finite element methods} (FEM).
This section describes the structure
of the equations that result upon application of FEM.
Readers wishing a proper treatment of FEM should
consult standard texts such as Becker, Carey, and Oden~\cite{becker81}
for an introduction, and Strang and Fix~\cite{strang73} for mathematical
analysis.

\begin{figure}
\begin{minipage}{3.45in}
\epsfxsize=3.4in \epsfbox{guitar2.pslg.eps}
\caption{\label{geometry}
Electric guitar.}
\end{minipage}
\begin{minipage}{3.45in}
\epsfxsize=3.4in \epsfbox{guitar2.1.ele.eps}
\caption{\label{mesh}
Finite element mesh.}
\end{minipage}
\end{figure}

Consider a heat conduction problem posed on the
two-dimensional domain of Figure \ref{geometry},
an electric guitar.  The problem is
to find the steady state temperature $u(x, y)$
of the guitar, given that the guitar is exposed to specified heat sources.
The physical behavior of this system is modeled
by a partial differential equation\footnote{
The steady-state heat equation is
$- \nabla\cdot\left[k(x,y)\nabla u(x,y)\right] = f(x,y)$, subject to
boundary conditions, where $f(x,y)$ represents internal heat sources,
and $k(x,y)$ is the thermal conductivity of the material.}.
We replace
this continuous problem
with a discrete approximation. Finite element
methods achieve this goal by partitioning the domain into
convex polygons or polyhedra, called {\em finite elements}.
Elements are
typically triangles or rectangles in two dimensions, and tetrahedra or
rectangular blocks in three dimensions.
The simplest
element is the linear triangle, which possesses three nodes, one
at each corner.  A {\em finite element mesh} is a collection of elements
(and nodes) covering a domain.
Figure \ref{mesh} shows a mesh of linear triangles covering
the guitar.

By considering the effects of the differential equation over each element,
we construct an approximate system of linear
equations of the form
\begin{equation}
\label{kuf}
{\bf K}{\bf u} = {\bf f}
\end{equation}
where ${\bf u} = [u_1, u_2, \ldots , u_n]^T$ denotes the temperatures at
each of the $n$ mesh nodes, the {\em force vector} ${\bf f}$ denotes the
heat sources at each mesh node, and the $n \times n$ {\em stiffness
matrix}\footnote{The
terminology for ${\bf K}$ and ${\bf f}$ comes from structural mechanics,
and has become widespread.}
${\bf K}$ describes the interactions among the nodes due to the principle
of conservation of energy, which governs heat conduction.
Once we have solved for ${\bf u}$, we can interpolate between
the nodal values to produce an approximate continuous solution
$\tilde{u}(x,y)$, which is piecewise linear on the mesh.

\begin{figure}
\centerline{\epsfbox{meshmatrix2.ps}}
\caption{\label{meshmatrix} A finite element mesh and its associated
stiffness matrix ${\bf K}$, finite element graph $G$, and element adjacency
graph $E$.
X's indicate nonzero matrix entries.}
\end{figure}

The structure of ${\bf K}$ is illustrated in
Figure \ref{meshmatrix}, which shows a simple mesh composed
of eight elements and nine nodes.  $K_{ij}$ is 
nonzero only if nodes $i$ and $j$ share some common element.  Hence,
$K_{24}$ is nonzero, $K_{37}$ is zero, and $K_{ii}$ is nonzero for any $i$.
For large meshes, ${\bf K}$ is very sparse.
The stiffness matrix ${\bf K}$ is computed by summing the action of
the differential equation over each individual element.  Hence, $K_{24}$
receives contributions from elements {\tt a} and {\tt b},
and $K_{66}$ from elements {\tt f} and {\tt g}.
The structure of ${\bf f}$ is composed similarly; for example,
$f_6$ receives contributions from elements {\tt f} and {\tt g}.

Given a finite element mesh, the associated {\em finite
element graph} $G$ is defined such that its vertices are the nodes of
the mesh, and there is an edge between any two vertices that belong to
a common element (see Figure \ref{meshmatrix}).  $G$ is isomorphic
to ${\bf K}$,
in the sense that $K_{uv} = 0$ if $(u, v) \not\in G$.
The associated {\em element adjacency graph} $E$ is defined such that its
vertices are the elements of the mesh, and there is an edge between any
two vertices whose corresponding elements share a common node.

Solving the linear equation (\ref{kuf}) accounts for most of the computation
and interprocessor communication of solving the PDE.  The linear system may
be solved either by a direct method such as Gaussian elimination~\cite{heath91}
or by an iterative method such as conjugate gradients~\cite{shewchuk94}.
In this
article, we consider only iterative methods, which are much easier to
parallelize than direct methods.

Most iterative methods have as their central
primitive the problem of computing the product of a large sparse matrix
${\bf K}$ and a vector ${\bf x}$.
For example, the iterative step of the conjugate
gradient method consists of several
vector operations and a single sparse matrix-vector multiplication.
Note that over the course of a simulation
the nonzero structure of the sparse matrix ${\bf K}$ remains fixed
(isomorphic to the finite element graph $G$), even though
the values of the nonzero entries of ${\bf K}$ may change over time (for
instance, in nonlinear problems).

The computation ${\bf y = Kx}$ involves the transfer of information along the
edges of $G$, because the value of ${\bf y}$
at node $i$ depends on the values of
${\bf x}$ at every node that is a neighbor of node $i$.  Nodes do not interact
directly with physically distant nodes; we will show how to
take advantage of this
locality to parallelize the computation.
