\section{Routines for barycentric coordinates}%
\label{S:bary_routines}

Operations on single elements are performed using barycentric
coordinates. In many applications, the world coordinates $x$ of the
local barycentric coordinates $\lambda$ have to be calculated (see
\secref{S:ass_tools}, e.g.). Some other applications will need
the calculation of barycentric coordinates for given world coordinates
(see \secref{S:traverse}, e.g.). Finally, derivatives of finite
element functions on elements involve the Jacobian of the barycentric
coordinates (see \secref{S:eval}, e.g.).

In case of a grid with parametric elements, these operations strongly
depend on the element parameterization and no general routines can be
supplied. For non-parametric simplices, \ALBERTA supplies
functions to perform these basic tasks:
%
\fdx{coord_to_world()@{\code{coord\_to\_world()}}}%
\idx{barycentric coordinates!coord_to_world()@{\code{coord\_to\_world()}}}%
\fdx{world_to_coord()@{\code{world\_to\_coord()}}}%
\idx{barycentric coordinates!world_to_coord()@{\code{world\_to\_coord()}}}%
\fdx{el_grd_lambda()@{\code{el\_grd\_lambda()}}}%
\idx{barycentric coordinates!el_grd_lambda()@{\code{el\_grd\_lambda()}}}%
\fdx{el_det()@{\code{el\_det()}}}%
\fdx{el_volume()@{\code{el\_volume()}}}%
\fdx{get_wall_normal()@{\code{get\_wall\_normal()}}}%
\bv\begin{lstlisting}
const REAL *coord_to_world(const EL_INFO *, const REAL *, REAL_D);
int world_to_coord(const EL_INFO *, const REAL *, REAL_B);
REAL el_grd_lambda(const EL_INFO *, REAL [N_LAMBDA][DIM_OF_WORLD]);
REAL el_det(const EL_INFO *);
REAL el_volume(const EL_INFO *);
REAL get_wall_normal(const EL_INFO *el_info, int i0, REAL *normal);
\end{lstlisting}\ev
Description:
\begin{descr}
\kitem{coord\_to\_world(el\_info, lambda, world)} returns a pointer
  to a vector, which contains the world coordinates of a point in
  barycentric coordinates \code{lambda} with respect to the element
  \code{el\_info->el};

  if \code{world} is not \nil the world coordinates are stored in this
  vector; otherwise the function itself provides memory for this
  vector; in this case the vector is overwritten during the next call
  of \code{coord\_to\_world()};

  \code{coord\_to\_world()} needs vertex coordinates information; the
  flag \code{FILL\_COORDS} has to be set during mesh traversal when
  calling this routine on elements.

\kitem{world\_to\_coord(el\_info, world, lambda)} calculates the
  barycentric coordinates with respect to the element
  \code{el\_info->el} of a point with world coordinates \code{world}
  and stores them in the vector given by \code{lambda}. The return
  value is \code{-1} when the point is inside the simplex (or on its
  boundary), otherwise the index of the barycentric coordinate with
  largest negative value (between 0 and $d$);

  \code{world\_to\_coord()} needs vertex coordinates information; the
  flag \code{FILL\_COORDS} has to be set during mesh traversal when
  calling this routine on elements.

  Note that -- with the exception of the 1d code -- this function is
  only implemented for the co-dimension $0$ case, i.e.
  \code{mesh->dim} and \DOW have to be equal, otherwise a call to this
  function will terminate the application with a corresponding
  error-message.

\kitem{el\_grd\_lambda(el\_info, Lambda)} calculates the Jacobian of
  the barycentric coordinates on \code{el\_info->el} and stores the
  matrix in \code{Lambda}; the return value of the function is the
  absolute value of the determinant of the affine linear
  parameterization's Jacobian. For $d<n$ the tangential gradient and
  the value of Gram's determinant are calculated.

  \code{el\_grd\_lambda()} needs vertex coordinates information; the
  flag \code{FILL\_COORDS} has to be set during mesh traversal when
  calling this routine on elements.

\kitem{el\_det(el\_info)} returns the the absolute value of the
  determinant of the affine linear parameterization's Jacobian, or
  Gram's determinant for $d<n$.

  \code{el\_det()} needs vertex coordinates information; the flag
  \code{FILL\_COORDS} has to be set during mesh traversal when calling
  this routine on elements.

\kitem{el\_volume(el\_info)} returns the the volume of the simplex;
  \code{el\_volume()} needs vertex coordinates information; the flag
  \code{FILL\_COORDS} has to be set during mesh traversal when calling
  this routine on elements.

\kitem{get\_wall\_normal(el\_info, wall, normal)} compute the outer
  unit normal of the face opposite to vertex \code{wall}. The result
  is stored in \code{normal}. The return value is the ``surface
  element'' of the given face, i.e. Gram's determinant of the
  transformation to the respective face of the reference element.
  \code{normal} may be \nil. In the case of non-zero co-dimension
  \code{normal} is contained in the sub-space spanned by the edges of
  the given simplex.
\end{descr}

\medskip

All functions described above also come with a \code{\dots\_Xd}
variant, e.g. \code{coord\_to\_world\_2d()}. For the case of $0$
co-dimension there are also wrapper functions \code{\dots\_0cd} which
call the appropriate \code{\dots\_Xd} variant with \code{X ==
  DIM\_OF\_WORLD}. The \code{\dots\_Xd} and \code{\dots\_0cd} variants
are very slightly faster because otherwise the dimension of the
underlying mesh has to be read out of the mesh structure -- e.g. via
\code{el\_info->mesh\_dim} -- and only then the functions branch to
the appropriate \code{\dots\_Xd} variant.

\section{Data structures for numerical quadrature}%
\label{S:quad_data}

For the numerical calculation of general integrals 
\[
\int_S f(x)\, dx
\]
we use quadrature formulas described in \ref{book:S:quadrature}.
\ALBERTA supports numerical integration in zero, one, two, and three dimensions
on the standard simplex $\Shat$ in barycentric coordinates.

\subsection{The \code{QUAD} data structure}%
\label{S:QUAD}

A quadrature formula is described by the following structure, which
is defined both as type \code{QUAD} and \code{QUADRATURE}:
\ddx{QUAD@{\code{QUAD}}}
\idx{numerical quadrature!QUAD@{\code{QUAD}}}
\ddx{QUADRATURE@{\code{QUADRATURE}}}
\idx{numerical quadrature!QUADRATURE@{\code{QUADRATURE}}}
\bv\begin{lstlisting}[name=QUAD,label=T:QUAD]
extern n_quad_points_max[DIM_MAX+1];
typedef struct quadrature   QUAD;
typedef struct quadrature   QUADRATURE;

struct quadrature
{
  char          *name;
  int           degree;

  int           dim;
  int           codim;
  int           subsplx;
  
  int           n_points;
  int           n_points_max;
  const REAL_B  *lambda;
  const REAL    *w;

  void          *metadata;

  INIT_ELEMENT_DECL;
};
\end{lstlisting}\ev
\label{T:QUADRATURE}
Description:
\begin{descr}
  \kitem{name} Textual description of the quadrature.
  %% 
  \kitem{degree} Quadrature is exact of degree \code{degree}.
  %% 
  \kitem{dim} Quadrature for dimension \code{dim}. The barycentric
  co-cordinates of the quadrature points always have \code{dim+1} valid
  components.
  %% 
  \kitem{codim} Co-dimension; \code{codim} is always \code{0} for
  quadratures returned by \code{get\_quadrature()}, and \code{1} for
  quadratures returned by \code{get\_wall\_quad()} and
  \code{get\_bndry\_quad()}.
  %%
  \kitem{subsplx} For \code{codim == 1} the number of the wall-simplex
  this quadrature can be used for; this implies that
  \code{lambda[iq][subsplx]} zero.
  %% 
  \kitem{n\_points} The number of quadrature points.
  %% 
  \kitem{n\_points\_max} The maximal number of quadrature points. The
  number of quadrature points can vary from simplex to simplex if
  \code{INIT\_ELEMENT\_METHOD(quad)} is not \nil.
  %% 
  \kitem{lambda} Vector
  $\mbox{\code{lambda[0]}},\dots,\mbox{\code{lambda[n\_points-1]}}$
  of quadrature points given in bary\-centric coordinates
  (thus having \code{N\_LAMBDA\_MAX} components).
  %% 
  \kitem{w} vector
  $\mbox{\code{w[0]}},\dots,\mbox{\code{w[n\_points-1]}}$ of
  quadrature weights.
  %%
  \kitem{metadata} Pointer to an internal data structure for
  per-element quadrature caches and the like, see e.g.
  \secref{S:fill_quad_el_cache}
  %%
  \kitem{INIT\_ELEMENT\_DECL} Function pointer to a per-element
  initializer. This pointer is always \nil for quadratures returned by
  \code{get\_quadrature()}, \code{get\_wall\_quad()} and
  \code{get\_bndry\_quad()}. External extension modules make use of it.
  See \secref{S:init_element}.
\end{descr}
Currently, numerical quadrature formulas exact up to degree 19 in one
(Gauss formulas), up to degree 17 in two, up to degree 7 in three
dimensions are implemented. We only use stable formulas; this results
in more quadrature points for some formulas (for example in 3d the
formula which is exact of degree 3). A compilation of quadrature
formulas on triangles and tetrahedra is given in
\cite{CoolsRabinowitz:93}. The implemented quadrature formulas are
taken from \cite{Dunavant:85,Gatermann:88,Kardestuncer:87,Stroud:71}.

Using a conical product rule it is possible to construct new
(non-symmetric) quadrature formulas from the existing ones, if that is
really needed.

Functions for numerical quadrature are
\fdx{get_quadrature()@{\code{get\_quadrature()}}}
\idx{numerical quadrature!get_quadrature()@{\code{get\_quadrature()}}}
\fdx{integrate_std_simp()@{\code{integrate\_std\_simp()}}}
\idx{numerical quadrature!integrate_std_simp()@{\code{integrate\_std\_simp()}}}
\fdx{get_product_quad()@{\code{get\_product\_quad()}}}
\idx{numerical quadrature!get_product_quad()@{\code{get\_product\_quad()}}}
\fdx{get_lumping_quadrature()@{\code{get\_lumping\_quadrature()}}}
\idx{numerical quadrature!get_lumping_quadrature()@{\code{get\_lumping\_quadrature()}}}
\fdx{register_quadrature()@{\code{register\_quadrature()}}}
\idx{numerical quadrature!register_quadrature()@{\code{register\_quadrature()}}}
\fdx{new_quadrature()@{\code{new\_quadrature()}}}
\idx{numerical quadrature!new_quadrature()@{\code{new\_quadrature()}}}
\bv\begin{lstlisting}
const QUAD *get_quadrature(int dim, int degree);
REAL  integrate_std_simp(const QUAD *quad, REAL (*f)(const REAL *));
const QUAD *get_product_quad(const QUAD *oq);
const QUAD *get_lumping_quadrature(int dim);
void register_quadrature(QUAD *quad);
bool new_quadrature(const QUAD *quad);
\end{lstlisting}\ev
Description:
\begin{descr}
  \kitem{get\_quadrature(dim, degree)} returns a pointer to a 
  \code{QUAD} structure for numerical integration in \code{dim}
  dimensions which is exact of degree 
  $\min(19,\code{degree})$ for \code{dim==1}, 
  $\min(17,\code{degree})$ for \code{dim==2},
  and $\min(7,\code{degree})$ for \code{dim==3}.

  It is possible to extend the maximal degrees by installing an
  application-defined quadrature rule via \code{new\_quadrature()}.
  %%
  \kitem{register\_quadrature(quad)} Equip an application-defined
  quadrature with internal used meta-data for quadrature caches; this
  function also updates \code{n\_quad\_points\_max[quad->dim]}. To
  install \code{quad} as a default quadrature which will returned on
  request by \code{get\_quadrature()} the functions
  \code{new\_quadrature()} has to be called additionally.
  \kitem{new\_quadrature(quad)} Install the given quadrature as new
  default quadrature for its dimension and polynomial degree; this
  means that \code{get\_quadrature()} will return a pointer to
  \code{quad} when called with \code{quad->dim} and \code{quad->degree}.
  %%
  \kitem{get\_product\_quad(quad)} Return a conical product quadrature
  rule. The returned quadrature formula is non-symmetric and works for
  one dimension higher than \code{quad} and is exact of the same
  degree as \code{quad}.  The 3D formula for degree $7$ found in
  \cite{Stroud:71} is of this type, for example.

  Note that \code{get\_product\_quad()} installs the new formula
  calling \code{new\_quadrature()}, so the formula will be available
  through \code{get\_quadrature()}.
  %%
  \kitem{get\_lumping\_quadrature(dim)} Returns a lumping quadrature
  with quadrature nodes at the vertices of the reference simplex.
  %%
  \kitem{integrate\_std\_simp(quad, f)} approximates an integral
  by the numerical quadrature described by \code{quad};

  \code{f} is a pointer to a function to be integrated, evaluated
  in barycentric coordinates; the return value is
  \[
  \sum_{\code{k = 0}}^{\mbox{\tiny\code{quad->n\_points-1}}} 
  \mbox{\code{quad->w[k] * (*f)(quad->lambda[k])}};
  \]
  for the approximation of $\int_S f$ we have to multiply this value
  with $d! |S|$ for a simplex $S$; for a parametric simplex,
  \code{f} should be a pointer to a function which calculates
  $f(\lambda) |\det DF_S(\hat x(\lambda))|$.
\end{descr}


The following functions initialize values and gradients of functions
at the quadrature nodes:
\fdx{f_at_qp()@{\code{f\_at\_qp()}}}%
\fdx{f_d_at_qp()@{\code{f\_d\_at\_qp()}}}%
\fdx{grd_f_at_qp()@{\code{grd\_f\_at\_qp()}}}%
\fdx{grd_f_d_at_qp()@{\code{grd\_f\_d\_at\_qp()}}}%
\fdx{f_loc_at_qp()@{\code{f\_loc\_at\_qp()}}}%
\fdx{f_loc_d_at_qp()@{\code{f\_loc\_d\_at\_qp()}}}%
\fdx{grd_f_loc_at_qp()@{\code{grd\_f\_loc\_at\_qp()}}}%
\fdx{param_grd_f_loc_at_qp()@{\code{param\_grd\_f\_loc\_at\_qp()}}}%
\fdx{grd_f_loc_d_at_qp()@{\code{grd\_f\_loc\_d\_at\_qp()}}}%
\fdx{param_grd_f_loc_d_at_qp()@{\code{param\_grd\_f\_loc\_d\_at\_qp()}}}%
\fdx{fx_at_qp()@{\code{fx\_at\_qp()}}}%
\fdx{fx_d_at_qp()@{\code{fx\_d\_at\_qp()}}}%
\fdx{grd_fx_at_qp()@{\code{grd\_fx\_at\_qp()}}}%
\fdx{grd_fx_d_at_qp()@{\code{grd\_fx\_d\_at\_qp()}}}%
\bv\begin{lstlisting}
REAL *f_at_qp(REAL vec[], const QUAD *quad, REAL (*f)(const REAL_B lambda));
REAL_D *grd_f_at_qp(REAL_D vec[], const QUAD *quad,
                    const REAL *(*f)(const REAL_B));
REAL_D *f_d_at_qp(REAL_D vec[], const QUAD *quad,
                  const REAL *(*f)(const REAL_B lambda));
REAL_DD *grd_f_d_at_qp(REAL_DD vec[], const QUAD *quad,
                       const REAL_D *(*f)(const REAL_B lambda));

REAL *f_loc_at_qp(REAL vec[], const EL_INFO *el_info, const QUAD *quad,
                  REAL (*f)(const EL_INFO *el_info, 
                            const QUAD *quad, int iq, void *ud),
                  void *ud);
REAL_D *grd_f_loc_at_qp(REAL_D vec[], const EL_INFO *el_info,
                        const QUAD *quad, const REAL_BD Lambda,
                        GRD_LOC_FCT_AT_QP grd_f, void *ud);
REAL_D *param_grd_f_loc_at_qp(REAL_D vec[], const EL_INFO *el_info,
                              const QUAD *quad, const REAL_BD Lambda[],
                              GRD_LOC_FCT_AT_QP grd_f, void *ud);
REAL_D *f_loc_d_at_qp(REAL_D vec[], const EL_INFO *el_info, const QUAD *quad,
                      const REAL *(*f)(REAL_D result, const EL_INFO *el_info,
                                       const QUAD *quad, int iq, void *ud),
                      void *ud);
REAL_DD *grd_f_loc_d_at_qp(REAL_DD vec[], const EL_INFO *el_info,
                           const QUAD *quad, const REAL_BD Lambda,
                           GRD_LOC_FCT_D_AT_QP grd_f, void *ud);
REAL_DD *param_grd_f_loc_d_at_qp(REAL_DD vec[], const EL_INFO *el_info,
                                 const QUAD *quad, const REAL_BD Lambda[],
                                 GRD_LOC_FCT_D_AT_QP grd_f, void *ud);

