\section{Finite element spaces and finite element discretization}%
\label{S:FEM}

\def\vI{\vec{I}}

In the sequel we assume that $\Omega \subset \R^d$ is a bounded
domain triangulated by $\tria$, i.e.
\[
\bar \Omega \; = \; \bigcup_{S \in \tria} S.
\]

The following considerations are also valid for a triangulation of an
immersed surface (with $n > d$). In this situation one has to exchange
derivatives (those with respect to $x$) by \emph{tangential}
derivatives (tangential to the actual element, derivatives are always
taken element--wise) and the determinant of the parameterization's
Jacobian has to be replaced by Gram's determinant of the
parameterization. But for the sake of clearness and simplicity we
restrict our considerations to the case $n = d$.

The values of a finite element function or the values of its
derivatives are uniquely defined by the values of its DOFs and the
values of the basis functions or the derivatives of the basis
functions connected with these DOFs. Usually, evaluation of those
values is performed element--wise. On a single element the value of a
finite element function at a point $x$ in this element is
determined by the DOFs associated with this specific element and the
values of the non vanishing basis functions at this point.

We follow the concept of finite elements which are  given on a
single element $S$ in local coordinates. We distinguish two
classes of finite elements: 

Finite element functions on an element $S$ defined by a finite
dimensional function space $\Pbar$ on a \emph{reference element}
$\Sbar$ and the (one to one) mapping $\lS$ from the reference element
$\Sbar$ to $S$. For this class the dependency on the actual element
$S$ is fully described by the mapping $\lS$. For example, all Lagrange
finite elements belong to this class.

Secondly, finite element functions depending on the actual element
$S$. Hence, the basis functions are not fully described by $\Pbar$ and
the one to one mapping $\lS$. But using an initialization of the
actual element (which defines a finite dimensional function space
$\Pbar$ with information about the actual element), we can implement
this class in the same way as the first one. This class is needed for
Hermite finite elements which are not affine equivalent to the
reference element. Examples in 2d are the \emph{Hsieh--Clough--Tocher}
or \emph{HCT element} or the \emph{Argyris element} where only
the normal derivative at the midpoint of edges are used
in the definition of finite element functions; both
elements lead to functions which are globally $C^1(\bar\Omega)$. The
concrete implementation for this class in \ALBERTA is future work.

All in all, for a very general situation, we only need a vector of
basis functions (and their derivatives) on $\Sbar$ and a function which
connects each of these basis functions with its degree of freedom on
the element. For the second class, we additionally need an
initialization routine for the actual element. By such information,
every finite element function is uniquely described on every element
of the grid.

\subsection{Barycentric coordinates}\label{S:bary_coord}%
\idx{barycentric coordinates|(}

For describing finite elements on simplicial grids, it is very
convenient to use $d+1$ barycentric coordinates as a local coordinate
system on an element of the triangulation. Using $d+1$ 
local coordinates, the \emph{reference simplex} $\Sbar$ is a subset of 
a hyper surface in $\R^{d+1}$:
\begin{equation*}\label{S:Sbar}
\Sbar := \big\{ (\lambda_0,\dots,\lambda_d)
\in \R^{d+1}; \lambda_k \ge 0,\; \sum\limits_{k = 0}^d \lambda_k = 1 \big\}
\end{equation*}
On the other hand, for numerical integration on an 
element it is much more convenient to use the standard element
$\Shat \in \R^d$ defined in \secref{S:ref.coarse} as
\[
\Shat = \mbox{conv hull}
\left\{\hat a_0 = 0, \hat a_1 = e_1, \dots, \hat a_d = e_d\right\}
\]
where $e_i$ are the unit vectors in $\R^d$; using $\Shat$ for 
the numerical integration, we only have to compute the determinant of
the parameterization's Jacobian and not Gram's determinant.

The relation between a given simplex $S$, the reference simplex $\Sbar$,
and the standard simplex $\Shat$ is now described in detail.

Let $S$ be an element of the triangulation with vertices
$\{a_0,\dots,a_d\}$; let $F_S\colon \Shat \to S$ be the diffeomorphic
parameterization of $S$ over $\Shat$ with regular Jacobian $DF_S$, such
that 
\[
  F_S(\hat a_k) = a_k, \qquad k = 0,\dots,d
\]
holds. For a point $x \in S$ we set
\begin{equation*}\label{E:xhat}
\hat x = F_S^{-1}(x) \in \Shat.
\end{equation*}
For a simplex $S$ the easiest choice for $F_S$ is the unique affine
mapping \mathref{E:FS} defined on page \pageref{E:FS}. For an affine mapping,
$DF_S$ is constant.  In the following, we assume that the
parameterization $F_S$ of a simplex $S$ is affine.

For a simplex $S$ the barycentric coordinates
\[
\lS(x) \; = \; (\lS_0,\dots,\lS_d)(x) \in \R^{d+1}
\]
of some point $x \in \R^d$ are (uniquely) determined by the 
$(d+1)$ equations
\[
\sum\limits_{k = 0}^d \lS_k(x)\; a_k = x 
\qquad\mbox{and}\qquad
\sum\limits_{k = 0}^d \lS_k(x) = 1.
\]
The following relation holds:
\[
x \in S \qquad \mbox{\iff} \qquad \lS_k(x) \in [0,1]\ \ \mbox{for all }
k = 0,\dots,d \qquad \mbox{iff} \qquad \lS \in \Sbar.
\]
On the other hand, each $\lambda \in \Sbar$ defines a unique point 
$\xS \in S$ by
\[
\xS(\lambda) = \sum\limits_{k = 0}^d \lambda_k\, a_k.
\]
Thus, $\xS\colon \Sbar \to S$ is an invertible mapping with inverse
$\lS\colon S \to \Sbar$. The barycentric coordinates of $x$ on $S$
are the same as those of $\hat x$ on $\Shat$, i.e. $\lS(x) =
\lShat(\hat x)$.

In the general situation, when $F_S$ may not be affine, i.e.\ 
we have a parametric element, the barycentric coordinates $\lS$ are 
given by the inverse of the parameterization $F_S$ and
the barycentric coordinates on $\Shat$:
\begin{equation*}\label{E:lS}
\lS(x) = \lShat(\hat x) = \lShat\left(F_S^{-1}(x)\right)
\end{equation*}
and the \emph{world coordinates} of a point $\xS \in S$ with barycentric
coordinates $\lambda$ are given by
\begin{equation*}\label{E:xS}
\xS(\lambda) \; = \; F_S\left(\sum_{k = 0}^d \lambda_k \hat a_k\right)
 \; = \; F_S\left(\xShat(\lambda)\right)
\end{equation*}
(see also Figure~\ref{F:S_Shat_Sb}).

\begin{figure}[htbp]
\centerline{\includegraphics[scale=1.0]{EPS/bary}}
\caption{Definition of $\lS \colon S \to \Sbar$ via $F_S^{-1}$ and $\lShat$,
and $\xS\colon \Sbar \to S$ via $\xShat$ and $F_S$}\label{F:S_Shat_Sb}
\end{figure}

Every function $f \colon S \to V$ defines (uniquely) two functions
\[
\begin{array}{rcl}
\bar f \colon\, \Sbar &\to & V\\
        \lambda & \mapsto & f(\xS(\lambda))
\end{array}
\qquad\mbox{and}\qquad
\begin{array}{rcl}
\hat f \colon\, \Shat & \to &V\\
        \xhat & \mapsto & f(F_S(\xhat)).
\end{array}
\]
Accordingly, $\hat f \colon\, \Shat \to V$ defines two functions
$f \colon\, S \to V$ and $\bar f \colon\, \Sbar \to V$, and 
$\bar f \colon\, \Sbar \to V$ defines $f \colon\, S \to V$ and 
$\hat f \colon\, \Shat \to V$.

Assuming that a function space $\Pbar \subset C^0(\Sbar)$ is given,
it uniquely defines function spaces $\Phat$ and $\PS$ by
\begin{equation}\label{Pbar_PS}
\Phat = \left\{\phat \in C^0(\Shat);\, \pbar \in \Pbar\right\}
\qquad \mbox{and} \qquad
\PS = \left\{\vphi \in C^0(S);\, \pbar \in \Pbar\right\}.
\end{equation}
We can also assume that the function space $\Phat$ is given and this
space uniquely defines $\Pbar$ and $\P_S$ in the same manner. In
\ALBERTA, we use the function space $\Pbar$ on $\Sbar$; the
implementation of a basis $\left\{\pbar^1,\dots,\pbar^m\right\}$ of
$\Pbar$ is much simpler than the implementation of a basis
$\left\{\phat^1,\dots,\phat^m\right\}$ of $\Phat$ as we are able to
use symmetry properties of the barycentric coordinates $\lambda$.

In the following we shall often drop the superscript $S$ of $\lS$
and $\xS$.  The mappings $\lambda(x) = \lS(x)$ and $x(\lambda) =
\xS(\lambda)$ are always defined with respect to the actual element 
$S \in \tria$.
\idx{barycentric coordinates|)}

