\def\jump#1{\left[\negthinspace\left[{#1}\right]\negthinspace\right]}
\section{Implementation of error estimators}%
\label{S:estimator}%
\idx{error estimators|(}

\subsection{Error estimator for elliptic problems}%
\label{S:ellipt_est}

\ALBERTA provides a residual type error estimator for non--linear
elliptic problems of the type
\begin{align*}
  -\nabla \cdot A \nabla u(x) + f\Big(x,u(x),\nabla u(x)\Big) &= 0 
&&x \in \Omega,\\
                        u(x)  &= g_d &&x \in \Gamma_D,\\
       \nu\cdot A\nabla u(x)  &= g_n &&x \in \Gamma_N,
\end{align*}
where $A \in \R^{n\times n}$ is a positive definite matrix and
$\partial\Omega=\Gamma_D\cup\Gamma_N$.
%%
\ALBERTA implements for this kind of equations the $L^2$ and $H^1$ 
per-element estimators $\eta_{S,0}$ and $\eta_{S,1}$ ($S\in\tria$)
\[
\begin{split}
  \eta_{S,0}^2&:=
  C_0^2\, h_S^4\, \|-\nabla \cdot A \nabla \uh + f(.,\uh,\nabla\uh)\|_{L^2(S)}^2\\
  &\qquad+ C_1^2\, \sum_{\Gamma\subset\partial S \cap\Omega}
  h_S^3\, \|\jump{A\nabla \uh}\|_{L^2(\Gamma)}^2
  + C_1^2\, \sum_{\Gamma\subset\partial S \cap\Gamma_N}
  h_S^3\, \|\nu\cdot A\nabla \uh-g_n\|_{L^2(\Gamma)}^2,
  \\
  \eta_{S,1}^2&:=
  C_0^2\,h_S^2\, \|-\nabla \cdot A \nabla \uh + f(.,\uh,\nabla \uh)\|_{L^2(S)}^2
  \\
  &\qquad + C_1^2\, \sum_{\Gamma\subset\partial S \cap \Omega} 
  h_S\, \|\jump{A\nabla \uh}\|_{L^2(\Gamma)}^2
  + C_1^2\, \sum_{\Gamma\subset\partial S \cap \Gamma_N} 
  h_S\, \|\nu\cdot A\nabla \uh-g_n\|_{L^2(\Gamma)}^2,
  \\
\end{split}
\]
where $\jump{.}$ denotes the jump of the normal component across an
interior co-dimension $1$ sub-simplex (vertex/edge/face) $\Gamma
\subset \partial S$.

Verf\"urth proved for $g_d\equiv 0$ and $g_n\equiv 0$ in
\cite{Verfuerth:94b} -- under suitable assumptions on $f$,
$u$ and $\uh$ in the non-linear case -- the estimate
\[
\|u - \uh\|_{H^1(\Omega)}^2\leq\sum_{S\in\tria} \eta_{S,1}^2,
\]
and B\"ansch and Siebert \cite{BaenschSiebert:95} proved a similar the
$L^2$-estimate for the semi--linear case $f = f(x,u)$ and $g_d\equiv
0$ and $\Gamma_N=\emptyset$:
\[
\|u - \uh\|_{L^2(\Omega)}^2\leq \sum_{S\in\tria}\eta_{S,0}^2.
\]
%%
The following functions implement above estimators for scalar and
vector-valued functions; the implementation works also for meshes with
non-zero co-dimension as well as for periodic meshes.
%%
\fdx{ellipt_est()@{\code{ellipt\_est()}}}%
\idx{error estimators!ellipt_est()@{\code{ellipt\_est()}}}%
\fdx{ellipt_est_dow()@{\code{ellipt\_est\_dow()}}}%
\idx{error estimators!ellipt_est_dow()@{\code{ellipt\_est\_dow()}}}%
\fdx{ellipt_est_d()@{\code{ellipt\_est\_d()}}}%
\idx{error estimators!ellipt_est_d()@{\code{ellipt\_est\_d()}}}%
\bv\begin{verbatim}
REAL ellipt_est(const DOF_REAL_VEC *uh, ADAPT_STAT *adapt,
                REAL *(*rw_est)(EL *), REAL *(*rw_estc)(EL *),
                int quad_deg,
                NORM norm, REAL C[3], const REAL_DD A,
                const BNDRY_FLAGS dirichlet_bndry,
                REAL (*f)(const EL_INFO *el_info,
                          const QUAD *quad, int qp,
                          REAL uh_qp, const REAL_D grd_uh_gp),
                FLAGS f_flags,
                REAL (*gn)(const EL_INFO *el_info,
                           const QUAD *quad, int qp,
                           REAL uh_qp, const REAL_D normal),
                FLAGS gn_flags);

REAL ellipt_est_dow(const DOF_REAL_VEC_D *uh, ADAPT_STAT *adapt,
                    REAL *(*rw_est)(EL *), REAL *(*rw_estc)(EL *),
                    int quad_deg,
                    NORM norm, REAL C[3],
                    const void *A, MATENT_TYPE A_type, MATENT_TYPE A_blocktype,
                    bool sym_grad,
                    const BNDRY_FLAGS dirichlet_bndry,
                    const REAL *(*f)(REAL_D result,
                                     const EL_INFO *el_info,
                                     const QUAD *quad, int qp,
                                     const REAL_D uh_qp,
                                     const REAL_DD grd_uh_gp),
                    FLAGS f_flags,
                    const REAL *(*gn)(REAL_D result,
                                      const EL_INFO *el_info,
                                      const QUAD *quad, int qp,
                                      const REAL_D uh_qp,
                                      const REAL_D normal),
                    FLAGS gn_flags);

REAL ellipt_est_d(const DOF_REAL_D_VEC *uh, ADAPT_STAT *adapt,
                  REAL *(*rw_est)(EL *), REAL *(*rw_estc)(EL *),
                  int quad_deg,
                  NORM norm, REAL C[3],
                  const void *A, MATENT_TYPE A_type, MATENT_TYPE A_blocktype,
                  bool sym_grad,
                  const BNDRY_FLAGS dirichlet_bndry,
                  const REAL *(*f)(REAL_D result,
                                   const EL_INFO *el_info,
                                   const QUAD *quad, int qp,
                                   const REAL_D uh_qp,
                                   const REAL_DD grd_uh_gp),
                  FLAGS f_flags,
                  const REAL *(*gn)(REAL_D result,
                                    const EL_INFO *el_info,
                                    const QUAD *quad, int qp,
                                    const REAL_D uh_qp,
                                    const REAL_D normal),
                  FLAGS gn_flags);
\end{verbatim}\ev
%%
Description:
\begin{descr}
\kitem{ellipt\_est(uh, adapt, rw\_est, rw\_estc, quad\_deg, norm, C,}
\kitem{~~~~~~~~~~~A, dirichlet\_bndry, f, f\_flags, gn, gn\_flags)}~\hfill

  computes an error estimate of the above type for the
  $H^1$ or $L^2$ norm; the return value is an
  approximation of the estimate $\|u - \uh\|$ by quadrature.

  \begin{descr}
    
  \kitem{uh} is a vector storing the coefficients of the
    discrete solution; if \code{uh} is a \nil pointer, nothing is
    done, the return value is \code{.0}.

  \kitem{adapt} is a pointer to an \code{ADAPT\_STAT} structure;
    if not \nil, the entries \code{adapt->p=2}, \code{err\_sum},
    and \code{err\_max} of \code{adapt} are set by
    \code{ellipt\_est()} (compare
    \secref{S:adapt_stat_in_ALBERTA}).

  \kitem{rw\_el\_est} is a function for writing the local error
    indicator for a single element (usually to some location
    inside \code{leaf\_data}, compare \secref{S:leaf_data_info});
    if this function is \nil, only the global estimate is
    computed, no local indicators are stored.
    \code{rw\_el\_est(el)} returns for each leaf element
    \code{el} a pointer to a \code{REAL} for storing the square
    of the element indicator, which can directly be used in the
    adaptive method, compare the \code{get\_el\_est()} function
    pointer in the \code{ADAPT\_STAT} structure (compare
    \secref{S:adapt_stat_in_ALBERTA}).

  \kitem{rw\_el\_estc} is a function for writing the local
    coarsening error indicator for a single element (usually to
    some location inside \code{leaf\_data}, compare
    \secref{S:leaf_data_info}); if this function is \nil, no
    coarsening error indicators are computed and stored;
    \code{rw\_el\_estc(el)} returns for each leaf element
    \code{el} a pointer to a \code{REAL} for storing the square
    of the element coarsening error indicator.

  \kitem{quad\_deg} is the degree of the quadrature that should
    be used for the approximation of the norms on the elements
    and edges/faces; if \code{degree} is less than zero a
    quadrature which is exact of degree
    \code{2*uh->fe\_space->bas\_fcts->degree} is used.

  \kitem{norm} can be either
    \code{H1\_NORM}\cdx{H1_NORM@{\code{H1\_NORM}}} or
    \code{L2\_NORM}\cdx{L2_NORM@{\code{L2\_NORM}}} (which are
    defined as symbolic constants in \code{alberta.h}) to
    indicate that the $H^1$ or $L^2$ error estimate has to be
    calculated.

  \kitem{C[0], C[1], C[2]} are the constants in
    front of the element residual, wall residual, and coarsening
    term respectively. If \code{C} is \nil, then all constants
    are set to $1.0$.

  \kitem{A} is the constant matrix of the second order term.

  \kitem{dirichlet\_bndry} A bit-mask marking those parts of the
    boundary which are subject to Dirichlet boundary conditions, see
    \secref{S:boundary}.

  \kitem{f} is a pointer to a function for the evaluation of the
    lower order terms at all quadrature nodes, i.e.
    $f(x(\lambda), u(\lambda), \nabla u(\lambda))$ ; if \code{f}
    is a \nil pointer, $f\equiv0$ is assumed;

    \code{f(el\_info, quad, qp, uh\_qp, grd\_uh\_qp)} returns the
    value of the lower oder terms on element \code{el\_info->el}
    at the quadrature node \code{quad->lambda[qp]}, where
    \code{uh\_qp} is the value and \code{grd\_uh\_qp} the
    gradient (with respect to the Cartesian coordinates) of the
    discrete solution at that quadrature point. See also
    \code{f\_flag} below:

  \kitem{f\_flag} specifies whether the function \code{f()}
    actually needs values of \code{uh\_qp} or \code{grd\_uh\_qp},
    \code{f\_flag} may be $0$ or
    \code{INIT\_UH}\cdx{INIT_UH@{\code{INIT\_UH}}} or
    \code{INIT\_GRD\_UH}\cdx{INIT_GRD_UH@{\code{INIT\_GRD\_UH}}}
    or their bitwise composition (\code{|}).  The arguments
    \code{uh\_qp} and \code{grd\_uh\_qp} of \code{f()} only hold
    valid information if the flags \code{INIT\_UH} respectively
    \code{INIT\_GRD\_UH} are set.

  \kitem{gn(el\_info, quad, qp, uh\_qp, normal)} is a pointer to a
    function for the evaluation of non-homogeneous Neumann boundary
    data.  \code{gn} may be \nil, in which case zero Neumann boundary
    conditions are assumed.  The argument \code{normal} always
    contains the normal of the Neumann boundary facet. In the case of
    non-vanishing co-dimension \code{normal} lies in the
    lower-dimensional space which is spanned by the mesh simplex
    defined by \code{el\_info}. \code{gn()} is evaluated on those
    parts of the boundary which are \emph{not} flagged as
    Dirichlet-boundaries by the argument \code{dirichlet\_bndry}.

  \kitem{gn\_flag} controls whether the argument \code{uh\_qp}
    of the function \code{gn()} actually contains the value of
    \code{uh} at the quadrature point \code{qp}. Note that the
    argument \code{normal} always contains valid data.
  \end{descr}
  
  The estimate is computed by traversing all leaf elements of
  \code{uh->fe\_space->mesh}, using the quadrature for the
  approximation of the residuals and storing the square of the
  element indicators on the elements (if \code{rw\_el\_est} and
  \code{rw\_el\_estc }are not \nil).

\kitem{ellipt\_est\_d(uh, adapt, rw\_est, rw\_estc, quad\_deg, norm, C,}
\kitem{~~~~~~~~~~~~~A, A\_type, A\_blocktype, sym\_grad,}
\kitem{~~~~~~~~~~~~~dirichlet\_bndry, f, f\_flags, gn, gn\_flags)}\hfill
\kitem{ellipt\_est\_dow(uh, adapt, rw\_est, rw\_estc, quad\_deg, norm, C,}
\kitem{~~~~~~~~~~~~~~~A, A\_type, A\_blocktype, sym\_grad,}
\kitem{~~~~~~~~~~~~~~~dirichlet\_bndry, f, f\_flags, gn, gn\_flags)}\hfill

Similar function for a (coupled) vector valued elliptic problem. We
document only the arguments which are different from the arguments of
\code{ellipt\_est()}:

\begin{descr}
  \kitem{A}
  now represents a tensor
  $(A_{ij}^{\mu\nu}\in\R^{n\times n,n\times n}$,
  $i,j,\mu,\nu=0,\dots,n-1$. The indexing is
  \begin{equation*}
    \code{A[i][j][mu][nu]} = A^{\code{mu},\code{nu}}_{\code{ij}},
  \end{equation*} 
  with \code{i,j,mu,nu==0,\dots,\code{DIM\_OF\_WORLD-1}}, see
  \secref{book:S:DisCoupled}. \code{A} describes the coefficients of the
  principal part of a coupled system of elliptical equations:
  \[
  -\sum_{\nu,i,j=0}^{n-1}\partial_i A_{ij}^{\mu\nu}\partial_j u^\nu +\text{lower order terms} = f^\mu\quad(\mu = 0,\,\dots,\,n-1).
  \]
  The \code{quasi-stokes.c} demo-program contains an example.

\kitem{A\_blocktype} must be one of \code{MATENT\_REAL},
  \code{MATENT\_REAL\_D} or \code{MATENT\_REAL\_DD}. It specifies the
  symmetry type for coupling of the PDE system. Note that the storage
  layout of \code{A} is determined by the argument
  \code{A\_blocktype}:
  \begin{descr}
    \kitem{MATENT\_REAL:~~~} \code{REAL A[DIM\_OF\_WORLD][DIM\_OF\_WORLD];}
    \kitem{MATENT\_REAL\_D:~} \code{REAL\_D A[DIM\_OF\_WORLD][DIM\_OF\_WORLD];}
    \kitem{MATENT\_REAL\_DD:} \code{REAL\_DD A[DIM\_OF\_WORLD][DIM\_OF\_WORLD];}
  \end{descr}
  \code{A\_blocktype == MATENT\_REAL} or
  \code{A\_blocktype == MATENT\_REAL\_D} means that the system is
  actually decoupled.

\kitem{A\_type} must be one of \code{MATENT\_REAL},
  \code{MATENT\_REAL\_D} or \code{MATENT\_REAL\_DD}. It specifies the
  symmetry type of \code{A} with respect to the first two indices. For
  a Laplacian, for example, one would use \code{DOWBM\_SCAL}. Note
  that the value of \code{A\_type} does \emph{not} change the storage
  layout of the array \code{A}.

  \kitem{sym\_grad} If set to \code{true} then it is assumed that the
  symmetric gradient has to be used for the computation of the jump-
  and Neumann-residuals. The demo-program \code{quasi-stokes.c} uses
  this feature to implement an error estimator for the Stokes equation
  with stress boundary conditions.

  \kitem{f} If the first argument of the function pointer
  \code{f(result,\dots)} is not \nil then the result \emph{must} be
  stored in the argument \code{result} and \code{f()} must return the
  base address of the array \code{result}. If \code{result} is \nil,
  then \code{f()} must store the result in a non-volatile storage area
  and return the address of that area.

  \kitem{dirichlet\_bndry} A bit-mask marking those parts of the
    boundary which are subject to Dirichlet boundary conditions, see
    \secref{S:boundary}.

  \kitem{f} is a pointer to a function for the evaluation of the
    lower order terms at all quadrature nodes, i.e.
    $f(x(\lambda), u(\lambda), \nabla u(\lambda))$ ; if \code{f}
    is a \nil pointer, $f\equiv0$ is assumed;

    \code{f(el\_info, quad, qp, uh\_qp, grd\_uh\_qp)} returns the
    value of the lower oder terms on element \code{el\_info->el}
    at the quadrature node \code{quad->lambda[qp]}, where
    \code{uh\_qp} is the value and \code{grd\_uh\_qp} the
    gradient (with respect to the Cartesian coordinates) of the
    discrete solution at that quadrature point. See also
    \code{f\_flag} below:

  \kitem{f\_flag} specifies whether the function \code{f()}
    actually needs values of \code{uh\_qp} or \code{grd\_uh\_qp},
    \code{f\_flag} may be $0$ or
    \code{INIT\_UH}\cdx{INIT_UH@{\code{INIT\_UH}}} or
    \code{INIT\_GRD\_UH}\cdx{INIT_GRD_UH@{\code{INIT\_GRD\_UH}}}
    or their bitwise composition (\code{|}).  The arguments
    \code{uh\_qp} and \code{grd\_uh\_qp} of \code{f()} only hold
    valid information if the flags \code{INIT\_UH} respectively
    \code{INIT\_GRD\_UH} are set.

  \kitem{gn(el\_info, quad, qp, uh\_qp, normal)} is a pointer to a
    function for the evaluation of non-homogeneous Neumann boundary
    data.  \code{gn} may be \nil, in which case zero Neumann boundary
    conditions are assumed.  The argument \code{normal} always
    contains the normal of the Neumann boundary facet. In the case of
    non-vanishing co-dimension \code{normal} lies in the
    lower-dimensional space which is spanned by the mesh simplex
    defined by \code{el\_info}. \code{gn()} is evaluated on those
    parts of the boundary which are \emph{not} flagged as
    Dirichlet-boundaries by the argument \code{dirichlet\_bndry}.

  \kitem{gn\_flag} controls whether the argument \code{uh\_qp}
    of the function \code{gn()} actually contains the value of
    \code{uh} at the quadrature point \code{qp}. Note that the
    argument \code{normal} always contains valid data.
\end{descr}
\end{descr}

\begin{example}[Linear problem]\label{E:est-impl}
Consider the scalar linear model problem \mathref{book:E:strong} with constant 
coefficients $A$, $b$, and $c$:
\begin{alignat*}{2}
-\nabla \cdot A \nabla u + b \cdot \nabla u + c\, u &= r \qquad & 
&\mbox{in } \Omega,\\
u &= 0        & &\mbox{on } \partial\Omega.
\end{alignat*}
Let \code{A} be a \code{REAL\_DD} matrix storing $A$, which is then the
eighth argument of \code{ellipt\_est()}. Assume that
\code{const REAL *b(const REAL\_D)} is a function returning a pointer
to a vector storing $b$, \code{REAL c(REAL\_D)} returns the
value of $c$ and \code{REAL r(const REAL\_D)} returns the value of the
right hand side $r$ of \mathref{book:E:strong} at some point in world
coordinates. The implementation of the function \code{f} is:
\bv\begin{verbatim}
static REAL f(const EL_INFO *el_info, const QUAD *quad, int iq, REAL uh_iq, 
              const REAL_D grd_uh_iq)
{
  FUNCNAME("f");
  const REAL  *bx, *x;
  extern const REAL b(const REAL_D);
  extern REAL       c(const REAL_D), r(const REAL_D);

  x = coord_to_world(el_info, quad->lambda[iq], nil);
  bx = b(x);

  return(SCP_DOW(bx, grd_uh_iq) + c(x)*uh_iq - r(x));
}
\end{verbatim}\ev
As both \code{uh\_iq} and \code{grd\_uh\_iq} are used, the estimator parameter
\code{f\_flag} must be given as \code{INIT\_UH|INIT\_GRD\_UH}.
\end{example}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection{Error estimator for parabolic problems}\label{S:para_est}

Similar to the stationary case, the \ALBERTA library provides
an error estimator for the non--linear parabolic problem
\begin{align*}
 \partial_t u -\nabla \cdot A \nabla u(x)
      + f\Big(x, t, u(x),\nabla u(x)\Big) &= 0    &&x \in \Omega, t>0,\\
                      u(x,t)  &= g_d   && x \in \Gamma_D, t>0,\\
     \nu\cdot A\nabla u(x,t)  &= g_n   && x \in \Gamma_N, t>0,\\
                      u(x,0)  &= u_0 && x \in \Omega,
\end{align*}
where $A \in \R^{d\times d}$ is a positive definite matrix and
$\partial\Omega=\Gamma_D\cup\Gamma_N$. The estimator is split in
several parts, where the initial error
\[
\eta_0 = \|u_0 - U_0\|_{L^2(\Omega)}
\]
can be approximated by the function \code{L2\_err()}, e.g. (compare 
\secref{S:error_cal}).

For the estimation of the spatial discretization error, the coarsening
error, and the time discretization error, the \ALBERTA estimator
implements the following (local) indicators
\begin{align*}
 \eta_S^2  &= 
      C_0^2\, h_S^4\, \left\| \frac{U_{n+1} - I_{n+1}U_{n}}{\tau_{n+1}}
      -\nabla \cdot A \nabla U_{n+1} + f(.,t_{n+1},U_{n+1},\nabla U_{n+1})
      \right\|_{L^2(S)}^2\\
      &\qquad + C_1^2\, \sum_{\Gamma\subset\partial S \cap \Omega} 
                h_S^3\, \|\jump{A\nabla U_{n+1}}\|_{L^2(\Gamma)}^2
          + C_1^2\, \sum_{\Gamma\subset\partial S \cap \Gamma_N} 
                h_S^3\, \|\nu\cdot A\nabla U_{n+1}-g_n\|_{L^2(\Gamma)}^2,
 \\[2mm]
 \eta_{S,c}^2  &= C_2^2\,
    h_S^3\, \|\jump{\nabla U_{n}}\|_{L^2(\Gamma_c)}^2
% \|U_{n} - I_{n+1}U_{n}\|_{L^2(S)}^2
 \\[2mm]
 \eta_\tau &= C_3 \|U_{n+1} - I_{n+1}U_{n}\|_{L^2(\Omega)}.
\end{align*}
The coarsening indicator is motivated by the fact that for piecewise
linear Lagrange finite element functions it holds $\|U_{n} -
I_{n+1}U_{n}\|_{L^2(S)}^2 = \eta_{S,c}^2$ with $C_2=C_2(d)$ and
$\Gamma_c$ the face that would be removed during a coarsening operation.
%%
The implementation is done by the functions
\fdx{heat_est()@{\code{heat\_est()}}}%
\fdx{heat_est_dow()@{\code{heat\_est\_dow()}}}%
\fdx{heat_est_d()@{\code{heat\_est\_d()}}}%
\idx{error estimators!heat_est()@{\code{heat\_est()}}}%
\idx{error estimators!heat_est_dow()@{\code{heat\_est\_dow()}}}%
\idx{error estimators!heat_est_d()@{\code{heat\_est\_d()}}}%
\bv\begin{verbatim}
REAL heat_est(const DOF_REAL_VEC *uh, ADAPT_INSTAT *adapt,
              REAL *(*rw_est)(EL *), REAL *(*rw_estc)(EL *),
              int quad_degree, REAL C[4], const DOF_REAL_VEC *uh_old,
              const REAL_DD A, const BNDRY_FLAGS dirichlet_bndry,
              REAL (*f)(const EL_INFO *el_info, const QUAD *quad, int qp,
                        REAL uh_qp, const REAL_D grd_uh_gp, REAL time),
              FLAGS f_flags,
              REAL (*gn)(const EL_INFO *el_info, const QUAD *quad, int qp,
                         REAL uh_qp, const REAL_D normal, REAL time),
              FLAGS gn_flags);
REAL heat_est_dow(const DOF_REAL_D_VEC *uh, ADAPT_INSTAT *adapt,
                  REAL *(*rw_est)(EL *), REAL *(*rw_estc)(EL *),
                  int quad_degree, REAL C[4], const DOF_REAL_D_VEC *uh_old,
                  const void *A, MATENT_TYPE A_type, MATENT_TYPE A_blocktype,
                  bool sym_grad,
                  BNDRY_FLAGS dirichlet_bndry,
                  const REAL *(*f)(REAL_D result,
                                   const EL_INFO *el_info,
                                   const QUAD *quad, int qp,
                                   const REAL_D uh_qp,
                                   const REAL_DD grd_uh_gp,
                                   REAL time),
                  FLAGS f_flags,
                  const REAL *(*gn)(REAL_D result,
                                    const EL_INFO *el_info,
                                    const QUAD *quad, int qp,
                                    const REAL_D uh_qp,
                                    const REAL_D normal,
                                    REAL time),
                  FLAGS gn_flags);
REAL heat_est_d(const DOF_REAL_D_VEC *uh,
		const DOF_REAL_D_VEC *uh_old,
		ADAPT_INSTAT *adapt,
		REAL *(*rw_est)(EL *),
		REAL *(*rw_estc)(EL *),
		int quad_degree,
		REAL C[4],
		const void *A,
		MATENT_TYPE A_type,
		MATENT_TYPE A_blocktype,
		bool sym_grad,
		const BNDRY_FLAGS dirichlet_bndry,
		const REAL *(*f)(REAL_D result,
				 const EL_INFO *el_info,
				 const QUAD *quad,
				 int qp,
				 const REAL_D uh_qp,
				 const REAL_DD grd_uh_gp,
				 REAL time),
		FLAGS f_flags,
		const REAL *(*gn)(REAL_D result,
				  const EL_INFO *el_info,
				  const QUAD *quad,
				  int qp,
				  const REAL_D uh_qp,
				  const REAL_D normal,
				  REAL time),
		FLAGS gn_flags);
\end{verbatim}\ev
Description:
\begin{descr}
\kitem{heat\_est(uh, adapt, rw\_el\_est, rw\_el\_estc, degree, C, uh\_old,}
\kitem{~~~~~~~~~A, dirichlet\_bndry, f, f\_flag, gn, gn\_flag)}\hfill

  computes an error estimate of the above type, the local and
  global space discretization estimators are stored in
  \code{adapt->adapt\_space} and via the \code{rw\_...} pointers;
  the return value is the time discretization indicator $\eta_\tau$.

  \begin{descr}
  \kitem{uh} is a vector storing the coefficients of the
    discrete solution $U_{n+1}$; if \code{uh} is a \nil pointer,
    nothing is done, the return value is \code{0.0}.
    
  \kitem{adapt} is a pointer to an \code{ADAPT\_INSTAT}
    structure; if it is not \nil, then the entries
    \code{adapt\_space->p=2}, \code{adapt\_space->err\_sum} and
    \code{adapt\_space->err\_max} of \code{adapt} are set by
    \code{heat\_est()} (compare
    \secref{S:adapt_stat_in_ALBERTA}).

  \kitem{rw\_el\_est} is a function for writing the local error
    indicator $\eta_S^2$ for a single element (usually to some
    location inside \code{leaf\_data}, compare
    \secref{S:leaf_data_info}); if this function is \nil, only
    the global estimate is computed, no local indicators are
    stored.  \code{rw\_el\_est(el)} returns for each leaf element
    \code{el} a pointer to a \code{REAL} for storing the square
    of the element indicator, which can directly be used in the
    adaptive method, compare the \code{get\_el\_est()} function
    pointer in the \code{ADAPT\_STAT} structure (compare
    \secref{S:adapt_stat_in_ALBERTA}).

  \kitem{rw\_el\_estc} is a function for writing the local
    coarsening error indicator $\eta_{S,c}^2$ for a single
    element (usually to some location inside \code{leaf\_data},
    compare \secref{S:leaf_data_info}); if this function is \nil,
    no coarsening error indicators are computed and stored;
    \code{rw\_el\_estc(el)} returns for each leaf element
    \code{el} a pointer to a \code{REAL} for storing the square
    of the element coarsening error indicator. The coarsening
    indicator is not used at the moment.

  \kitem{degree} is the degree of the quadrature that should be
    used for the approximation of the norms on the elements and
    edges/faces; if \code{degree} is less than zero a quadrature
    which is exact of degree
    \code{2*uh->fe\_space->bas\_fcts->degree} is used.

  \kitem{C[0]}, \code{C[1]}, \code{C[2]}, \code{C[3]} are the
    constants in front of the element residual, wall residual,
    coarsening term, and time residual, respectively. If \code{C}
    is \nil, then all constants are set to $1.0$.

  \kitem{uh\_old} is a vector storing the coefficients of the
    discrete solution $U_{n}$ from previous time step; if
    \code{uh\_old} is a \nil pointer, nothing is done, the return
    value is \code{0.0}.

  \kitem{A} is the constant matrix of the second order term.

  \kitem{dirichlet\_bndry} A bit mask marking those parts of the
    boundary which are subject to Dirichlet boundary conditions. See
    \secref{S:boundary}.

  \kitem{f} is a pointer to a function for the evaluation of the
    lower order terms at all quadrature nodes, i.e.
    $f(x(\lambda), t, u(\lambda), \nabla u(\lambda))$ ; if
    \code{f} is a \nil pointer, $f\equiv0$ is assumed;

    \code{f(el\_info, quad, iq, t, uh\_iq, grd\_uh\_iq)} returns
    the value of the lower oder terms on element
    \code{el\_info->el} at the quadrature node
    \code{quad->lambda[iq]}, where \code{uh\_iq} is the value and
    \code{grd\_uh\_iq} the gradient (with respect to the world
    coordinates) of the discrete solution at that quadrature
    node.

  \kitem{f\_flag} specifies whether the function \code{f()}
    actually needs values of \code{uh\_iq} or \code{grd\_uh\_iq}.
    This flag may hold zero, the predefined values
    \code{INIT\_UH} or \code{INIT\_GRD\_UH}, or their composition
    \code{INIT\_UH|INIT\_GRD\_UH}; the arguments \code{uh\_iq}
    and \code{grd\_uh\_iq} of \code{f()} only hold valid
    information, if the flags \code{INIT\_UH} respectively
    \code{INIT\_GRD\_UH} are set.

  \kitem{gn(el\_info, quad, qp, uh\_qp, normal)} is a pointer to
    a function for the evaluation of non-homogeneous Neumann
    boundary data.  \code{gn} may be \nil, in which case zero
    Neumann boundary conditions are assumed.  The argument
    \code{normal} always contains the normal of the Neumann
    boundary facet. In the case of non-vanishing co-dimension
    \code{normal} lies in the lower-dimensional space which is
    spanned by the mesh simplex defined by \code{el\_info}.

  \kitem{gn\_flag} controls whether the argument \code{uh\_qp}
    of the function \code{gn()} actually contains the value of
    \code{uh} at the quadrature point \code{qp}. Note that the
    argument \code{normal} always contains valid data.
  \end{descr}
  \smallskip

  The estimate is computed by traversing all leaf elements of
  \code{uh->fe\_space->mesh}, using the quadrature for the
  approximation of the residuals and storing the square of the
  element indicators on the elements (if \code{rw\_el\_est} and
  \code{rw\_el\_estc }are not \nil).

\kitem{heat\_est\_d(uh, adapt, rw, rwc, deg, C, uh\_old,}
\kitem{~~~~~~~~~~A, A\_type, A\_blocktype, sym\_grad,}
\kitem{~~~~~~~~~~dirichlet\_bndry, f, f\_flag)}\hfill
\kitem{heat\_est\_dow(uh, adapt, rw, rwc, deg, C, uh\_old,}
\kitem{~~~~~~~~~~~~A, A\_type, A\_blocktype, sym\_grad,}
\kitem{~~~~~~~~~~~~dirichlet\_bndry, f, f\_flag)}\hfill

  Coupled vector valued version. See \code{ellipt\_est\_dow()} above.
\end{descr}

There are also some less high-level support functions which allow for
custom contributions to the per-element error estimates. We will not
document this in detail, but rather refer the reader to the
\code{stokes.c} and \code{quasi-stokes.c} demo-programs.
%%
\fdx{ellipt_est_init()@{\code{ellipt\_est\_init()}}}
\fdx{heat_est_init()@{\code{heat\_est\_init()}}}
\fdx{element_est()@{\code{element\_est()}}}
\fdx{element_est_finish()@{\code{element\_est\_finish()}}}
\fdx{element_est_uh()@{\code{element\_est\_uh()}}}
\fdx{element_est_grd_uh()@{\code{element\_est\_grd\_uh()}}}
\fdx{ellipt_est_finish()@{\code{ellipt\_est\_finish()}}}
\fdx{heat_est_finish()@{\code{heat\_est\_finish()}}}
%%
\fdx{ellipt_est_dow_init()@{\code{ellipt\_est\_dow\_init()}}}
\fdx{heat_est_dow_init()@{\code{heat\_est\_dow\_init()}}}
\fdx{element_est_dow()@{\code{element\_est\_dow()}}}
\fdx{element_est_dow_finish()@{\code{element\_est\_dow\_finish()}}}
\fdx{element_est_uh_dow()@{\code{element\_est\_uh\_dow()}}}
\fdx{element_est_grd_uh_dow()@{\code{element\_est\_grd\_uh\_dow()}}}
\fdx{ellipt_est_dow_finish()@{\code{ellipt\_est\_dow\_finish()}}}
\fdx{heat_est_dow_finish()@{\code{heat\_est\_dow\_finish()}}}
%%
\bv\begin{lstlisting}
const void *ellipt_est_init(const DOF_REAL_VEC *uh,
                            ADAPT_STAT *adapt,
                            REAL *(*rw_est)(EL *),
                            REAL *(*rw_estc)(EL *),
                            const QUAD *quad,
                            const WALL_QUAD *wall_quad,
                            NORM norm,
                            REAL C[3],
                            const REAL_DD A,
                            const BNDRY_FLAGS dirichlet_bndry,
                            REAL (*f)(const EL_INFO *el_info,
                                      const QUAD *quad,
                                      int qp,
                                      REAL uh_qp,
                                      const REAL_D grd_uh_gp),
                            FLAGS f_flags,
                            REAL (*gn)(const EL_INFO *el_info,
                                       const QUAD *quad,
                                       int qp,
                                       REAL uh_qp,
                                       const REAL_D normal),
                            FLAGS gn_flags);
const void *heat_est_init(const DOF_REAL_VEC *uh,
                          const DOF_REAL_VEC *uh_old,
                          ADAPT_INSTAT *adapt,
                          REAL *(*rw_est)(EL *),
                          REAL *(*rw_estc)(EL *),
                          const QUAD *quad,
                          const WALL_QUAD *wall_quad,
                          REAL C[4],
                          const REAL_DD A,
                          const BNDRY_FLAGS dirichlet_bndry,
                          REAL (*f)(const EL_INFO *el_info,
                                    const QUAD *quad,
                                    int qp,
                                    REAL uh_qp,
                                    const REAL_D grd_uh_gp,
                                    REAL time),
                          FLAGS f_flags,
                          REAL (*gn)(const EL_INFO *el_info,
                                     const QUAD *quad,
                                     int qp,
                                     REAL uh_qp,
                                     const REAL_D normal,
                                     REAL time),
                          FLAGS gn_flags);

REAL element_est(const EL_INFO *el_info, const void *est_handle);
void element_est_finish(const EL_INFO *el_info,
                        REAL est_el, const void *est_handle);
const REAL *element_est_uh(const void *est_handle);
const REAL_D *element_est_grd_uh(const void *est_handle);
REAL ellipt_est_finish(ADAPT_STAT *adapt, const void *est_handle);
REAL heat_est_finish(ADAPT_INSTAT *adapt, const void *est_handle);
\end{lstlisting}\ev
%%
There are similar proto-types for the vector-valued case. Now, what
are these functions good for? The \code{stokes.c} program makes use of
this framework to add a contribution concerning the divergence
constraint. Of course, this is an ad-hoc error indicator, and only
meant to demonstrate the programming frame-work. The functions
\code{element\_est\_uh[\_dow]()} and
\code{element\_est\_grd\_uh[\_dow]()} give the application access to
the values of the discrete solution at the quadrature points
(respectively to its Jaocbians). Otherwise, the general layout is like
follows:
%%
\bv\begin{lstlisting}
void *est_handle = ellipt_est_init(...);
TRAVERSE_FIRST(mesh, -1, <suitable fill-flags>) {
  REAL est_el = element_est(el_info, est_handle);

  ... /* add whatever you like to est_el */

  element_est_finish(el_info, est_el, est_handle);
} TRAVERSE_NEXT();
REAL est = ellipt_est_finish(adapt, est_handle);
\end{lstlisting}\ev
%%
The relevant excerpt from \code{stokes.c} reads as follows:
%%
\bv\begin{lstlisting}
  est_handle = ellipt_est_dow_init(u_h, adapt, rw_el_est, NULL /* rw_estc */,
				 quad, NULL /* wall_quad */,
				 H1_NORM, C,
				 A, MATENT_REAL, MATENT_REAL,
				 false /* !sym_grad */,
				 dirichlet_mask,
				 r, INIT_GRD_UH,
				 NULL /* inhomog. Neumann res. */, 0);

  fill_flags = FILL_NEIGH|FILL_COORDS|FILL_OPP_COORDS|FILL_BOUND|CALL_LEAF_EL;
  fill_flags |= u_fe_space->bas_fcts->fill_flags;
  fill_flags |= p_fe_space->bas_fcts->fill_flags;  
  TRAVERSE_FIRST(mesh, -1, fill_flags) {
    const EL_GEOM_CACHE *elgc;
    const QUAD_EL_CACHE *qelc;
    REAL est_el;

    est_el = element_est_dow(el_info, est_handle);

    if (C[3]) {
      REAL div_uh_el, div_uh_qp;
      const REAL_DD *grd_uh_qp;
      int qp, i;

      grd_uh_qp = element_est_grd_uh_d(est_handle);
      div_uh_el = 0.0;
      if (!(el_info->fill_flag & FILL_COORDS)) {
	qelc = fill_quad_el_cache(el_info, quad, FILL_EL_QUAD_DET);

	for (qp = 0; qp < quad->n_points; qp++) {
	  div_uh_qp = 0;
	  for (i = 0; i < DIM_OF_WORLD; i++) {
	    div_uh_qp += grd_uh_qp[qp][i][i];
	  }
	  div_uh_el += qelc->param.det[qp]*quad->w[qp]*SQR(div_uh_qp);
	}
      } else {
	elgc = fill_el_geom_cache(el_info, FILL_EL_DET);

	for (qp = 0; qp < quad->n_points; qp++) {
	  div_uh_qp = 0;
	  for (i = 0; i < DIM_OF_WORLD; i++) {
	    div_uh_qp += grd_uh_qp[qp][i][i];
	  }
	  div_uh_el += quad->w[qp]*SQR(div_uh_qp);
	}
	div_uh_el *= elgc->det;
      }

      est_el += C[3] * div_uh_el;
    }

    element_est_dow_finish(el_info, est_el, est_handle);    
  } TRAVERSE_NEXT();
  est = ellipt_est_dow_finish(adapt, est_handle);
\end{lstlisting}\ev

\idx{error estimators|)}


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