REAL *fx_at_qp(REAL vec[], const EL_INFO *el_info, const QUAD *quad,
               FCT_AT_X f);
REAL_D *grd_fx_at_qp(REAL_D vec[], const EL_INFO *el_info, const QUAD *quad,
                     GRD_FCT_AT_X grd_f);
REAL_D *fx_d_at_qp(REAL_D vec[], const EL_INFO *el_info, const QUAD *quad,
                   FCT_D_AT_X f);
REAL_DD *grd_fx_d_at_qp(REAL_DD vec[], const EL_INFO *el_info,
                        const QUAD *quad, GRD_FCT_D_AT_X grd_f);


\end{lstlisting}\ev
Description:
\def\DOW{\code{DIM\_OF\_WORLD}\xspace}
\def\div{\mathrm{div\,}}
\begin{descr}
\kitem{f\_at\_qp(vec, quad, f)} returns a pointer \code{ptr} to a
       vector \code{quad->n\_points} storing the values of a \code{REAL}
       valued function at all quadrature points of \code{quad};
       \code{f} is a pointer to that function, evaluated in
       barycentric coordinates; if \code{vec} is not \nil, the values
       are stored in this vector, otherwise the values are stored in
       some static local vector, which is overwritten on the next call;

       \code{ptr[i]}$ = $\code{(*f)(quad->lambda[i])} for
       $0\leq\code{i}<\mbox{\code{quad->n\_points}}$.
\kitem{grd\_f\_at\_qp(vec, quad, grd\_f)} returns a pointer \code{ptr} to a
       vector \code{quad->n\_points} storing the gradient (with
       respect to world coordinates) of a \code{REAL} valued function
       at all quadrature points of \code{quad};

       \code{grd\_f} is a pointer to a function, evaluated in
       barycentric coordinates and returning a pointer to a vector of
       length \DOW storing the gradient;

       if \code{vec} is not \nil, the values are stored in this
       vector, otherwise the values are stored in some local static
       vector, which is overwritten on the next call;

       \code{ptr[i][j]}$ = $\code{(*grd\_f)(quad->lambda[i])[j]}, for 
       $0\leq\code{j}<\mbox{\DOW}$ and
       $0\leq\code{i}<\mbox{\code{quad->n\_points}}$,
\kitem{f\_d\_at\_qp(vec, quad, fd)} returns a pointer \code{ptr} to a
       vector \code{quad->n\_points} storing the values of a
       \code{REAL\_D} valued function at all quadrature points of
       \code{quad}; 

       \code{fd} is a pointer to that function, evaluated
       in barycentric coordinates and returning a pointer to a vector
       of length \DOW storing all components; if the second argument
       \code{val} of \code{(*fd)(lambda, val)} is not \nil, the values
       have to be stored at \code{val}, otherwise \code{fd} has to
       provide memory for the vector which may be overwritten on the
       next call; 

       if \code{vec} is not \nil, the values are stored in
       this vector, otherwise the values are stored in some static
       local vector, which is overwritten on the next call;

       \code{ptr[i][j]}$ = $\code{(*fd)(quad->lambda[i],val)[j]}, for 
       $0\leq\code{j}<\mbox{\DOW}$ and
       $0\leq\code{i}<\mbox{\code{quad->n\_points}}$.
\kitem{grd\_f\_d\_at\_qp(vec, quad, grd\_fd)} returns a pointer 
       \code{ptr} to a vector \code{quad->n\_points} storing the
       Jacobian (with respect to world coordinates) of a
       \code{REAL\_D} valued function at all quadrature points of
       \code{quad}; 

       \code{grd\_fd} is a pointer to a function, evaluated in
       barycentric coordinates and returning a pointer to a matrix of
       size \DOW$\times$\DOW storing the Jacobian; if the second
       argument \code{val} of \code{(*grd\_fd)(x, val)} is not \nil, the
       Jacobian has to be stored at \code{val}, otherwise \code{grd\_fd} has
       to provide memory for the matrix which may be overwritten on
       the next call;

       if \code{vec} is not \nil, the values are stored in this
       vector, otherwise the values are stored in some static local
       vector, which is overwritten on the next call;

       \code{ptr[i][j][k]}$ = $\code{(*grd\_fd)(quad->lambda[i],val)[j][k]},
       for $0\leq\code{j},\code{k}<\mbox{\DOW}$ and
       $0\leq\code{i}<\mbox{\code{quad->n\_points}}$, 

\kitem{f\_loc\_at\_qp(vec, el\_info, quad, f, ud)}
\kitem{grd\_f\_loc\_at\_qp(vec, el\_info, quad, Lambda, grd\_f, ud)}
\kitem{param\_grd\_f\_loc\_at\_qp(vec, el\_info, quad, Lambda, grd\_f, ud)}
\kitem{f\_loc\_d\_at\_qp(vec, el\_info, quad, fd, ud)}
\kitem{grd\_f\_loc\_d\_at\_qp(vec, el\_info, quad, Lambda, grd\_fd, ud)}
\kitem{param\_grd\_f\_loc\_d\_at\_qp(vec, el\_info, quad, Lambda, grd\_fd, ud)}

\kitem{fx\_at\_qp(vec, el\_info, quad, f)}
\kitem{grd\_fx\_at\_qp(vec, el\_info, quad, grd\_f)}
\kitem{fx\_d\_at\_qp(vec, el\_info, quad, fd)}
\kitem{grd\_fx\_d\_at\_qp(vec, el\_info, quad, grd\_fd)}
\end{descr}

\subsection{The \code{QUAD\_FAST} data structure}
\label{S:QUAD_FAST}

Often numerical integration involves basis functions, such as the
assembling of the system matrix and right hand side, or the
integration of finite element functions.  Since numerical quadrature
involves only the values at the quadrature points and the values of
basis functions and its derivatives (with respect to barycentric
coordinates) are the same at these points for
all elements of the grid, such routines can be much more efficient, if
they can use pre--computed values of the basis functions at the
quadrature points. In this case the basis functions do not have to be
evaluated for each quadrature point newly on every element.

Information that should be pre--computed can be specified by the
following symbolic constants:
\cdx{INIT_PHI@{\code{INIT\_PHI}}}%
\idx{numerical quadrature!INIT_PHI@{\code{INIT\_PHI}}}%
\cdx{INIT_GRD_PHI@{\code{INIT\_GRD\_PHI}}}%
\idx{numerical quadrature!INIT_GRD_PHI@{\code{INIT\_GRD\_PHI}}}%
\cdx{INIT_D2_PHI@{\code{INIT\_D2\_PHI}}}%
\idx{numerical quadrature!INIT_D2_PHI@{\code{INIT\_D2\_PHI}}}%
\bv\begin{lstlisting}[label=M:QFAST_FLAGS]
#define INIT_PHI        0x01
#define INIT_GRD_PHI    0x02
#define INIT_D2_PHI     0x04
#define INIT_D3_PHI     0x08
#define INIT_D4_PHI     0x10
#define INIT_TANGENTIAL 0x80
\end{lstlisting}\ev
Description:
\begin{descr}
\kitem{INIT\_PHI} pre--compute the values of all basis functions at
       all quadrature nodes;
\kitem{INIT\_GRD\_PHI} pre--compute the gradients (with respect to
       the barycentric coordinates) of all basis functions at
       all quadrature nodes;
\kitem{INIT\_D2\_PHI}pre--compute all 2nd derivatives (with respect to
       the barycentric coordinates) of all basis functions at
       all quadrature nodes.
\end{descr}
In order to store such information for one set of basis functions
we define the data structure
\ddx{QUAD_FAST@{\code{QUAD\_FAST}}}%
\idx{numerical quadrature!QUAD_FAST@{\code{QUAD\_FAST}}}%
\bv\begin{lstlisting}[name=QUAD_FAST,label=T:QUAD_FAST]
typedef struct quad_fast QUAD_FAST;

struct quad_fast
{
  const QUAD      *quad;
  const BAS_FCTS  *bas_fcts;

  FLAGS           init_flag;

  int             dim;
  int             n_points;
  int             n_bas_fcts;
  int             n_points_max;
  int             n_bas_fcts_max;
  const REAL      *w;               /* shallow copy of quad->w */
  const REAL      (*const*phi);     /* [qp][bf] */
  const REAL_B    (*const*grd_phi);
  const REAL_BB   (*const*D2_phi);
  const REAL_BBB  (*const*D3_phi);
  const REAL_BBBB (*const*D4_phi);

  /* For vector valued basis functions with a p.w. constant
   * directional derivative we cache that direction and make it
   * available for applications. The component is initialized by the
   * INIT_ELEMENT() method.
   *
   * So: phi_d[i] gives the value of the directional factor for the
   * i-th basis function. If (!bas_fcts->dir_pw_const), then phi_d is
   * NULL.
   */
  const REAL_D *phi_d;

  /* chain to next structure, if bas_fcts->chain is non-empty */
  DBL_LIST_NODE  chain;

  /* a clone of this structure, but as single item. */
  const QUAD_FAST *unchained;

  INIT_ELEMENT_DECL;

  void *internal;
};
\end{lstlisting}\ev
The entries yield following information:
\begin{descr}
  \kitem{quad} Values stored for numerical quadrature \code{quad}.
  %% 
  \kitem{bas\_fcts} Values stored for basis functions \code{bas\_fcts}.
  %% 
  \kitem{dim} Clone of \code{quad->dim}.
  %% 
  \kitem{init\_flag} Indicates which information is initialized; may be
  one of, or a bitwise \textsf{OR} of several of \code{INIT\_PHI},
  \code{INIT\_GRD\_PHI}, \code{INIT\_D2\_PHI}, \code{INIT\_D3\_PHI} or
  \code{INIT\_D4\_PHI}. Not all basis functions have support for higher
  derivatives. There is one additional fill-flag,
  \code{INIT\_TANGENTIAL} with the meaning that only the tangential
  derivatives of the basis functions will be computed if \code{quad} is
  a co-dimension $1$ quadrature rule.
  %% 
  \kitem{n\_points} The number of quadrature points; equals
  \code{quad->n\_points}.
  %% 
  \kitem{n\_bas\_fcts} number of basis functions; equals
  \code{bas\_fcts->n\_bas\_fcts}.
  %% 
  \kitem{n\_points\_max} The maximum number of quadrature points. If
  \code{quad->init\_element()} is non-\nil, then the number of basis
  functions can vary on a per-element basis.
  %% 
  \kitem{n\_bas\_fcts\_max} The maximum
  number of basis functions. If \code{bas\_fcts->init\_element} is
  non-\nil, then the number of basis functions can vary on a per-element
  basis.
  %% 
  \kitem{w} Vector of quadrature weights; \code{w = quad->w}.
  %% 
  \kitem{phi} Matrix storing function values if the flag 
  \code{INIT\_PHI} is set.

  \code{phi[i][j]} stores the value
  \code{bas\_fcts->phi[j](quad->lambda[i])},
  $0\leq\code{j}<\mbox{\code{n\_bas\_fcts}}$ and
  $0\leq\code{i}<\mbox{\code{n\_points}}$;
  %%
  \kitem{grd\_phi} Matrix storing all gradients (with respect to the
  barycentric coordinates) if the flag
  \code{INIT\_GRD\_PHI} is set;
  
  \code{grd\_phi[i][j][k]} Stores the value 
  \code{bas\_fcts->grd\_phi[j](quad->lambda[i])[k]} for
  $0\leq\code{j}<\mbox{\code{n\_bas\_fcts}}$,
  $0\leq\code{i}<\mbox{\code{n\_points}}$, and
  $0\leq\code{k}\leq d$;
  %% 
  \kitem{D2\_phi} Matrix storing all second derivatives (with respect to the
  barycentric coordinates) if the flag
  \code{INIT\_D2\_PHI} is set; 

  \code{D2\_phi[i][j][k][l]}  Stores the value 
  \code{bas\_fcts->D2\_phi[j](quad->lambda[i])[k][l]} for
  $0\leq\code{j}<\mbox{\code{n\_bas\_fcts}}$,
  $0\leq\code{i}<\mbox{\code{n\_points}}$, and
  $0\leq\code{k,l}\leq d$.
  %%
  \kitem{D3\_phi} Matrix storing all third derivatives (with respect
  to the barycentric coordinates) if the flag \code{INIT\_D3\_PHI} is
  set;

  \code{D3\_phi[i][j][k][l]}  Stores the value 
  \code{bas\_fcts->D3\_phi[j](quad->lambda[i])[k][l][m]} for
  $0\leq\code{j}<\mbox{\code{n\_bas\_fcts}}$,
  $0\leq\code{i}<\mbox{\code{n\_points}}$, and
  $0\leq\code{k,l,m}\leq d$.
  %%
  %%
  \kitem{D4\_phi} Matrix storing all fourth derivatives (with respect
  to the barycentric coordinates) if the flag \code{INIT\_D4\_PHI} is
  set;

  \code{D4\_phi[i][j][k][l]}  Stores the value 
  \code{bas\_fcts->D4\_phi[j](quad->lambda[i])[k][l][m][n]} for
  $0\leq\code{j}<\mbox{\code{n\_bas\_fcts}}$,
  $0\leq\code{i}<\mbox{\code{n\_points}}$, and
  $0\leq\code{k,l,m,n}\leq d$.
  %%
  \kitem{phi\_d} The directional part of vector-valued basis
  functions, if that is constant on each element. This means, if
  \code{bas\_fcts->rdim == \DOW} and \code{bas\_fcts->dir\_pw\_const},
  then \code{phi\_d} contains valid data, probably after calling
  \code{QUAD\_FAST.init\_element()} with the current element and the
  instance of the \code{QUAD\_FAST} structure in question. See also
  \secref{S:vector_bfcts}.
  %%
  \kitem{chain} If \code{bas\_fcts} forms part of a chain of basis
  functions because the corresponding finite element space is a direct
  sum, then this code{get\_quad\_fast()} will also generate a chain of
  \code{QUAD\_FAST}-structures, one for each component. The chain
  forms a doubly linked list, and the \code{chain}-component is the
  list node. See also \secref{S:bfcts_chains} and
  \secref{S:chain_impl}.
  %%
  \kitem{unchained} A clone of the current structure, but as single
  element. Points back to the structure itself if the underlying basis
  functions do not form part of chain of basis function sets.  See
  \secref{S:bfcts_chains} and \secref{S:chain_impl}.
  %%
  \kitem{INIT\_ELEMENT\_DECL} Per element initializer, see
  \secref{S:init_element}.
  %%
  \kitem{internal} Pointer to internal meta-data stuff.
\end{descr}
%
A filled structure can be accessed by a call of 
\fdx{get_quad_fast()@{\code{get\_quad\_fast()}}}%
\idx{numerical quadrature!get_quad_fast()@{\code{get\_quad\_fast()}}}%
\bv\begin{lstlisting}
const QUAD_FAST *get_quad_fast(const BAS_FCTS *, const QUAD *, U_CHAR);
\end{lstlisting}\ev
Description:
\begin{descr}
\kitem{get\_quad\_fast(bas\_fcts, quad, init\_flag)} \code{bas\_fcts}
       is a pointer to a filled \code{BAS\_FCTS} structure,
       \code{quad} a pointer to some quadrature (accessed by
       \code{get\_quadrature()}, e.g.)  and \code{init\_flag}
       indicates which information should be filled into the
       \code{QUAD\_FAST} structure; it may be one of, or a bitwise
       \textsf{OR} of several of \code{INIT\_PHI}, \code{INIT\_GRD\_PHI},
       \code{INIT\_D2\_PHI}; the function returns a pointer to a
       filled \code{QUAD\_FAST} structure where all demanded information
       is computed and stored.

       All used \code{QUAD\_FAST} structures are stored in a linked
       list and are identified uniquely by the members \code{quad} and
       \code{bas\_fcts}; first, \code{get\_quad\_fast()} looks for a
       matching structure in the linked list; if no structure is
       found, a new structure is generated and linked to the list;
       thus for one combination \code{bas\_fcts} and \code{quad} only
       one \code{QUAD\_FAST} structure is created.

       Then \code{get\_quad\_fast()} allocates memory for all
       information demanded by \code{init\_flag} and which is not yet
       initialized for this structure; only such information
       is then computed and stored; on the first call for 
       \code{bas\_fcts} and \code{quad}, all information demanded
       \code{init\_flag} is generated, on a subsequent call only
       missing information is generated.

       \code{get\_quad\_fast()} will return a \nil pointer, if
       \code{INIT\_PHI} flag is set and \code{bas\_fcts->phi} is \nil,
       \code{INIT\_GRD\_PHI} flag is set and
       \code{bas\_fcts->grd\_phi} is \nil, and \code{INIT\_D2\_PHI}
       flag is set and \code{bas\_fcts->D2\_phi} is \nil.

       There may be several \code{QUAD\_FAST} structures in the list for
       the same set of basis functions for different quadratures,
       and there may be several \code{QUAD\_FAST} structures for one 
       quadrature for different sets of basis functions.

       The function \code{get\_quad\_fast()} should not be called on
       each element during mesh traversal, because it has to look in a
       list for an existing entry for a set of basis functions and a
       quadrature; a pointer to the \code{QUAD\_FAST} structure should
       be accessed before mesh traversal and passed to the 
       element routine.