\subsection{Finite element spaces}%
\label{S:FES}%
\idx{finite element spaces|(}

\ALBERTA supports scalar and vector valued finite element functions.
The basis functions are always real valued; thus, the coefficients
of a finite element function in a representation by the basis functions
are either scalars or vectors of length $n$.

For a given function space $\Pbar$ and some given real valued function space
$C$ on $\Omega$, a finite element space $X_h$ is defined by
\begin{equation*}\label{E:Xh}
X_h = X_h(\tria, \Pbar, C) = 
\left\{\vphi \in C; \; \vphi_{|S} \in \P_S \mbox{ for all } S \in \tria\right\}
\end{equation*}
for scalar finite elements. For vector valued finite elements, $X_h$ is
given by
\begin{align*}\label{E:Xhv}
X_h &= X_h(\tria, \Pbar, C)
%\notag\\ & 
= \left\{
\vphi = (\vphi_1,\dots,\vphi_n) \in C^n;
\, \vphi_{i\,|S} \in \P_S \mbox{ for all } i=1,\dots,n,\, S \in \tria\right\}.
\end{align*}
The spaces $\P_S$ are defined by $\Pbar$ via \mathref{Pbar_PS}.

For conforming finite element discretizations, $C$ is the continuous
space $X$ (for 2nd order problems, $C = X = H^1(\Omega)$). For
non--conforming finite element discretizations, $C$ may control the
\emph{non conformity} of the ansatz space which has to be controlled in
order to obtain optimal error estimates (for example, in the
discretization of the incompressible Navier--Stokes equation by the
non--conforming Crouzeix--Raviart element, the finite element
functions are continuous only in the midpoints of the edges).

The abstract definition of finite element spaces as a triple (mesh,
local basis functions, and DOFs) now matches the mathematical
definition as $X_h = X_h(\tria,\Pbar,C)$ in the following way: The
mesh is the underlying triangulation $\tria$, the local basis
functions are given by the function space $\Pbar$, and together with
the relation between global and local degrees of freedom every finite
element function satisfies the global continuity constraint
$C$. This relation between global and local DOFs is now described in
detail.\idx{finite element spaces|)}

\subsection{Evaluation of finite element functions}%
\label{S:eval_fe}%
\idx{DOFs!relation global and local DOFs}%
\idx{evaluation of finite element functions}

Let $\left\{\pbar^1,\dots,\pbar^m\right\}$ be a basis of $\Pbar$ and
let $\left\{\vphi_1,\dots,\vphi_N\right\}$ be a basis of $X_h$, $N =
\dim X_h$, such that for every $S \in \tria$ and for all 
basis functions $\vphi_j$ which do not vanish on $S$
\[
\vphi_j|_S(x(\lambda)) = \pbar^i(\lambda) \qquad \mbox{for all } \lambda
\in \Sbar
\]
holds with some $i \in \{1,\dots,m\}$ depending on $j$ and $S$.  Thus, the
numbering of the basis functions in $\Pbar$ and the mapping $\xS$
induces a local numbering of the non vanishing basis functions on
$S$. Denoting by $J_S$ the index set of all non vanishing basis
functions on $S$, the mapping $i_S: J_S \to \{1,\dots,m\}$ is one to one
and uniquely connects the degrees of freedom of a finite element
function on $S$ with the local numbering of basis functions.
If $j_S: \{1,\dots,m\} \to J_S$ denotes the inverse mapping of $i_S$,
the connection between global and local basis functions 
is uniquely determined on each element $S$ by
\begin{subequations}\label{E:global-lokal}
\begin{alignat}{2}
\vphi_j(x(\lambda)) &= \pbar^{i_S(j)}(\lambda),
&\qquad & \mbox{for all } \lambda \in \Sbar,\,j \in J_S,
\label{E:global-lokal.a}\\
\vphi_{j_S(i)}(x(\lambda)) &= \pbar^{i} (\lambda), 
&\qquad & \mbox{for all } \lambda \in \Sbar,\,i \in \{1,\dots,m\}.
\label{E:global-lokal.b}
\end{alignat}
\end{subequations}

For $\uh \in X_h$ denote by $\left(u_1,\dots,u_N\right)$ the 
\emph{global coefficient vector} of the basis $\left\{\vphi_j\right\}$ with 
$u_j \in \R$ for scalar finite elements, $u_j \in \R^n$ for
vector valued finite elements, i.e.
\[
\uh(x) = \sum\limits_{j = 1}^N u_j \vphi_j(x) \qquad
\mbox{for all } x \in \Omega,
\]
and the \emph{local coefficient vector}
\begin{equation*}\label{E:local_coeff}
\left(u_S^1,\dots,u_S^m\right) = 
\left(u_{j_S(1)},\dots,u_{j_S(m)}\right)
\end{equation*}
of $\uh$ on $S$ with respect to the local numbering of the non
vanishing basis functions (local numbering is denoted by a superscript
index and global numbering is denoted by a subscript index). Using the
local coefficient vector and the local basis functions we obtain the
\emph{local representation} of $u_h$ on $S$:
\begin{equation*}\label{E:uh_on_Sx}
\uh(x) = \sum\limits_{i = 1}^m u^i_S\, \pbar^i(\lambda(x))
\qquad \mbox{for all } x \in S.
\end{equation*}
In finite element codes, finite element functions are \emph{not}
evaluated at \emph{world coordinates} $x$ as in \mathref{E:uh_on_Sx},
but they are evaluated on a single element $S$ at barycentric coordinates
$\lambda$ on $S$ giving the value at the world coordinates $x(\lambda)$:
\begin{equation*}\label{E:uh_on_S}
\uh(x(\lambda)) = \sum\limits_{i = 1}^m u^i_S\, \pbar^i(\lambda)
\qquad \mbox{for all } \lambda \in \Sbar.
\end{equation*}
The mapping $j_S$, which allows the access to the local
coefficient vector from the global one, is the relation between
local DOFs and global DOFs. Global DOFs for a finite element 
space are stored on the mesh elements at positions which are known
to the DOF administration of the finite element space. Thus,
the corresponding DOF administration will provide information
for the implementation of the function $j_S$ and therefore
information for reading/writing local coefficients from/to global 
coefficient vectors (compare Section~\ref{man:S:basfct_data}).

\subsubsection{Evaluating derivatives of finite element functions}%
\label{S:eval_Dfe}%
\idx{evaluation of derivatives}

Derivatives of finite element functions are also evaluated on single
elements $S$ in barycentric coordinates. In the calculation of
derivatives in terms of the basis functions $\pbar^i$, the Jacobian
$\Lambda = \Lambda_S \in \R^{d\times\code{DIM\_OF\_WORLD}}$ 
of the barycentric coordinates on $S$ is involved (we consider here 
only the case $d=\code{DIM\_OF\_WORLD}=n$):
\begin{equation*}\label{E:LambdaS}
\Lambda(x) =
\left(
\begin{matrix}
\lambda_{0,x_1}(x) &\lambda_{0,x_2}(x) & \cdots & \lambda_{0,x_n} (x)\\
\vdots      & \vdots     &        &  \vdots\\
\lambda_{d,x_1}(x) &\lambda_{d,x_2}(x) & \cdots & \lambda_{d,x_n} (x)\\
\end{matrix}
\right)
= 
\left(
\begin{matrix}
\mbox{--}&\nabla \lambda_0(x)^t&\mbox{--}\\
& \vdots &\\
\mbox{--}&\nabla \lambda_d(x)^t&\mbox{--}
\end{matrix}
\right).
\end{equation*}
Now, using the chain rule we obtain for every function $\vphi \in \PS$
\begin{equation*}\label{E:grad_phi}
\nabla \vphi(x) = \nabla \left(\pbar(\lambda(x))\right)
= \sum_{k = 0}^d  \pbar_{,\lambda_k}(\lambda(x)) \nabla \lambda_k(x)
= \Lambda^t(x) \nablal \pbar(\lambda(x)), \qquad x \in S
\end{equation*}
and
\begin{equation*}\label{E:D2_phi}
D^2 \vphi(x) = \Lambda^t(x) D^2_\lambda \pbar(\lambda(x)) \Lambda(x)
 + \sum_{k = 0}^d  D^2 \lambda_k(x) \, \pbar_{,\lambda_k}(\lambda(x)),
\qquad x \in S.
\end{equation*}
For a simplex $S$ with an affine parameterization $F_S$,
$\nabla \lambda_k$ is constant on $S$ and we get
\[
\nabla \vphi(x) = \Lambda^t \nablal \pbar(\lambda(x)) 
\quad \mbox{and} \quad 
D^2 \vphi(x) = \Lambda^t D^2_\lambda \pbar(\lambda(x)) \Lambda,
\qquad x \in S.
\]
Using these equations, we immediately obtain
\begin{equation*}\label{E:grd_uh_on_Sx}
\nabla u_h(x) = \Lambda^t(x) \sum_{i = 1}^m u_S^i \, \nablal 
\pbar^i(\lambda(x)), \qquad x \in S
\end{equation*}
and
\begin{equation*}\label{E:D2_uh_on_Sx}
D^2 u_h(x) = \Lambda^t(x) \sum_{i = 1}^m u_S^i \,
D^2_\lambda \pbar^i(\lambda(x))\Lambda(x)
 +  
\sum_{k = 0}^d  D^2 \lambda_k(x) \sum_{i = 1}^m u_S^i \,
 \pbar^i_{,\lambda_k} (\lambda(x)), \qquad x \in S.
\end{equation*}
Since the evaluation is actually done in barycentric coordinates,
this turns on $S$ into
\begin{equation*}\label{E:grd_uh_on_S}
\nabla u_h(x(\lambda)) = \Lambda^t(x(\lambda)) \sum_{i = 1}^m u_S^i \, \nablal 
\pbar^i(\lambda), \qquad \lambda \in \Sbar
\end{equation*}
and
\begin{equation*}\label{E:D2_uh_on_S}
D^2 u_h(x(\lambda)) = \Lambda^t(x(\lambda)) \sum_{i = 1}^m u_S^i \,
D^2_\lambda \pbar^i(\lambda)\Lambda(x(\lambda))
 + 
\sum_{k = 0}^d  D^2 \lambda_k(x(\lambda)) \sum_{i = 1}^m u_S^i \,
 \pbar^i_{,\lambda_k} (\lambda), \qquad \lambda \in \Sbar.
\end{equation*}

Once the values of the basis functions, its derivatives, and the local
coefficient vector $(u_S^1, \dots, u_S^{m})$ are known, the
evaluation of $u_h$ and its derivatives depends only on $\Lambda$
and can be performed by some general evaluation routine (compare
Section~\ref{man:S:eval}).