\end{descr}
Many functions using the \code{QUAD\_FAST} structure need vectors
for storing values at all quadrature points; for these functions
it can be of interest to get the count of the maximal number of
quadrature nodes used by the all initialized \code{quad\_fast}
structures in order to avoid several memory reallocations. This
count can be accessed by the function
\fdx{max_quad_points()@{\code{max\_quad\_points()}}}%
\idx{numerical quadrature!max_quad_points()@{\code{max\_quad\_points()}}}%
\bv\begin{lstlisting}
int max_quad_points(void);
\end{lstlisting}\ev
Description:
\begin{descr}
\kitem{max\_quad\_points()} returns the maximal number of quadrature
       points for all yet initialized \code{quad\_fast} structures;
       this value may change after a new initialization of a 
       \code{quad\_fast} structures;

       this count is \emph{not} the maximal number of quadrature
       points of all used \code{QUAD} structures, since new quadratures
       can be used at any time without an initialization.
\end{descr}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Integration over subsimplices (walls)}%
\label{S:subsimplices}

The weak formulation of non-homogeneous Neumann or Robin boundary
values needs integration over $d-1$ dimensional boundary
simplices of $d$ dimensional mesh elements (compare \secref{book:S:Submeshes}), 
and the evaluation of jump residuals for error estimators (compare Sections
\ref{book:S:adaptive_methods}, \ref{S:estimator}) needs integration over
all interior $d-1$ dimensional sub--simplices.  The quadrature
formulas and data structures described above are available for any $d$
dimensional simplex, $d=0,1,2,3$.  The above task can therefore be
accomplished by using a $d-1$ dimensional quadrature formula
and augmenting the corresponding $d$ dimensional barycentric
coordinates of quadrature points on edges/faces to $d+1$
dimensional coordinates on adjacent mesh elements.

When an integral over an edge/face involves values from both adjacent
elements (in the computation of jump residuals e.\,g.)  it is
necessary to have a common orientation of the edge/face from both
elements. Only a common orientation of the edges/faces ensures that
augmenting $d$ dimensional barycentric coordinates of
quadrature points on the edge/face to $d+1$ dimensional
barycentric coordinates on the adjacent mesh elements results in the
same points from both sides.

This augmentation process, taking the relative orientation of
neighboring simplices into account, is taken care of by dedicated
co-dimension $1$ quadrature rules, see \secref{S:WALL_QUAD} and
\ref{S:WALL_QUAD_FAST}. Additionally, the calculation of Gram's
determinant for the $d-1$ dimensional transformation as well as
vertex/edge/face normals is needed. See \secref{S:bary_routines}
above.

Low-level access to the relative orientation of neighboring simplices
is provided through the routines and look-up tables
%%
\fdx{wall_orientation()@{\code{wall\_orientation()}}}%
\fdx{wall_rel_orientation()@{\code{wall\_rel\_orientation()}}}%
\ddx{sorted_wall_vertices_1d@{\code{sorted\_wall\_vertices\_1d}}}
\ddx{sorted_wall_vertices_2d@{\code{sorted\_wall\_vertices\_2d}}}
\ddx{sorted_wall_vertices_3d@{\code{sorted\_wall\_vertices\_3d}}}
\bv\begin{lstlisting}
int wall_orientation(int dim, const EL *el, int wall, int **vec);
int wall_rel_orientation(
  int dim, const EL *el, const EL *neigh, int wall, int oppv);

const int sorted_wall_vertices_1d[N_WALLS_1D][DIM_FAC_1D][2*N_VERTICES_0D-1];
const int sorted_wall_vertices_2d[N_WALLS_2D][DIM_FAC_2D][2*N_VERTICES_1D-1];
const int sorted_wall_vertices_3d[N_WALLS_3D][DIM_FAC_3D][2*N_VERTICES_2D-1];
\end{lstlisting}\ev
Description:
\begin{descr}
  \kitem{wall\_orientation(dim, el, wall, vec)} can be used to match
  the local enumeration of the vertices of faces separating
  neighbouring simplexes. The return value is a unique number between
  \code{1} and \code{dim!}. On return \code{vec} -- if non-\nil --
  contains a permutation of the local numbering of the vertices of
  face number \code{wall} on \code{el}. If \code{neigh\_vec} is the
  corresponding permutation for the neighbour element, then
  \code{(*vec)[i]} and \code{(*neigh\_vec)[i]} refer to the same
  vertex, e.g.  \code{el\_info->coord[(*vec)[i]]} is the same as
  \code{neigh\_info->coord[(*neigh\_vec)[i]]}.

  Actually, the return value of \code{wall\_orientation()} is just the
  index into the look-up tables
  \code{sorted\_wall\_vertices\_Xd[][][]}, such that \code{vec}, if
  non-\nil point upon return to
  \code{sorted\_wall\_vertices\_Xd[wall][retval]}.

  The principal purpose of this function is to match quadrature points
  during the numerical integration of jumps of derivatives of finite
  element function across the faces of the triangulation, see
  \secref{S:subsimplices}.

  \kitem{wall\_rel\_orientation(dim, el, neigh, wall, oppv)} can be
  used to compute a relative orientation of a given wall separating
  two elements with respect to both elements. The return value
  \bv\begin{lstlisting}
    perm = wall_rel_orientation(dim, el, neigh, wall, oppv);
  \end{lstlisting}\ev
  can be used as an offset into \code{sorted\_wall\_vertices\_Xd} in
  the sense that
  %% 
  \bv\begin{lstlisting}
    nv = sorted_wall_vertices_Xd[oppv][perm][i];
  \end{lstlisting}\ev
  %% 
  matches
  %% 
  \bv\begin{lstlisting}
    v = vertex_of_wall_Xd[wall][i];
  \end{lstlisting}\ev
  %% 
  So it holds \code{el->dof[v][0] == neigh->dof[nv][0]}.
  %% 
\end{descr}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{The \code{WALL\_QUAD} data structure}%
\label{S:WALL_QUAD}

A collection of quadrature rules for the integration over walls (3d:
faces, 2d: edges) of a simplex. The quadrature points of these rules
are given in barycentric coordinates with \code{dim+1} valid
components; the component corresponding to the respective wall will be
set to zero.

Each of the quadrature rules \code{WALL\_QUAD\:\:quad[wall]} may have
its own \code{INIT\_ELEMENT} method.  \code{INIT\_ELEMENT(el\_info,
  WALL\_QUAD)} may or may not be called: it is legal to only call
\code{INIT\_ELEMENT(el\_info, WALL\_QUAD\:\:quad[wall])} individually.
If \code{INIT\_ELEMENT(el\_info, WALL\_QUAD)} is called, then it has
to initialize all quadrature rules for all walls, so the sub-ordinate
initializers need not be called in this case.
%%%
\ddx{WALL_QUAD@{\code{WALL\_QUAD}}}%
\idx{numerical quadrature!WALL_QUAD@{\code{WALL\_QUAD}}}
\bv\begin{lstlisting}[name=WALL_QUAD,label=T:WALL_QUAD]
typedef struct wall_quadrature   WALL_QUAD;

struct wall_quadrature
{
  const char *name;
  int  degree;
  int  dim;
  int  n_points_max;
  QUAD quad[N_WALLS_MAX];

  INIT_ELEMENT_DECL;

  void *metadata;
};
\end{lstlisting}\ev
Description:
\begin{descr}
  \kitem{name} Textual description of the quadrature.
  %% 
  \kitem{degree} Quadrature is exact of degree \code{degree}.
  %% 
  \kitem{dim} Quadrature for dimension \code{dim}; the barycentric
  coordinates of the quadrature points have \code{dim+1} valid
  components.
  %% 
  \kitem{n\_points\_max} The maximal number of quadrature points.
  %% 
  \kitem{quad} Quadrature rules for each wall. These are co-dimension
  $1$ rules.
  %%
  \kitem{INIT\_ELEMENT\_DECL} Function pointer to a per-element
  initializer. This pointer is always \nil for quadratures returned by
  \code{get\_wall\_quad()}. External extension modules may make use of
  it. See \secref{S:init_element}.
  %% 
  \kitem{metadata} Pointer to an internal data structure for
  per-element quadrature caches and the like.
\end{descr}
Functions for numerical quadrature are:
\fdx{get_wall_quad()@{\code{get\_wall\_quad()}}}%
\idx{numerical quadrature!get_wall_quad()@{\code{get\_wall\_quad()}}}%
\fdx{register_wall_quadrature()@{\code{register\_wall\_quadrature()}}}%
\idx{numerical quadrature!register_wall_quadrature()@{\code{register\_wall\_quadrature()}}}%
\bv\begin{lstlisting}
const WALL_QUAD *get_wall_quad(int dim, int degree);
void register_wall_quadrature(WALL_QUAD *wall_quad);
const QUAD *get_neigh_quad(const EL_INFO *el_info, const WALL_QUAD *wall_quad, int neigh);
\end{lstlisting}\ev

Description:
\begin{descr}
  \kitem{get\_wall\_quad(dim, degree)} returns a pointer to a 
  \code{WALL\_QUAD} structure for numerical integration in \code{dim}
  dimensions.
%%
  \kitem{register\_wall\_quadrature(wall\_quad)} initializes the meta-data
  for the given \code{WALL\_QUAD}, no need to call this if the
  \code{WALL\_QUAD} has been acquired by \code{get\_wall\_quad()},
  only needed for externally defined extension quadrature rules.
\end{descr}

\subsection{The \code{WALL\_QUAD\_FAST} data structure}
\label{S:WALL_QUAD_FAST}

Convenience structure for \code{WALL\_QUAD}: its is legal to call
\code{get\_quad\_fast(bas\_fcts, WALL\_QUAD::quad[wall], ...)}
(see \secref{S:QUAD_FAST} ``The \code{QUAD\_FAST} data structure'')
individually, however \code{get\_wall\_quad\_fast()} does this
in a single run. If \code{INIT\_ELEMENT(el\_info, WALL\_QUAD\_FAST)}
is called, then the sub-ordinate initializers
\code{INIT\_ELEMENT(el\_info,WALL\_QUAD\_FAST::quad\_fast[wall])}
need not be called.

\ddx{WALL_QUAD_FAST@{\code{WALL\_QUAD\_FAST}}}
\idx{numerical quadrature!WALL_QUAD_FAST@{\code{WALL\_QUAD\_FAST}}}
%%
\bv\begin{lstlisting}[label=T:WALL_QUAD_FAST]
typedef struct wall_quad_fast  WALL_QUAD_FAST;

struct wall_quad_fast
{
  const WALL_QUAD *wall_quad;
  const BAS_FCTS  *bas_fcts;

  FLAGS           init_flag;
  const QUAD_FAST *quad_fast[N_WALLS_MAX];

  INIT_ELEMENT_DECL;
};
\end{lstlisting}\ev
The entries yield following information:
\begin{descr}
\kitem{wall\_quad} values stored for numerical quadrature \code{quad};
\kitem{bas\_fcts} values stored for basis functions \code{bas\_fcts};
\kitem{init\_flag} indicates which information is initialized; may be
  one of, or a bitwise \textsf{OR} of several of \code{INIT\_PHI},
  \code{INIT\_GRD\_PHI}, \code{INIT\_D2\_PHI};
\kitem{\*quad\_fast[N\_WALLS\_MAX]} Pointer to \code{N\_WALLS\_MAX}
  \code{quad\_fast} structures.
\kitem{INIT\_ELEMENT\_DECL} Function pointer to for a per-element
  initialiser. This pointer is always \nil for quadratures returned by
  \code{get\_quadrature()}, \code{get\_wall\_quad()} and
  \code{get\_bndry\_quad()}. External extension modules make use of it.
  See \secref{S:init_element}.
\end{descr}

\fdx{get_wall_quad_fast()@{\code{get\_wall\_quad\_fast()}}}%
\idx{numerical quadrature!get_wall_quad_fast()@{\code{get\_wall\_quad\_fast()}}}%
\fdx{get_neigh_quad()@{\code{get\_neigh\_quad()}}}%
\idx{numerical quadrature!get_neigh_quad()@{\code{get\_neigh\_quad()}}}%
\fdx{get_neigh_quad_fast()@{\code{get\_neigh\_quad\_fast()}}}%
\idx{numerical quadrature!get_neigh_quad_fast()@{\code{get\_neigh\_quad\_fast()}}}%

\bv\begin{lstlisting}
const WALL_QUAD *get_wall_quad_fast(const BAS_FCTS *, const WALL_QUAD *, FLAGS init_flag);
QUAD_FAST *get_neigh_quad_fast(const EL_INFO *el_info,
				     const WALL_QUAD_FAST *wqfast,
				     int neigh);
\end{lstlisting}\ev
Description:
\begin{descr}
  \kitem{get\_wall\_quad\_fast(bas\_fcts, wall\_quad, init\_flag)}
  \code{bas\_fcts}
  is a pointer to a filled \code{BAS\_FCTS} structure,
  \code{wall\_quad} a pointer to some quadrature (accessed by
  \code{get\_wall\_quad()}, e.g.)  and \code{init\_flag}
  indicates which information should be filled into the
  \code{QUAD\_FAST} structure. The function returns a pointer to a
  filled \code{QUAD\_FAST} structure where all demanded information
  is computed and stored.
  %% 
  \kitem{get\_neigh\_quad(el\_info, wall\_quad, neigh)} returns a
  suitable quadrature for integrating over the given wall (neigh number),
  but the barycentric co-ordinates of \code{QUAD->lambda} are relative
  to the neighbour element. 

  \kitem{get\_neigh\_quad\_fast(el\_info, wall\_quad, neigh)} returns a
  suitable \code{QUAD\_FAST} structure for integrating over the given
  wall, but relative to the neighbour element. If the returned
  \code{QUAD\_FAST} object has a per-element initializer, then it must be
  called with an \code{EL\_INFO} structure for the neighbour element.
  It is also legal to just call
  \code{get\_quad\_fast(bas\_fcts, get\_neigh\_quad(el\_info, wall\_quad, neigh), ...)} but \code{get\_neigh\_quad\_fast()} is slightly more efficient.
\end{descr}

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


\subsection{Caching of geometric quantities on quadrature
  nodes}
\label{S:fill_quad_el_cache}

Like for geometric quantities which are constant on a given
mesh-element it is useful to share geometric data attached to
quadrature nodes between different places of program code, see also
\secref{S:fill_el_geom_cache}. For this purpose there is a
per-quadrature-per-element cache, called \code{QUAD\_EL\_CACHE}, which
can be filled and accessed through calls to the function
\code{fill\_quad\_el\_cache()}. This is in particular useful for
higher-order parametric meshes, where for example the transformation
to the reference element is no longer piece-wise constant on each
element. Internally, the \code{QUAD\_EL\_CACHE} is maintained as part
of the ``metadata'' attached to each quadrature rule, see
\ref{T:QUAD}. The per-quadrature node cache and the related
definitions and proto-types are as follows:
%%
\bv\begin{lstlisting}[name=QUAD_EL_CACHE,label=T:QUAD_EL_CACHE]
typedef struct quad_el_cache QUAD_EL_CACHE;

struct quad_el_cache
{
  EL      *current_el;
  FLAGS   fill_flag;
  REAL_D  *world;
  struct {
    REAL      *det;
    REAL_BD   *Lambda;
    REAL_BDD  *DLambda;
    REAL_BD   *grd_world;
    REAL_BDB  *D2_world;
    REAL_BDBB *D3_world;
    REAL      *wall_det;    /* for co-dim 1 */
    REAL_D    *wall_normal; /* for co-dim 1 */
    REAL_DB   *grd_normal;  /* for co-dim 1 */
    REAL_DBB  *D2_normal;   /* for co-dim 1 */
  } param;
};

#define FILL_EL_QUAD_WORLD       0x0001
#define FILL_EL_QUAD_DET         0x0002
#define FILL_EL_QUAD_LAMBDA      0x0004
#define FILL_EL_QUAD_DLAMBDA     0x0008
#define FILL_EL_QUAD_GRD_WORLD   0x0010
#define FILL_EL_QUAD_D2_WORLD    0x0020
#define FILL_EL_QUAD_D3_WORLD    0x0040
#define FILL_EL_QUAD_WALL_DET    0x0100
#define FILL_EL_QUAD_WALL_NORMAL 0x0200
#define FILL_EL_QUAD_GRD_NORMAL  0x0400
#define FILL_EL_QUAD_D2_NORMAL   0x0800