\subsection{Interpolation and restriction during refinement and coarsening}%
\label{S:inter_restrict}%
\idx{interpolation and restriction of DOF vectors|(}

We assume the following situation: Let $S$ be a (non--parametric)
simplex with children $S_0$ and $S_1$ generated by the bisection of
$S$ (compare \algorithmref{A:recursive_refine}). Let $X_S$,
$X_{S_0,S_1}$ be the finite element spaces restricted to $S$ and $S_0
\cup S_1$ respectively.

Throughout this section we denote by $\big\{\vphi^i\big\}_{i = 1, \dots, m}$
the basis of the coarse grid space $X_S$ and by
$\big\{\psi^j\}_{j = 1, \dots, k}$ the basis functions of
$X_{S_0 \cup S_1}$. For a function $\uh \in X_S$ we denote by 
$\vec{u}_\vphi = (u^1_\vphi, \dots, u^m_\vphi)^t$ 
the coefficient vector with respect to the basis $\big\{\vphi^i\big\}$
and for a function $\vh \in X_{S_0 \cup S_1}$ by
$\vec{v}_\psi = (v^1_\psi, \dots, v^k_\psi)^t$ the
coefficient vector with respect to $\big\{\psi^j\big\}$. 

We now derive equations for the transformation of local coefficient
vectors for finite element functions that are interpolated during
refinement and coarsening, and for vectors storing values of a
linear functional applied to the basis functions on the fine grid
which are restricted to the coarse functions during coarsening.

Let
\[
\I^S_{S_0,S_1} \colon\, X_S \to X_{S_0 \cup S_1}
\]
be an interpolation operator. For nested finite element spaces, i.e.
$X_S \subset X_{S_0 \cup S_1}$, every coarse grid function $\uh \in
X_S$ belongs also to $X_{S_0 \cup S_1}$, so the natural choice is
$\I^S_{S_0,S_1} = id$ on $X_S$ (for example, Lagrange finite elements
are nested). The interpolants $\I^S_{S_0,S_1} \vphi^i$ can be written
in terms of the fine grid basis functions
\[
\I^S_{S_0,S_1} \vphi^i = \sum_{j = 1}^k a_{ij} \psi^j
\]
defining the $(m \times k)$--matrix 
\begin{equation}\label{E:a_matrix}
A = (a_{ij})_{\atop{i = 1,\dots,m}{j = 1,\dots,k}}.
\end{equation}
This matrix $A$ is involved in the interpolation during refinement and the
transformation of a linear functional during coarsening.

For the interpolation of functions during coarsening, we need an
interpolation operator $\I^{S_0,S_1}_S \colon\, X_{S_0 \cup S_1} \to X_S$.
The interpolants $\I^{S_0,S_1}_S \psi^j$ of the fine grid
functions can now be represented by the coarse grid basis
\[
  \I^{S_0,S_1}_S \psi^j = \sum_{i = 1}^m b_{ij}\, \vphi^i
\]
defining the $(m \times k)$--matrix
\begin{equation}\label{E:b_matrix}
B = (b_{ij})_{\atop{i = 1,\dots,m}{j = 1,\dots,k}}.
\end{equation}
This matrix $B$ is used for the interpolation during coarsening.

Both matrices depend only on the set of local basis functions on
parent and children. Thus, they depend on the reference element
$\Sbar$ and one single bisection of the reference element into
$\Sbar_0$, $\Sbar_1$.  The matrices do depend on the local numbering
of the basis functions on the children with respect to the
parent. Thus, in 3d the matrices depend on the element type of $S$
also (for an element of type \code{0} the numbering of basis functions
on $\Sbar_1$ differs from the numbering on $\Sbar_1$ for an element of
type \code{1}, \code{2}). But all matrices can be calculated by the
local set of basis functions on the reference element.

DOFs can be shared by several elements, compare Section~\ref{S:DOFs1}.
Every DOF is connected to a basis function which has a support on all
elements sharing this DOF\@. Each DOF refers to one coefficient of a
finite element function, and this coefficient has to be calculated
only once during interpolation. During the restriction of a linear
functional, contributions from several basis functions are added to
the coefficient of another basis function. Here we have to control
that for two DOFs, both shared by several elements, the contribution
of the basis function at one DOF is only added once to the other DOF
and vice versa. This can only be done by performing the interpolation
and restriction on the whole refinement/coarsening patch at the same
time.

\subsubsection{Interpolation during refinement}%
\idx{refinement!interpolation of DOF vectors}%
\idx{interpolation of DOF vectors}

Let $\vec{u}_\vphi = (u^1_\vphi, \dots, u^m_\vphi)^t$ be the coefficient
vector of  a finite element function $\uh \in X_S$ with respect to 
$\big\{\vphi^i\big\}$,
and let $\vec{u}_\psi = (u^1_\psi, \dots, u^k_\psi)^t$ the
coefficient vector of $\I^S_{S_0,S_1} \uh$  with respect to
$\big\{\psi^j\big\}$. Using matrix $A$ defined in \mathref{E:a_matrix}
we conclude
\[
\I^S_{S_0,S_1} \uh = \sum_{i = 1}^m  u^i_\vphi\,\I^S_{S_0,S_1} \vphi^i 
=    \sum_{i = 1}^m u^i_\vphi \, \sum_{j = 1}^k a_{ij} \, \psi^j
  =  \sum_{j = 1}^k \left(A^t \vec{u}_\vphi \right)_j \psi^j,
\]
or equivalently
\[
\vec{u}_\psi  =  A^t \vec{u}_\vphi.
\]
A subroutine which interpolates a finite element function during
refinement is an efficient implementation of this matrix--vector
multiplication.

\subsubsection{Restriction during coarsening}%
\idx{coarsening!restriction of DOF vectors}%
\idx{restriction of DOF vectors}

In an (Euler, e.g.) discretization of a time dependent problem, the
term $(\uh^\mathrm{old}, \vphi_i)_{L^2(\Omega)}$ appears on the right
hand side of the discrete system, where $\uh^\mathrm{old}$ is the
solution from the last time step. Such an expression can be calculated
exactly, if the grid does not change from one time step to another.
Assuming that the finite element spaces are nested, it is also
possible to calculate this expression exactly when the grid was
refined, since $\uh^\mathrm{old}$ belongs to the fine grid finite
element space also. Usually, during coarsening information is lost,
since we can not represent $\uh^\mathrm{old}$ exactly on a coarser
grid. But we can calculate $(\uh^\mathrm{old}, \psi_i)_{L^2(\Omega)}$
exactly on the fine grid; using the representation of the
coarse grid basis functions $\vphi_i$ by the fine grid functions
$\psi_i$, we can transform data during coarsening such that
$(\uh^\mathrm{old}, \vphi_i)_{L^2(\Omega)}$ is calculated exactly
for the coarse grid functions too.

More general, assume that the finite element spaces are nested and
that we can evaluate a linear functional $f$ exactly for all basis
functions of the fine grid. Knowing the values 
$\vec{f}_\psi = (\langle f,\,\psi^1\rangle, \dots, 
\langle f,\,\psi^k\rangle)^t$
for the fine grid functions, we obtain with matrix $A$ 
from \mathref{E:a_matrix} for the values
$\vec{f}_\vphi = (\langle f,\,\vphi^1\rangle, \dots, \langle
f,\,\vphi^m\rangle)^t $ on the coarse grid
\[
\vec{f}_\vphi  =  A \vec{f}_\psi
\]
since
\[
  \langle f, \, \vphi^i \rangle 
 =  
  \langle f,\,\sum_{j = 1}^k a_{ij}\psi^j \rangle 
 =  
  \sum_{j = 1}^k a_{ij}  \langle f,\,\psi^j \rangle
\]
holds (here we used the fact, that $\I^S_{S_0,S_1} = id$ on $X_S$
since the spaces are nested).

Thus, given a functional $f$ which we can evaluate exactly
for all basis functions of a grid $\tilde \tria$ and its refinements,
we can also calculate $\langle f, \, \vphi^i \rangle$ exactly for all basis
functions $\vphi^i$ of a grid  $\tria$ obtained by 
refinement and coarsening of $\tilde \tria$ in the following way:
First refine all elements of the grid that have to be refined;
calculate $\langle f, \, \vphi \rangle$ for all basis functions $\vphi$
of this intermediate grid; in the last step coarsen all elements
that may be coarsened and restrict this vector during each coarsening
step as described above. 

In \ALBERTA the assemblage of the discrete system inside the adaptive
method can be split into three steps: one initialization step before
refinement, the second step between refinement and coarsening, and the
last, and usually most important, step after coarsening, when all grid
modifications are completed, see \secref{man:S:adapt_stat_in_ALBERTA}.
The second assemblage step can be used for an exact computation of a
functional $\langle f, \, \vphi \rangle$ as described above.


\subsubsection{Interpolation during coarsening}%
\idx{coarsening!interpolation of DOF vectors}%
\idx{interpolation of DOF vectors}

Finally, we can interpolate a finite element function
during coarsening. The matrix for transforming the coefficient vector
$\vec{u}_\psi = (u^1_\psi, \dots, u^k_\psi)^t$ of a fine grid function $\uh$
to the coefficient vector $\vec{u}_\vphi = (u^1_\vphi, \dots, u^m_\vphi)^t$
of the interpolant on the coarse grid is given by matrix $B$ defined
in \mathref{E:b_matrix}:
\begin{align*}
\I^{S_0,S_1}_S \uh & = \I^{S_0,S_1}_S  \sum_{j = 1}^k u^j_\psi\, \psi^j
 =  \sum_{j = 1}^k u^j_\psi\, \I^{S_0,S_1}_S \psi^j\\
& = \sum_{j = 1}^k u^j_\psi\, \sum_{i = 1}^m b_{ij}\, \vphi^i
 =  \sum_{i = 1}^m \left(\sum_{j = 1}^k b_{ij}\,  u^j_\psi\right) \vphi^i.
\end{align*}
Thus we have the following equation for the coefficient vector of 
the interpolant of $\uh$:
\[
\vec{u}_\vphi  =  B\, \vec{u}_\psi.
\]
In contrast to the interpolation during refinement and the above
described transformation of a linear functional, information is
lost during an interpolation to the coarser grid.

\begin{example}[Lagrange elements]
Lagrange finite elements are connected to Lagrange nodes $x^i$. For
linear elements, these nodes are the vertices of the triangulation,
and for quadratic elements, the vertices and the edge--midpoints.
The Lagrange basis functions $\big\{\vphi^i\big\}$ satisfy
\[
\vphi_i(x_j) = \delta_{ij} \qquad \mbox{for } i,j = 1,\dots, \dim X_h.
\]
Consider the situation of a simplex $S$ with children $S_0$, $S_1$.
Let $\big\{\vphi^i\big\}_{i = 1,\dots,m}$ be the Lagrange basis functions
of $X_S$ with Lagrange nodes $\big\{x^i_\vphi\big\}_{i = 1,\dots,m}$ on $S$
and $\big\{\psi^j\}_{j = 1,\dots,k}$ be the Lagrange basis functions of
$X_{S_0 \cup S_1}$ with Lagrange nodes
$\big\{x^j_\psi\}_{j = 1,\dots,k}$ on $S_0 \cup S_1$. The Matrix $A$ is then
given by
\[
 a_{ij} = \vphi^i(x^j_\psi), \qquad i = 1,\dots,m,\; j = 1,\dots,k
\]
and matrix $B$ is given by
\[
 b_{ij} = \psi^j(x^i_\vphi), \qquad i = 1,\dots,m,\; j = 1,\dots,k.
\]
\end{example}
\idx{interpolation and restriction of DOF vectors|)}