static inline const QUAD_EL_CACHE *fill_quad_el_cache(const EL_INFO *el_info,
						      const QUAD *quad,
						      FLAGS fill);
\end{lstlisting}\ev

The quadrature cache can be obtained and filled by calls to
\code{fill\_quad\_el\_cache()}, see also below
\exampleref{E:fill_quad_el_cache}. The members of
\code{QUAD\_EL\_CACHE} have the following meaning:
\begin{descr}
  \kitem{current\_el} For internal use only.
  %%
  \kitem{fill\_flag} A bit-mask, bit-wise or of the fill flags
  listed above (\ref{T:QUAD_EL_CACHE}).
  %% 
  \kitem{world} The world co-ordinates of the quadrature points,
  filled by \code{fill\_quad\_el\_cache(..., FILL\_EL\_QUAD\_WORLD)}.
  %% 
  \kitem{param} A cache for geometric quantities which are constant on
  each element for affine-linear meshes, but vary between quadrature
  points for higher-order parametric meshes.
  \begin{descr}
    \kitem{det} The determinant of the transformation to the reference
    element, filled by 
    filled by \code{fill\_quad\_el\_cache(..., FILL\_EL\_QUAD\_DET)}.
    %%
    \kitem{Lambda} The derivative of the barycentric coordinates
    w.r.t.  the Cartesian coordinates, filled by
    \code{fill\_quad\_el\_cache(..., FILL\_EL\_QUAD\_LAMBDA)}.
    %%
    \kitem{DLambda} The second derivatives of the barycentric
    coordinates w.r.t. the Cartesian coordinates, filled by
    \code{fill\_quad\_el\_cache(..., FILL\_EL\_QUAD\_DLAMBDA)}.
    %%
    \kitem{grd\_world} The first derivatives of the Cartesian
    coordinates w.r.t. the barycentric coordinates, filled by
    \code{fill\_quad\_el\_cache(..., FILL\_EL\_QUAD\_GRD\_WORLD)}.
    %%
    \kitem{D2\_world} The second derivatives of the Cartesian
    coordinates w.r.t. the barycentric coordinates, filled by
    \code{fill\_quad\_el\_cache(..., FILL\_EL\_QUAD\_D2\_WORLD)}.
    %%
    \kitem{D3\_world} The third derivatives of the Cartesian
    coordinates w.r.t. the barycentric coordinates, filled by
    \code{fill\_quad\_el\_cache(..., FILL\_EL\_QUAD\_D3\_WORLD)}.
    %%
    \kitem{wall\_det} The determinant of the transformation of the
    walls to the reference element's walls. This can be filled only
    for co-dimension $1$ quadratures by
    \code{fill\_quad\_el\_cache(..., FILL\_EL\_QUAD\_WALL\_DET)}.
    %%
    \kitem{wall\_normal} The outer wall-normal. This can be filled
    only for co-dimension $1$ quadratures by
    \code{fill\_quad\_el\_cache(..., FILL\_EL\_QUAD\_WALL\_NORMAL)}.
    %%
    \kitem{grd\_normal} The first derivative of the outer normal-field
    with respect to the barycentric coordinates.  This can be filled
    only for co-dimension $1$ quadratures by
    \code{fill\_quad\_el\_cache(..., FILL\_EL\_QUAD\_GRD\_NORMAL)}.
    %%
    \kitem{D2\_normal} The second derivative of the outer normal-field
    with respect to the barycentric coordinates.  This can be filled
    only for co-dimension $1$ quadratures by
    \code{fill\_quad\_el\_cache(..., FILL\_EL\_QUAD\_D2\_NORMAL)}.
  \end{descr}  
\end{descr}

\begin{example}
  \label{E:fill_quad_el_cache}
  A simple example which computes the measure of the region occupied
  by the mesh (of course, this can be achieved more efficiently by
  computing a boundary integral \dots). This example is, of course,
  quite artificial -- and in this context it would be more efficient
  \emph{not} to read through the per-element caches.
  %%
  \bv\begin{lstlisting}[name=fill_quad_el_cache,label=C:quad_el_cache]
    const PARAMERIC *param = mesh->parametric;
    const QUAD *quad = get_quadrature(mesh->dim, 3 /* degree */);
    REAL meas = 0.0;
    TRAVERSE_FIRST(mesh, -1, CALL_LEAF_EL|FILL_COORDS) {
      int iq;
      if (param->init_element(el_info, param)) {
        const QUAD_EL_CACHE *qelc = fill_quad_el_cache(el_info, quad, FILL_EL_QUAD_DET);
        for (iq = 0; iq < quad->n_points; iq++) {
          meas += quad->w[iq] * qelc->param.det[iq];
        }
      } else {
        const EL_GEOM_CACHE *elgc = fill_el_geom_cache(el_info, FILL_EL_DET);
        meas += elgc->det / (REAL)DIM_FAC(mesh->dim);
      }
    } TRAVERSE_NEXT()
  \end{lstlisting}\ev
\end{example}


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

\section{Functions for the evaluation of finite elements}
\label{S:eval}

Finite element functions are evaluated locally on single elements
using barycentric coordinates (compare \secref{book:S:eval_fe}). \ALBERTA
supplies several functions for calculating values and first and second
derivatives of finite element functions on single elements. Functions
for the calculation of derivatives are currently only implemented
for (non--parametric) simplices.

Recalling \mathref{book:E:uh_on_S} on page \pageref{book:E:uh_on_S} we obtain
for the value of a finite element function $\uh$ on an element $S$
\[
\uh(x(\lambda)) = \sum\limits_{i = 1}^m u^i_S\, \pbar^i(\lambda)
\qquad \mbox{for all } \lambda \in \Sbar,
\]
where $\left(\pbar^1,\dots,\pbar^m\right)$ is a basis of $\Pbar$ and
$\left(u_S^1,\dots,u_S^m\right)$ the local coefficient vector
of $\uh$ on $S$. Derivatives are evaluated on $S$ by
\[
\nabla u_h(x(\lambda)) = \Lambda^t
 \sum_{i = 1}^m u_S^i \, \nablal \pbar^i(\lambda), \qquad \lambda \in \Sbar
\]
and
\[
D^2 u_h(x(\lambda)) = \Lambda^t \sum_{i = 1}^m u_S^i \,
D^2_\lambda \pbar^i(\lambda)\Lambda, \qquad \lambda \in \Sbar,
\]
where $\Lambda$ is the Jacobian of the barycentric coordinates,
compare \secref{book:S:eval_Dfe}.

These formulas are used for all evaluation routines. Information
about values of basis functions and their derivatives can be 
calculated via function pointers in the \code{BAS\_FCTS} structure.
Additionally, the local coefficient vector and the Jacobian
of the barycentric coordinates are needed (for the calculation of
derivatives).

The following routines calculate values of a finite element function
at a single point, given in barycentric coordinates:
\fdx{eval_uh()@{\code{eval\_uh()}}}%
\idx{evaluation of finite element functions!eval_uh()@{\code{eval\_uh()}}}
\fdx{eval_grd_uh()@{\code{eval\_grd\_uh()}}}%
\idx{evaluation of finite element functions!eval_grd_uh()@{\code{eval\_grd\_uh()}}}%
\fdx{eval_D2_uh()@{\code{eval\_D2\_uh()}}}%
\idx{evaluation of finite element functions!eval_D2_uh()@{\code{eval\_D2\_uh()}}}%
\fdx{eval_uh_d()@{\code{eval\_uh\_d()}}}%
\idx{evaluation of finite element functions!eval_uh_d()@{\code{eval\_uh\_d()}}}
\fdx{eval_grd_uh_d()@{\code{eval\_grd\_uh\_d()}}}%
\idx{evaluation of finite element functions!eval_grd_uh_d()@{\code{eval\_grd\_uh\_d()}}}%
\fdx{eval_div_uh_d()@{\code{eval\_div\_uh\_d()}}}%
\idx{evaluation of finite element functions!eval_div_uh_d()@{\code{eval\_div\_uh\_d()}}}%
\fdx{eval_D2_uh_d()@{\code{eval\_D2\_uh\_d()}}}%
\idx{evaluation of finite element functions!eval_D2_uh_d()@{\code{eval\_D2\_uh\_d()}}}%
\fdx{eval_uh_dow()@{\code{eval\_uh\_dow()}}}%
\idx{evaluation of finite element functions!eval_uh_dow()@{\code{eval\_uh\_dow()}}}
\fdx{eval_grd_uh_dow()@{\code{eval\_grd\_uh\_dow()}}}%
\idx{evaluation of finite element functions!eval_grd_uh_dow()@{\code{eval\_grd\_uh\_dow()}}}%
\fdx{eval_div_uh_dow()@{\code{eval\_div\_uh\_dow()}}}%
\idx{evaluation of finite element functions!eval_div_uh_dow()@{\code{eval\_div\_uh\_dow()}}}%
\fdx{eval_D2_uh_dow()@{\code{eval\_D2\_uh\_dow()}}}%
\idx{evaluation of finite element functions!eval_D2_uh_dow()@{\code{eval\_D2\_uh\_dow()}}}%
\bv\begin{lstlisting}
REAL eval_uh(const REAL_B lambda, const EL_REAL_VEC *uh_loc,
             const BAS_FCTS *bfcts);
REAL *eval_grd_uh(REAL_D result, const REAL_B lambda, const REAL_BD Lambda,
                  const EL_REAL_VEC *uh_loc, const BAS_FCTS *bfcts);
REAL_D *eval_D2_uh(REAL_DD result, const REAL_B lambda, const REAL_BD Lambda,
                   const EL_REAL_VEC *uh_loc, const BAS_FCTS *bfcts);

REAL *eval_uh_d(REAL_D result, const REAL_B lambda,
                const EL_REAL_D_VEC *uh_loc, const BAS_FCTS *bfcts);
REAL_D *eval_grd_uh_d(REAL_DD result, const REAL_B lambda,
                      const REAL_BD Lambda, const EL_REAL_D_VEC *uh_loc,
                      const BAS_FCTS *bfcts);
REAL eval_div_uh_d(const REAL_B lambda, const REAL_BD Lambda,
                   const EL_REAL_D_VEC *uh_loc, const BAS_FCTS *bfcts);
REAL_DD *eval_D2_uh_d(REAL_DDD result, const REAL_B lambda,
                      const REAL_BD Lambda, const EL_REAL_D_VEC *uh_loc,
                      const BAS_FCTS *bfcts);

REAL *eval_uh_dow(REAL_D result, const REAL_B lambda,
                  const EL_REAL_VEC_D *uh_loc, const BAS_FCTS *bfcts);
REAL_D *eval_grd_uh_dow(REAL_DD result, const REAL_B lambda,
                        const REAL_BD Lambda, const EL_REAL_VEC_D *uh_loc,
                        const BAS_FCTS *bfcts);
REAL eval_div_uh_dow(const REAL_B lambda, const REAL_BD Lambda,
                     const EL_REAL_VEC_D *uh_loc, const BAS_FCTS *bfcts);
REAL_DD *eval_D2_uh_dow(REAL_DDD result, const REAL_B lambda,
                        const REAL_BD Lambda, const EL_REAL_VEC_D *uh_loc,
                        const BAS_FCTS *bfcts);
\end{lstlisting}\ev
Description:

In the following $\code{lambda} = \lambda$ are the barycentric
coordinates at which the function is evaluated, $\code{Lambda} =
\Lambda$ is the Jacobian of the barycentric coordinates, \code{uh} the
local coefficient vector $\left(u_S^0,\dots,u_S^{\code{m}-1}\right)$
(where $u_S^i$ is a \code{REAL} or a \code{REAL\_D}), and
\code{bas\_fcts} is a pointer to a \code{BAS\_FCTS} structure, storing
information about the set of local basis functions
$\left(\pbar^0,\dots,\pbar^{\code{m}-1}\right)$.

All functions returning a pointer to a vector or matrix provide memory
for the vector or matrix in a statically allocated memory area. This
area is overwritten during the next call. If the first argument of such
a function is not \nil, then it is a pointer to a storage area where
the results are stored. This memory area must be of correct size, no
check is performed.

\begin{compatibility}
  \label{compat:resultspace}
  Former versions of \ALBERTA expected the argument providing optional
  storage for the result at the last place in the parameter list. In
  the current version of the library, storage for the result is still
  optional, but generally passed as first argument to the respective
  function.
\end{compatibility}

The functions for \DOW-valued discrete functions come in two variants,
one for discrete functions based on scalar-valued local basis function
sets, where the coefficients are \DOW-valued, and one for discrete
functions which may be based on either scalar-valued or \DOW-valued
local basis functions, modeled by \code{DOF\_REAL\_VEC\_D} -- and
locally by \code{EL\_REAL\_VEC\_D} -- objects. The names for the
latter functions have a \code{\dots\_dow} suffix, the others a
\code{\dots\_d} suffix. Besides the slightly differing argument types
the calling conventions for both variants are the same, so they are
documented together in the descriptions following below.

\begin{descr}
\kitem{eval\_uh(lambda, uh\_loc, bas\_fcts)} the function returns
  $\uh(\lambda)$.
  %% 
\kitem{eval\_grd\_uh(result, lambda, Lambda, uh\_loc, bas\_fcts)} the function 
  returns a pointer \code{ptr} to a vector of length \DOW storing
  $\nabla \uh(\lambda)$, i.e.
  \[
  \mbox{\code{ptr[i]}} = {\uh}_{,x_\code{i}}(\lambda),
  \qquad \code{i}=0,\dots,\mbox{\DOW}-1;
  \]
  \code{result} is optional and provides storage for the result if non-\nil.
  See \compatref{compat:resultspace}.
  %% 
\kitem{eval\_D2\_uh(result, lambda, Lambda, uh\_loc, bas\_fcts)} the function 
  returns a pointer \code{ptr} to a matrix of size
  $(\DOW\times\DOW)$ storing $D^2 \uh(\lambda)$, i.e.
  \[
  \mbox{\code{ptr[i][j]}} = {\uh}_{,x_\code{i}x_\code{j}}(\lambda),
  \qquad \code{i},\code{j}=0,\dots,\mbox{\DOW}-1;
  \]
  \code{result} is optional and provides storage for the result if non-\nil.
  See \compatref{compat:resultspace}.
  %% 
\kitem{eval\_uh\_[d|dow](result, lambda, uh\_loc, bas\_fcts)}
  the function 
  returns a pointer \code{ptr} to a vector of length \DOW storing
  $\uh(\lambda)$, i.e.
  \[
  \mbox{\code{ptr[k]}} = {\uh}_\code{k}(\lambda),
  \qquad \code{k}=0,\dots,\mbox{\DOW}-1;
  \]
  \code{result} is optional and provides storage for the result if non-\nil.
  See \compatref{compat:resultspace}.
 %% 
\kitem{eval\_grd\_uh\_[d|dow](result, lambda, Lambda, uh\_loc, bas\_fcts)}
  the function returns a pointer \code{ptr} to a vector of \DOW
  vectors of length \DOW storing $\nabla {\uh}(\lambda)$, i.e.
  \[
  \mbox{\code{ptr[k][i]}} = {\uh}_{\code{k},x_\code{i}}(\lambda),
  \qquad \code{k},\code{i}=0,\dots,\mbox{\DOW}-1;
  \]
  \code{result} is optional and provides storage for the result if non-\nil.
  See \compatref{compat:resultspace}.
 %% 
\kitem{eval\_div\_uh\_[d|dow](lambda, Lambda, uh\_loc, bas\_fcts)} the function 
  returns \, $\div \uh(\lambda)$.
  %% 
\kitem{eval\_D2\_uh\_[d|dow](result, lambda, Lambda, uh\_loc, bas\_fcts)} the function 
  returns a pointer \code{ptr} to a vector of $(\DOW\times\DOW)$
  matrices of length \DOW storing $D^2 {\uh}(\lambda)$, i.e.
  \[
  \mbox{\code{ptr[k][i][j]}} = {\uh}_{\code{k},x_\code{i}x_\code{j}}(\lambda),
  \qquad \code{k},\code{i},\code{j}=0,\dots,\mbox{\DOW}-1;
  \]
  \code{result} is optional and provides storage for the result if non-\nil.
  See \compatref{compat:resultspace}.
 %% 
\end{descr}

Using pre--computed values of basis functions at the evaluation point,
these routines can be implemented more efficiently.
\fdx{eval_uh_fast()@{\code{eval\_uh\_fast()}}}
\idx{evaluation of finite element functions!eval_uh_fast()@{\code{eval\_uh\_fast()}}}
\fdx{eval_grd_uh_fast()@{\code{eval\_grd\_uh\_fast()}}}
\idx{evaluation of finite element functions!eval_grd_uh_fast()@{\code{eval\_grd\_uh\_fast()}}}
\fdx{eval_D2_uh_fast()@{\code{eval\_D2\_uh\_fast()}}}
\idx{evaluation of finite element functions!eval_D2_uh_fast()@{\code{eval\_D2\_uh\_fast()}}}
\fdx{eval_uh_d_fast()@{\code{eval\_uh\_d\_fast()}}}
\idx{evaluation of finite element functions!eval_uh_d_fast()@{\code{eval\_uh\_d\_fast()}}}
\fdx{eval_grd_uh_d_fast()@{\code{eval\_grd\_uh\_d\_fast()}}}
\idx{evaluation of finite element functions!eval_grd_uh_d_fast()@{\code{eval\_grd\_uh\_d\_fast()}}}
\fdx{eval_div_uh_d_fast()@{\code{eval\_div\_uh\_d\_fast()}}}
\idx{evaluation of finite element functions!eval_div_uh_d_fast()@{\code{eval\_div\_uh\_d\_fast()}}}
\fdx{eval_D2_uh_d_fast()@{\code{eval\_D2\_uh\_d\_fast()}}}
\idx{evaluation of finite element functions!eval_D2_uh_d_fast()@{\code{eval\_D2\_uh\_d\_fast()}}}
\fdx{eval_uh_dow_fast()@{\code{eval\_uh\_dow\_fast()}}}
\idx{evaluation of finite element functions!eval_uh_dow_fast()@{\code{eval\_uh\_dow\_fast()}}}
\fdx{eval_grd_uh_dow_fast()@{\code{eval\_grd\_uh\_dow\_fast()}}}
\idx{evaluation of finite element functions!eval_grd_uh_dow_fast()@{\code{eval\_grd\_uh\_dow\_fast()}}}
\fdx{eval_div_uh_dow_fast()@{\code{eval\_div\_uh\_dow\_fast()}}}
\idx{evaluation of finite element functions!eval_div_uh_dow_fast()@{\code{eval\_div\_uh\_dow\_fast()}}}
\fdx{eval_D2_uh_dow_fast()@{\code{eval\_D2\_uh\_dow\_fast()}}}
\idx{evaluation of finite element functions!eval_D2_uh_dow_fast()@{\code{eval\_D2\_uh\_dow\_fast()}}}
\bv\begin{lstlisting}
REAL eval_uh_fast(const EL_REAL_VEC *uh_loc, const QUAD_FAST *qfast, int iq);
const REAL *eval_grd_uh_fast(REAL_D grd_uh, const REAL_BD Lambda,
                             const EL_REAL_VEC *uh_loc,
                             const QUAD_FAST *qfast, int iq);
const REAL_D *eval_D2_uh_fast(REAL_DD result, const REAL_BD Lambda,
                              const EL_REAL_VEC *uh_loc,
                              const QUAD_FAST *qfast, int iq);

const REAL *eval_uh_d_fast(REAL_D result, const EL_REAL_D_VEC *uh_loc,
                           const QUAD_FAST *qfast, int iq);
const REAL_D *eval_grd_uh_d_fast(REAL_DD result, const REAL_BD Lambda,
                                 const EL_REAL_D_VEC *uh_loc,
                                 const QUAD_FAST *qfast, int iq);
REAL eval_div_uh_d_fast(const REAL_BD Lambda, const EL_REAL_D_VEC *uh_loc,
                        const QUAD_FAST *qfast, int iq);
const REAL_DD *eval_D2_uh_d_fast(REAL_DDD result, const REAL_BD Lambda,
                                 const EL_REAL_D_VEC *uh_loc,
                                 const QUAD_FAST *qfast, int iq);

const REAL *eval_uh_dow_fast(REAL_D result, const EL_REAL_VEC_D *uh_loc,
                             const QUAD_FAST *qfast, int iq);
const REAL_D *eval_grd_uh_dow_fast(REAL_DD result, const REAL_BD Lambda,
                                   const EL_REAL_VEC_D *uh_loc,
                                   const QUAD_FAST *qfast, int iq);
REAL eval_div_uh_dow_fast(const REAL_BD Lambda, const EL_REAL_VEC_D *uh_loc,
                          const QUAD_FAST *qfast, int iq);
const REAL_DD *eval_D2_uh_dow_fast(REAL_DDD result, const REAL_BD Lambda,
                                   const EL_REAL_VEC_D *uh_loc,
                                   const QUAD_FAST *qfast, int iq);
\end{lstlisting}\ev
\begin{compatibility}
  \label{compat:evalcallconv}
  Former versions of \ALBERTA didn't expect the arguments
  \begin{lstlisting}
    ..., const QUAD_FAST *qfast, int iq, ...    
  \end{lstlisting}
  -- meaning the \hyperref[S:QUAD_FAST]{quadrature cache} and the
  index of the quadrature point -- but instead expected the actual
  cached-values to be passed, i.e. for the computation of the gradient
  \begin{lstlisting}
    ..., qfast->grd_phi[iq], qfast->n_bas_fcts, ...
  \end{lstlisting}
  There is some potential for confusion, in particular because the
  proto-types listed in the old documentation often omit the parameter
  name and only give the parameter type. In the new version,
  \lstinline!.., int iq, ...! denotes the index of the quadrature
  point. The number of basis functions on the reference element is not
  needed, because the evaluation functions fetch this quantity
  themselves from the \hyperref[S:QUAD_FAST]{\code{QUAD\_FAST}} data
  structure.
\end{compatibility}
Description: 
In the following $\code{Lambda} = \Lambda$ denotes the Jacobian of the
barycentric coordinates, \code{uh\_loc} the local coefficient vector (of
type \hyperref[T:EL_REAL_VEC]{\code{EL\_REAL\_VEC}},
\hyperref[T:EL_REAL_D_VEC]{\code{EL\_REAL\_D\_VEC}} etc.) on an element.
\begin{descr}
\kitem{eval\_uh\_fast(uh\_loc, qfast, iq)} the function returns $\uh(\lambda)$;

 \code{qfast} is a quadrature cache storing the values 
  $\pbar^0(\lambda),\dots,\pbar^{\code{m}-1}(\lambda)$.
  %% 
\kitem{eval\_grd\_uh\_fast(grd, Lambda, uh\_loc, qfast, iq)} 
  the function returns a pointer \code{ptr} to a vector of length
  \DOW storing $\nabla \uh(\lambda)$, i.e.
  \[
  \mbox{\code{ptr[i]}} = {\uh}_{,x_\code{i}}(\lambda),
  \qquad \code{i}=0,\dots,\mbox{\DOW}-1;
  \]
  \code{grd} is optional and provides storage for the result if non-\nil.
  See \compatref{compat:resultspace}.

  \code{qfast} is a quadrature cache storing 
  $\nablal \pbar^0(\lambda),\dots,\nablal \pbar^{\code{m}-1}(\lambda)$;
 %% 
\kitem{eval\_D2\_uh\_fast(D2, Lambda, uh\_loc, qfast, iq)} the function 
  returns a pointer \code{ptr} to a matrix of size
  $(\DOW\times\DOW)$ storing $D^2 \uh(\lambda)$, i.e.
  \[
  \mbox{\code{ptr[i][j]}} = {\uh}_{,x_\code{i}x_\code{j}}(\lambda),
  \qquad \code{i},\code{j}=0,\dots,\mbox{\DOW}-1;
  \]
  \code{D2} is optional and provides storage for the result if non-\nil.
  See \compatref{compat:resultspace}.

  \code{qfast} is a quadrature cache storing 
  $D^2_\lambda \pbar^0(\lambda),\dots,D2_\lambda\pbar^{\code{m}-1}
  (\lambda)$.
 %% 
\kitem{eval\_uh\_[d|dow]\_fast(result, uh\_loc, qfast, iq)} the function 
  returns a pointer \code{ptr} to a vector of \DOW vectors of
  length \DOW storing $\nabla {\uh}(\lambda)$, i.e. 
  \[
  \mbox{\code{ptr[k][i]}} = {\uh}_{\code{k},x_\code{i}}(\lambda),
  \qquad \code{k},\code{i}=0,\dots,\mbox{\DOW}-1;
  \]
  \code{qfast} is a quadrature cache storing the values 
  $\pbar^0(\lambda),\dots,\pbar^{\code{m}-1}(\lambda)$;

  \code{result} is optional and provides storage for the resulty if non-\nil.
  See \compatref{compat:resultspace}.
  %% 
\kitem{eval\_grd\_uh\_[d|dow]\_fast(grd, Lambda, uh\_loc, qfast, iq)} 
  the function returns a pointer \code{ptr} to a vector of \DOW
  vectors of length \DOW storing $\nabla {\uh}(\lambda)$, i.e.
  \[
  \mbox{\code{ptr[k][i]}} = {\uh}_{\code{k},x_\code{i}}(\lambda),
  \qquad \code{k},\code{i}=0,\dots,\mbox{\DOW}-1;
  \]
  \code{qfast} is a quadrature cache storing
  $\nablal \pbar^0(\lambda),\dots,\nablal
  \pbar^{\code{m}-1}(\lambda)$;
  \code{grd} is optional storage for the result if non-\nil.
  See \compatref{compat:resultspace}.
  %% 
\kitem{eval\_div\_uh\_[d|dow]\_fast(Lambda, uh\_loc, qfast, iq)} the function 
  returns $\div \uh(\lambda)$;

  \code{qfast} is a quadrature cache storing 
  $\nablal \pbar^0(\lambda),\dots,\nablal \pbar^{\code{m}-1}(\lambda)$.
  Unused entries must be set to \code{0.0}.
  %% 
\kitem{eval\_D2\_uh\_[d|dow]\_fast(D2, Lambda, uh\_loc, qfast, iq)}
  the function returns a pointer \code{ptr} to a vector of
  $(\DOW\times\DOW)$ matrices of length \DOW storing
  $D^2 {\uh}(\lambda)$, i.e.
  \[
  \mbox{\code{ptr[k][i][j]}} = {\uh}_{\code{k},x_\code{i}x_\code{j}}(\lambda),
  \qquad \code{k},\code{i},\code{j}=0,\dots,\mbox{\DOW}-1;
  \]
  \code{qfast} is a quadrature cache storing 
  $D^2_\lambda \pbar^0(\lambda),\dots,D2_\lambda\pbar^{\code{m}-1}
  (\lambda)$;

  \code{D2} is optional storage for the result if non-\nil.
  See \compatref{compat:resultspace}.
\end{descr}

One important task is the evaluation of finite element functions at
all quadrature nodes for a given quadrature formula. Using the
\code{QUAD\_FAST} data structures, the values of the basis functions
are known at the quadrature nodes which results in an efficient
calculation of values and derivatives of finite element functions
at these quadrature points.
\fdx{uh_at_qp()@{\code{uh\_at\_qp()}}}%
\idx{evaluation of finite element functions!uh_at_qp()@{\code{uh\_at\_qp()}}}%
\idx{numerical quadrature!uh_at_qp()@{\code{uh\_at\_qp()}}}%
\fdx{grd_uh_at_qp()@{\code{[param\_]grd\_uh\_at\_qp()}}}%
\idx{evaluation of finite element functions!grd_uh_at_qp()@{\code{[param\_]grd\_uh\_at\_qp()}}}%
\idx{numerical quadrature!grd_uh_at_qp()@{\code{[param\_]grd\_uh\_at\_qp()}}}%
\fdx{D2_uh_at_qp()@{\code{[param\_]D2\_uh\_at\_qp()}}}
\idx{evaluation of finite element functions!D2_uh_at_qp()@{\code{[param\_]D2\_uh\_at\_qp()}}}
\idx{numerical quadrature!D2_uh_at_qp()@{\code{[param\_]D2\_uh\_at\_qp()}}}
\fdx{uh_d_at_qp()@{\code{uh\_d\_at\_qp()}}}%
\idx{evaluation of finite element functions!uh_d_at_qp()@{\code{uh\_d\_at\_qp()}}}%
\idx{numerical quadrature!uh_d_at_qp()@{\code{uh\_d\_at\_qp()}}}%
\fdx{grd_uh_d_at_qp()@{\code{[param\_]grd\_uh\_d\_at\_qp()}}}%
\idx{evaluation of finite element functions!grd_uh_d_at_qp()@{\code{[param\_]grd\_uh\_d\_at\_qp()}}}%
\idx{numerical quadrature!grd_uh_d_at_qp()@{\code{[param\_]grd\_uh\_d\_at\_qp()}}}%
\fdx{div_uh_d_at_qp()@{\code{[param\_]div\_uh\_d\_at\_qp()}}}%
\idx{evaluation of finite element functions!div_uh_d_at_qp()@{\code{[param\_]div\_uh\_d\_at\_qp()}}}%
\idx{numerical quadrature!div_uh_d_at_qp()@{\code{[param\_]div\_uh\_d\_at\_qp()}}}%
\fdx{D2_uh_d_at_qp()@{\code{[param\_]D2\_uh\_d\_at\_qp()}}}%
\idx{evaluation of finite element functions!D2_uh_d_at_qp()@{\code{[param\_]D2\_uh\_d\_at\_qp()}}}%
\idx{numerical quadrature!D2_uh_d_at_qp()@{\code{[param\_]D2\_uh\_d\_at\_qp()}}}%
\fdx{uh_dow_at_qp()@{\code{uh\_dow\_at\_qp()}}}%
\idx{evaluation of finite element functions!uh_dow_at_qp()@{\code{uh\_dow\_at\_qp()}}}%
\idx{numerical quadrature!uh_dow_at_qp()@{\code{uh\_dow\_at\_qp()}}}%
\fdx{grd_uh_dow_at_qp()@{\code{[param\_]grd\_uh\_dow\_at\_qp()}}}%
\idx{evaluation of finite element functions!grd_uh_dow_at_qp()@{\code{[param\_]grd\_uh\_dow\_at\_qp()}}}%
\idx{numerical quadrature!grd_uh_dow_at_qp()@{\code{[param\_]grd\_uh\_dow\_at\_qp()}}}%
\fdx{div_uh_dow_at_qp()@{\code{[param\_]div\_uh\_dow\_at\_qp()}}}%
\idx{evaluation of finite element functions!div_uh_dow_at_qp()@{\code{[param\_]div\_uh\_dow\_at\_qp()}}}%
\idx{numerical quadrature!div_uh_dow_at_qp()@{\code{[param\_]div\_uh\_dow\_at\_qp()}}}%
\fdx{D2_uh_dow_at_qp()@{\code{[param\_]D2\_uh\_dow\_at\_qp()}}}%
\idx{evaluation of finite element functions!D2_uh_dow_at_qp()@{\code{[param\_]D2\_uh\_dow\_at\_qp()}}}%
\idx{numerical quadrature!D2_uh_dow_at_qp()@{\code{[param\_]D2\_uh\_dow\_at\_qp()}}}%
\bv\begin{lstlisting}
REAL *uh_at_qp(REAL *result, const QUAD_FAST *qfast,
               const EL_REAL_VEC *uh_loc);
REAL_D *grd_uh_at_qp(REAL_D *result, const QUAD_FAST *qfast,
                     const REAL_BD Lambda, const EL_REAL_VEC *uh_loc);
REAL_DD *D2_uh_at_qp(REAL_DD *result, const QUAD_FAST *qfast,
                     const REAL_BD Lambda, const EL_REAL_VEC *uh_loc);

REAL_D *param_grd_uh_at_qp(REAL_D vec[], const QUAD_FAST *qfast,
                           const REAL_BD Lambda[], const EL_REAL_VEC *uh_loc);
REAL_DD *param_D2_uh_at_qp(REAL_DD *result,  const QUAD_FAST *qfast,
                           const REAL_BD Lambda[], const REAL_BDD DLambda[],
                           const EL_REAL_VEC *uh_loc);

REAL_D *uh_d_at_qp(REAL_D *result, const QUAD_FAST *qfast,
                   const EL_REAL_D_VEC *uh_loc);
REAL_DD *grd_uh_d_at_qp(REAL_DD *result, const QUAD_FAST *qfast,
                        const REAL_BD Lambda, const EL_REAL_D_VEC *uh_loc);
REAL *div_uh_d_at_qp(REAL *result, const QUAD_FAST *qfast,
                     const REAL_BD Lambda, const EL_REAL_D_VEC *uh_loc);
REAL_DDD *D2_uh_d_at_qp(REAL_DDD vec[], const QUAD_FAST *qfast,
                        const REAL_BD Lambda, const EL_REAL_D_VEC *uh_loc);

REAL_DD *param_grd_uh_d_at_qp(REAL_DD vec[], const QUAD_FAST *qfast,
                              const REAL_BD Lambda[],
                              const EL_REAL_D_VEC *uh_loc);
REAL *param_div_uh_d_at_qp(REAL vec[], const QUAD_FAST *qfast,
                           const REAL_BD Lambda[],
                           const EL_REAL_D_VEC *uh_loc);
REAL_DDD *param_D2_uh_d_at_qp(REAL_DDD vec[], const QUAD_FAST *qfast,
                              const REAL_BD grd_lam[],
                              const REAL_BDD DLambda[],
                              const EL_REAL_D_VEC *uh_loc);

REAL_D *uh_dow_at_qp(REAL_D *result, const QUAD_FAST *qfast,
                     const EL_REAL_VEC_D *uh_loc);
REAL_DD *grd_uh_dow_at_qp(REAL_DD *result, const QUAD_FAST *qfast,
                          const REAL_BD Lambda, const EL_REAL_VEC_D *uh_loc);
REAL *div_uh_dow_at_qp(REAL *result, const QUAD_FAST *qfast,
                       const REAL_BD Lambda,  const EL_REAL_VEC_D *uh_loc);
REAL_DDD *D2_uh_dow_at_qp(REAL_DDD vec[], const QUAD_FAST *qfast,
                          const REAL_BD Lambda, const EL_REAL_VEC_D *uh_loc);

REAL_DD *param_grd_uh_dow_at_qp(REAL_DD vec[], const QUAD_FAST *qfast,
                                const REAL_BD *Lambda,
                                const EL_REAL_VEC_D *uh_loc);
REAL *param_div_uh_dow_at_qp(REAL vec[], const QUAD_FAST *qfast,
                             const REAL_BD *Lambda,
                             const EL_REAL_VEC_D *uh_loc);
REAL_DDD *param_D2_uh_dow_at_qp(REAL_DDD *result, const QUAD_FAST *qfast,
                                const REAL_BD *Lambda, const REAL_BDD *DLambda,
                                const EL_REAL_VEC_D *uh_loc);
\end{lstlisting}\ev
Description:
In the following \code{uh\_loc} denotes the local coefficient vector (of
type \hyperref[T:EL_REAL_VEC]{\code{EL\_REAL\_VEC}},
\hyperref[T:EL_REAL_D_VEC]{\code{EL\_REAL\_D\_VEC}} etc.) on an element.
\begin{descr}
\kitem{uh\_at\_qp(result, qfast, uh\_loc)} the function returns a pointer
  \code{ptr} to a vector of length \code{qfast->n\_points} storing the
  values of $\uh$ at all quadrature points of \code{qfast->quad}, i.e.
  \[
  \mbox{\code{ptr[l]}} = \uh(\mbox{\code{qfast->quad->lambda[l]}})
  \]
  where $\code{l}=0,\dots,\mbox{\code{qfast->quad->n\_points}}-1$;

  the \code{INIT\_PHI} flag must be set in \code{qfast->init\_flag};

  \code{result} is optional and provides storage for the result if non-\nil.
  See \compatref{compat:resultspace}.
 %% 
\kitem{grd\_uh\_at\_qp(result, qfast, Lambda, uh\_loc)} the function returns
  a pointer \code{ptr} to a vector of length \code{qfast->n\_points}
  of \DOW vectors storing $\nabla \uh$ at all quadrature points of
  \code{qfast->quad}, i.e.
  \[
  \mbox{\code{ptr[l][i]}} = 
  {\uh}_{,x_\code{i}}(\mbox{\code{qfast->quad->lambda[l]}})
  \]
  where $\code{l}=0,\dots,\mbox{\code{qfast->quad->n\_points}}-1$, and
  $\code{i}=0,\dots,\mbox{\DOW}-1$;

  the \code{INIT\_GRD\_PHI} flag must be set in
  \code{qfast->init\_flag};

  \code{result} is optional and provides storage for the result if non-\nil.
  See \compatref{compat:resultspace}.

  %% 
\kitem{D2\_uh\_at\_qp(result, qfast, Lambda, uh\_loc)}

  %% 
\kitem{param\_grd\_uh\_at\_qp(result, qfast, Lambdas, uh\_loc)} version for
  parametric meshes; must be passed a vector storing the gradients of
  the barycentric coordinates at each quadrature point. The same holds
  for the other \code{param\_}-prefixed routines.

  %% 
\kitem{[param\_]D2\_uh\_at\_qp(result, qfast, Lambda[s], uh\_loc, D2)} The
  function returns a pointer \code{ptr} to a vector of length
  \code{qfast->n\_points} of $(\DOW\times\DOW)$ matrices storing $D^2
  \uh$ at all quadrature points of \code{qfast->quad}, i.e.
  \[
  \mbox{\code{ptr[l][i][j]}} =
  {\uh}_{,x_\code{i}x_\code{j}}(\mbox{\code{qfast->quad->lambda[l]}})
  \]
  where $\code{l}=0,\dots,\mbox{\code{qfast->quad->n\_points}}-1$, and
  $\code{i},\code{j}=0,\dots,\mbox{\DOW}-1$;
  
  the \code{INIT\_D2\_PHI} flag must be set in
  \code{qfast->init\_flag};

  \code{result} is optional and provides storage for the result if non-\nil.
  See \compatref{compat:resultspace}.

  %% 
\kitem{uh\_[d|dow]\_at\_qp(result, qfast, uh\_loc)} The function returns a pointer
  \code{ptr} to a vector of length \code{qfast->n\_points} of \DOW
  vectors storing the values of $\uh$ at all quadrature points of
  \code{qfast->quad}, i.e.
  \[
  \mbox{\code{ptr[l][k]}} = 
  {\uh}_\code{k}(\mbox{\code{qfast->quad->lambda[l]}})
  \]
  where $\code{l}=0,\dots,\mbox{\code{qfast->quad->n\_points}}-1$, and
  $\code{k}=0,\dots,\mbox{\DOW}-1$;

  the \code{INIT\_PHI} flag must be set in \code{qfast->init\_flag};

  \code{result} is optional and provides storage for the result if non-\nil.
  See \compatref{compat:resultspace}.

  %% 
\kitem{grd\_uh\_[d|dow]\_at\_qp(result, qfast, Lambda, uh\_loc)}
  %% 
\kitem{div\_uh\_[d|dow]\_at\_qp(result, qfast, Lambda, uh\_loc)}
  %% 
\kitem{D2\_uh\_[d|dow]\_at\_qp(result, qfast, Lambda[], uh\_loc)} The function
  returns a pointer \code{ptr} to a vector of length
  \code{qfast->n\_points} of $(\DOW\times\DOW\times\DOW)$ tensors
  storing $D^2 \uh$ at all quadrature points \code{qfast->quad}, i.e.
  \[
  \mbox{\code{ptr[l][k][i][j]}} =
  {\uh}_{\code{k},x_\code{i}x_\code{j}}(\mbox{\code{qfast->quad->lambda[l]}})
  \]
  where $\code{l}=0,\dots,\mbox{\code{qfast->quad->n\_points}}-1$, and
  $\code{k},\code{i},\code{j}=0,\dots,\mbox{\DOW}-1$;

  the \code{INIT\_D2\_PHI} flag must be set in
  \code{qfast->init\_flag};

  \code{result} is optional and provides storage for the result if non-\nil.
  See \compatref{compat:resultspace}.

  %% 
\kitem{param\_grd\_uh\_[d|dow]\_at\_qp(vec[], qfast, Lambda[], uh\_loc)} The
  function returns a pointer \code{ptr} to a vector of length
  \code{qfast->n\_points} of $(\DOW\times\DOW)$ matrices storing
  $\nabla \uh$ at all quadrature points of \code{qfast->quad}, i.e.
  \[
  \mbox{\code{ptr[l][k][i]}} = 
  {\uh}_{\code{k},x_\code{i}}(\mbox{\code{qfast->quad->lambda[l]}})
  \]
  where $\code{l}=0,\dots,\mbox{\code{qfast->quad->n\_points}}-1$, and
  $\code{k},\code{i}=0,\dots,\mbox{\DOW}-1$;

  the \code{INIT\_GRD\_PHI} flag must be set in
  \code{qfast->init\_flag};

  \code{vec} is optional and provides storage for the result if non-\nil.
  See \compatref{compat:resultspace}.

  %% 
\kitem{param\_div\_uh\_[d|dow]\_at\_qp(result[], qfast, Lambda[], uh\_loc)}

  %% 
\kitem{param\_D2\_uh\_[d|dow]\_at\_qp(result, qfast, Lambda, DLambda, uh\_loc)}
  Second derivatives for parametric meshes.  Note that one needs the
  second derivatives \code{DLambda} of the barycentric co-ordinates
  with respect to the cartesian co-ordiantes for this function. Also
  note that -- in the case of non-zero co-dimension -- the matrix
  $(\ul\nabla(\ul\nabla u)_i)_j$ built from the components of the
  second tangential derivatives is \emph{not} symmetric in general.

  \code{vec} is optional and provides storage for the result if non-\nil.
  See \compatref{compat:resultspace}.
\end{descr}

\fdx{eval_bar_grd_uh()@{\code{eval\_bar\_grd\_uh()}}}
\fdx{eval_bar_grd_uh_d()@{\code{eval\_bar\_grd\_uh\_d()}}}
\fdx{eval_bar_grd_uh_dow()@{\code{eval\_bar\_grd\_uh\_dow()}}}
\fdx{eval_bar_grd_uh_fast()@{\code{eval\_bar\_grd\_uh\_fast()}}}
\fdx{eval_bar_grd_uh_d_fast()@{\code{eval\_bar\_grd\_uh\_d\_fast()}}}
\fdx{eval_bar_grd_uh_dow_fast()@{\code{eval\_bar\_grd\_uh\_dow\_fast()}}}
\fdx{bar_grd_uh_at_qp()@{\code{bar\_grd\_uh\_at\_qp()}}}
\fdx{bar_grd_uh_d_at_qp()@{\code{bar\_grd\_uh\_d\_at\_qp()}}}
\fdx{bar_grd_uh_dow_at_qp()@{\code{bar\_grd\_uh\_dow\_at\_qp()}}}
\fdx{eval_bar_D2_uh()@{\code{eval\_bar\_D2\_uh()}}}
\fdx{eval_bar_D2_uh_d()@{\code{eval\_bar\_D2\_uh\_d()}}}
\fdx{eval_bar_D2_uh_dow()@{\code{eval\_bar\_D2\_uh\_dow()}}}
\fdx{eval_bar_D2_uh_fast()@{\code{eval\_bar\_D2\_uh\_fast()}}}
\fdx{eval_bar_D2_uh_d_fast()@{\code{eval\_bar\_D2\_uh\_d\_fast()}}}
\fdx{eval_bar_D2_uh_dow_fast()@{\code{eval\_bar\_D2\_uh\_dow\_fast()}}}
\fdx{bar_D2_uh_at_qp()@{\code{bar\_D2\_uh\_at\_qp()}}}
\fdx{bar_D2_uh_d_at_qp()@{\code{bar\_D2\_uh\_d\_at\_qp()}}}
\fdx{bar_D2_uh_dow_at_qp()@{\code{bar\_D2\_uh\_dow\_at\_qp()}}}

\bv\begin{lstlisting}
REAL *eval_bar_grd_uh(REAL_B result, const REAL_B lambda,
                      const EL_REAL_VEC *uh_loc, const BAS_FCTS *bfcts);
REAL_B *eval_bar_grd_uh_d(REAL_DB result, const REAL_B lambda,
                          const EL_REAL_D_VEC *uh_loc, const BAS_FCTS *bfcts);
REAL_B *eval_bar_grd_uh_dow(REAL_DB result, const REAL_B lambda,
                            const EL_REAL_VEC_D *uh_loc,
                            const BAS_FCTS *bfcts);

REAL *eval_bar_grd_uh_fast(REAL_B result, const EL_REAL_VEC *uh_loc,
                           const QUAD_FAST *qfast, int iq);
REAL_B *eval_bar_grd_uh_d_fast(REAL_DB result, const EL_REAL_D_VEC *uh_loc,
                               const QUAD_FAST *qfast, int iq);
REAL_B *eval_bar_grd_uh_dow_fast(REAL_DB result, const EL_REAL_VEC_D *uh_loc,
                                 const QUAD_FAST *qfast, int iq);

REAL_B *bar_grd_uh_at_qp(REAL_B *result, const QUAD_FAST *qfast,
                         const EL_REAL_VEC *uh_loc);
REAL_DB *bar_grd_uh_d_at_qp(REAL_DB *result, const QUAD_FAST *qfast,
                            const EL_REAL_D_VEC *uh_loc);
REAL_DB *bar_grd_uh_dow_at_qp(REAL_DB *result, const QUAD_FAST *fast,
                              const EL_REAL_VEC_D *uh_loc);

REAL_B *eval_bar_D2_uh(REAL_BB result, const REAL_B lambda,
                       const EL_REAL_VEC *uh_loc, const BAS_FCTS *bfcts);
REAL_BB *eval_bar_D2_uh_d(REAL_DBB result, const REAL_B lambda,
                          const EL_REAL_D_VEC *uh_loc, const BAS_FCTS *bfcts);
REAL_BB *eval_bar_D2_uh_dow(REAL_DBB result, const REAL_B lambda,
                            const EL_REAL_VEC_D *uh_loc,
                            const BAS_FCTS *bfcts);

REAL_B *eval_bar_D2_uh_fast(REAL_BB result, const EL_REAL_VEC *uh_loc,
                            const QUAD_FAST *qfast, int iq);
REAL_BB *eval_bar_D2_uh_d_fast(REAL_DBB result, const EL_REAL_D_VEC *uh_loc,
                               const QUAD_FAST *qfast, int iq, bool update)
REAL_BB *eval_bar_D2_uh_dow_fast(REAL_DBB result, const EL_REAL_VEC_D *uh_loc,
                                 const QUAD_FAST *qfast, int iq);

REAL_BB *bar_D2_uh_at_qp(REAL_BB *result, const QUAD_FAST *qfast,
                         const EL_REAL_VEC *uh_loc);
REAL_DBB *bar_D2_uh_d_at_qp(REAL_DBB vec[], const QUAD_FAST *qfast,
                            const EL_REAL_D_VEC *uh_loc);
REAL_DBB *bar_D2_uh_dow_at_qp(REAL_DBB vec[], const QUAD_FAST *qfast,
                              const EL_REAL_VEC_D *uh_loc);
\end{lstlisting}\ev

Description: These functions compute the respective derivatives with
respect to barycentric co-ordinates. Otherwise they are functionally
equivalent to the functions without the \code{bar\_}-prefix.

%\begin{descr}
%\kitem{eval\_bar\_grd\_uh(result, lambda, uh\_loc, bfcts)}
%\kitem{eval\_bar\_grd\_uh\_d(result, lambda, uh\_loc, bfcts)}
%\kitem{eval\_bar\_grd\_uh\_dow(result, lambda, uh\_loc, bfcts)}

%\kitem{eval\_bar\_grd\_uh\_fast(result, uh\_loc, qfast, iq)}
%\kitem{eval\_bar\_grd\_uh\_d\_fast(result, uh\_loc, qfast, iq)}
%\kitem{eval\_bar\_grd\_uh\_dow\_fast(result, uh\_loc, qfast, iq)}

%\kitem{bar\_grd\_uh\_at\_qp(result, qfast, uh\_loc)}
%\kitem{bar\_grd\_uh\_d\_at\_qp(result, qfast, uh\_loc)}
%\kitem{bar\_grd\_uh\_dow\_at\_qp(result, qfast, uh\_loc)}

%\kitem{eval\_bar\_D2\_uh(result, lambda, uh\_loc, bfcts)}
%\kitem{eval\_bar\_D2\_uh\_d(result, lambda, uh\_loc, bfcts)}
%\kitem{eval\_bar\_D2\_uh\_dow(result, lambda, uh\_loc, bfcts)}

%\kitem{eval\_bar\_D2\_uh\_fast(result, uh\_loc, qfast, iq)}
%\kitem{eval\_bar\_D2\_uh\_d\_fast(result, uh\_loc, qfast, iq)}
%\kitem{eval\_bar\_D2\_uh\_dow\_fast(result, uh\_loc, qfast, iq)}

%\kitem{bar\_D2\_uh\_at\_qp(result, qfast, uh\_loc)}
%\kitem{bar\_D2\_uh\_d\_at\_qp(result, qfast, uh\_loc)}
%\kitem{bar\_D2\_uh\_dow\_at\_qp(result, qfast, uh\_loc)}
%\end{descr}

\section{Calculation of norms for finite element functions}\label{S:eval_norm}

\ALBERTA supplies functions for the calculation of the $L^2$ norm and
$H^1$ semi--norm of a given scalar or vector valued finite element
function.
%% 
\fdx{H1_norm_uh()@{\code{H1\_norm\_uh()}}}
\idx{evaluation tools!H1_norm_uh()@{\code{H1\_norm\_uh()}}}
\fdx{L2_norm_uh()@{\code{L2\_norm\_uh()}}}
\idx{evaluation tools!L2_norm_uh()@{\code{L2\_norm\_uh()}}}
\fdx{H1_norm_uh_d()@{\code{H1\_norm\_uh\_d()}}}
\idx{evaluation tools!H1_norm_uh_d()@{\code{H1\_norm\_uh\_d()}}}
\fdx{L2_norm_uh_d()@{\code{L2\_norm\_uh\_d()}}}
\idx{evaluation tools!L2_norm_uh_d()@{\code{L2\_norm\_uh\_d()}}}
\fdx{H1_norm_uh_dow()@{\code{H1\_norm\_uh\_dow()}}}
\idx{evaluation tools!H1_norm_uh_dow()@{\code{H1\_norm\_uh\_dow()}}}
\fdx{L2_norm_uh_dow()@{\code{L2\_norm\_uh\_dow()}}}
\idx{evaluation tools!L2_norm_uh_dow()@{\code{L2\_norm\_uh\_dow()}}}
\bv\begin{lstlisting}[label=P:DISC_FCT_NORMS]
REAL H1_norm_uh(const QUAD *, const DOF_REAL_VEC *);
REAL L2_norm_uh(const QUAD *, const DOF_REAL_VEC *);
REAL H1_norm_uh_d(const QUAD *, const DOF_REAL_D_VEC *);
REAL L2_norm_uh_d(const QUAD *, const DOF_REAL_D_VEC *);
REAL H1_norm_uh_dow(const QUAD *, const DOF_REAL_VEC_D *);
REAL L2_norm_uh_dow(const QUAD *, const DOF_REAL_VEC_D *);
\end{lstlisting}\ev
\paragraph{Descriptions}
\begin{descr}
\kitem{H1\_norm\_uh(quad, uh)}
%%
\fdx{H1_norm_uh()@{\code{H1\_norm\_uh()}}}
\idx{evaluation tools!H1_norm_uh()@{\code{H1\_norm\_uh()}}}
%%
returns an approximation to the $H^1$ semi norm $(\int_\Omega |\nabla
u_h|^2)^{1/2}$ of a finite element function; the coefficient vector of
the vector is stored in \code{uh}; the domain is given by
\code{uh->fe\_space->mesh}; the element integrals are approximated by
the numerical quadrature \code{quad}, if \code{quad} is not \nil;
otherwise a quadrature which is exact of degree
\code{2*uh->fe\_space->bas\_fcts->degree-2} is used.
%%
\kitem{L2\_norm\_uh(quad, uh)}
%%
\fdx{L2_norm_uh()@{\code{L2\_norm\_uh()}}}
\idx{evaluation tools!L2_norm_uh()@{\code{L2\_norm\_uh()}}}
%%
returns an approximation to the $L^2$ norm $(\int_\Omega
|u_h|^2)^{1/2}$ of a finite element function; the coefficient vector
of the vector is stored in \code{uh}; the domain is given by
\code{uh->fe\_space->mesh}; the element integrals are approximated by
the numerical quadrature \code{quad}, if \code{quad} is not \nil;
otherwise a quadrature which is exact of degree
\code{2*uh->fe\_space->bas\_fcts->degree} is used.
%%
\kitem{H1\_norm\_uh\_[d|dow](quad, uh\_d)}
%%
\fdx{H1_norm_uh_d()@{\code{H1\_norm\_uh\_d()}}}
\idx{evaluation tools!H1_norm_uh_d()@{\code{H1\_norm\_uh\_d()}}}
\fdx{H1_norm_uh_dow()@{\code{H1\_norm\_uh\_dow()}}}
\idx{evaluation tools!H1_norm_uh_dow()@{\code{H1\_norm\_uh\_dow()}}}
%%
returns an approximation to the $H^1$ semi norm of a vector valued
finite element function; the coefficient vector of the vector is
stored in \code{uh\_d}; the domain is given by
\code{uh\_d->fe\_space->mesh}; the element integrals are approximated
by the numerical quadrature \code{quad}, if \code{quad} is not \nil;
otherwise a quadrature which is exact of degree
\code{2*uh\_d->fe\_space->bas\_fcts->degree-2} is used.
%%
\kitem{L2\_norm\_uh\_[d|dow](quad, uh\_d)}
%%
\fdx{L2_norm_uh_d()@{\code{L2\_norm\_uh\_d()}}}
\idx{evaluation tools!L2_norm_uh_d()@{\code{L2\_norm\_uh\_d()}}}
\fdx{L2_norm_uh_dow()@{\code{L2\_norm\_uh\_dow()}}}
\idx{evaluation tools!L2_norm_uh_dow()@{\code{L2\_norm\_uh\_dow()}}}
%%
returns an approximation to the $L^2$ norm of a vector valued finite
element function; the coefficient vector of the vector is stored in
\code{uh\_d}; the domain is given by \code{uh\_d->fe\_space->mesh};
the element integrals are approximated by the numerical quadrature
\code{quad}, if \code{quad} is not \nil; otherwise a quadrature which
is exact of degree \code{2*uh\_d->fe\_space->bas\_fcts->degree} is
used.
\end{descr}

\section{Interface for application provided functions}
\label{S:user_fcts} %% oh, what a user ...

Often the library function in the \ALBERTA package require certain
application provided functions, e.g. for assembling the ``right hand
side'', for computations of the ``true'' error, for inhomogeneous
boundary conditions or for interpolation of (non-discrete) functions
onto finite element spaces. This section defines some basis calling
conventions concerning these application provided functions.

Most of these function must conform to one of the following proto-types:
%%
\ddx{FCT_AT_X@{\code{FCT\_AT\_X}}}
\idx{application function types!FCT_AT_X@{\code{FCT\_AT\_X}}}
\ddx{GRD_FCT_AT_X@{\code{GRD\_FCT\_AT\_X}}}
\idx{application function types!GRD\_FCT_AT_X@{\code{GRD\_FCT\_AT\_X}}}
\ddx{D2_FCT_AT_X@{\code{D2\_FCT\_AT\_X}}}
\idx{application function types!D2\_FCT_AT_X@{\code{D2\_FCT\_AT\_X}}}
\ddx{FCT_D_AT_X@{\code{FCT\_D\_AT\_X}}}
\idx{application function types!FCT_D_AT_X@{\code{FCT\_D\_AT\_X}}}
\ddx{GRD_FCT_D_AT_X@{\code{GRD\_FCT\_D\_AT\_X}}}
\idx{application function types!GRD\_FCT_D_AT_X@{\code{GRD\_FCT\_D\_AT\_X}}}
\ddx{D2_FCT_D_AT_X@{\code{D2\_FCT\_D\_AT\_X}}}
\idx{application function types!D2\_FCT_D_AT_X@{\code{D2\_FCT\_D\_AT\_X}}}
\ddx{LOC_FCT_AT_QP@{\code{LOC\_FCT\_AT\_QP}}}
\idx{application function types!LOC_FCT_AT_QP@{\code{LOC\_FCT\_AT\_QP}}}
\ddx{GRD_LOC_FCT_AT_QP@{\code{GRD\_LOC\_FCT\_AT\_QP}}}
\idx{application function types!GRD_LOC_FCT_AT_QP@{\code{GRD\_LOC\_FCT\_AT\_QP}}}
\ddx{LOC_FCT_D_AT_QP@{\code{LOC\_FCT\_D\_AT\_QP}}}
\idx{application function types!LOC_FCT_D_AT_QP@{\code{LOC\_FCT\_D\_AT\_QP}}}
\ddx{GRD_LOC_FCT_D_AT_QP@{\code{GRD\_LOC\_FCT\_D\_AT\_QP}}}
\idx{application function types!GRD_LOC_FCT_D_AT_QP@{\code{GRD\_LOC\_FCT\_D\_AT\_QP}}}
%%
\bv\begin{lstlisting}
typedef REAL (*FCT_AT_X)(const REAL_D x);
typedef const REAL *(*GRD_FCT_AT_X)(const REAL_D x, REAL_D result);
typedef const REAL_D *(*D2_FCT_AT_X)(const REAL_D x, REAL_DD result);

typedef const REAL    *(*FCT_D_AT_X)(const REAL_D x, REAL_D result);
typedef const REAL_D  *(*GRD_FCT_D_AT_X)(const REAL_D x, REAL_DD result);
typedef const REAL_DD *(*D2_FCT_D_AT_X)(const REAL_D x, REAL_DDD result);

typedef REAL (*LOC_FCT_AT_QP)(const EL_INFO *el_info,
			      const QUAD *quad, int iq,
			      void *ud);
typedef const REAL *(*LOC_FCT_D_AT_QP)(REAL_D result,
				       const EL_INFO *el_info,
				       const QUAD *quad, int iq,
				       void *ud);
typedef const REAL *(*GRD_LOC_FCT_AT_QP)(REAL_D res,
					 const EL_INFO *el_info,
					 const REAL_BD Lambda,
					 const QUAD *quad, int iq,
					 void *ud);
typedef const REAL_D *(*GRD_LOC_FCT_D_AT_QP)(REAL_DD res,
					     const EL_INFO *el_info,
					     const REAL_BD Lambda,
					     const QUAD *quad, int iq,
					     void *ud);
\end{lstlisting}\ev

\begin{datatype}{FCT\_AT\_X}
\label{S:FCT_AT_X_fptr}
\ddx{FCT_AT_X@{\code{FCT\_AT\_X}}}
\idx{application function types!FCT_AT_X@{\code{FCT\_AT\_X}}}

\item[Prototype]~\hfill

\begin{lstlisting}
typedef REAL (*FCT_AT_X)(const REAL_D x);
\end{lstlisting}

\item[Synopsis]~\hfill

\begin{lstlisting}
FCT_AT_X fptr;

result = fptr(x);
\end{lstlisting}

\item[Description]~\hfill

  Evaluate at the point \code{x} and return a
  scalar value. This is the simplest function-type.
  %% 

\item[Parameters]~\hfill

  \begin{descr}
  \kitem{x} The point of evaluation.
  \end{descr}

\item[Return Value]~\hfill

  The function value.

\end{datatype}

\begin{datatype}{GRD\_FCT\_AT\_X}
\label{S:GRD_FCT_AT_X_fptr}
\ddx{GRD_FCT_AT_X@{\code{GRD\_FCT\_AT\_X}}}
\idx{application function types!GRD\_FCT_AT_X@{\code{GRD\_FCT\_AT\_X}}}

\item[Prototype]~\hfill

\begin{lstlisting}
const REAL *GRD_FCT_AT_X(const REAL_D x, REAL\_D result);
\end{lstlisting}

\item[Synopsis]~\hfill

\begin{lstlisting}
GRD_FCT_AT_X fptr;

result = fptr(x, result);
result = fptr(x, NULL);
\end{lstlisting}

\item[Description]~\hfill

  Evaluate the first derivative at the point \code{x}.
\item[Parameters]~\hfill
  \begin{descr}
  \kitem{x} The point of evaluation, in Cartesian coordinates.
  \kitem{result} Storage for the result, or \nil.
  \end{descr}
\item[Return Value] The address of \code{result}, if \code{result !=
    \nil}, otherwise a pointer to a statically allocated storage
  area, see \exampleref{E:grd_user_fct} below.
%%
\end{datatype}

\begin{example}
  \label{E:grd_user_fct}
  \bv\begin{lstlisting}[name=GRD_FCT_AT_X_ABI,label=C:GRD_FCT_AT_X_ABI]
      const REAL *grd_g_implementation(const REAL_D x, REAL_D result) {
        static REAL_D storage; /* mind the "static" key-word!!! */

        if (result == NULL) {
          result = storage;
        }
        
        ... /* mighty complicated computations for "result" */
        
        return result;
      }
\end{lstlisting}\ev
\end{example}

\hrulefill

\begin{descr}

  \kitem{const REAL\_D *D2\_FCT\_AT\_X(const REAL\_D x, REAL\_DD result)}
  Evaluate the second derivative at the point \code{x}.
  \begin{description}
  \item[Parameters]~\hfill
    \begin{descr}
      \kitem{x} The point of evaluation, in Cartesian coordinates.
      \kitem{result} Storage for the result, or \nil.
    \end{descr}
  \item[Return Value] The address of \code{result}, if \code{result !=
      \nil}, otherwise a pointer to a statically allocated storage
    area, see \exampleref{E:D2_user_fct} below.
  \end{description}
  \begin{example}
    \label{E:D2_user_fct}
    \bv\begin{lstlisting}[name=D2_FCT_AT_X_ABI,label=C:D2_FCT_AT_X_ABI]
      const REAL_D *D2_g_implementation(const REAL_D x, REAL_DD result) {
        static REAL_DD storage; /* mind the "static" key-word!!! */

        if (result == NULL) {
          result = storage;
        }
        
        ... /* mighty complicated computations for "result" */
        
        return (const REAL_D *)result;
      }
    \end{lstlisting}\ev
  \end{example}
  %%

  \hrulefill

  \kitem{const REAL *FCT\_D\_AT\_X(const REAL\_D x, REAL\_D result)}
  \kitem{const REAL\_D *GRD\_FCT\_D\_AT\_X(const REAL\_D x, REAL\_DD result)}
  \kitem{const REAL\_DD *D2\_FCT\_D\_AT\_X(const REAL\_D x, REAL\_DDD result)}
  %%
  ~\hfill
  Evaluate a vector valued function at the point \code{x}. There is, of
  course, no difference between the \code{GRD\_FCT\_AT\_X} and the
  \code{FCT\_D\_AT\_X} function pointers.
  \begin{description}
  \item[Parameters]~\hfill
    \begin{descr}
      \kitem{x} The point of evaluation, in Cartesian coordinates.
      \kitem{result} Storage for the result, or \nil.
    \end{descr}
  \item[Return Value] The address of \code{result}, if \code{result !=
      \nil}, otherwise a pointer to a statically allocated storage
    area, see \exampleref{E:user_fct_d} below.
  \end{description}
  \begin{example}
    \label{E:user_fct_d}
    \bv\begin{lstlisting}[name=FCT_D_AT_X_ABI,label=C:FCT_D_AT_X_ABI]
      const REAL *g_implementation(const REAL_D x, REAL_D result) {
        static REAL_D storage; /* mind the "static" key-word!!! */

        if (result == NULL) {
          result = storage;
        }
        
        ... /* mighty complicated computations for "result" */
        
        return result;
      }
    \end{lstlisting}\ev
  \end{example}
  
\pagebreak[2]
  \kitem{REAL LOC\_FCT\_AT\_QP(}
  \kitem{           const EL\_INFO *el\_info, const QUAD *quad, int iq, void *ud)}
  ~\hfill

  %% 
  Evaluate the function at \code{quad->lambda[iq]}. This looks
  slightly more complicated than the simple \code{FCT\_AT\_X} types,
  but passing the \code{EL\_INFO} descriptor along with quadrature
  rule opens the door to implement even complicated functions in an
  efficient and simpler way than is possible with the simple
  \code{FCT\_AT\_X} types. See also \exampleref{E:loc_fct_at_qp}.
  \begin{description}
  \item[Parameters]~\hfill
    \begin{descr}
      \kitem{el\_info} The current \code{EL\_INFO} descriptor.
      \kitem{quad} The quadrature rule storing the evaluation points.
      \kitem{iq} The number of the evaluation point.
      \kitem{ud} Application data pointer.
    \end{descr}
  \item[Return Value] The function value.
  \end{description}

  \hrulefill

  \kitem{const REAL *GRD\_LOC\_FCT\_AT\_QP(}
  \kitem{              REAL\_D res, const EL\_INFO *el\_info, const QUAD *quad, int iq, void *ud)}
  \kitem{const REAL *LOC\_FCT\_D\_AT\_QP(}
  \kitem{              REAL\_D res, const EL\_INFO *el\_info, const QUAD *quad, int iq, void *ud)}
  \kitem{const REAL\_D *GRD\_LOC\_FCT\_D\_AT\_QP(}
  \kitem{              REAL\_DD res, const EL\_INFO *el\_info, const QUAD *quad, int iq, void *ud)}
  %%
  More or less self-explanatory, the convention for the \code{res}
  argument are the same as for the \code{FCT\_AT\_X} types: a \nil
  pointer must be accepted, and then a pointer to a statically
  allocated storage area has to be returned, otherwise the result has
  to be stored in \code{res}, and the return value must be \code{res},
  too.

\end{descr}


\section{Calculation of errors of finite element approximations}
\label{S:error_cal}

For test purposes it is convenient to calculate the ``exact error''
between a finite element approximation and the exact solution.
\ALBERTA supplies functions to calculate the error in several norms.
For test purposes, the integral error routines may be used as ``error
estimators'' in an adaptive method. The local element error $\int_S
|\nabla (u-u_h)|^2$ or $\int_S |u-u_h|^2$ can be used as an error
indicator and can be stored on the element leaf data, for example.
\ALBERTA provides also functions for the computation of the mean value
of a given function, respectively the mean value difference of a given
non-discrete function and a discrete function.

Not all variants of the functions listed below will be explained in
detail further below. For the calling conventions for the application
supplied function pointer we refer the reader to \secref{S:user_fcts}
on page \pageref{S:user_fcts}.
%%
\fdx{max_err_at_qp()@{\code{max\_err\_at\_qp()}}}
\idx{evaluation tools!max_err_at_qp()@{\code{max\_err\_at\_qp()}}}
\fdx{max_err_at_qp_loc()@{\code{max\_err\_at\_qp\_loc()}}}
\idx{evaluation tools!max_err_at_qp_loc()@{\code{max\_err\_at\_qp\_loc()}}}
\fdx{max_err_dow_at_qp()@{\code{max\_err\_dow\_at\_qp()}}}
\idx{evaluation tools!max_err_dow_at_qp()@{\code{max\_err\_dow\_at\_qp()}}}
\fdx{max_err_dow_at_qp_loc()@{\code{max\_err\_dow\_at\_qp\_loc()}}}
\idx{evaluation tools!max_err_dow_at_qp_loc()@{\code{max\_err\_dow\_at\_qp\_loc()}}}
%%
\fdx{max_err_at_vert()@{\code{max\_err\_at\_vert()}}}
\idx{evaluation tools!max_err_at_vert()@{\code{max\_err\_at\_vert()}}}
\fdx{max_err_at_vert_loc()@{\code{max\_err\_at\_vert\_loc()}}}
\idx{evaluation tools!max_err_at_vert_loc()@{\code{max\_err\_at\_vert\_loc()}}}
\fdx{max_err_dow_at_vert()@{\code{max\_err\_dow\_at\_vert()}}}
\idx{evaluation tools!max_err_dow_at_vert()@{\code{max\_err\_dow\_at\_vert()}}}
\fdx{max_err_dow_at_vert_loc()@{\code{max\_err\_dow\_at\_vert\_loc()}}}
\idx{evaluation tools!max_err_dow_at_vert_loc()@{\code{max\_err\_dow\_at\_vert\_loc()}}}
%%
\fdx{L2_err()@{\code{L2\_err()}}}
\idx{evaluation tools!L2_err()@{\code{L2\_err()}}}
\fdx{L2_err_loc()@{\code{L2\_err\_loc()}}}
\idx{evaluation tools!L2_err_loc()@{\code{L2\_err\_loc()}}}
\fdx{L2_err_weighted()@{\code{L2\_err\_weighted()}}}
\idx{evaluation tools!L2_err_weighted()@{\code{L2\_err\_weighted()}}}
\fdx{L2_err_dow()@{\code{L2\_err\_dow()}}}
\idx{evaluation tools!L2_err_dow()@{\code{L2\_err\_dow()}}}
\fdx{L2_err_loc_dow()@{\code{L2\_err\_loc\_dow()}}}
\idx{evaluation tools!L2_err_loc_dow()@{\code{L2\_err\_loc\_dow()}}}
\fdx{L2_err_dow_weighted()@{\code{L2\_err\_dow\_weighted()}}}
\idx{evaluation tools!L2_err_dow_weighted()@{\code{L2\_err\_dow\_weighted()}}}
%%
\fdx{H1_err()@{\code{H1\_err()}}}
\idx{evaluation tools!H1_err()@{\code{H1\_err()}}}
\fdx{H1_err_loc()@{\code{H1\_err\_loc()}}}
\idx{evaluation tools!H1_err_loc()@{\code{H1\_err\_loc()}}}
\fdx{H1_err_weighted()@{\code{H1\_err\_weighted()}}}
\idx{evaluation tools!H1_err_weighted()@{\code{H1\_err\_weighted()}}}
\fdx{H1_err_dow()@{\code{H1\_err\_dow()}}}
\idx{evaluation tools!H1_err_dow()@{\code{H1\_err\_dow()}}}
\fdx{H1_err_loc_dow()@{\code{H1\_err\_loc\_dow()}}}
\idx{evaluation tools!H1_err_loc_dow()@{\code{H1\_err\_loc\_dow()}}}
\fdx{H1_err_dow_weighted()@{\code{H1\_err\_dow\_weighted()}}}
\idx{evaluation tools!H1_err_dow_weighted()@{\code{H1\_err\_dow\_weighted()}}}
%%
\fdx{mean_value()@{\code{mean\_value()}}}
\idx{evaluation tools!mean_value()@{\code{mean\_value()}}}
\fdx{mean_value_dow()@{\code{mean\_value\_dow()}}}
\idx{evaluation tools!mean_value_dow()@{\code{mean\_value\_dow()}}}
\fdx{mean_value_loc()@{\code{mean\_value\_loc()}}}
\idx{evaluation tools!mean_value_loc()@{\code{mean\_value\_loc()}}}
\fdx{mean_value_loc_dow()@{\code{mean\_value\_loc\_dow()}}}
\idx{evaluation tools!mean_value_loc_dow()@{\code{mean\_value\_loc\_dow()}}}
\bv\begin{lstlisting}
REAL max_err_at_qp(FCT_AT_X u, const DOF_REAL_VEC *uh, const QUAD *quad);
REAL max_err_at_qp_loc(LOC_FCT_AT_QP u_loc, void *ud, FLAGS fill_flag,
		       const DOF_REAL_VEC *uh,
		       const QUAD *quad);
REAL max_err_dow_at_qp(FCT_D_AT_X u,
		       const DOF_REAL_VEC_D *uh,
		       const QUAD *quad);
REAL max_err_dow_at_qp_loc(LOC_FCT_D_AT_QP u_loc,
			   void *ud, FLAGS fill_flag,
			   const DOF_REAL_VEC_D *uh,
			   const QUAD *quad);

REAL max_err_at_vert(FCT_AT_X u, const DOF_REAL_VEC *uh);
REAL max_err_at_vert_loc(LOC_FCT_AT_QP u_at_qp,
			 void *ud, FLAGS fill_flag,
			 const DOF_REAL_VEC *uh);
REAL max_err_dow_at_vert(FCT_D_AT_X u, const DOF_REAL_VEC_D *uh);
REAL max_err_dow_at_vert_loc(LOC_FCT_D_AT_QP u_at_qp,
			     void *ud, FLAGS fill_flag,
			     const DOF_REAL_VEC_D *uh);

REAL L2_err(FCT_AT_X u, const DOF_REAL_VEC *uh,
	    const QUAD *quad,
	    bool rel_err, bool mean_value_adjust,
	    REAL *(*rw_err_el)(EL *el), REAL *max_l2_err2);
REAL L2_err_loc(LOC_FCT_AT_QP u_loc, void *ud, FLAGS fill_flag,
		const DOF_REAL_VEC *uh,
		const QUAD *quad,
		bool rel_err, bool mean_value_adjust,
		REAL *(*rw_err_el)(EL *el), REAL *max_l2_err2);
REAL L2_err_weighted(FCT_AT_X weight, FCT_AT_X u, const DOF_REAL_VEC *uh,
		     const QUAD *quad,
		     bool rel_err, bool mean_value_adjust,
		     REAL *(*rw_err_el)(EL *el), REAL *max_l2_err2);
REAL L2_err_dow(FCT_D_AT_X u,
		const DOF_REAL_VEC_D *uh,
		const QUAD *quad,
		bool rel_err, bool mean_value_adjust,
		REAL *(*rw_err_el)(EL *el), REAL *max_l2_err2);
REAL L2_err_loc_dow(LOC_FCT_D_AT_QP u_loc,
		    void *ud, FLAGS fill_flag,
		    const DOF_REAL_VEC_D *uh,
		    const QUAD *quad,
		    bool rel_err, bool mean_value_adjust,
		    REAL *(*rw_err_el)(EL *el), REAL *max_l2_err2);
REAL L2_err_dow_weighted(FCT_AT_X weight, FCT_D_AT_X u,
			 const DOF_REAL_VEC_D *uh,
			 const QUAD *quad,
			 bool rel_err, bool mean_value_adjust,
			 REAL *(*rw_err_el)(EL *el), REAL *max_l2_err2);

REAL H1_err(GRD_FCT_AT_X grd_u, const DOF_REAL_VEC *uh, 
	    const QUAD *quad, bool rel_err, REAL *(*rw_err_el)(EL *),
	    REAL *max_el_err2);
REAL H1_err_loc(GRD_LOC_FCT_AT_QP grd_u_loc,
		void *ud, FLAGS fill_flag,
		const DOF_REAL_VEC *uh, 
		const QUAD *quad, bool rel_err, REAL *(*rw_err_el)(EL *),
		REAL *max_el_err2);
REAL H1_err_weighted(FCT_AT_X weight, GRD_FCT_AT_X grd_u,
		     const DOF_REAL_VEC *uh, const QUAD *quad,
		     bool rel_err, REAL *(*rw_err_el)(EL *),
		     REAL *max_el_err2);
REAL H1_err_dow(GRD_FCT_D_AT_X grd_u, const DOF_REAL_VEC_D *uh,
		const QUAD *quad,
		bool rel_err, REAL *(*rw_err_el)(EL *), REAL *max_el_err2);
REAL H1_err_loc_dow(GRD_LOC_FCT_D_AT_QP grd_u_loc, void *ud, FLAGS fill_flag,
		    const DOF_REAL_VEC_D *uh, const QUAD *quad,
		    bool rel_err, 
		    REAL *(*rw_err_el)(EL *), REAL *max_el_err2);
REAL H1_err_dow_weighted(FCT_AT_X weight, GRD_FCT_D_AT_X grd_u,
			 const DOF_REAL_VEC_D *uh,
			 const QUAD *quad,
			 bool rel_err, REAL *(*rw_err_el)(EL *),
			 REAL *max_el_err2);

REAL mean_value(MESH *mesh, REAL (*f)(const REAL_D), const DOF_REAL_VEC *fh,
                const QUAD *quad);
REAL mean_value_loc(MESH *mesh, LOC_FCT_AT_QP f_at_qp,
                    void *ud, FLAGS fill_flags,
                    const DOF_REAL_VEC *fh, const QUAD *quad);
const REAL *mean_value_dow(MESH *mesh, FCT_D_AT_X f, const DOF_REAL_VEC_D *fh,
                           const QUAD *quad, REAL_D mean);
const REAL *mean_value_loc_dow(REAL_D mean, MESH *mesh,
                               LOC_FCT_D_AT_QP f_at_qp,
                               void *ud, FLAGS fill_flag,
                               const DOF_REAL_VEC_D *fh,
                               const QUAD *quad);
\end{lstlisting}\ev
\paragraph{Descriptions}
\begin{descr}
\kitem{max\_err\_at\_qp(u, uh, quad)}
%%
the function returns the maximal error, $\max |u-u_h|$, between the
true solution and the approximation at all quadrature nodes on all
elements of a mesh; \code{u} is a pointer to a function for the
evaluation of the true solution, \code{uh} stores the coefficients of
the approximation, \code{uh->fe\_space->mesh} is the underlying mesh,
and \code{quad} is the quadrature which gives the quadrature nodes; if
\code{quad} is \nil, a quadrature which is exact of degree
\code{2*uh->fe\_space->bas\_fcts->degree-2} is used.
%%
\kitem{H1\_err(grd\_u, uh, quad, rel\_err, rw\_el\_err, max)}
%%
the function returns an approximation to the absolute error
$(\int_\Omega |\nabla (u-u_h)|^2)^{1/2}$ (\code{rel\_err == 0}) or
relative error $(\int_\Omega|\nabla(u-u_h)|^2/\int_\Omega|\nabla
u|^2)^{1/2}$ (\code{rel\_err == 1}) between the true solution and the
approximation in the $H^1$ semi norm;
       
\code{grd\_u} is a pointer to a function for the evaluation of the
gradient of the true solution returning a \DOW vector storing this
gradient, \code{uh} stores the coefficients of the approximation,
\code{uh->fe\_space->mesh} is the underlying mesh, and \code{quad} is
the quadrature for the approximation of the element integrals; if
\code{quad} is \nil, a quadrature which is exact of degree
\code{2*uh->fe\_space->bas\_fcts->degree-2} is used;

if \code{rw\_el\_err} is not \nil, the return value of
\code{(*rw\_el\_err)(el)} provides for each mesh element \code{el} an
address where the local error is stored; if \code{max} is not \nil,
\code{*max} is the maximal local error on an element on output.
%%
\kitem{L2\_err(u, uh, quad, rel\_err, rw\_el\_err, max)}
%%
the function returns an approximation to the absolute error
$(\int_\Omega |u-u_h|^2)^{1/2}$ (\code{rel\_err == 0}) or the relative
error $(\int_\Omega|u-u_h|^2/\int_\Omega|u|^2)^{1/2}$ (\code{rel\_err
  == 1}) between the true solution and the approximation in the $L^2$
norm,
       
\code{u} is a pointer to a function for the evaluation of the true
solution, \code{uh} stores the coefficients of the approximation,
\code{uh->fe\_space->mesh} is the underlying mesh, and \code{quad} is
the quadrature for the approximation of the element integrals; if
\code{quad} is \nil, a quadrature which is exact of degree
\code{2*uh->fe\_space->bas\_fcts->degree-2} is used;

if \code{rw\_el\_err} is not \nil, the return value of
\code{(*rw\_el\_err)(el)} provides for each mesh element \code{el} an
address where the local error is stored; if \code{max} is not \nil,
\code{*max} is the maximal local error on an element on output.
%%
\kitem{max\_err\_at\_qp\_[d|dow](u\_d, uh\_d, quad)}
%%
the function returns the maximal error between the true solution and
the approximation at all quadrature nodes on all elements of a mesh;
\code{u\_d} is a pointer to a function for the evaluation of the true
solution returning a \DOW vector storing the value of the function,
\code{uh\_d} stores the coefficients of the approximation,
\code{uh\_d->fe\_space->mesh} is the underlying mesh, and \code{quad}
is the quadrature which gives the quadrature nodes; if \code{quad} is
\nil, a quadrature which is exact of degree
\code{2*uh\_d->fe\_space->bas\_fcts->degree-2} is used.
%%
\kitem{H1\_err2\_[d|dow](grd\_u\_d, uh\_d, quad, rel\_err, rw\_el\_err, max)} 
%%
the function returns an approximation to the absolute error
(\code{rel\_err == 0}) or relative error (\code{rel\_err == 1})
between the true solution and the approximation in the $H^1$ semi
norm;

\code{grd\_u\_d} is a pointer to a function for the evaluation of the
Jacobian of the true solution returning a $\DOW\times\DOW$ matrix
storing this Jacobian, \code{uh\_d} stores the coefficients of the
approximation, \code{uh\_d->fe\_space->mesh} is the underlying mesh,
and \code{quad} is the quadrature for the approximation of the element
integrals; if \code{quad} is \nil, a quadrature which is exact of
degree \code{2*uh\_d->fe\_space->bas\_fcts->degree-2} is used;

if \code{rw\_el\_err} is not \nil, the return value of
\code{(*rw\_el\_err)(el)} provides for each mesh element \code{el} an
address where the local error is stored; if \code{max} is not \nil,
\code{*max} is the maximal local error on an element on output.

\kitem{L2\_err2\_[d|dow](u\_d, uh\_d, quad, rel\_err, rw\_el\_err, max)} 
%%
the function returns an approximation to the absolute error
(\code{rel\_err == 0}) or relative error (\code{rel\_err == 1})
between the true solution and the approximation in the $L^2$ norm;

\code{u\_d} is a pointer to a function for the evaluation of the true
solution returning a \DOW vector storing the value of the function,
\code{uh\_d} stores the coefficients of the approximation,
\code{uh\_d->fe\_space->mesh} is the underlying mesh, and \code{quad}
is the quadrature for the approximation of the element integrals; if
\code{quad} is \nil, a quadrature which is exact of degree
\code{2*uh\_d->fe\_space->bas\_fcts->degree-2} is used;

if \code{rw\_el\_err} is not \nil, the return value of
\code{(*rw\_el\_err)(el)} provides for each mesh element \code{el} an
address where the local error is stored; if \code{max} is not \nil,
\code{*max} is the maximal local error on an element on output.
%%
\kitem{mean\_value(mesh, f, fh, quad)}
%%
\fdx{mean_value()@{\code{mean\_value()}}}
\idx{evaluation tools!mean_value()@{\code{mean\_value()}}}
%%
\kitem{mean\_value\_[d|dow](mesh, f, fh, quad, mean)}
%%
\fdx{mean_value_dow()@{\code{mean\_value\_dow()}}}
\idx{evaluation tools!mean_value_dow()@{\code{mean\_value\_dow()}}}
%%
\kitem{mean\_value\_loc(mesh, f\_at\_qp, ud, fill\_flags, fh, quad)}
%%
\fdx{mean_value_loc()@{\code{mean\_value\_loc()}}}
\idx{evaluation tools!mean_value_loc()@{\code{mean\_value\_loc()}}}
%%
\kitem{mean\_value\_loc\_[d|dow](mean, mesh, f\_at\_qp, ud, fill\_flags, fh, quad)}
%%
\fdx{mean_value_loc_dow()@{\code{mean\_value\_loc\_dow()}}}
\idx{evaluation tools!mean_value_loc_dow()@{\code{mean\_value\_loc\_dow()}}}
%%
compute the mean value of either a finite element function or a
non-discrete function. If both are given return the difference of
their mean values \code{(f-fh)}.
\end{descr}

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