\subsection{Discretization of 2nd order problems}%
\label{S:Dis2ndorder}%
\idx{finite element discretization|(}

In this section we describe the assembling of the discrete
system in detail. 
We consider the following second order differential equation in
divergence form:
\begin{subequations}\label{E:strong}
\begin{alignat}{2}
L u := -\nabla \cdot A \nabla u + b \cdot \nabla u + c\, u 
       &= f \qquad & &\mbox{in } \Omega,\label{E:strong.a}\\
     u &= g        & &\mbox{on } \GD,\label{E:strong.b}\\
\nO \cdot A \nabla u & = 0 &&\mbox{on } \GN\label{E:strong.c},
\end{alignat}
\end{subequations}
where $A \in L^\infty(\Omega;\R^{n \times n})$, 
$b \in L^\infty(\Omega;\R^n)$, $c \in L^\infty(\Omega)$, and
$f \in L^2(\Omega)$.
By $\GD \subset \partial \Omega$ (with $|\GD| \ne 0$) we denote the
Dirichlet boundary and we assume that the Dirichlet boundary values
$g \colon \GD \to \R$ have an extension to some function
$g \in H^1(\Omega)$.

By $\GN = \partial \Omega \backslash \GD$ we denote the 
Neumann boundary, and by $\nO$ we denote the outer unit normal vector
on $\partial \Omega$. The boundary condition \mathref{E:strong.c}
is a so called \emph{natural} Neumann condition.

Equations \mathref{E:strong} describe not only a simple model problem.
The same kind of equations result from a linearization of nonlinear
elliptic problems (for example by a Newton method) as well as from a time
discretization scheme for \hbox{(non--)} linear parabolic problems.

Setting
\begin{equation*}\label{E:X-X0}
X =  H^{1}(\Omega) \qquad \mbox{and} \qquad
\Xc = \left\{v \in H^{1}(\Omega);\, v = 0 \mbox{ on } \GD\right\}
\end{equation*}
this equation has the following weak formulation: We are looking
for a solution $u \in X$, such that $u \in g+\Xc$ and
\begin{equation}\label{E:weak}
\int_\Omega (\nabla \varphi(x)) \cdot A(x) \nabla u(x) 
+ \varphi(x)\, b(x) \cdot  \nabla u(x)
+ c(x)\,\varphi(x) \, u(x)   \, dx = 
\int_\Omega f(x)\, \varphi(x)\, dx
\end{equation}
for all $\vphi \in \Xc$

Denoting by $\Xc^*$ the dual space of $\Xc$ we 
identify the differential operator $L$ with the linear operator
$L \in \Lin(X,\Xc^*)$ defined by
\begin{equation*}\label{E:Lu}
\ldual{L v}{\varphi}{\Xc^* \times \Xc} := 
\int_\Omega \nabla \varphi \cdot A \nabla v 
 +  \int_\Omega \varphi \, b \cdot  \nabla v
 +  \int_\Omega c\,\varphi \, v
\qquad \mbox{for all } v, \varphi \in \Xc
\end{equation*}
and the right hand side $f$ with
the linear functional $f \in \Xc^*$ defined by
\begin{equation*}
\ldual{F}{\varphi}{\Xc^* \times \Xc} :=
\int_\Omega f\, \varphi \qquad \mbox{for all } \varphi \in \Xc.
\end{equation*}

With these identifications we use the following reformulation of
\mathref{E:weak}: Find $u \in X$ such that
\begin{equation}\label{E:weak2}
u \in g + \Xc: \qquad L\, u  =  f \qquad\mbox{in }\Xc^*
\end{equation}
holds.

Suitable assumptions on the coefficients imply that $L$ is elliptic,
i.e.\ there is a constant $C = C_{A,b,c,\Omega}$ such that
\[
 \ldual{L\,\vphi}{\vphi}{\Xc^* \times \Xc}  \ge  C\,\lnorm{\vphi}{X}^2
 \qquad\mbox{for all }\vphi\in\Xc.
\]
The existence of a unique solution $u \in X$ of \mathref{E:weak2}
is then a direct consequence of the Lax--Milgram--Theorem.

We consider a finite dimensional subspace
$X_h \subset X$ for the discretization of \mathref{E:weak2} with 
$N = \dim\ X_h$. We set $\Xc_h = X_h \cap \Xc$ with
$\Nc = \dim\ \Xc_h$. Let $\gh \in X_h$ be an approximation of $g \in X$.
A discrete solution of \mathref{E:weak2} is then given by: Find
$\uh \in X_h$ such that
%
\begin{equation}\label{E:discrete}
\uh \in \gh + \Xc_h: \qquad L\, \uh  = f \qquad\mbox{in } \Xc_h^*,
\end{equation}
i.e.
\[
\uh \in \gh + \Xc_h: \qquad
\ldual{L \uh}{\ph}{\Xc_h^* \times \Xc_h}
 = 
\ldual{f}{\ph}{\Xc_h^* \times \Xc_h} \quad
\mbox{for all } \ph \in \Xc_h
\]
holds. If $L$ is elliptic, we have a unique discrete solution 
$\uh \in X_h$ of \mathref{E:discrete}, again using the 
Lax--Milgram--Theorem.

\def\vv{\vec{v}}
\def\vu{\vec{u}}
\def\vL{\vec{L}}
\def\vf{\vec{f}}

Choose a basis $\left\{\vphi_1,\dots,\vphi_{N\vphantom{\Nc}}\right\}$
of $X_h$ such that $\left\{\vphi_1,\dots,\vphi_{\Nc}\right\}$ is a
basis of $\Xc_h$.  For a function $\vh \in X_h$ we denote by $\vv =
(v_1,\dots,v_N)$ the coefficient vector of $\vh$ with respect to the
basis $\left\{\vphi_1,\dots,\vphi_{N\vphantom{\Nc}}\right\}$, i.e.
\[
\vh  =  \sum_{j = 1}^N v_j \vphi_j.
\]
Using \mathref{E:discrete} with test functions $\vphi_i$, $i =
1,\dots,\Nc$, we get the following $N$ equations for the coefficient
vector $\vu = (u_1,\dots,u_N)$ of $\uh$:
%\begin{subequations*}\label{E:discrete.system}
\begin{alignat*}{2}
\sum\limits_{j = 1}^N u_j \,\ldual{L \vphi_j}{\vphi_i}{\Xc_h^* \times \Xc_h}
&= \ldual{f}{\vphi_i}{\Xc_h^* \times \Xc_h} 
\qquad & &\mbox{for } i = 1,\dots,\Nc, \\
u_i & = g_i & &\mbox{for } i = \Nc+1,\dots,N.
\end{alignat*}
%\end{subequations*}
Defining the \emph{system matrix} $\vec{L}$ by
\begin{equation}\label{E:matrix}
\vL := 
\left[
\begin{matrix}
\ldual{L \vphi_1}{\vphi_1}{} & \dots & \ldual{L \vphi_{\Nc}}{\vphi_1}{} &
\ldual{L \vphi_{\Nc+1}}{\vphi_1}{} & \dots & \ldual{L \vphi_{N\vphantom{\Nc}}}{\vphi_1}{}\\
\vdots & \ddots & \vdots & \vdots & \ddots & \vdots\\
\ldual{L \vphi_1}{\vphi_{\Nc}}{} & \dots & \ldual{L \vphi_{\Nc}}{\vphi_{\Nc}}{} &
\ldual{L \vphi_{\Nc+1}}{\vphi_{\Nc}}{} & \dots & \ldual{L \vphi_{N\vphantom{\Nc}}}{\vphi_{\Nc}}{}\\
0 & \dots & 0 & \multicolumn{3}{c}{1 \hfil 0 \hfil \dots \hfil 0}\\
0 & \dots & 0 & \multicolumn{3}{c}{0 \hfil 1 \hfil \dots \hfil 0}\\
\vdots & \ddots &0&\multicolumn{3}{c}{0 \hfil 0 \hfil \ddots \hfil \vdots\,}\\
0 & \dots & 0 & \multicolumn{3}{c}{0 \hfil 0 \hfil \dots \hfil 1}\\
\end{matrix}
\right]
\end{equation}
and the \emph{right hand side vector} or \emph{load vector} $\vf$ by
\begin{equation}\label{E:Fright}
\vf :=
\left[
\begin{matrix}
\ldual{f}{\vphi_1}{}\\
\vdots\\
\ldual{f}{\vphi_{\Nc}}{}\\
g_{\Nc+1}\\
\vdots\\
g_{N\vphantom{\Nc}}
\end{matrix}
\right],
\end{equation}
we can write the discrete equations as the linear $N \times N$ system
\begin{equation}\label{E:LuF}
\vL \, \vu = \vf,
\end{equation}
which has to be assembled and solved numerically.

\subsection{Discretization of coupled vector valued problems}%
\label{S:DisCoupled}%
\idx{coupled vector valued problems}

This section describes the discretization and assembling of coupled 
vector valued problems. Consider the following artificial coupled Poisson 
problem:

Let $C\in\R^{n\times n}$ be a regular coupling matrix, $f \in
L^2(\Omega;\R^n)$ a given right hand side, and
$g\colon\partial\Omega\to\R^n$ suitable boundary values. Find
$u:\Omega\to\R^n$ with

\begin{subequations}\label{E:coupled_example}
\begin{alignat}{2}
  -\sum_{\nu=1}^n C_{\mu\nu} \Delta u_\nu &= f_\mu \qquad\text{in }\Omega\text{
  for }\mu=1,\dots,n\\ u &= g \qquad\text{on }\partial\Omega,
\end{alignat}
\end{subequations}

By a left multiplication with $C^{-1}$ this problem decouples 
into a set of independent scalar Poisson problems, for which we could 
apply the same existence and uniqueness theory as above. However, we
will refrain from doing this in order to illustrate the concepts of
this section. Generally, the weak form of a coupled system of linear second 
order equations can be written as follows:

Define vector valued spaces $X=H^1(\Omega;\R^n)$,
$\Xc=H_0^1(\Omega;\R^n)$. Find a solution $u \in X$, such that $u \in
g+\Xc$ and
\begin{equation}\label{E:coupled_weak}
 \ldual{L u}{\vphi}{\Xc^* \times \Xc} := 
 \sum_{\mu,\nu=1}^n\int_\Omega \nabla \vphi_\mu \cdot A^{\mu\nu} \nabla u_\nu 
  + \vphi_\mu\, b^{\mu\nu} \cdot \nabla u_\nu
  + c^{\mu\nu}\,\vphi_\mu \, u_\nu \, dx = 
  \int_\Omega f \cdot \vphi\, dx
\end{equation}
for all $\vphi \in \Xc$.

To obtain the weak form of problem \mathref{E:coupled_example} for
example, we would set $b:=0$, $c:=0$ and
$A_{ij}^{\mu\nu}:=\delta_{ij}C_{\mu\nu}$. The next step is to derive a
suitable linear system for a discretization as in the last section.

As mentioned in \secref{S:FES}, basis functions are always
scalar-valued. To gain a vector valued finite element space $X_h$, we
use vector valued coefficients. Choose a set of scalar basis functions
$\left\{\vphi_1,\dots,\vphi_{N\vphantom{\Nc}}\right\}$ as above. For a
function $\vh \in X_h$ we denote by $\vv = (v_1,\dots,v_N)$ the
coefficient vector of $\vh$. Each $v_j$ is now itself a vector
$v_j=(v_{j\mu})_{\mu=1}^n$. Thus, we have the following decomposition:
\[
\vh  =  \sum_{j = 1}^N v_j \vphi_j.
\]
Define $\vphi_j^\mu:=(\delta_{\mu\nu}\vphi_j)_{\nu=1}^n$. The discrete problem
can now be written as a set of linear equations for the coefficients
$u_{j\mu}$:
\begin{alignat*}{2}
  \sum_{j = 1}^N \sum_{\mu=1}^n u_{j\mu} \,\ldual{L
    \vphi_j^\mu}{\vphi_i^\nu}{\Xc_h^* \times \Xc_h} &=
  \ldual{f}{\vphi_i^\nu}{\Xc_h^* \times \Xc_h} \qquad & &\text{for } i =
  1,\dots,\Nc; \nu = 1,\dots,n,\\ u_{i\nu} & = g_{i\nu} & &\text{for } i =
  \Nc+1,\dots,N; \nu=1,\dots,n.
\end{alignat*}
The corresponding system matrix $\vL$ is defined by
\begin{equation}\label{E:coupled_matrix}
\vL := 
\left[
\begin{matrix}
  \vL^{11} & \dots & \vL^{1\Nc} & \vL^{1,\Nc+1} & \dots & \vL^{1N}\\
  \vdots & \ddots & \vdots & \vdots & \ddots & \vdots\\
  \vL^{\Nc1} & \dots & \vL^{\Nc\Nc} & \vL^{\Nc,\Nc+1} & \dots & \vL^{\Nc N}\\
  0 & \dots & 0 & \multicolumn{3}{c}{\vI \hfil 0 \hfil \dots \hfil 0}\\
  0 & \dots & 0 & \multicolumn{3}{c}{0 \hfil \vI \hfil \dots \hfil 0}\\
  \vdots & \ddots &0&\multicolumn{3}{c}{0 \hfil 0 \hfil \ddots \hfil\vdots\,}\\
  0 & \dots & 0 & \multicolumn{3}{c}{0 \hfil 0 \hfil \dots \hfil \vI}\\
  \end{matrix}
\right]
\end{equation}
with $\vI\in\R^{n\times n}$ an identity matrix and 
\[
\vL^{ij} := (\ldual{L\vphi_j^\mu}{\vphi_i^\nu}{\Xc_h^* \times \Xc_h})_{\mu,\nu=1}^n.
\]
The load vector $\vf$ is defined by
\begin{equation}\label{E:coupled_Fright}
\vf :=
\left[
\begin{matrix}
\ldual{f}{\vphi_1^1}{}\\
\vdots\\
\ldual{f}{\vphi_1^n}{}\\
\vdots\\
\ldual{f}{\vphi_{\Nc}^1}{}\\
\vdots\\
\ldual{f}{\vphi_{\Nc}^n}{}\\
g_{\Nc+1,1}\\
\vdots\\
g_{Nd}
\end{matrix}
\right].
\end{equation}
The problem can now be written as the linear $Nd \times Nd$ system
\begin{equation}\label{E:coupled_LuF}
  \vL \, \vu = \vf,
\end{equation}
which has to be assembled and solved. The organization of vectors and matrices
using small $n$-size blocks as components was chosen with the goal of
efficient cache usage during matrix-vector multiplication.

\subsection{Numerical quadrature}%
\label{S:quadrature}%
\idx{numerical quadrature}

For the assemblage of the system matrix and right hand side
vector of the linear system \mathref{E:LuF}, we have to compute
integrals, for example
\[
\int_\Omega f(x) \vphi_i(x)\,dx.
\]
For general data $A$, $b$, $c$, and $f$, these integrals can not be
calculated exactly. Quadrature formulas have to be used in order to
calculate the integrals approximately. Numerical integration in finite
element methods is done by looping over all grid elements and using a
quadrature formula on each element.

\begin{definition}[Numerical quadrature]\label{D:numer-quad}
\idx{numerical quadrature}
A \emph{numerical quadrature} $\hat Q$ on $\Shat$ is a set
$\{(w_k,\lambda_k) \in \R \times \R^{d+1}; k = 0, \dots, n_Q-1\}$
of weights $w_k$ and quadrature points $\lambda_k \in \bar S$ 
(i.e.\ given in barycentric coordinates) such that
\[
 \int_{\Shat} f(\hat x)\, d\hat{x}  \approx  \hat Q(f)  :=  
 \sum_{k = 0}^{n_Q-1}  w_k f(\hat{x}(\lambda_k)).
\]
It is called \emph{exact of degree $p$} for some $p \in \N$ if 
\[
 \int_{\Shat} q(\hat x)\, d\hat{x}  =  \hat Q(q)
 \qquad\mbox{for all }q \in \P_p(\Shat).
\]
It is called \emph{stable} if
\[
 w_k > 0 \qquad\mbox{for all } k = 0, \dots, n_Q-1.
\]
\end{definition}

\begin{remark}\label{R:numerical_int} A given numerical quadrature $\hat Q$ on 
$\Shat$ defines for each element $S$ a numerical quadrature $Q_S$.
Using the transformation rule we define 
$Q_S$ on an element $S$ which is parameterized by 
$F_S \colon\, \Shat \to S$ and a function $f : S \to \R$:
\begin{equation*}\label{E:numer-quad-S}
\int_{S} f(x)\, dx  \approx  Q_S(f) :=  
\hat Q((f \circ F_S) |\det DF_S|)  =
\sum_{k = 0}^{n_Q-1}  w_k f(x(\lambda_k)) |\det DF_S(\hat x(\lambda_k)|.
\end{equation*}
For a simplex $S$ this results in
\begin{equation*}\label{E:numer-quad-sim}
Q_S(f) = d!\, |S| \sum_{k = 0}^{n_Q-1}  w_k f(x(\lambda_k)).
\end{equation*}
\end{remark}

\subsection{Finite element discretization of 2nd order problems}%
\label{S:FE-dis-2nd}%
\idx{assemblage of discrete system|(}

Let $\Pbar$ be a finite dimensional function space on $\Sbar$ with
basis $\{\pbar^1,\dots,\pbar^m\}$. For a conforming finite element
discretization of \mathref{E:weak} we use the finite element space
$X_h = X_h(\tria,\Pbar,X)$. For this space $\Xc_h$ is
given by $X_h(\tria,\Pbar,\Xc)$.

\idx{assemblage of discrete system!load vector}
By the relation \mathref{E:global-lokal.a} for
global and local basis functions, we obtain for the $j$th component of
the right hand side vector $\vf$
\begin{align*}
\ldual{f}{\vphi_j}{} &= 
\int_\Omega f(x)\, \vphi_j(x)\,dx
=
\sum_{S \in \tria} \int_S f(x)\, \vphi_j(x)\,dx
=
\sum_{\atop{S \in \tria}{S \subset \supp(\vphi_j)}} 
   \int_S f(x) \pbar^{i_S(j)}(\lambda(x))\,dx\\
&=
\sum_{\atop{S \in \tria}{S \subset \supp(\vphi_j)}}
   \int_{\Shat} f(F_S(\xhat)) \pbar^{i_S(j)}(\lambda(\xhat)) 
       |\det DF_S(\xhat)|\,d\xhat,
\end{align*}
where $S$ is parameterized by $F_S\colon\Shat\to S$.  
The above sum is reduced to a sum over all $S \subset \supp(\vphi_j)$
which are only few elements due to the small support of $\vphi_j$.

The right hand side vector can be assembled in the following way:
First, the right hand side vector is initialized with zeros. For each
element $S$ of $\tria$ we calculate the \emph{element load vector}
$\vf_S = (f_S^1,\dots,f_S^m)^t$, where
\begin{equation}\label{E:int_S-f-phi}
f_S^i =    \int_{\Shat} f(F_S(\xhat)) \pbar^{i}(\lambda(\xhat)) 
       |\det DF_S(\xhat)|\,d\xhat, \qquad i = 1,\dots,m.
\end{equation}
Denoting again by $j_S: \{1,\dots,m\} \to J_S$ the function which connects
the local DOFs with the global DOFs (defined in \mathref{E:global-lokal.b}),
the values $f_S^i$ are then added to the $j_S(i)$th 
component of the right hand side vector $\vf$, $i = 1,\dots,m$.

For general $f$, the integrals in \mathref{E:int_S-f-phi} can not be
calculated exactly and we have to use a quadrature formula for the
approximation of the integrals (compare Section~\ref{S:quadrature}).
Given a numerical quadrature $\hat Q$ on $\Shat$ we approximate
\begin{align}\label{EX:quad_right_param}
f_S^i &\approx \hat Q\left(
         (f \circ F_S)\,(\pbar^{i}\circ \lambda) |\det DF_S| \right)
      = \sum_{k = 0}^{n_Q-1}  w_k f(x(\lambda_k))\,\pbar^{i}(\lambda_k)
 |\det DF_S(\hat x(\lambda_k)|.
\end{align}
For a simplex $S$ this is simplified to
\begin{equation*}\label{E:quad_right}
f_S^i \approx
d!\, |S| \,\sum_{k = 0}^{n_Q-1}  w_k f(x(\lambda_k))\,\pbar^{i}(\lambda_k).
\end{equation*}

In \ALBERTA, information about values of basis functions and its
derivatives as well as information about the connection of global and
local DOFs (i.e.\ information about $j_S$) is stored in special data
structures for local basis functions (compare
Section~\ref{man:S:basfct_impl}).  By such information, the element
load vector can be assembled by a general routine if a function for
the evaluation of the right hand side is supplied. For parametric
elements, a function for evaluating $|\det DF_S|$ is additionally
needed. The assemblage into the global load vector $\vf$ can again be
performed automatically, using information about the connection of
global and local DOFs.

\idx{assemblage of discrete system!system matrix}
The calculation of the system matrix is also done element--wise.
Additionally, we have to handle derivatives of basis functions.
Looking first at the second order term we obtain by the chain
rule \mathref{E:grad_phi} and the relation \mathref{E:global-lokal} for
global and local basis functions
\begin{align*}
\int_S \nabla \vphi_i(x) \cdot A(x) \nabla \vphi_j(x)\,dx &=
\int_S \nabla (\pbar^{i_S(i)} \circ \lambda)(x) \cdot A(x)
       \nabla (\pbar^{i_S(j)} \circ \lambda)(x) \,dx\\
&=
\int_S 
 \nablal \pbar^{i_S(i)}(\lambda(x)) \cdot (\Lambda(x)\, A(x)\, \Lambda^t(x))
    \nablal \pbar^{i_S(j)}(\lambda(x))\,dx,
\end{align*}
where $\Lambda$ is the Jacobian of the barycentric coordinates $\lambda$
on $S$. In the same manner we obtain for the first and zero order terms
\[
 \int_S \vphi_i(x) \, b(x) \cdot  \nabla \vphi_j(x) \,dx =
\int_S \pbar^{i_S(i)}(\lambda(x))\,(\Lambda(x)\, b(x)) \cdot 
       \nablal \pbar^{i_S(j)}(\lambda(x)) \,dx
\]
and
\[
 \int_S c(x)\,\vphi_i(x) \, \vphi_j(x)\,dx =
 \int_S c(x)\, \pbar^{i_S(i)}(\lambda(x)) \, \pbar^{i_S(j)}(\lambda(x))\,dx.
\]
Using on $S$ the abbreviations
%\begin{subequations*}\label{E:bar_A_b_c}
\begin{align*}
\bar A(\lambda) &:=
\left(\bar a_{kl}(\lambda)\right)_{k,l = 0,\ldots,d} 
  := |\det DF_S(\xhat(\lambda))| \,
\Lambda(x(\lambda)) \, A(x(\lambda)) \, \Lambda^t(x(\lambda)),\\
\bar b(\lambda) & :=
\left(\bar b_{l}(\lambda)\right)_{l = 0,\ldots,d} 
 := |\det DF_S(\xhat(\lambda))| \,
\Lambda(x(\lambda)) \, b(x(\lambda)), \quad \mbox{and}\\
\bar c(\lambda) & := |\det DF_S(\xhat(\lambda))| \, c(x(\lambda))
\end{align*}
%\end{subequations*}
and transforming the integrals to the reference simplex, we can write
the \emph{element matrix} $\vL_S = (L_S^{ij})_{i,j=1,\dots,m}$ as
\begin{align}\label{E:L-phii-phij}
L_S^{ij} &= 
\int_{\Shat}
 \nablal \pbar^{i}(\lambda(\xhat)) \cdot \bar A(\lambda(\xhat))\,
    \nablal \pbar^{j}(\lambda(\xhat))\,d\xhat
+
\int_{\Shat} \pbar^{i}(\lambda(\xhat))\;
  \bar b(\lambda(\xhat)) \cdot 
       \nablal \pbar^{j}(\lambda(\xhat)) \,d\xhat  \notag\\
&\qquad +
\int_{\Shat} \bar c(\lambda(\xhat))\, \pbar^{i}(\lambda(\xhat)) \,
                             \pbar^{j}(\lambda(\xhat))\,d\xhat,
\end{align}
or writing the matrix--vector and vector--vector products explicitly
\begin{align*}\label{E:L-phii-phij'}
L_{S}^{ij} &=
\sum_{k,l = 0}^d \int_{\Shat}
\bar a_{kl}(\lambda(\xhat)) \, \pbar_{,\lambda_k}^i(\lambda(\xhat)) \, 
                \pbar_{,\lambda_l}^j(\lambda(\xhat)) \, d\xhat
+ \sum_{l = 0}^d \int_{\Shat}
\bar b_{l}(\lambda(\xhat)) \, \pbar^i(\lambda(\xhat)) \,
                \pbar_{,\lambda_l}^j(\lambda(\xhat))\, d\xhat
%\notag
\\
%\tag{\ref{E:L-phii-phij}'}
&\qquad +  \int_{\Shat}
\bar c(\lambda(\xhat)) \, \pbar^i(\lambda(\xhat))\,\pbar^j(\lambda(\xhat))\,
 d\xhat,
\end{align*}
$i,j = 1,\dots,m$. Using quadrature formulas $\hat Q_2$, $\hat Q_1$, and
$\hat Q_0$ for the second, first and zero order term we approximate
the element matrix
\begin{equation*}\label{E:quad_L}
L_S^{ij}  \approx 
\hat Q_2\left(\sum_{k,l = 0}^d (\bar a_{kl} \pbar^i_{,\lambda_k} \, 
    \pbar^j_{,\lambda_l}) \circ\lambda\right)
  + 
\hat Q_1\left(\sum_{l = 0}^d (
\bar b_{l} \, \pbar^i \, \pbar^j_{,\lambda_l})\circ\lambda\right)
  +  \hat Q_0\Big(
 (\bar c \, \pbar^i \,\pbar^j)\circ\lambda\Big),
\end{equation*}
$i,j=1,\dots,m$. Having access to functions for the evaluation of
\begin{equation*}\label{E:LALt-Lb-c}
\bar a_{kl}(\lambda_q), \quad \bar b_{l}(\lambda_q),\quad
\bar c(\lambda_q)
\end{equation*}
at all quadrature points $\lambda_q$ on $S$, $\vec{L}_S$ can
be computed by some general routine. The assemblage into the
system matrix can also be done automatically (compare the assemblage
of the load vector).

\begin{remark}
Due to the small support of the global basis
function, the system matrix is a sparse matrix, i.e.\ the maximal
number of entries in all matrix rows is much smaller than the
size of the matrix. Special data structures are needed for an
efficient storage of sparse matrices and they are described in
Section~\ref{man:S:DOF_MATRIX}.
\end{remark}

\begin{remark}
The calculation of the gradient of the barycentric coordinates
$\Lambda(x(\lambda))$ usually involves the
determinant of the parameterization's Jacobian
$|\det DF_S(\xhat(\lambda))|$. Thus, a calculation of
$|\det DF_S(\xhat(\lambda))| \,\Lambda(x(\lambda)) \, A(x(\lambda))
\, \Lambda^t(x(\lambda))$  may be much faster than the calculation of
$\Lambda(x(\lambda)) \, A(x(\lambda)) \, \Lambda^t(x(\lambda))$ only;
the same holds for the first order term.
\end{remark}

Assuming that the coefficients $A$, $b$, and $c$ are piecewise
constant on a non--parametric triangulation, $\bar A(\lambda)$,
$\bar b(\lambda)$, and $\bar c(\lambda)$ are constant on each simplex $S$
and thus simplified to
\begin{equation*}\label{E:bar_A_b_cS}
\bar A_S = \left(\bar a_{kl}\right)_{k,l = 0,\ldots,d} 
 =  d!|S| \, \Lambda \, A_{|S} \, \Lambda^t, \quad
\bar b_S = \left(\bar b_{l}\right)_{l = 0,\ldots,d} 
 =  d!|S| \,\Lambda \, b_{|S}, \quad
\bar c_S  =  d!|S| \, c_{|S}.
\end{equation*}
For the approximation of the element matrix by quadrature we then
obtain
\begin{equation}\label{E:quad_LS}
\vec{L}_S^{ij}  \approx 
\sum_{k,l = 0}^d \bar a_{kl} 
\hat Q_2\Big((\pbar^i_{,\lambda_k} \, \pbar^j_{,\lambda_l})\circ\lambda\Big)
  + 
\sum_{l = 0}^d \bar b_{l} \, 
\hat Q_1\Big((\pbar^i \, \pbar^j_{,\lambda_l})\circ\lambda\Big)
  +  \bar c_S\, \hat Q_0\Big(
 (\pbar^i \,\pbar^j)\circ\lambda\Big)
\end{equation}
$i,j=1,\dots,m$. Here, the numerical quadrature is only applied
for approximating integrals of the basis
functions on the standard simplex. Theses values can be computed
only once, and can then be used on each simplex $S$. This will speed
up the assembling of the system matrix. Additionally, for polynomial
basis functions we can use quadrature formulas which integrate these
integrals exactly.
\idx{assemblage of discrete system|)}
\idx{assemblage of discrete system!coupled case|(}
So far we have only considered the case of scalar problems. The
transition to (coupled) vector valued problems is straight forward and
simply involves two more indices. The entries of the element matrix
are now $d\times d$ matrices themselves:
\begin{align*}
L_{S,\mu\nu}^{ij} &:= 
 \int_S \nabla \vphi_i(x) \cdot A^{\mu\nu}(x) \nabla 
 \vphi_j(x) + \vphi_i(x)\, b^{\mu\nu}(x) \cdot \nabla \vphi_j(x)
  + c^{\mu\nu}(x)\,\vphi_i(x) \, \vphi_j(x) \, dx\\ 
  &= 
  \int_{\Shat}
  \nablal \pbar^{i}(\lambda(\xhat)) \cdot \bar A^{\mu\nu}(\lambda(\xhat))\,
  \nablal \pbar^{j}(\lambda(\xhat))\,d\xhat
  +
  \int_{\Shat} \pbar^{i}(\lambda(\xhat))\;
  \bar b^{\mu\nu}(\lambda(\xhat)) \cdot 
  \nablal \pbar^{j}(\lambda(\xhat)) \,d\xhat  \notag\\
  &\qquad +
  \int_{\Shat} \bar c^{\mu\nu}(\lambda(\xhat))\, \pbar^{i}(\lambda(\xhat)) \,
  \pbar^{j}(\lambda(\xhat))\,d\xhat,
\end{align*}
with
\begin{align*}
  \bar A^{\mu\nu}(\lambda) &:= \left(\bar
  a_{kl}^{\mu\nu}(\lambda)\right)_{k,l = 0,\ldots,d} := |\det
  DF_S(\xhat(\lambda))| \, \Lambda(x(\lambda)) \,
  A^{\mu\nu}(x(\lambda)) \, \Lambda^t(x(\lambda)),\\
  \bar b^{\mu\nu}(\lambda) & := \left(\bar
  b_{l}^{\mu\nu}(\lambda)\right)_{l = 0,\ldots,d} := |\det
  DF_S(\xhat(\lambda))| \, \Lambda(x(\lambda)) \, b^{\mu\nu}(x(\lambda)),
  \quad\text{and}\\
  \bar c^{\mu\nu}(\lambda) & := |\det
  DF_S(\xhat(\lambda))| \, c^{\mu\nu}(x(\lambda))
\end{align*}
for $\mu,\nu=1,\dots,d$. The approximation of the integrals
using quadratures is done analogously to the scalar case.

As a result, using information about values of basis functions and
their derivatives, and information about the connection of global and
local DOFs, the linear system can be assembled automatically by some
general routines. Only functions for the evaluation of given data have
to be provided for special applications. The general assemble
routines are described in Section~\ref{man:S:ass_tools}.
\idx{assemblage of discrete system!coupled case|)}%
\idx{finite element discretization|)}

%%% Local Variables: 
%%% mode: latex
%%% TeX-master: "alberta-book"
%%% End: