File: assemble-tools.tex

package info (click to toggle)
alberta 3.1.1-3
links: PTS, VCS
area: main
in suites: forky, sid, trixie
size: 19,176 kB
sloc: ansic: 135,836; cpp: 6,601; makefile: 2,801; sh: 333; fortran: 180; lisp: 177; xml: 30
file content (4076 lines) | stat: -rw-r--r-- 187,839 bytes
\def\vL{\vec{L}}

\section{Tools for the assemblage of linear systems}%
\label{S:ass_tools}%
\idx{assemblage tools|(}

This section describes data structures and subroutines for matrix and
vector assembly. Section \ref{S:update} presents basic routines for
the update of global matrices and vectors by adding contributions from
one single element. Data structures and routines for global matrix
assembly are described in Section \ref{S:matrix_assemblage}. This
includes library routines for the efficient implementation of a
general second order linear elliptic operator. Section
\ref{S:ass_info} presents data structures and routines for the
handling of pre--computed integrals, which are used to speed up
calculations in the case of problems with constant coefficients. The
assembly of (right hand side) vectors is described in Section
\ref{S:vector_update}. The incorporation of Dirichlet boundary values into
the right hand side is presented in Section \ref{S:dirichlet_bound}.
Finally, routines for generation of interpolation coefficients are
described in Section \ref{S:I_FES}.

\subsection{Element matrices and vectors}%
\label{S:update}

The usual way to assemble the system matrix and the load vector is to
loop over all (leaf) elements, calculate the local element
contributions and add these to the global system matrix and the global
load vector. The updating of the load vector is rather easy.  The
contribution of a local degree of freedom is added to the value of the
corresponding global degree of freedom. Here we have to use the
function $j_S$ defined on each element $S$ in \mathref{book:E:global-lokal} 
on page \pageref{book:E:global-lokal}. It combines uniquely the
local DOFs with the global ones. The basis functions provide in the
\code{BAS\_FCTS} structure the entry \code{get\_dof\_indices()}
which is an implementation of $j_S$, see \secref{S:basfct_data}.

The updating of the system matrix is not that easy. As mentioned
in \secref{book:S:FE-dis-2nd}, the system matrix is usually sparse and we use
special data structures for storing these matrices, compare
\secref{S:DOF_MATRIX}. For sparse matrices we do not have
for each DOF a matrix row storing values for all other DOFs; only the
values for pairs of DOFs are stored, where the corresponding
\emph{global} basis functions have a common support.  Usually, the
exact number of entries in one row of a sparse matrix is not know a
priori and can change during grid modifications.

Thus, we use the following concept: A call of
\code{clear\_dof\_matrix()} will not set all matrix entries to zero,
but will remove all matrix rows from the matrix, compare the
description of this function on page \pageref{Func:clear_dof_matrix}.
During the updating of a matrix for the value corresponding to a pair
of local DOFs $(i,j)$, we look in the $j_S(i)$th row of the
matrix for a column $j_S(j)$ (the \code{col} member of 
\code{matrix\_row}); if such an entry exists, we add
the current contribution; if this entry does not yet exist we will
create a new entry, set the current value and column number. This
creation may include an enlargement of the row, by linking a new
matrix row to the list of matrix rows, if no space for a new entry is
left. After the assemblage we then have a sparse matrix, storing all
values for pairs of global basis functions with common support.

The functions which we describe now allows also to handle matrices
where the DOFs indexing the rows can differ from the DOFs indexing the
columns; this makes the combination of DOFs from different finite
element spaces possible.

\begin{compatibility}
  \label{compat:DOWB_MATRIX}
  Previous versions of \ALBERTA defined extra-types for vector-valued
  problems, like \code{DOF\_DOWB\_MATRIX}, \code{DOWB\_OPERATOR\_INFO}
  etc. The ``\code{DOWB}'' (``DimOfWorldBlocks'') variants, however,
  already incorporated all the functionality of the ordinary
  scalar-only versions. Therefore the scalar-ony versions of most
  data-structures have been abandoned and were replaced by the
  ``\code{DOWB}'' variants, which in turn were renamed to use the
  scalar-only names. For example, in the current implementation a
  \code{DOF\_MATRIX} is in fact what older versions called a
  \code{DOF\_DOWB\_MATRIX}; and implements the scalar-only case as
  well as the block-matrix case.
\end{compatibility}

\subsubsection{Element matrix and vector structures}
\label{S:elvecmat}

\paragraph{Block-matrix types}
\idx{assemblage tools!MATENT_TYPE@{\code{MATENT\_TYPE}}}
\ddx{MATENT_TYPE@{\code{MATENT\_TYPE}}}
\bv\begin{lstlisting}[label=T:MATENT_TYPE,name=MATENT_TYPE,caption={[data-type: \code{MATENT\_TYPE}]}]
typedef enum matent_type {
  MATENT_NONE    = -1,
  MATENT_REAL    =  0,
  MATENT_REAL_D  =  1,
  MATENT_REAL_DD =  2
} MATENT_TYPE;
\end{lstlisting}\ev

\noindent
Description: This enumeration type defines symbolic types for
block-matrix entries. \code{MATENT\_REAL} means scalar blocks,
\code{MATENT\_REAL\_D} stands for diagonal blocks, and
\code{MATENT\_REAL\_DD} is a code for full matrix blocks. In general,
data-structures make use of these types to store the matrix blocks in
an efficient way.

\paragraph{Structure for element matrices}
\idx{assemblage tools!EL_MATRIX@{\code{EL\_MATRIX}}}
\ddx{EL_MATRIX@{\code{EL\_MATRIX}}}
\bv\begin{lstlisting}[label=T:EL_MATRIX,name=EL_MATRIX,caption={[data-type: \code{EL\_MATRIX}]}]
typedef struct el_matrix EL_MATRIX;

struct el_matrix 
{
  MATENT_TYPE type;
  int n_row, n_col;
  int n_row_max, n_col_max;
  union {
    REAL    *const*real;
    REAL_D  *const*real_d;
    REAL_DD *const*real_dd;
  } data;
  DBL_LIST_NODE row_chain;
  DBL_LIST_NODE col_chain;
};
\end{lstlisting}\ev

\noindent
Description: A data structure to store per-element contributions
during the assembling of discrete systems. There is some limited
support for the operation of element-matrices on element-vectors and
global DOF-vectors, see \secref{S:elementblas}.
\begin{descr}
  \kitem{type} One out of \code{MATENT\_REAL}, \code{MATENT\_REAL\_D}
  or \code{MATENT\_REAL\_DD}. The entries stored in
  \code{EL\_MATRIX->data} have to be interpreted accordingly. See
  \code{MATENT\_TYPE} on page \pageref{T:MATENT_TYPE}.
%%
 \kitem{n\_row} is the number of rows of the element matrix
%%
 \kitem{n\_col} is the number of columns of the element matrix
%%
 \kitem{n\_row\_max} is the maximal number of rows. The number of rows
 can vary from element to element if the underlying basis functions
 have a per-element initializer.
%%
 \kitem{n\_col\_max} is the maximal number of columns.
%
 \kitem{data, data.real, data.real\_d, data.real\_dd}
 \code{EL\_MATRIX->data} is a union, its components should be accessed
 according to the symmetry type indicated by \code{EL\_MATRIX->type}.
 %%
 \kitem{row\_chain, col\_chain} If the underlying finite element
 spaces are a direct sum of function spaces, then the resulting
 element matrices have a block-layout. The link to the other parts of
 the resulting block-matrix is implemented using cyclic doubly linked
 lists, \code{row\_chain} and \code{col\_chain} are the corresponding
 list-nodes. There is a separate section explaining how to handle such
 chains of objects, see \secref{S:chain_impl}.
\end{descr}

\paragraph{Structures for element vectors}
\ddx{EL_INT_VEC@{\code{EL\_INT\_VEC}}}
\idx{assemblage tools!EL_INT_VEC@{\code{EL\_INT\_VEC}}}
\ddx{EL_DOF_VEC@{\code{EL\_DOF\_VEC}}}
\idx{assemblage tools!EL_DOF_VEC@{\code{EL\_DOF\_VEC}}}
\ddx{EL_UCHAR_VEC@{\code{EL\_UCHAR\_VEC}}}
\idx{assemblage tools!EL_UCHAR_VEC@{\code{EL\_UCHAR\_VEC}}}
\ddx{EL_SCHAR_VEC@{\code{EL\_SCHAR\_VEC}}}
\idx{assemblage tools!EL_SCHAR_VEC@{\code{EL\_SCHAR\_VEC}}}
\ddx{EL_BNDRY_VEC@{\code{EL\_BNDRY\_VEC}}}
\idx{assemblage tools!EL_BNDRY_VEC@{\code{EL\_BNDRY\_VEC}}}
\ddx{EL_PTR_VEC@{\code{EL\_PTR\_VEC}}}
\idx{assemblage tools!EL_PTR_VEC@{\code{EL\_PTR\_VEC}}}
\ddx{EL_REAL_VEC@{\code{EL\_REAL\_VEC}}}
\idx{assemblage tools!EL_REAL_VEC@{\code{EL\_REAL\_VEC}}}
\ddx{EL_REAL_VEC_D@{\code{EL\_REAL\_VEC\_D}}}
\idx{assemblage tools!EL_REAL_VEC_D@{\code{EL\_REAL\_VEC\_D}}}
\ddx{EL_REAL_D_VEC@{\code{EL\_REAL\_D\_VEC}}}
\idx{assemblage tools!EL_REAL_D_VEC@{\code{EL\_REAL\_D\_VEC}}}
\bv\begin{lstlisting}[label=T:EL_INT_VEC,name=EL_XXX_VEC,caption={[data-types: \code{EL\_XXX\_VEC}]}] 
typedef struct el_int_vec       EL_INT_VEC;
typedef struct el_dof_vec       EL_DOF_VEC;
typedef struct el_uchar_vec     EL_UCHAR_VEC;
typedef struct el_schar_vec     EL_SCHAR_VEC;
typedef struct el_bndry_vec     EL_BNDRY_VEC;
typedef struct el_ptr_vec       EL_PTR_VEC;
typedef struct el_real_vec      EL_REAL_VEC;
typedef struct el_real_vec_d    EL_REAL_VEC_D;
typedef struct el_real_d_vec    EL_REAL_D_VEC;
\end{lstlisting}\ev
\label{T:EL_DOF_VEC}
\label{T:EL_UCHAR_VEC}
\label{T:EL_SCHAR_VEC}
\label{T:EL_BNDRY_VEC}
\label{T:EL_PTR_VEC}
\label{T:EL_REAL_D_VEC}
The \code{el\_*\_vec} structures are declared similary, the only
difference between them is the type of the structure entry vec. Below,
the \code{EL\_REAL\_VEC} structure is given:
%%
\bv\begin{lstlisting}[label=T:EL_REAL_VEC,name=EL_REAL_VEC,caption={data-type: \code{EL\_REAL\_VEC}}] 
struct el_real_vec
{
  int           n_components;
  int           n_components_max;
  DBL_LIST_NODE chain;
  int           reserved;
  REAL          vec[1]; /* different type in EL_INT_VEC, ... */
}; 
\end{lstlisting}\ev
%%
and the \code{EL\_REAL\_VEC\_D} structure is described in detail:
%%
\bv\begin{lstlisting}[label=T:EL_REAL_VEC_D,name=EL_REAL_VEC_D]
struct el_real_vec_d
{
  int           n_components;
  int           n_components_max;
  DBL_LIST_NODE chain;
  int           stride; /* either 1 or DIM_OF_WORLD */
  REAL          vec[1];
}; 
\end{lstlisting}\ev
%%
Description:
\begin{descr}
  \kitem{n\_components} The actual number of components available in
  and following \code{EL\_XXX\_VEC->vec}. Note that the actual number
  of components is -- of course -- larger than $1$ in general
  \code{get\_el\_XXX\_vec(bas\_fcts)} takes care of allocating enough
  space.
%
  \kitem{n\_components\_max} Behing \code{EL\_XXX\_VEC->vec[0]} may
  actually be more space available than the number of currently valid
  entries as indicated by \code{EL\_XXX\_VEC->n\_components}; this is
  the maximum size to access without crossing the bounds of the data
  segment allocated for this element vector.
%
  \kitem{chain} If the underlying basis-function implementation is
  part of a chain of sets of basis functions, then this structure is
  inherited also by the element vectors: they are chained using a
  doubly linked list, \code{chain} is the corresponding list-node.
  There is a separate section about such chained objects, see
  \secref{S:chain_impl}.
%
 \kitem{stride, reserved} For element vectors other than an
 \code{EL\_REAL\_VEC\_D} this is a reserved value and actually tied to
 the constant value $1$ with the exception of a
 \code{EL\_REAL\_D\_VEC} were \code{reserved} is fixed at
 \code{DIM\_OF\_WORLD}. For \code{EL\_REAL\_VEC\_D} this varies,
 based on the dimension of the range of the underlying basis function
 implementation. For vector-valued basis functions
 \code{EL\_REAL\_VEC\_D->stride} is again tied to $1$, for
 scalar-valued basis functions \code{EL\_REAL\_VEC\_D->stride} is
 fixed at \code{DIM\_OF\_WORLD}, in both cases it gives the number of
 \code{REAL}'s belonging to a single \code{DOF}. See also
 \code{DOF\_REAL\_VEC\_D} on page \pageref{T:DOF_REAL_VEC_D}.
%
 \kitem{vec[1]} Start of the data-segment,
 \code{EL\_XXX\_VEC->n\_components} items contain valid data,
 \code{EL\_XXX\_VEC->n\_components\_max} items are allocated. Note
 that for a \code{EL\_REAL\_VEC\_D} vector the numbers have to be
 multiplied by \code{EL\_REAL\_VEC\_D->stride} to get the actual
 number of \code{REAL}'s allocated.
\end{descr}

\subsubsection{Accumulating per-element contributions}
The following functions can be used on elements for updating matrices
and vectors. 
\fdx{add_element_matrix()@{\code{add\_element\_matrix()}}}
\idx{assemblage tools!add_element_matrix()@{\code{add\_element\_matrix()}}}
\fdx{add_element_vec()@{\code{add\_element\_vec()}}}
\idx{assemblage tools!add_element_vec()@{\code{add\_element\_vec()}}}
\fdx{add_element_d_vec()@{\code{add\_element\_d\_vec()}}}
\idx{assemblage tools!add_element_d_vec()@{\code{add\_element\_d\_vec()}}}
\fdx{add_element_vec_dow()@{\code{add\_element\_vec\_dow()}}}
\idx{assemblage tools!add_element_vec_dow()@{\code{add\_element\_vec\_dow()}}}
\bv\begin{lstlisting}[name={proto-type add_element_matrix()},label=C:add_element_matrix]
void add_element_matrix(DOF_MATRIX *matrix, REAL factor, 
                        const EL_MATRIX *el_matrix, MatrixTranspose transpose,
                        const EL_DOF_VEC *row_dof, const EL_DOF_VEC *col_dof, 
                        const EL_SCHAR_VEC *bound);
void add_element_vec(DOF_REAL_VEC *drv, REAL factor, const EL_REAL_VEC *el_vec,
                     const EL_DOF_VEC *dof, const EL_SCHAR_VEC *bound);
void add_element_d_vec(DOF_REAL_D_VEC *drdv, REAL factor, 
                       const EL_REAL_D_VEC *el_vec, const EL_DOF_VEC *dof,
                       const EL_SCHAR_VEC *bound);
void add_element_vec_dow(DOF_REAL_VEC_D *drdv, REAL factor, 
                         const EL_REAL_VEC_D *el_vec, const EL_DOF_VEC *dof,
                         const EL_SCHAR_VEC *bound);
\end{lstlisting}\ev
\paragraph{Descriptions}
\begin{descr}
\kitem{add\_element\_matrix(mat,
\fdx{add_element_matrix()@{\code{add\_element\_matrix()}}}
\idx{assemblage tools!add_element_matrix()@{\code{add\_element\_matrix()}}}
  factor, el\_mat, transpose, row\_dof, col\_dof, bound)}\hfill

Updates the global \code{DOF\_MATRIX} \code{mat} by adding element
contributions. If \code{row\_dof} equals \code{col\_dof}, the diagonal
element is \emph{always} the first entry in a matrix row; this makes
the access to the diagonal element easy for a diagonal preconditioner,
for example. In general, \code{add\_element\_matrix()} does the
following: for all \code{i} the values
\code{fac*el\_mat->data.\{REAL,REAL\_D,REAL\_DD\}[i][j]} are added to
the entries at the position (\code{row\_dof->vec[i],col\_dof->vec[j]})
in the global matrix \code{mat}
($\code{0}\leq\code{i}<\code{el\_mat->n\_row}$,
$\code{0}\leq\code{j}<\code{el\_mat->n\_col}$). If such an entry
exists in the row number \code{row\_dof->vec[i]} the global matrix
\code{mat} the value is simply added. Otherwise a new entry is created
in the row, the value is set and the column number is set to
\code{col\_dof[j]}. This may imply an enlargement of the row by adding
a new \code{MATRIX\_ROW} structure to the list of matrix rows.

Note that the first element matrix added to \code{mat} after calling
\code{clear\_dof\_matrix()} determines the block-type of the global
matrix \code{mat}. It is possible to add element-matrices with higher
block-symmetry to global \code{DOF\_MATRIX}es with lower
block-symmetry, for example it is allowed to add \code{el\_mat} to
\code{mat} if \code{el\_mat->type == MATENT\_REAL} and \code{mat->type
  == MATENT\_REAL\_DD}.

\begin{description}
\item[Parameters]\hfill
  \begin{description}
  \item[\code{mat}] the global \code{DOF\_MATRIX}.
  \item[\code{factor}] is a multiplier for the element contributions;
    usually \code{factor} is \code{1} or \code{-1};
    
  \item[\code{el\_mat}] is a matrix of size
    $\code{n\_row}\times\code{n\_col}$ storing the element
    contributions;
    
  \item[\code{transpose}] the original matrix is used if
    \code{transpose} == \code{NoTranspose} (= 0) and the transposed
    matrix if \code{transpose} == \code{Transpose} (= 1);
    
  \item[\code{row\_dof}] is a vector of length
    \code{row\_dof->n\_components} storing the global row indices;

  \item[\code{col\_dof}] is a vector of length
    \code{col\_dof->n\_components} storing the global column indices,
    \code{col\_dof} may be a \nil pointer if the DOFs indexing the
    columns are the same as the DOFs indexing the rows; in this case
    \code{col\_dof = row\_dof} is used;
    
  \item[\code{bound}] is either \code{NULL} or an
    \code{EL\_SCHAR\_Vec} stucture storing a vector of length
    \code{bound->n\_components}.  In this case
    \code{bound->n\_components} must match either
    \code{row\_dof->n\_components} or \code{col\_dof->n\_components},
    depending on the value of \code{transpose}.

    If \code{bound->vec[i] >= DIRICHLET}, then the following happens:
    \begin{description}
    \item[\code{row\_dof == col\_dof}] In the global \code{matrix} the
      row \code{row\_dof->vec[i]} is cleared to zero, with the
      exception of the diagonal entry, which is set to \code{1.0}.
    \item[\code{row\_dof != col\_dof}] In the global \code{matrix} the
      row \code{row\_dof->vec[i]} is cleared to zero.
    \end{description}

    All other contributions of \code{el\_mat} are added to
    \code{matrix} as usual. This allows for a convenient way to
    implement inhomogeneous Dirichlet boundary conditions, without
    having to modify the right-hand-side of the discrete systems
    explicitly.
  \end{description}
\end{description}

 \hrulefill
\fdx{add_element_vec()@{\code{add\_element\_vec()}}}
\idx{assemblage tools!add_element_vec()@{\code{add\_element\_vec()}}}
\fdx{add_element_d_vec()@{\code{add\_element\_d\_vec()}}}
\idx{assemblage tools!add_element_d_vec()@{\code{add\_element\_d\_vec()}}}
\fdx{add_element_vec_dow()@{\code{add\_element\_vec\_dow()}}}
\idx{assemblage tools!add_element_vec_dow()@{\code{add\_element\_vec\_dow()}}}
 \kitem{add\_element\_vec(drv, factor, el\_vec, dof, bound)}
 \kitem{add\_element\_d\_vec(drv, factor, el\_vec, dof, bound)}
 \kitem{add\_element\_vec\_d(drv, factor, el\_vec, dof, bound)}\hfill

 These do similar things as \code{add\_element\_matrix()}, but with
 element vectors. \secref{S:elementblas} also lists other routines
 which might be helpful in this context.
\end{descr}

\subsubsection{Allocation and filling of element vectors}
\label{S:fillgetelvec}

\paragraph{Prototypes}
\fdx{get_dof_indices()@{\code{get\_dof\_indices()}}}
\fdx{get_bound()@{\code{get\_bound()}}}
\fdx{el_interpol()@{\code{el\_interpol()}}}
\fdx{el_interpol_dow()@{\code{el\_interpol\_dow()}}}
\fdx{dirichlet_map()@{\code{dirichlet\_map()}}}
%%
%\fdx{init_el_int_vec()@{\code{init\_el\_int\_vec()}}}
%\fdx{init_el_dof_vec()@{\code{init\_el\_dof\_vec()}}}
%\fdx{init_el_real_vec()@{\code{init\_el\_real\_vec()}}}
%\fdx{init_el_real_d_vec()@{\code{init\_el\_real\_d\_vec()}}}
%\fdx{init_el_uchar_vec()@{\code{init\_el\_uchar\_vec()}}}
%\fdx{init_el_schar_vec()@{\code{init\_el\_schar\_vec()}}}
%\fdx{init_el_ptr_vec()@{\code{init\_el\_ptr\_vec()}}}
%%
\fdx{fill_el_int_vec()@{\code{fill\_el\_int\_vec()}}}
\fdx{fill_el_real_vec()@{\code{fill\_el\_real\_vec()}}}
\fdx{fill_el_real_d_vec()@{\code{fill\_el\_real\_d\_vec()}}}
\fdx{fill_el_real_vec_d()@{\code{fill\_el\_real\_vec\_d()}}}
\fdx{fill_el_uchar_vec()@{\code{fill\_el\_uchar\_vec()}}}
\fdx{fill_el_schar_vec()@{\code{fill\_el\_schar\_vec()}}}
%%
\fdx{get_el_int_vec()@{\code{get\_el\_int\_vec()}}}
\fdx{get_el_dof_vec()@{\code{get\_el\_dof\_vec()}}}
\fdx{get_el_uchar_vec()@{\code{get\_el\_uchar\_vec()}}}
\fdx{get_el_schar_vec()@{\code{get\_el\_schar\_vec()}}}
\fdx{get_el_bndry_vec()@{\code{get\_el\_bndry\_vec()}}}
\fdx{get_el_ptr_vec()@{\code{get\_el\_ptr\_vec()}}}
\fdx{get_el_real_vec()@{\code{get\_el\_real\_vec()}}}
\fdx{get_el_real_d_vec()@{\code{get\_el\_real\_d\_vec()}}}
\fdx{get_el_real_vec_d()@{\code{get\_el\_real\_vec\_d()}}}
%%
\fdx{free_el_int_vec()@{\code{free\_el\_int\_vec()}}}
\fdx{free_el_dof_vec()@{\code{free\_el\_dof\_vec()}}}
\fdx{free_el_uchar_vec()@{\code{free\_el\_uchar\_vec()}}}
\fdx{free_el_schar_vec()@{\code{free\_el\_schar\_vec()}}}
\fdx{free_el_bndry_vec()@{\code{free\_el\_bndry\_vec()}}}
\fdx{free_el_ptr_vec()@{\code{free\_el\_ptr\_vec()}}}
\fdx{free_el_real_vec()@{\code{free\_el\_real\_vec()}}}
\fdx{free_el_real_d_vec()@{\code{free\_el\_real\_d\_vec()}}}
\fdx{free_el_real_vec_d()@{\code{free\_el\_real\_vec\_d()}}}
\bv\begin{lstlisting}[label=P:ELVEC_PROTOS,name=ELVEC_PROTOS,caption={[proto-types: filling element vectors]}]
EL_DOF_VEC *get_dof_indices(EL_DOF_VEC *dofs, const FE_SPACE *fe_space,
                            const EL *el);
EL_BNDRY_VEC *get_bound(EL_BNDRY_VEC *bndry, const BAS_FCTS *bas_fcts,
                        const EL_INFO *el_info);
void el_interpol(EL_REAL_VEC *coeff, const EL_INFO *el_info, int wall,
                 const EL_INT_VEC *indices, LOC_FCT_AT_QP f, void *ud,
                 const BAS_FCTS *bas_fcts);
void el_interpol_dow(EL_REAL_VEC_D *coeff, const EL_INFO *el_info, int wall,
                     const EL_INT_VEC *indices, LOC_FCT_D_AT_QP f,
                     void *f_data, const BAS_FCTS *bas_fcts);
void dirichlet_map(EL_SCHAR_VEC *bound, const EL_BNDRY_VEC *bndry_bits,
                   const BNDRY_FLAGS mask);

const EL_INT_VEC *
fill_el_int_vec(EL_INT_VEC *el_vec, EL *el, const DOF_INT_VEC *dof_vec);
const EL_REAL_VEC *
fill_el_real_vec(EL_REAL_VEC *el_vec, EL *el, const DOF_REAL_VEC *dof_vec);
const EL_REAL_D_VEC *
fill_el_real_d_vec(EL_REAL_D_VEC *el_vec, EL *el, const DOF_REAL_D_VEC *dof_vec);
const EL_REAL_VEC_D *
fill_el_real_vec_d(EL_REAL_VEC_D *el_vec, EL *el, const DOF_REAL_VEC_D *dof_vec);
const EL_UCHAR_VEC *
fill_el_uchar_vec(EL_UCHAR_VEC *el_vec, EL *el, const DOF_UCHAR_VEC *dof_vec);
const EL_SCHAR_VEC *
fill_el_schar_vec(EL_SCHAR_VEC *el_vec, EL *el, const DOF_SCHAR_VEC *dof_vec);

EL_INT_VEC *get_el_int_vec(const BAS_FCTS *bas_fcts);
EL_DOF_VEC *get_el_dof_vec(const BAS_FCTS *bas_fcts);
EL_UCHAR_VEC *get_el_uchar_vec(const BAS_FCTS *bas_fcts);
EL_SCHAR_VEC *get_el_schar_vec(const BAS_FCTS *bas_fcts);
EL_BNDRY_VEC *get_el_bndry_vec(const BAS_FCTS *bas_fcts);
EL_PTR_VEC *get_el_ptr_vec(const BAS_FCTS *bas_fcts);
EL_REAL_VEC *get_el_real_vec(const BAS_FCTS *bas_fcts);
EL_REAL_D_VEC *get_el_real_d_vec(const BAS_FCTS *bas_fcts);
EL_REAL_VEC_D *get_el_real_vec_d(const BAS_FCTS *bas_fcts);

void free_el_int_vec(EL_INT_VEC *el_vec);
void free_el_dof_vec(EL_DOF_VEC *el_vec);
void free_el_uchar_vec(EL_UCHAR_VEC *el_vec);
void free_el_schar_vec(EL_SCHAR_VEC *el_vec);
void free_el_bndry_vec(EL_BNDRY_VEC *el_vec);
void free_el_ptr_vec(EL_PTR_VEC *el_vec);
void free_el_real_vec(EL_REAL_VEC *el_vec);
void free_el_real_d_vec(EL_REAL_D_VEC *el_vec);
void free_el_real_vec_d(EL_REAL_VEC_D *el_vec);

DEF_EL_VEC_VAR(VECNAME, name, _size, _size_max, _init);
DEF_EL_VEC_CONST(VECNAME, name, _size, _size_max);
ALLOC_EL_VEC(VECNAME, _size, _size_max);
\end{lstlisting}\ev

\paragraph{Descriptions}
\begin{descr}
  \kitem{get\_dof\_indices(dofs, fe\_space, el)}
  
  Compute the mapping
  \fdx{get_dof_indices()@{\code{get\_dof\_indices()}}}
  between the local \code{DOF}-indices on \code{el} and the global
  \code{DOF}-indices according to \code{fe\_space->admin}.
 \begin{description}
 \item[Parameters]~
   \begin{description}
   \item[\code{dofs}] Storage for the result or \code{NULL}. In the
     latter case the mapping is returned in a statically allocated
     \code{EL\_DOF\_VEC}. \emph{Note: this storage area will be
       overwritten on the next call to this function, even if the
       \code{fe\_space} argument differs.}
   \item[\code{fe\_space}] The finite element space to compute the mapping for.
   \item[\code{el}] The current mesh element (\emph{not} the current
     \code{EL\_INFO} pointer, use \code{EL\_INFO->el}).
   \end{description}
 \item[return] Either again the argument \code{dofs} or -- if
   \code{dofs == NULL} -- a pointer to a statically allocated
   \code{EL\_DOF\_VEC}.
 \item[examples]
\begin{samepage}
With pre-allocated \code{EL\_DOF\_VEC}:
\bv\begin{lstlisting}[caption={[Example for \code{get\_dof\_indices()}]},label=C:get_dof_indices_example]
EL_DOF_VEC *dofs = get_el_dof_vec(fe_space->bas_fcts);
TRAVERSE_FIRST(mesh, -1, CALL_LEAF_EL) {
  int i;

  get_dof_indices(dofs, fe_space, el_info->el);
  for (i = 0; i < bas_fcts->n_bas_fcts; i++) {
    MSG("dofs[%d] = %d\n", dofs->vec[i]);
  }
} TRAVERSE_NEXT();
free_el_dof_vec(dofs);
\end{lstlisting}\ev
\end{samepage}

\begin{samepage}
Without pre-allocated \code{EL\_DOF\_VEC}:
\bv\begin{lstlisting}[caption={[Example for \code{get\_dof\_indices()}]},label=C:get_dof_indices_example2]
TRAVERSE_FIRST(mesh, -1, CALL_LEAF_EL) {
  int i;
  EL_DOF_VEC *dofs = get_dof_indices(NULL, fe_space, el_info->el);

  for (i = 0; i < bas_fcts->n_bas_fcts; i++) {
    MSG("dofs[%d] = %d\n", dofs->vec[i]);
  }
} TRAVERSE_NEXT();
\end{lstlisting}\ev
\end{samepage}
%%\item[see also]
 \end{description}
 %%
 \hrulefill
 %%
 \kitem{get\_bound(bndry, bas\_fcts, el\_info)}
 \fdx{get_bound()@{\code{get\_bound()}}} Extract the boundary types of
 the local DOFs of \code{bas\_fcts}.  The boundary types are returned
 in form of a bit-mask. If bit \code{j} in the bit-mask
 \code{bndry[i]} is set, then the local DOF number \code{i} belongs to
 the boundary segment which has been assigned the number \code{j} in
 the macro-triangulation. Boundary types range from $1$ to $255$.
 \begin{description}
 \item[Parameters]~
   \begin{description}
   \item[\code{EL\_BNDRY\_VEC *bndry}] Storage for the result or
     \code{NULL}. In the latter case the data is returned in a
     statically allocated \code{EL\_BNDRY\_VEC}.
   \item[\code{BAS\_FCTS *bas\_fcts}] The local basis functions.
   \item[\code{const EL\_INFO *el\_info}] The current mesh element
     info structure.  (\emph{not} the current \code{EL\_INFO} pointer.
   \end{description}
 \item[return] Either again the argument \code{bndry} or -- if
   \code{bndry == NULL} -- a pointer to a statically allocated
   \code{EL\_BNDRY\_VEC}.
 \item[examples]
\begin{samepage}
With pre-allocated \code{EL\_BNDRY\_VEC}:
   \bv\begin{lstlisting}
EL_BNDRY_VEC *bndry = get_el_bndry_vec(bas_fcts);
TRAVERSE_FIRST(mesh, -1, CALL_LEAF_EL|FILL_BOUND) {
  int i, j;
  get_bound(bndry, bas_fcts, el_info);
  for (i = 0; i < bas_fcts->n_bas_fcts; i++) {
    for (j = 1; j < N_BNDRY_TYPES; j++) {
      if (BNNDRY_FLAGS_IS_INTERIOR(bndry->vec[i])) {
        MSG("Local dof %d is an interior DOF\n");
      } else if (BNDRY_FLAGS_IS_AT_BNDRY(bndry->vec[i], j)) {
        MSG("Local dof %d belongs to boundary segment %d\n", i, j);
      }
    }
  }
} TRAVERSE_NEXT();
free_el_bndry_vec(bndry);
   \end{lstlisting}\ev
\end{samepage}

\begin{samepage}
Without pre-allocated \code{EL\_BNDRY\_VEC}:
   \bv\begin{lstlisting}
TRAVERSE_FIRST(mesh, -1, CALL_LEAF_EL) {
  int i,j;
  EL_BNDRY_VEC *bndry = get_bound(NULL, bas_fcts, el_info);

  for (i = 0; i < bas_fcts->n_bas_fcts; i++) {
    for (j = 1; j < N_BNDRY_TYPES; j++) {
      if (BNNDRY_FLAGS_IS_INTERIOR(bndry->vec[i])) {
        MSG("Local dof %d is an interior DOF\n");
      } else if (BNDRY_FLAGS_IS_AT_BNDRY(bndry->vec[i], j)) {
        MSG("Local dof %d belongs to boundary segment %d\n", i, j);
      }
    }
  }
} TRAVERSE_NEXT();
   \end{lstlisting}\ev
\end{samepage}
 \end{description}
 %%
 \hrulefill
 %%
 \kitem{fill\_el\_int\_vec(el\_vec, el, dof\_vec)}
 \kitem{fill\_el\_real\_vec(el\_vec, el, dof\_vec)}
 \kitem{fill\_el\_real\_d\_vec(el\_vec, el, dof\_vec)}
 \kitem{fill\_el\_real\_vec\_d(el\_vec, el, dof\_vec)}
 \kitem{fill\_el\_uchar\_vec(el\_vec, el, dof\_vec)}
 \kitem{fill\_el\_schar\_vec(el\_vec, el, dof\_vec)}
 \fdx{fill_el_int_vec()@{\code{fill\_el\_int\_vec()}}}
 \fdx{fill_el_real_vec()@{\code{fill\_el\_real\_vec()}}}
 \fdx{fill_el_real_d_vec()@{\code{fill\_el\_real\_d\_vec()}}}
 \fdx{fill_el_real_vec_d()@{\code{fill\_el\_real\_vec\_d()}}}
 \fdx{fill_el_uchar_vec()@{\code{fill\_el\_uchar\_vec()}}}
 \fdx{fill_el_schar_vec()@{\code{fill\_el\_schar\_vec()}}}~\hfill
 %%

 Fill the respective element vector with data. The description below
 is for \code{fill\_el\_real\_vec()}, the other versions work similar.
 \begin{description}
 \item[Parameters]\hfill
   \begin{description}
   \item[\code{EL\_REAL\_VEC *el\_vec}] Storage for the result or
     \code{NULL}. In the latter case the return value is
     \hyperlink{DOF_REAL_VEC:vec_loc}{\code{DOF\_REAL\_VEC->vec\_loc}};
     the data will be overwritten on the next call to
     \code{fill\_el\_real\_vec()} with the same \code{dof\_vec}
     argument. Calling \code{fill\_el\_real\_vec()} with \emph{other}
     DOF-vectors will \emph{not} invalidate the data.
   \item[\code{EL *el}] The current mesh element (\emph{not} the
     current \code{EL\_INFO} pointer, use \code{EL\_INFO->el}).
   \item[\code{DOF\_REAL\_VEC *dof\_vec}] The global \code{DOF}-vector
     to extract the data from.
   \end{description}
 \item[return] Either again a the pointer \code{el\_vec} or --
   if \code{el\_vec == NULL} a pointer to a statically allocated
   result space which will be overwritten on the next call to
   \code{fill\_el\_real\_vec()}. \emph{Warning:} see ``bugs'' below.
 \end{description}
 %%
 \hrulefill
 %% 
 \kitem{get\_el\_int\_vec(bas\_fcts)}
 \kitem{get\_el\_dof\_vec(bas\_fcts)}
 \kitem{get\_el\_uchar\_vec(bas\_fcts)}
 \kitem{get\_el\_schar\_vec(bas\_fcts)}
 \kitem{get\_el\_bndry\_vec(bas\_fcts)}
 \kitem{get\_el\_ptr\_vec(bas\_fcts)}
 \kitem{get\_el\_real\_vec(bas\_fcts)}
 \kitem{get\_el\_real\_d\_vec(bas\_fcts)}
 \kitem{get\_el\_real\_vec\_d(bas\_fcts)}\hfill
\fdx{get_el_int_vec()@{\code{get\_el\_int\_vec()}}}
\fdx{get_el_dof_vec()@{\code{get\_el\_dof\_vec()}}}
\fdx{get_el_uchar_vec()@{\code{get\_el\_uchar\_vec()}}}
\fdx{get_el_schar_vec()@{\code{get\_el\_schar\_vec()}}}
\fdx{get_el_bndry_vec()@{\code{get\_el\_bndry\_vec()}}}
\fdx{get_el_ptr_vec()@{\code{get\_el\_ptr\_vec()}}}
\fdx{get_el_real_vec()@{\code{get\_el\_real\_vec()}}}
\fdx{get_el_real_d_vec()@{\code{get\_el\_real\_d\_vec()}}}
\fdx{get_el_real_vec_d()@{\code{get\_el\_real\_vec\_d()}}}

The \code{get\_el\_*\_vec()} routines automatically allocates enough 
memory for the element data vector \code{vec} as indicated by
\code{bas\_fcts->n\_bas\_fcts}.
\begin{description}
\item[Parameters]
  \begin{description}
  \item[\code{const BAS\_FCTS *bas\_fcts}] 
  \end{description}
\item[return] A pointer to a dynamically allocated element
  vector of the respective type.
\item[examples] See the first example for the
  \code{fill\_el\_real\_vec()} function.
\end{description}
 %%
 \hrulefill
 %% 
 \kitem{free\_el\_int\_vec(el\_vec)}
 \kitem{free\_el\_dof\_vec(el\_vec)}
 \kitem{free\_el\_uchar\_vec(el\_vec)}
 \kitem{free\_el\_schar\_vec(el\_vec)}
 \kitem{free\_el\_bndry\_vec(el\_vec)}
 \kitem{free\_el\_ptr\_vec(el\_fcts)}
 \kitem{free\_el\_real\_vec(bas\_fcts)}
 \kitem{free\_el\_real\_d\_vec(bas\_fcts)}
 \kitem{free\_el\_real\_vec\_d(bas\_fcts)}
\fdx{free_el_int_vec()@{\code{free\_el\_int\_vec()}}}
\fdx{free_el_dof_vec()@{\code{free\_el\_dof\_vec()}}}
\fdx{free_el_uchar_vec()@{\code{free\_el\_uchar\_vec()}}}
\fdx{free_el_schar_vec()@{\code{free\_el\_schar\_vec()}}}
\fdx{free_el_bndry_vec()@{\code{free\_el\_bndry\_vec()}}}
\fdx{free_el_ptr_vec()@{\code{free\_el\_ptr\_vec()}}}
\fdx{free_el_real_vec()@{\code{free\_el\_real\_vec()}}}
\fdx{free_el_real_d_vec()@{\code{free\_el\_real\_d\_vec()}}}
\fdx{free_el_real_vec_d()@{\code{free\_el\_real\_vec\_d()}}}
\hfill

The \code{free\_el\_XXX\_vec()} routines free all previously allocated 
storage for \code{el\_XXX\_vec} data.
\begin{description}
\item[Parameters]
  \begin{description}
  \item[\code{const BAS\_FCTS *bas\_fcts}] 
  \end{description}
\item[return] \code{void}
\item[examples] See the first example for the
  \code{fill\_el\_real\_vec()} function.
\end{description}
 %%
 \hrulefill
 %%
\kitem{DEF\_EL\_VEC\_VAR(VECNAME, name, size, size\_max, init)}
\fdx{DEF_EL_VEC_VAR()@{\code{DEF\_EL\_VEC\_VAR()}}}
\mdx{DEF_EL_VEC_VAR()@{\code{DEF\_EL\_VEC\_VAR()}}}\hfill

This is a macro which defines a (local) variable with id \code{name},
pointing to an \code{EL\_VECNAME\_VEC} of size \code{size}, holding a
maximal number of elements \code{max\_size}, which is initialised if
\code{init} is \code{true}. \code{size} and \code{size\_max} may be
variables.

 \hrulefill
 %%
\kitem{DEF\_EL\_VEC\_CONST(VECNAME, name, size, size\_max)}\hfill
\fdx{DEF_EL_VEC_CONST()@{\code{DEF\_EL\_VEC\_CONST()}}}
\mdx{DEF_EL_VEC_CONST()@{\code{DEF\_EL\_VEC\_CONST()}}}

This is a macro which defines a (local) variable with id \code{name},
pointing to an \code{EL\_VECNAME\_VEC} of size \code{size}, holding a
maximal number of elements \code{max\_size}. \code{size} and
\code{size\_max} must be constant values.

 \hrulefill
 %%
\kitem{ALLOC\_EL\_VEC(VECNAME, size, size\_max)}
\fdx{ALLOC_EL_VEC()@{\code{ALLOC\_EL\_VEC()}}}
\mdx{ALLOC_EL_VEC()@{\code{ALLOC\_EL\_VEC()}}}

This macro allocates a \code{EL\_VECNAME\_VEC} with enough storage to
hold \code{size\_max} elements; the \code{n\_components} component of
the element vector structure is set to \code{size}.

\hrulefill
 %%
 \kitem{el\_interpol(coeff, el\_info, wall, indices, f, ud, bas\_fcts)}
 \fdx{el_interpol()@{\code{el\_interpol()}}}
 %%
 \kitem{el\_interpol\_dow(coeff, el\_info, wall, indices, f, f\_data, ud, bas\_fcts)}
 \fdx{el_interpol_dow()@{\code{el\_interpol\_dow()}}}
 %%
 \kitem{dirichlet\_map(bound, bndry\_bits, mask)}
 \fdx{dirichlet_map()@{\code{dirichlet\_map()}}}\hfill

\end{descr}

\subsubsection{BLAS-like Element-matrix and -vector operations}
\label{S:elementblas}

The source code listing below lists the proto-types, refer to
\tableref{tab:elementblas1} and \tableref{tab:elementblas2} for a
description of the respective operations. The routines in
\tableref{tab:elementblas2} take an argument
%%
\bv
\begin{lstlisting}
const EL_SCHAR_VEC *bound.
\end{lstlisting}\ev
%%
In this case the operations will act only on the rows $r$ which are
not masked-out by \code{bound->vec[r] >= DIRICHLET}. The \code{bound}
argument maybe \nil.
%%
\fdx{el_bi_mat_vec()@{\code{el\_bi\_mat\_vec()}}}
\idx{BLAS for element vectors and matrices!el_bi_mat_vec()@{\code{el\_bi\_mat\_vec()}}}
\fdx{el_bi_mat_vec_d()@{\code{el\_bi\_mat\_vec\_d()}}}
\idx{BLAS for element vectors and matrices!el_bi_mat_vec_d()@{\code{el\_bi\_mat\_vec\_d()}}}
\fdx{el_bi_mat_vec_dow()@{\code{el\_bi\_mat\_vec\_dow()}}}
\idx{BLAS for element vectors and matrices!el_bi_mat_vec_dow()@{\code{el\_bi\_mat\_vec\_dow()}}}
\fdx{el_bi_mat_vec_rrd()@{\code{el\_bi\_mat\_vec\_rrd()}}}
\idx{BLAS for element vectors and matrices!el_bi_mat_vec_rrd()@{\code{el\_bi\_mat\_vec\_rrd()}}}
\fdx{el_bi_mat_vec_scl_dow()@{\code{el\_bi\_mat\_vec\_scl\_dow()}}}
\idx{BLAS for element vectors and matrices!el_bi_mat_vec_scl_dow()@{\code{el\_bi\_mat\_vec\_scl\_dow()}}}
\fdx{el_bi_mat_vec_rdr()@{\code{el\_bi\_mat\_vec\_rdr()}}}
\idx{BLAS for element vectors and matrices!el_bi_mat_vec_rdr()@{\code{el\_bi\_mat\_vec\_rdr()}}}
\fdx{el_bi_mat_vec_dow_scl()@{\code{el\_bi\_mat\_vec\_dow\_scl()}}}
\idx{BLAS for element vectors and matrices!el_bi_mat_vec_dow_scl()@{\code{el\_bi\_mat\_vec\_dow\_scl()}}}
%%
\fdx{el_gen_mat_vec()@{\code{el\_gen\_mat\_vec()}}}
\idx{BLAS for element vectors and matrices!el_gen_mat_vec()@{\code{el\_gen\_mat\_vec()}}}
\fdx{el_gen_mat_vec_d()@{\code{el\_gen\_mat\_vec\_d()}}}
\idx{BLAS for element vectors and matrices!el_gen_mat_vec_d()@{\code{el\_gen\_mat\_vec\_d()}}}
\fdx{el_gen_mat_vec_dow()@{\code{el\_gen\_mat\_vec\_dow()}}}
\idx{BLAS for element vectors and matrices!el_gen_mat_vec_dow()@{\code{el\_gen\_mat\_vec\_dow()}}}
\fdx{el_gen_mat_vec_rrd()@{\code{el\_gen\_mat\_vec\_rrd()}}}
\idx{BLAS for element vectors and matrices!el_gen_mat_vec_rrd()@{\code{el\_gen\_mat\_vec\_rrd()}}}
\fdx{el_gen_mat_vec_scl_dow()@{\code{el\_gen\_mat\_vec\_scl\_dow()}}}
\idx{BLAS for element vectors and matrices!el_gen_mat_vec_scl_dow()@{\code{el\_gen\_mat\_vec\_scl\_dow()}}}
\fdx{el_gen_mat_vec_rdr()@{\code{el\_gen\_mat\_vec\_rdr()}}}
\idx{BLAS for element vectors and matrices!el_gen_mat_vec_rdr()@{\code{el\_gen\_mat\_vec\_rdr()}}}
\fdx{el_gen_mat_vec_dow_scl()@{\code{el\_gen\_mat\_vec\_dow\_scl()}}}
\idx{BLAS for element vectors and matrices!el_gen_mat_vec_dow_scl()@{\code{el\_gen\_mat\_vec\_dow\_scl()}}}
%%
\fdx{el_mat_vec()@{\code{el\_mat\_vec()}}}
\idx{BLAS for element vectors and matrices!el_mat_vec()@{\code{el\_mat\_vec()}}}
\fdx{el_mat_vec_d()@{\code{el\_mat\_vec\_d()}}}
\idx{BLAS for element vectors and matrices!el_mat_vec_d()@{\code{el\_mat\_vec\_d()}}}
\fdx{el_mat_vec_dow()@{\code{el\_mat\_vec\_dow()}}}
\idx{BLAS for element vectors and matrices!el_mat_vec_dow()@{\code{el\_mat\_vec\_dow()}}}
\fdx{el_mat_vec_rrd()@{\code{el\_mat\_vec\_rrd()}}}
\idx{BLAS for element vectors and matrices!el_mat_vec_rrd()@{\code{el\_mat\_vec\_rrd()}}}
\fdx{el_mat_vec_scl_dow()@{\code{el\_mat\_vec\_scl\_dow()}}}
\idx{BLAS for element vectors and matrices!el_mat_vec_scl_dow()@{\code{el\_mat\_vec\_scl\_dow()}}}
\fdx{el_mat_vec_rdr()@{\code{el\_mat\_vec\_rdr()}}}
\idx{BLAS for element vectors and matrices!el_mat_vec_rdr()@{\code{el\_mat\_vec\_rdr()}}}
\fdx{el_mat_vec_dow_scl()@{\code{el\_mat\_vec\_dow\_scl()}}}
\idx{BLAS for element vectors and matrices!el_mat_vec_dow_scl()@{\code{el\_mat\_vec\_dow\_scl()}}}
%%
\fdx{bi_mat_el_vec()@{\code{bi\_mat\_el\_vec()}}}
\idx{BLAS for element vectors and matrices!bi_mat_el_vec()@{\code{bi\_mat\_el\_vec()}}}
\fdx{bi_mat_el_vec_d()@{\code{bi\_mat\_el\_vec\_d()}}}
\idx{BLAS for element vectors and matrices!bi_mat_el_vec_d()@{\code{bi\_mat\_el\_vec\_d()}}}
\fdx{bi_mat_el_vec_dow()@{\code{bi\_mat\_el\_vec\_dow()}}}
\idx{BLAS for element vectors and matrices!bi_mat_el_vec_dow()@{\code{bi\_mat\_el\_vec\_dow()}}}
\fdx{bi_mat_el_vec_rrd()@{\code{bi\_mat\_el\_vec\_rrd()}}}
\idx{BLAS for element vectors and matrices!bi_mat_el_vec_rrd()@{\code{bi\_mat\_el\_vec\_rrd()}}}
\fdx{bi_mat_el_vec_scl_dow()@{\code{bi\_mat\_el\_vec\_scl\_dow()}}}
\idx{BLAS for element vectors and matrices!bi_mat_el_vec_scl_dow()@{\code{bi\_mat\_el\_vec\_scl\_dow()}}}
\fdx{bi_mat_el_vec_rdr()@{\code{bi\_mat\_el\_vec\_rdr()}}}
\idx{BLAS for element vectors and matrices!bi_mat_el_vec_rdr()@{\code{bi\_mat\_el\_vec\_rdr()}}}
\fdx{bi_mat_el_vec_dow_scl()@{\code{bi\_mat\_el\_vec\_dow\_scl()}}}
\idx{BLAS for element vectors and matrices!bi_mat_el_vec_dow_scl()@{\code{bi\_mat\_el\_vec\_dow\_scl()}}}
%%
\fdx{gen_mat_el_vec()@{\code{gen\_mat\_el\_vec()}}}
\idx{BLAS for element vectors and matrices!gen_mat_el_vec()@{\code{gen\_mat\_el\_vec()}}}
\fdx{gen_mat_el_vec_d()@{\code{gen\_mat\_el\_vec\_d()}}}
\idx{BLAS for element vectors and matrices!gen_mat_el_vec_d()@{\code{gen\_mat\_el\_vec\_d()}}}
\fdx{gen_mat_el_vec_dow()@{\code{gen\_mat\_el\_vec\_dow()}}}
\idx{BLAS for element vectors and matrices!gen_mat_el_vec_dow()@{\code{gen\_mat\_el\_vec\_dow()}}}
\fdx{gen_mat_el_vec_rrd()@{\code{gen\_mat\_el\_vec\_rrd()}}}
\idx{BLAS for element vectors and matrices!gen_mat_el_vec_rrd()@{\code{gen\_mat\_el\_vec\_rrd()}}}
\fdx{gen_mat_el_vec_scl_dow()@{\code{gen\_mat\_el\_vec\_scl\_dow()}}}
\idx{BLAS for element vectors and matrices!gen_mat_el_vec_scl_dow()@{\code{gen\_mat\_el\_vec\_scl\_dow()}}}
\fdx{gen_mat_el_vec_rdr()@{\code{gen\_mat\_el\_vec\_rdr()}}}
\idx{BLAS for element vectors and matrices!gen_mat_el_vec_rdr()@{\code{gen\_mat\_el\_vec\_rdr()}}}
\fdx{gen_mat_el_vec_dow_scl()@{\code{gen\_mat\_el\_vec\_dow\_scl()}}}
\idx{BLAS for element vectors and matrices!gen_mat_el_vec_dow_scl()@{\code{gen\_mat\_el\_vec\_dow\_scl()}}}
%%
\fdx{mat_el_vec()@{\code{mat\_el\_vec()}}}
\idx{BLAS for element vectors and matrices!mat_el_vec()@{\code{mat\_el\_vec()}}}
\fdx{mat_el_vec_d()@{\code{mat\_el\_vec\_d()}}}
\idx{BLAS for element vectors and matrices!mat_el_vec_d()@{\code{mat\_el\_vec\_d()}}}
\fdx{mat_el_vec_dow()@{\code{mat\_el\_vec\_dow()}}}
\idx{BLAS for element vectors and matrices!mat_el_vec_dow()@{\code{mat\_el\_vec\_dow()}}}
\fdx{mat_el_vec_rrd()@{\code{mat\_el\_vec\_rrd()}}}
\idx{BLAS for element vectors and matrices!mat_el_vec_rrd()@{\code{mat\_el\_vec\_rrd()}}}
\fdx{mat_el_vec_scl_dow()@{\code{mat\_el\_vec\_scl\_dow()}}}
\idx{BLAS for element vectors and matrices!mat_el_vec_scl_dow()@{\code{mat\_el\_vec\_scl\_dow()}}}
\fdx{mat_el_vec_rdr()@{\code{mat\_el\_vec\_rdr()}}}
\idx{BLAS for element vectors and matrices!mat_el_vec_rdr()@{\code{mat\_el\_vec\_rdr()}}}
\fdx{mat_el_vec_dow_scl()@{\code{mat\_el\_vec\_dow\_scl()}}}
\idx{BLAS for element vectors and matrices!mat_el_vec_dow_scl()@{\code{mat\_el\_vec\_dow\_scl()}}}
%%
\fdx{el_mat_set()@{\code{el\_mat\_set()}}}
\idx{BLAS for element vectors and matrices!el_mat_set()@{\code{el\_mat\_set()}}}
\fdx{el_mat_axey()@{\code{el\_mat\_axey()}}}
\idx{BLAS for element vectors and matrices!el_mat_axey()@{\code{el\_mat\_axey()}}}
\fdx{el_mat_axpy()@{\code{el\_mat\_axpy()}}}
\idx{BLAS for element vectors and matrices!el_mat_axpy()@{\code{el\_mat\_axpy()}}}
\fdx{el_mat_axpby()@{\code{el\_mat\_axpy()}}}
\idx{BLAS for element vectors and matrices!el_mat_axpby()@{\code{el\_mat\_axpby()}}}
\bv\begin{lstlisting}
EL_REAL_VEC *el_bi_mat_vec(REAL a, const EL_MATRIX *A,
                           REAL b, const EL_MATRIX *B,
                           const EL_REAL_VEC *u_h,
                           REAL c, EL_REAL_VEC *f_h);
EL_REAL_D_VEC *el_bi_mat_vec_d(REAL a, const EL_MATRIX *A,
                               REAL b, const EL_MATRIX *B,
                               const EL_REAL_D_VEC *u_h,
                               REAL c, EL_REAL_D_VEC *f_h);
EL_REAL_VEC_D *el_bi_mat_vec_dow(REAL a, const EL_MATRIX *A,
                                 REAL b, const EL_MATRIX *B,
                                 const EL_REAL_VEC_D *u_h,
                                 REAL c, EL_REAL_VEC_D *f_h);
EL_REAL_VEC *el_bi_mat_vec_rrd(REAL a, const EL_MATRIX *A,
                               REAL b, const EL_MATRIX *B,
                               const EL_REAL_D_VEC *u_h,
                               REAL c, EL_REAL_VEC *f_h);
EL_REAL_VEC *el_bi_mat_vec_scl_dow(REAL a, const EL_MATRIX *A,
                                   REAL b, const EL_MATRIX *B,
                                   const EL_REAL_VEC_D *u_h,
                                   REAL c, EL_REAL_VEC *f_h);
EL_REAL_D_VEC *el_bi_mat_vec_rdr(REAL a, const EL_MATRIX *A,
                                 REAL b, const EL_MATRIX *B,
                                 const EL_REAL_VEC *u_h,
                                 REAL c, EL_REAL_D_VEC *f_h);
EL_REAL_VEC_D *el_bi_mat_vec_dow_scl(REAL a, const EL_MATRIX *A,
                                     REAL b, const EL_MATRIX *B,
                                     const EL_REAL_VEC *u_h,
                                     REAL c, EL_REAL_VEC_D *f_h);

EL_REAL_VEC *el_gen_mat_vec(REAL a, const EL_MATRIX *A,
                            const EL_REAL_VEC *u_h,
                             REAL b, EL_REAL_VEC *f_h);
EL_REAL_D_VEC *el_gen_mat_vec_d(REAL a, const EL_MATRIX *A,
                                const EL_REAL_D_VEC *u_h,
                                REAL b, EL_REAL_D_VEC *f_h);
EL_REAL_VEC_D *el_gen_mat_vec_dow(REAL a, const EL_MATRIX *A,
                                  const EL_REAL_VEC_D *u_h,
                                  REAL b, EL_REAL_VEC_D *f_h);
EL_REAL_VEC *el_gen_mat_vec_rrd(REAL a, const EL_MATRIX *A,
                                const EL_REAL_D_VEC *u_h,
                                REAL b, EL_REAL_VEC *f_h);
EL_REAL_VEC *el_gen_mat_vec_scl_dow(REAL a, const EL_MATRIX *A,
                                    const EL_REAL_VEC_D *u_h,
                                    REAL b, EL_REAL_VEC *f_h);
EL_REAL_D_VEC *el_gen_mat_vec_rdr(REAL a, const EL_MATRIX *A,
                                  const EL_REAL_VEC *u_h,
                                  REAL b, EL_REAL_D_VEC *f_h);
EL_REAL_VEC_D *el_gen_mat_vec_dow_scl(REAL a, const EL_MATRIX *A,
                                      const EL_REAL_VEC *u_h,
                                      REAL b, EL_REAL_VEC_D *f_h);

EL_REAL_VEC *el_mat_vec(REAL a, const EL_MATRIX *A,
                        const EL_REAL_VEC *u_h, EL_REAL_VEC *f_h);
EL_REAL_D_VEC *el_mat_vec_d(REAL a, const EL_MATRIX *A,
                            const EL_REAL_D_VEC *u_h, EL_REAL_D_VEC *f_h);
EL_REAL_VEC_D *el_mat_vec_dow(REAL a, const EL_MATRIX *A,
                              const EL_REAL_VEC_D *u_h, EL_REAL_VEC_D *f_h);
EL_REAL_VEC *el_mat_vec_rrd(REAL a, const EL_MATRIX *A,
                            const EL_REAL_D_VEC *u_h, EL_REAL_VEC *f_h);
EL_REAL_VEC *el_mat_vec_scl_dow(REAL a, const EL_MATRIX *A,
                                const EL_REAL_VEC_D *u_h, EL_REAL_VEC *f_h);
EL_REAL_D_VEC *el_mat_vec_rdr(REAL a, const EL_MATRIX *A,
                              const EL_REAL_VEC *u_h, EL_REAL_D_VEC *f_h);
EL_REAL_VEC_D *el_mat_vec_dow_scl(REAL a, const EL_MATRIX *A,
                                  const EL_REAL_VEC *u_h, EL_REAL_VEC_D *f_h);

void bi_mat_el_vec(REAL a, const EL_MATRIX *A,
                   REAL b, const EL_MATRIX *B,
                   const EL_REAL_VEC *u_h, REAL c, DOF_REAL_VEC *f_h,
                   const EL_DOF_VEC *dof, const EL_SCHAR_VEC *bound);
void bi_mat_el_vec_d(REAL a, const EL_MATRIX *A,
                     REAL b, const EL_MATRIX *B,
                     const EL_REAL_D_VEC *u_h, REAL c, DOF_REAL_D_VEC *f_h,
                     const EL_DOF_VEC *dof, const EL_SCHAR_VEC *bound);
void bi_mat_el_vec_dow(REAL a, const EL_MATRIX *A,
                       REAL b, const EL_MATRIX *B,
                       const EL_REAL_VEC_D *u_h, REAL c, DOF_REAL_VEC_D *f_h,
                       const EL_DOF_VEC *dof, const EL_SCHAR_VEC *bound);
void bi_mat_el_vec_rrd(REAL a, const EL_MATRIX *A,
                       REAL b, const EL_MATRIX *B,
                       const EL_REAL_D_VEC *u_h, REAL c, DOF_REAL_VEC *f_h,
                       const EL_DOF_VEC *dof, const EL_SCHAR_VEC *bound);
void bi_mat_el_vec_scl_dow(REAL a, const EL_MATRIX *A,
                           REAL b, const EL_MATRIX *B,
                           const EL_REAL_VEC_D *u_h, REAL c, DOF_REAL_VEC *f_h,
                           const EL_DOF_VEC *dof, const EL_SCHAR_VEC *bound);
void bi_mat_el_vec_rdr(REAL a, const EL_MATRIX *A,
                       REAL b, const EL_MATRIX *B,
                       const EL_REAL_VEC *u_h, REAL c, DOF_REAL_D_VEC *f_h,
                       const EL_DOF_VEC *dof, const EL_SCHAR_VEC *bound);
void bi_mat_el_vec_dow_scl(REAL a, const EL_MATRIX *A,
                           REAL b, const EL_MATRIX *B,
                           const EL_REAL_VEC *u_h, REAL c, DOF_REAL_VEC_D *f_h,
                           const EL_DOF_VEC *dof, const EL_SCHAR_VEC *bound);

void gen_mat_el_vec(REAL a, const EL_MATRIX *A,
                    const EL_REAL_VEC *u_h, REAL b, DOF_REAL_VEC *f_h,
                    const EL_DOF_VEC *dof, const EL_SCHAR_VEC *bound);
void gen_mat_el_vec_d(REAL a, const EL_MATRIX *A,
                      const EL_REAL_D_VEC *u_h, REAL b, DOF_REAL_D_VEC *f_h,
                      const EL_DOF_VEC *dof, const EL_SCHAR_VEC *bound);
void gen_mat_el_vec_dow(REAL a, const EL_MATRIX *A,
                        const EL_REAL_VEC_D *u_h, REAL b, DOF_REAL_VEC_D *f_h,
                        const EL_DOF_VEC *dof, const EL_SCHAR_VEC *bound);
void gen_mat_el_vec_rrd(REAL a, const EL_MATRIX *A,
                        const EL_REAL_D_VEC *u_h, REAL b, DOF_REAL_VEC *f_h,
                        const EL_DOF_VEC *dof, const EL_SCHAR_VEC *bound);
void gen_mat_el_vec_scl_dow(REAL a, const EL_MATRIX *A,
                            const EL_REAL_VEC_D *u_h, REAL b, DOF_REAL_VEC *f_h,
                            const EL_DOF_VEC *dof, const EL_SCHAR_VEC *bound);
void gen_mat_el_vec_rdr(REAL a, const EL_MATRIX *A,
                        const EL_REAL_VEC *u_h, REAL b, DOF_REAL_D_VEC *f_h,
                        const EL_DOF_VEC *dof, const EL_SCHAR_VEC *bound);
void gen_mat_el_vec_dow_scl(REAL a, const EL_MATRIX *A,
                            const EL_REAL_VEC *u_h, REAL b, DOF_REAL_VEC_D *f_h,
                            const EL_DOF_VEC *dof, const EL_SCHAR_VEC *bound);

void mat_el_vec(REAL a, const EL_MATRIX *A,
                const EL_REAL_VEC *u_h, DOF_REAL_VEC *f_h,
                const EL_DOF_VEC *dof, const EL_SCHAR_VEC *bound);
void mat_el_vec_d(REAL a, const EL_MATRIX *A,
                  const EL_REAL_D_VEC *u_h, DOF_REAL_D_VEC *f_h,
                  const EL_DOF_VEC *dof, const EL_SCHAR_VEC *bound);
void mat_el_vec_dow(REAL a, const EL_MATRIX *A,
                    const EL_REAL_VEC_D *u_h, DOF_REAL_VEC_D *f_h,
                    const EL_DOF_VEC *dof, const EL_SCHAR_VEC *bound);
void mat_el_vec_rrd(REAL a, const EL_MATRIX *A,
                    const EL_REAL_D_VEC *u_h, DOF_REAL_VEC *f_h,
                    const EL_DOF_VEC *dof, const EL_SCHAR_VEC *bound);
void mat_el_vec_scl_dow(REAL a, const EL_MATRIX *A,
                        const EL_REAL_VEC_D *u_h, DOF_REAL_VEC *f_h,
                        const EL_DOF_VEC *dof, const EL_SCHAR_VEC *bound);
void mat_el_vec_rdr(REAL a, const EL_MATRIX *A,
                    const EL_REAL_VEC *u_h, DOF_REAL_D_VEC *f_h,
                    const EL_DOF_VEC *dof, const EL_SCHAR_VEC *bound);
void mat_el_vec_dow_scl(REAL a, const EL_MATRIX *A,
                        const EL_REAL_VEC *u_h, DOF_REAL_VEC_D *f_h,
                        const EL_DOF_VEC *dof, const EL_SCHAR_VEC *bound);

EL_MATRIX *el_mat_set(REAL a, EL_MATRIX *result);
EL_MATRIX *el_mat_axey(REAL a, const EL_MATRIX *A, EL_MATRIX *result);
EL_MATRIX *el_mat_axpy(REAL a, const EL_MATRIX *A, EL_MATRIX *result);
EL_MATRIX *el_mat_axpby(REAL a, const EL_MATRIX *A,
                        REAL b, const EL_MATRIX *B, EL_MATRIX *result);
\end{lstlisting}\ev

\begin{table}[htbp]
\idx{BLAS for element vectors and matrices}
\fdx{el_bi_mat_vec()@{\code{el\_bi\_mat\_vec()}}}
\idx{BLAS for element vectors and matrices!el_bi_mat_vec()@{\code{el\_bi\_mat\_vec()}}}
\fdx{el_bi_mat_vec_d()@{\code{el\_bi\_mat\_vec\_d()}}}
\idx{BLAS for element vectors and matrices!el_bi_mat_vec_d()@{\code{el\_bi\_mat\_vec\_d()}}}
\fdx{el_bi_mat_vec_dow()@{\code{el\_bi\_mat\_vec\_dow()}}}
\idx{BLAS for element vectors and matrices!el_bi_mat_vec_dow()@{\code{el\_bi\_mat\_vec\_dow()}}}
\fdx{el_bi_mat_vec_rrd()@{\code{el\_bi\_mat\_vec\_rrd()}}}
\idx{BLAS for element vectors and matrices!el_bi_mat_vec_rrd()@{\code{el\_bi\_mat\_vec\_rrd()}}}
\fdx{el_bi_mat_vec_scl_dow()@{\code{el\_bi\_mat\_vec\_scl\_dow()}}}
\idx{BLAS for element vectors and matrices!el_bi_mat_vec_scl_dow()@{\code{el\_bi\_mat\_vec\_scl\_dow()}}}
\fdx{el_bi_mat_vec_rdr()@{\code{el\_bi\_mat\_vec\_rdr()}}}
\idx{BLAS for element vectors and matrices!el_bi_mat_vec_rdr()@{\code{el\_bi\_mat\_vec\_rdr()}}}
\fdx{el_bi_mat_vec_dow_scl()@{\code{el\_bi\_mat\_vec\_dow\_scl()}}}
\idx{BLAS for element vectors and matrices!el_bi_mat_vec_dow_scl()@{\code{el\_bi\_mat\_vec\_dow\_scl()}}}
%%
\fdx{el_gen_mat_vec()@{\code{el\_gen\_mat\_vec()}}}
\idx{BLAS for element vectors and matrices!el_gen_mat_vec()@{\code{el\_gen\_mat\_vec()}}}
\fdx{el_gen_mat_vec_d()@{\code{el\_gen\_mat\_vec\_d()}}}
\idx{BLAS for element vectors and matrices!el_gen_mat_vec_d()@{\code{el\_gen\_mat\_vec\_d()}}}
\fdx{el_gen_mat_vec_dow()@{\code{el\_gen\_mat\_vec\_dow()}}}
\idx{BLAS for element vectors and matrices!el_gen_mat_vec_dow()@{\code{el\_gen\_mat\_vec\_dow()}}}
\fdx{el_gen_mat_vec_rrd()@{\code{el\_gen\_mat\_vec\_rrd()}}}
\idx{BLAS for element vectors and matrices!el_gen_mat_vec_rrd()@{\code{el\_gen\_mat\_vec\_rrd()}}}
\fdx{el_gen_mat_vec_scl_dow()@{\code{el\_gen\_mat\_vec\_scl\_dow()}}}
\idx{BLAS for element vectors and matrices!el_gen_mat_vec_scl_dow()@{\code{el\_gen\_mat\_vec\_scl\_dow()}}}
\fdx{el_gen_mat_vec_rdr()@{\code{el\_gen\_mat\_vec\_rdr()}}}
\idx{BLAS for element vectors and matrices!el_gen_mat_vec_rdr()@{\code{el\_gen\_mat\_vec\_rdr()}}}
\fdx{el_gen_mat_vec_dow_scl()@{\code{el\_gen\_mat\_vec\_dow\_scl()}}}
\idx{BLAS for element vectors and matrices!el_gen_mat_vec_dow_scl()@{\code{el\_gen\_mat\_vec\_dow\_scl()}}}
%%
\fdx{el_mat_vec()@{\code{el\_mat\_vec()}}}
\idx{BLAS for element vectors and matrices!el_mat_vec()@{\code{el\_mat\_vec()}}}
\fdx{el_mat_vec_d()@{\code{el\_mat\_vec\_d()}}}
\idx{BLAS for element vectors and matrices!el_mat_vec_d()@{\code{el\_mat\_vec\_d()}}}
\fdx{el_mat_vec_dow()@{\code{el\_mat\_vec\_dow()}}}
\idx{BLAS for element vectors and matrices!el_mat_vec_dow()@{\code{el\_mat\_vec\_dow()}}}
\fdx{el_mat_vec_rrd()@{\code{el\_mat\_vec\_rrd()}}}
\idx{BLAS for element vectors and matrices!el_mat_vec_rrd()@{\code{el\_mat\_vec\_rrd()}}}
\fdx{el_mat_vec_scl_dow()@{\code{el\_mat\_vec\_scl\_dow()}}}
\idx{BLAS for element vectors and matrices!el_mat_vec_scl_dow()@{\code{el\_mat\_vec\_scl\_dow()}}}
\fdx{el_mat_vec_rdr()@{\code{el\_mat\_vec\_rdr()}}}
\idx{BLAS for element vectors and matrices!el_mat_vec_rdr()@{\code{el\_mat\_vec\_rdr()}}}
\fdx{el_mat_vec_dow_scl()@{\code{el\_mat\_vec\_dow\_scl()}}}
\idx{BLAS for element vectors and matrices!el_mat_vec_dow_scl()@{\code{el\_mat\_vec\_dow\_scl()}}}
%%
\fdx{el_mat_set()@{\code{el\_mat\_set()}}}
\idx{BLAS for element vectors and matrices!el_mat_set()@{\code{el\_mat\_set()}}}
\fdx{el_mat_axey()@{\code{el\_mat\_axey()}}}
\idx{BLAS for element vectors and matrices!el_mat_axey()@{\code{el\_mat\_axey()}}}
\fdx{el_mat_axpy()@{\code{el\_mat\_axpy()}}}
\idx{BLAS for element vectors and matrices!el_mat_axpy()@{\code{el\_mat\_axpy()}}}
\fdx{el_mat_axpby()@{\code{el\_mat\_axpby()}}}
\idx{BLAS for element vectors and matrices!el_mat_axpby()@{\code{el\_mat\_axpby()}}}
\begin{center}{\small
    \begin{tabular}{|l|l|}
      \hline
      \Strut\verb|f = el_mat_vec(a, A, u, f)| &
      $f_i \leftarrow (a\,A\,u)_i$ \\
      \Strut\verb|f = el_mat_vec_d(a, A, u, f)| & \\
      \Strut\verb|f = el_mat_vec_dow(a, A, u, f)| & \\
      \Strut\verb|f = el_mat_vec_rrd(a, A, u, f)| & \\
      \Strut\verb|f = el_mat_vec_scl_dow(a, A, u, f)| &\\
      \Strut\verb|f = el_mat_vec_rdr(a, A, u, f)| &\\
      \Strut\verb|f = el_mat_vec_dow_scl(a, A, u, f)|&\\
      \hline
      \Strut\verb|f = el_gen_mat_vec(a, A, u, b, f)| & 
      $f_i \leftarrow (a\,A\,u + b\,f)_i$ \\
      \Strut\verb|f = el_gen_mat_vec_d(a, A, u, b, f)| & \\
      \Strut\verb|f = el_gen_mat_vec_dow(a, A, u, b, f)| & \\
      \Strut\verb|f = el_gen_mat_vec_rrd(a, A, u, b, f)| & \\
      \Strut\verb|f = el_gen_mat_vec_scl_dow(a, A, u, b, f)| &\\
      \Strut\verb|f = el_gen_mat_vec_rdr(a, A, u, b, f)| &\\
      \Strut\verb|f = el_gen_mat_vec_dow_scl(a, A, u, b, f)|&\\
      \hline
      \Strut\verb|f = el_bi_mat_vec(a, A, b, B, u, c, f)| & 
      $f_i \leftarrow ((a\,A + b\,B)\,u + c\,f)_i$ \\
      \Strut\verb|f = el_bi_mat_vec_d(a, A, b, B, u, c, f)| & \\
      \Strut\verb|f = el_bi_mat_vec_dow(a, A, b, B, u, c, f)| & \\
      \Strut\verb|f = el_bi_mat_vec_rrd(a, A, b, B, u, c, f)| & \\
      \Strut\verb|f = el_bi_mat_vec_scl_dow(a, A, b, B, u, c, f)| &\\
      \Strut\verb|f = el_bi_mat_vec_rdr(a, A, b, B, u, c, f)| &\\
      \Strut\verb|f = el_bi_mat_vec_dow_scl(a, A, b, B, u, c, f)|&\\
      \hline
      \Strut\verb|A = el_mat_set(a, A)| &
      $A_{ij} \leftarrow a$\\
      \Strut\verb|B = el_mat_axey(a, A, B)| &
      $B_{ij} \leftarrow a\,A_{ij}$\\
      \Strut\verb|B = el_mat_axpy(a, A, B)| &
      $B_{ij} \leftarrow a\,A_{ij}+B_{ij}$\\
      \Strut\verb|C = el_mat_axpby(a, A, b, B, C)| &
      $C_{ij} \leftarrow a\,A_{ij}+b\,B_{ij}$\\
      \hline
    \end{tabular}}\end{center}
\caption[BLAS-operations for element-vectors and -matrices]
{BLAS-operations for element-vectors and -matrices. $A$ and $B$ denote element matrices, $u$ and $f$ element vectors, $a$, $b$, $c$ are numbers.}
\label{tab:elementblas1}
\end{table}

\begin{table}[htbp]
\idx{BLAS for element vectors and matrices}
\fdx{bi_mat_el_vec()@{\code{bi\_mat\_el\_vec()}}}
\idx{BLAS for element vectors and matrices!bi_mat_el_vec()@{\code{bi\_mat\_el\_vec()}}}
\fdx{bi_mat_el_vec_d()@{\code{bi\_mat\_el\_vec\_d()}}}
\idx{BLAS for element vectors and matrices!bi_mat_el_vec_d()@{\code{bi\_mat\_el\_vec\_d()}}}
\fdx{bi_mat_el_vec_dow()@{\code{bi\_mat\_el\_vec\_dow()}}}
\idx{BLAS for element vectors and matrices!bi_mat_el_vec_dow()@{\code{bi\_mat\_el\_vec\_dow()}}}
\fdx{bi_mat_el_vec_rrd()@{\code{bi\_mat\_el\_vec\_rrd()}}}
\idx{BLAS for element vectors and matrices!bi_mat_el_vec_rrd()@{\code{bi\_mat\_el\_vec\_rrd()}}}
\fdx{bi_mat_el_vec_scl_dow()@{\code{bi\_mat\_el\_vec\_scl\_dow()}}}
\idx{BLAS for element vectors and matrices!bi_mat_el_vec_scl_dow()@{\code{bi\_mat\_el\_vec\_scl\_dow()}}}
\fdx{bi_mat_el_vec_rdr()@{\code{bi\_mat\_el\_vec\_rdr()}}}
\idx{BLAS for element vectors and matrices!bi_mat_el_vec_rdr()@{\code{bi\_mat\_el\_vec\_rdr()}}}
\fdx{bi_mat_el_vec_dow_scl()@{\code{bi\_mat\_el\_vec\_dow\_scl()}}}
\idx{BLAS for element vectors and matrices!bi_mat_el_vec_dow_scl()@{\code{bi\_mat\_el\_vec\_dow\_scl()}}}
%%
\fdx{gen_mat_el_vec()@{\code{gen\_mat\_el\_vec()}}}
\idx{BLAS for element vectors and matrices!gen_mat_el_vec()@{\code{gen\_mat\_el\_vec()}}}
\fdx{gen_mat_el_vec_d()@{\code{gen\_mat\_el\_vec\_d()}}}
\idx{BLAS for element vectors and matrices!gen_mat_el_vec_d()@{\code{gen\_mat\_el\_vec\_d()}}}
\fdx{gen_mat_el_vec_dow()@{\code{gen\_mat\_el\_vec\_dow()}}}
\idx{BLAS for element vectors and matrices!gen_mat_el_vec_dow()@{\code{gen\_mat\_el\_vec\_dow()}}}
\fdx{gen_mat_el_vec_rrd()@{\code{gen\_mat\_el\_vec\_rrd()}}}
\idx{BLAS for element vectors and matrices!gen_mat_el_vec_rrd()@{\code{gen\_mat\_el\_vec\_rrd()}}}
\fdx{gen_mat_el_vec_scl_dow()@{\code{gen\_mat\_el\_vec\_scl\_dow()}}}
\idx{BLAS for element vectors and matrices!gen_mat_el_vec_scl_dow()@{\code{gen\_mat\_el\_vec\_scl\_dow()}}}
\fdx{gen_mat_el_vec_rdr()@{\code{gen\_mat\_el\_vec\_rdr()}}}
\idx{BLAS for element vectors and matrices!gen_mat_el_vec_rdr()@{\code{gen\_mat\_el\_vec\_rdr()}}}
\fdx{gen_mat_el_vec_dow_scl()@{\code{gen\_mat\_el\_vec\_dow\_scl()}}}
\idx{BLAS for element vectors and matrices!gen_mat_el_vec_dow_scl()@{\code{gen\_mat\_el\_vec\_dow\_scl()}}}
%%
\fdx{mat_el_vec()@{\code{mat\_el\_vec()}}}
\idx{BLAS for element vectors and matrices!mat_el_vec()@{\code{mat\_el\_vec()}}}
\fdx{mat_el_vec_d()@{\code{mat\_el\_vec\_d()}}}
\idx{BLAS for element vectors and matrices!mat_el_vec_d()@{\code{mat\_el\_vec\_d()}}}
\fdx{mat_el_vec_dow()@{\code{mat\_el\_vec\_dow()}}}
\idx{BLAS for element vectors and matrices!mat_el_vec_dow()@{\code{mat\_el\_vec\_dow()}}}
\fdx{mat_el_vec_rrd()@{\code{mat\_el\_vec\_rrd()}}}
\idx{BLAS for element vectors and matrices!mat_el_vec_rrd()@{\code{mat\_el\_vec\_rrd()}}}
\fdx{mat_el_vec_scl_dow()@{\code{mat\_el\_vec\_scl\_dow()}}}
\idx{BLAS for element vectors and matrices!mat_el_vec_scl_dow()@{\code{mat\_el\_vec\_scl\_dow()}}}
\fdx{mat_el_vec_rdr()@{\code{mat\_el\_vec\_rdr()}}}
\idx{BLAS for element vectors and matrices!mat_el_vec_rdr()@{\code{mat\_el\_vec\_rdr()}}}
\fdx{mat_el_vec_dow_scl()@{\code{mat\_el\_vec\_dow\_scl()}}}
\idx{BLAS for element vectors and matrices!mat_el_vec_dow_scl()@{\code{mat\_el\_vec\_dow\_scl()}}}
\begin{center}{\small
    \begin{tabular}{|l|l|}
      \hline
      \Strut\verb|f = mat_el_vec(a, A, u, f, dof, mask)| &
      $f[dof[i]] \leftarrow (a\,A\,u)_i$ \\
      \Strut\verb|mat_el_vec_d(a, A, u, f, dof, mask)| &
      if \code{mask[i] != DIRICHLET}\\
      \Strut\verb|mat_el_vec_dow(a, A, u, f, dof, mask)| &
      or \code{mask == \nil}\\
      \Strut\verb|mat_el_vec_rrd(a, A, u, f, dof, mask)| & \\
      \Strut\verb|mat_el_vec_scl_dow(a, A, u, f, dof, mask)| &\\
      \Strut\verb|mat_el_vec_rdr(a, A, u, f, dof, mask)| &\\
      \Strut\verb|mat_el_vec_dow_scl(a, A, u, f, dof, mask)|&\\
      \hline
      \Strut\verb|gen_mat_el_vec(a, A, u, b, f, dof, mask)| & 
      $f[dof[i]] \leftarrow (a\,A\,u)_i + b\,f[dof[i]]$ \\
      \Strut\verb|gen_mat_el_vec_d(a, A, u, b, f, dof, mask)| &
      if \code{mask[i] != DIRICHLET}\\
      \Strut\verb|gen_mat_el_vec_dow(a, A, u, b, f, dof, mask)| &
      or \code{mask == \nil}\\
      \Strut\verb|gen_mat_el_vec_rrd(a, A, u, b, f, dof, mask)| & \\
      \Strut\verb|gen_mat_el_vec_scl_dow(a, A, u, b, f, dof, mask)| &\\
      \Strut\verb|gen_mat_el_vec_rdr(a, A, u, b, f, dof, mask)| &\\
      \Strut\verb|gen_mat_el_vec_dow_scl(a, A, u, b, f, dof, mask)|&\\
      \hline
      \Strut\verb|bi_mat_el_vec(a, A, b, B, u, c, f, dof, mask)| & 
      $f[dof[i]]$ \\
      \Strut\verb|bi_mat_el_vec_d(a, A, b, B, u, c, f, dof, mask)| &
      $\leftarrow ((a\,A + b\,B)\,u)_i + c\,f[dof[i]]$ \\
      \Strut\verb|bi_mat_el_vec_dow(a, A, b, B, u, c, f, dof, mask)| &
      if \code{mask[i] != DIRICHLET}\\
      \Strut\verb|bi_mat_el_vec_rrd(a, A, b, B, u, c, f, dof, mask)| &
      or \code{mask == \nil}\\
      \Strut\verb|bi_mat_el_vec_scl_dow(a, A, b, B, u, c, f, dof, mask)| &\\
      \Strut\verb|bi_mat_el_vec_rdr(a, A, b, B, u, c, f, dof, mask)| &\\
      \Strut\verb|bi_mat_el_vec_dow_scl(a, A, b, B, u, c, f, dof, mask)|&\\
      \hline
\end{tabular}}\end{center}
\caption[BLAS-operations for element-vectors and -matrices]
{BLAS-operations for element-vectors and -matrices. $A$ and $B$ denote element matrices, $u$ an element vector. $f$ is a global \code{DOF}-vector. \code{mask} is an \code{EL\_SCHAR\_VEC} masking out certain \emph{local} \code{DOF}s. \code{mask} may be \nil. \code{dof} is \code{EL\_DOF\_VEC} mapping local to global \code{DOF}s. $a$, $b$, $c$ are numbers.}
\label{tab:elementblas2}
\end{table}

\subsection{Data structures and functions for matrix assemblage}%
\label{S:matrix_assemblage}%

The following structure holds full information for the assembling
of scalar element matrices. This structure is used by the function
\code{update\_matrix()} described below.
\ddx{EL_MATRIX_INFO@{\code{EL\_MATRIX\_INFO}}}
\idx{assemblage tools!EL_MATRIX_INFO@{\code{EL\_MATRIX\_INFO}}}
%%
\newcommand{\ELMATRIXINFO}{\hyperref[T:EL_MATRIX_INFO]{\code{EL\_MATRIX\_INFO}}\xspace}
%%
\bv\begin{lstlisting}[label=T:EL_MATRIX_INFO,name=EL_MATRIX_INFO]
typedef struct el_matrix_info  EL_MATRIX_INFO;

struct el_matrix_info
{
  const FE_SPACE *row_fe_space;
  const FE_SPACE *col_fe_space;

  MATENT_TYPE    krn_blk_type;

  BNDRY_FLAGS    dirichlet_bndry;
  REAL           factor;

  EL_MATRIX_FCT  el_matrix_fct;
  void           *fill_info;

  const EL_MATRIX_FCT *neigh_el_mat_fcts;
  void           *neigh_fill_info;

  FLAGS          fill_flag;
};
\end{lstlisting}\ev
Description:
\begin{descr}
  \kitem{row\_fe\_space} pointer to a finite element space connected
  to the row DOFs of the resulting matrix.
  %% 
  \kitem{col\_fe\_space} pointer to a finite element space connected 
  to the columns DOFs of the resulting matrix.
  %% 
  \kitem{krn\_blk\_type} defines the block-matrix type of matrix entries
  %% 
  \kitem{dirichlet\_boundary} bndry-type bit-mask for Dirichlet-boundary 
  conditions built into the matrix
  %% 
  \kitem{factor} is a multiplier for the element contributions; usually
  \code{factor} is \code{1} or \code{-1}.
  %% 
  \kitem{el\_matrix\_fct}is a pointer to a function for the computation
  of the element matrix; \code{el\_matrix\_fct(el\_info, fill\_info)}
  returns a pointer to a matrix of size 
  $\code{n\_row}\times\code{n\_col}$ storing the element matrix
  on element \code{el\_info->el}; \code{fill\_info} is a pointer
  to data needed by \code{el\_matrix\_fct()}; the function has
  to provide memory for storing the element matrix, which
  can be overwritten on the next call.
  %% 
  \kitem{fill\_info} pointer to data needed by \code{el\_matrix\_fct()};
  will be given as second argument to this function.
  %% 
  \kitem{neigh\_el\_mat\_fcts} If the \code{BNDRY\_OPERATOR\_INFO}
  (code-listing \ref{T:BNDRY_OPERATOR_INFO}) structure passed to
  \code{fill\_matrix\_info()} was flagged as discontinuous, then this
  is the base-address of an array storing \code{N\_NEIGH(mesh->dim)}
  many element-matrix functions which pair the local basis-function
  set with the local basis function set on the neighbor.
  Intentionally, this is meant to support assembling linear systems in
  the context of DG-methods. The idea is that
  %%
  \code{EL\_MATRIX\_INFO.neigh\_el\_mat\_fcts[neigh\_nr](el\_info, EL\_MATRIX\_INFO.neigh\_fill\_info)}
  %%
  assembles a jump-term where the local basis functions on the element
  described by \code{el\_info} are used as test-functions
  (corresponding to the rows of the element matrix) and the local
  basis function set on the neighbour element defines the local space
  of ansatz-functions (column-space).
  %% 
  \kitem{neigh\_fill\_info} Data pointer passed to the element-matrix
  functions stored in \code{neigh\_el\_mat\_fcts}.
  %% 
  \kitem{fill\_flag}the flag for the mesh traversal for assembling the
  matrix.
\end{descr}
The following function updates a matrix by assembling element
contributions during mesh traversal; information for computing
the element matrices is provided in an \code{EL\_MATRIX\_INFO} structure:
\fdx{update_matrix()@{\code{update\_matrix()}}}
\idx{assemblage tools!update_matrix()@{\code{update\_matrix(}}}
\bv\begin{lstlisting}[name=update_matrix(),label=T:update_matrix_proto]
void update_matrix(DOF_MATRIX *dof_matrix, const EL_MATRIX_INFO *minfo,
		   MatrixTranspose transpose);;
\end{lstlisting}\ev
Description:
\begin{descr}
  \kitem{update\_matrix(matrix, info, transpose)} updates the matrix
  \code{matrix} by traversing the underlying mesh and assembling the
  element contributions into the matrix; information about the
  computation of element matrices and connection of local and global
  DOFs is stored in \code{info}.
       
  The flags for the mesh traversal of the mesh
  \code{matrix->fe\_space->mesh} are stored at \code{info->fill\_flag}
  which specifies the elements to be visited and information that
  should be present on the elements for the calculation of the element
  matrices and boundary information (if \code{info->get\_bound} is not
  \nil).
       
  On the elements, information about the row DOFs is accessed by
  \code{info->get\_row\_dof} using \code{info->row\_admin}; this
  vector is also used for the column DOFs if \code{info->n\_col} is
  less or equal to zero, or \code{info->get\_col\_admin} or
  \code{info->get\_col\_dof} is a \nil pointer; when row and column
  DOFs are the same, the boundary type of the DOFs is accessed by
  \code{info->get\_bound} if \code{info->get\_bound} is not a \nil
  pointer; then the element matrix is computed by
  \code{info->el\_matrix\_fct(el\_info, info->fill\_info)}; these
  contributions, multiplied by \code{info->factor}, are eventually
  added to \code{matrix} by a call of \code{add\_element\_matrix()}
  with all information about row and column DOFs, the element matrix,
  and boundary types, if available.

  \code{update\_matrix()} acts \emph{additive}, the
  element-contributions are added to the data already present in
  \code{dof\_matrix}. This makes several calls for the assemblage of
  one matrix possible. \code{clear\_dof\_matrix()} can be used to
  erase the contents of \code{dof\_matrix} prior to calling
  \code{update\_matrix()}.
\begin{description}
\item[Parameters]\hfill
  \begin{descr}
    \kitem{dof\_matrix} The global \code{DOF\_MATRIX} to add data to.
    %%
    \kitem{minfo} The element-matrix handle, as returned by
    \code{fill\_matrix\_info()} or 
    %%
    \kitem{transpose}
  \end{descr}
\end{description}
\end{descr}

\subsection{Matrix assemblage for second order problems}%
\label{S:matrix_assemblage_scalar}%

Now we want to describe some tools which enable an easy assemblage of
the system matrix in the case of scalar elliptic problems. For this we
have to provide a function for the calculation of the element
matrix. For a general scalar problem the element matrix
$\vL_S = (L_S^{ij})_{i,j=1,\dots,m}$ is given by (recall
\mathref{book:E:L-phii-phij} on page
\pageref{book:E:L-phii-phij})
%in \secref{book:S:Dis2ndorder}
\begin{align*}
L_S^{ij} &= 
\int_{\Shat}
 \nablal \pbar^{i}(\lambda(\xhat)) \cdot \bar A(\lambda(\xhat))\,
    \nablal \pbar^{j}(\lambda(\xhat))\,d\xhat
+
\int_{\Shat} \pbar^{i}(\lambda(\xhat))\;
  \bar b(\lambda(\xhat)) \cdot 
       \nablal \pbar^{j}(\lambda(\xhat)) \,d\xhat\\
&\qquad +
\int_{\Shat} \bar c(\lambda(\xhat))\, \pbar^{i}(\lambda(\xhat)) \,
                             \pbar^{j}(\lambda(\xhat))\,d\xhat,
\end{align*}
where $\bar A$, $\bar b$, and $\bar c$ are functions depending
on given data and on the actual element, namely
\begin{align*}
\bar A(\lambda) &:=
\left(\bar a_{kl}(\lambda)\right)_{k,l = 0,\ldots,d} 
  := |\det DF_S(\xhat(\lambda))| \,
\Lambda(x(\lambda)) \, A(x(\lambda)) \, \Lambda^t(x(\lambda)),\\
\bar b(\lambda) & :=
\left(\bar b_{l}(\lambda)\right)_{l = 0,\ldots,d} 
 := |\det DF_S(\xhat(\lambda))| \,
\Lambda(x(\lambda)) \, b(x(\lambda)), \quad \mbox{and}\\
\bar c(\lambda) & := |\det DF_S(\xhat(\lambda))| \, c(x(\lambda)).
\end{align*}
Having access to functions for the evaluation of $\bar A$, $\bar b$,
and $\bar c$ at given quadrature nodes, the above integrals can be
computed by some general routine for any set of local basis functions
using quadrature. Additionally, if a coefficient is piecewise constant
on the mesh, only an integration of basis functions has to be done
(compare \mathref{book:E:quad_LS} on page \pageref{book:E:quad_LS})
for this term. Here we can use pre--computed integrals of the basis
functions on the standard element and transform th-em to the actual
element. Such a computation is usually much faster than using
quadrature on each single element. Data structures for storing such
pre--computed values are described in \secref{S:ass_info}.

For the assemblage routines which we will describe now, we use
the following slight generalization: In the discretization of
the first order term, sometimes integration by parts is used
too. For a divergence free vector field
$b$ and purely Dirichlet boundary values this leads for instance to 
\[
\int_\Omega \varphi(x)\, b(x) \cdot  \nabla u(x)\,dx =
\frac12 \left(\int_\Omega \varphi(x)\, b(x) \cdot  \nabla u(x)\,dx
- \int_\Omega \nabla \varphi(x) \cdot b(x) \,  u(x)\,dx\right)
\]
yielding a modified first order term for the element matrix
\[
\int_{\Shat} \pbar^{i}(\lambda(\xhat))\;
  \frac12 \bar b(\lambda(\xhat)) \cdot 
       \nablal \pbar^{j}(\lambda(\xhat)) \,d\xhat
- \int_{\Shat} \nablal \pbar^{i}(\lambda(\xhat))\cdot 
   \frac12 \bar b(\lambda(\xhat)) \; \pbar^{j}(\lambda(\xhat)) \,d\xhat.
\]
Secondly, we allow that we have two finite element spaces with local
basis functions $\{\bar\psi_i\}_{i=1,\dots,n}$ and
$\{\bar\vphi_i\}_{i=1,\dots,m}$.

In general the following contributions of the element matrix 
$\vL_S=(L_S^{ij})_{\substack{i=1,\dots,n\\j=1,\dots,m}}$ have to be 
computed:
\begin{align*}
&\int_{\Shat}
 \nablal \bar\psi^{i}(\lambda(\xhat)) \cdot \bar A(\lambda(\xhat))\,
    \nablal \pbar^{j}(\lambda(\xhat))\,d\xhat & &\mbox{second order term,}\\
&\begin{array}{l}%
\ds \int_{\Shat} \bar\psi^{i}(\lambda(\xhat))\;
  \bar b^0(\lambda(\xhat)) \cdot \nablal \pbar^{j}(\lambda(\xhat)) \,d\xhat
\\[3mm]
\ds \int_{\Shat}^{} \nablal \bar\psi^{i}(\lambda(\xhat)) \cdot\;
  \bar b^1(\lambda(\xhat))\, \pbar^{j}(\lambda(\xhat)) \,d\xhat\\
\end{array}
& &\mbox{first order terms,}\\
&
\int_{\Shat} \bar c(\lambda(\xhat))\, \bar\psi^{i}(\lambda(\xhat)) \,
                             \pbar^{j}(\lambda(\xhat))\,d\xhat
& &\mbox{zero order term,}
\end{align*}
where for instance $\bar b^0 = \bar b$ and  $\bar b^1 = 0$, or using
integration by parts $\bar b^0 = \frac12\bar b$ and 
$\bar b^1=-\frac12\bar b$.

In order to store information about the finite element spaces, the
problem dependent functions $\bar A$, $\bar b^0$, $\bar b^1$, $\bar c$
and the quadrature that should be used for the numerical
integration of the element matrix, we define the following data
structure:
%%
\ddx{OPERATOR_INFO@{\code{OPERATOR\_INFO}}}
\idx{assemblage tools!OPERATOR_INFO@{\code{OPERATOR\_INFO}}}
%%
\newcommand{\OPERATORINFO}{\hyperref[T:OPERATOR_INFO]{\code{OPERATOR\_INFO}}\xspace}
%%
\bv\begin{lstlisting}[label=T:OPERATOR_INFO,name=OPERATOR_INFO]
typedef struct operator_info  OPERATOR_INFO;

struct operator_info
{
  const FE_SPACE *row_fe_space; /* range fe-space */
  const FE_SPACE *col_fe_space; /* domain fe-space */

  const QUAD     *quad[3];

  bool           (*init_element)(const EL_INFO *el_info,
				 const QUAD *quad[3], void *apd);
  
  union {
    const REAL_B *(*real)(const EL_INFO *el_info,
			  const QUAD *quad, int iq, void *apd);
    const REAL_BD *(*real_d)(const EL_INFO *el_info,
			     const QUAD *quad, int iq, void *apd);
    const REAL_BDD *(*real_dd)(const EL_INFO *el_info,
			       const QUAD *quad, int iq, void *apd);
  } LALt;
  MATENT_TYPE    LALt_type; /* MATENT_REAL, _REAL_D or _REAL_DD */
  bool           LALt_pw_const;
  bool           LALt_symmetric;
  int            LALt_degree;

  union {
    const REAL *(*real)(const EL_INFO *el_info,
			const QUAD *quad, int iq, void *apd);
    const REAL_D *(*real_d)(const EL_INFO *el_info,
			    const QUAD *quad, int iq, void *apd);
    const REAL_DD *(*real_dd)(const EL_INFO *el_info,
			      const QUAD *quad, int iq, void *apd);
  } Lb0;
  bool           Lb0_pw_const;
  union {
    const REAL *(*real)(const EL_INFO *el_info,
			const QUAD *quad, int iq, void *apd);
    const REAL_D *(*real_d)(const EL_INFO *el_info,
			    const QUAD *quad, int iq, void *apd);
    const REAL_DD *(*real_dd)(const EL_INFO *el_info,
			      const QUAD *quad, int iq, void *apd);
  } Lb1;
  bool           Lb1_pw_const;
  MATENT_TYPE    Lb_type; /* MATENT_REAL, _REAL_D or _REAL_DD */
  bool           Lb0_Lb1_anti_symmetric;
  int            Lb_degree;

  union {
    REAL (*real)(const EL_INFO *el_info,
		 const QUAD *quad, int iq, void *apd);
    const REAL *(*real_d)(const EL_INFO *el_info,
			  const QUAD *quad, int iq, void *apd);
    const REAL_D *(*real_dd)(const EL_INFO *el_info,
			     const QUAD *quad, int iq, void *apd);
  } c;
  bool           c_pw_const;
  MATENT_TYPE    c_type; /* MATENT_REAL, _REAL_D or _REAL_DD */
  int            c_degree;
  
  BNDRY_FLAGS    dirichlet_bndry; /* bndry-type bit-mask for
				   * Dirichlet-boundary conditions
				   * built into the matrix
				   */
  FLAGS          fill_flag;
  void           *user_data; /* application data, passed to init_element */
};
\end{lstlisting}\ev

\begin{compatibility}
  Former versions of the \ALBERTA toolkit had special
  ``\code{DOWB\_OPERATOR\_INFO}'' and \code{DOF\_DOWB\_MATRIX}''
  definitions to model block-matrix structures with $\DOW\times\DOW$
  blocks, $1\times\DOW$ and $\DOW\times 1$ blocks and $1\times 1$
  blocks (i.e. not-blocked). Because those structures included
  the scalar case as well, the ordinary scalar-only
  \code{OPERATOR\_INFO} and \code{DOF\_MATRIX} structures have been
  abandoned altogether, and the \code{\dots\_DOWB\_\dots} versions
  were renamed, dropping the bizarre \code{DOWB} component of their
  names.
\end{compatibility}

Description of the \code{OPERATOR\_INFO} structure:
\begin{descr}
\hyperitem{OPERATOR_INFO:row_fe_space}{row\_fe\_space} pointer to a
  finite element space connected to the row DOFs of the resulting
  matrix.
  %% 
\hyperitem{OPERATOR_INFO:col\_fe\_space}{col\_fe\_space} pointer to a
  finite element space connected to the column DOFs of the resulting
  matrix.
  %% 
\hyperitem{OPERATOR_INFO:quad}{quad} vector with pointers to
  quadratures; \code{quad[0]} is used for the integration of the zero
  order term, \code{quad[1]} for the first order term(s), and
  \code{quad[2]} for the second order term.
  %% 
  \idx{Per-element initializers!OPERATOR_INFO@{\code{OPERATOR\_INFO}}}
  \idx{Per-element initializers!BNDRY_OPERATOR_INFO@{\code{BNDRY\_OPERATOR\_INFO}}}
  \idx{init_element()@{\code{init\_element()}}!OPERATOR_INFO@{\code{OPERATOR\_INFO}}}
  \idx{init_element()@{\code{init\_element()}}!BNDRY_OPERATOR_INFO@{\code{BNDRY\_OPERATOR\_INFO}}}
\hyperitem{OPERATOR_INFO:init_element}{init\_element}\idx{init_element()@{\code{init\_element()}}}
  pointer to a function for doing an initialization step on each
  element; \code{init\_element} may be a \nil pointer;
  
  if \code{init\_element} is not \nil, \code{init\_element(el\_info,
    quad, user\_data)} is the first statement executed on each element
  \code{el\_info->el} and may initialize data which is used by the
  functions \code{LALt()}, \code{Lb0()}, \code{Lb1()}, and/or
  \code{c()} (calculate the Jacobian of the barycentric coordinates in
  the 1st and 2nd order terms or the element volume for all order
  terms, e.g.); \code{quad} is a pointer to a vector of quadratures
  which is actually used for the integration of the various order
  terms and \code{user\_data} may hold a pointer to user data, filled
  by \code{init\_element()}, e.g.; the return value is of interest in
  the case of parametric meshes and should be \true if the element is
  a curved element and \false otherwise.
  %% 
\hyperitem{OPERATOR_INFO:LALt}{LALt, LALt.real, LALt.real\_d,
    LALt.real\_dd}\idx{LALt()} is a pointer to a function for the
  evaluation of $\bar A$ at quadrature nodes on the element;
  \code{LALt} may be a \nil pointer, if no second order term has to be
  integrated.

  if \code{LALt} is not \nil, \code{LALt(el\_info, quad, iq,
    user\_data)} returns a pointer to a matrix of size
  ${\code{N\_LAMBDA}}\times{\code{N\_LAMBDA}}$ storing the value of
  $\bar A$ at \code{quad->lambda[iq]}; \code{quad} is the quadrature
  for the second order term and \code{user\_data} is a pointer to user
  data and \code{EL\_INFO} the current element descriptor.

  The element-type of the returned matrix is determined by
  \code{LALt\_type}, i.e. the actual return type is either
  \code{REAL\_BB} for \code{MATENT\_REAL}, \code{REAL\_BBD} for
  \code{MATENT\_REAL\_D} or \code{REAL\_BBDD} for
  \code{MATENT\_REAL\_DD}. Note that one of the \code{B}'s is missing
  in the structure definition above because the \code{LALt} is
  supposed to return the address of the first element of the matrix.
  %% 
\hyperitem{OPERATOR_INFO:LALt_type}{LALt\_type} codes the
  block-matrix type, see \code{MATENT\_TYPE} on page
  \pageref{T:MATENT_TYPE}.
  %% 
\hyperitem{OPERATOR_INFO:LALt_pw_const}{LALt\_pw\_const} should be
  \code{true} if $\bar A$ is piecewise constant on the mesh (constant
  matrix $A$ on a non--parametric mesh, e.g.); thus integration of the
  second order term can use pre--computed integrals of the basis
  functions on the standard element; otherwise integration is done by
  using quadrature on each element; this entry also influences the
  assembly on parametric meshes with \code{strategy>0}, see
  \secref{S:access_param_mesh}: \ALBERTA will assume a constant value
  of $\bar A$ for non--curved elements on a parametric mesh and
  optimize by only calling \code{LALt} once with \code{iq==0};
  
\hyperitem{OPERATOR_INFO:LALt_symmetric}{LALt\_symmetric} should be
  \code{true} if $\bar A$ is a symmetric matrix; if the finite element
  spaces for rows and columns are the same, only the diagonal and the
  upper part of the element matrix for the second order term have to
  be computed; elements of the lower part can then be set using the
  symmetry; otherwise the complete element matrix has to be
  calculated;

\hyperitem{OPERATOR_INFO:LALt_degree}{LALt\_degree} If
  \code{LALt\_pw\_const == false}, the \code{LALt\_degree} gives a
  hint about which quadrature rule should be used to integrate the
  second order term. This has only an effect if
  \hyperlink{OPERATOR_INFO:quad}{\code{quad[2]} == \nil}. In that
  case, \ALBERTA takes \code{LALt\_degree} and the row- and column
  finite element spaces into account to select a suitable quadrature
  formula.

\hyperitem{OPERATOR_INFO:Lb0}{Lb0, Lb0.real, Lb0.real\_d,
    Lb0.real\_dd} is a pointer to a function for the evaluation of
  $\bar b^0$, at quadrature nodes on the element; \code{Lb0} may be a
  \nil pointer, if this first order term has not to be integrated;

  if \code{Lb0} is not \nil, \code{Lb0(el\_info, quad, iq,
    user\_data)} returns a pointer to a vector of length
  \code{N\_LAMBDA} storing the value of $\bar b^0$ at
  \code{quad->lambda[iq]}; \code{quad} is the quadrature for the first
  order term and \code{user\_data} is a pointer to user data;

\hyperitem{OPERATOR_INFO:Lb0_pw_const}{Lb0\_pw\_const} should be
  \code{true} if $\bar b^0$ is piecewise constant on the mesh
  (constant vector $b$ on a non--parametric mesh, e.g.); thus
  integration of the first order term can use pre--computed integrals
  of the basis functions on the standard element; otherwise
  integration is done by using quadrature on each element; for
  parametric meshes the same remarks as for \code{LALt\_symmetric}
  above hold;

\hyperitem{OPERATOR_INFO:Lb1}{Lb1, Lb1.real, Lb1.real\_d,
    Lb1.real\_dd} is a pointer to a function for the evaluation of
  $\bar b^1$, at quadrature nodes on the element; \code{Lb1} may be a
  \nil pointer, if this first order term has not to be integrated;

  if \code{Lb1} is not \nil, \code{Lb1(el\_info, quad, iq,
    user\_data)} returns a pointer to a vector of length
  \code{N\_LAMBDA} storing the value of $\bar b^1$ at
  \code{quad->lambda[iq]}; \code{quad} is the quadrature for the first
  order term and \code{user\_data} is a pointer to user data;
  
\hyperitem{OPERATOR_INFO:Lb1_pw_const}{Lb1\_pw\_const} should be
  \code{true} if $\bar b^1$ is piecewise constant on the mesh
  (constant vector $b$ on a non--parametric mesh, e.g.); thus
  integration of the first order term can use pre--computed integrals
  of the basis functions on the standard element; otherwise
  integration is done by using quadrature on each element;

\hyperitem{OPERATOR_INFO:Lb_type}{Lb\_type} see \code{LALt\_type}.

\hyperitem{OPERATOR_INFO:Lb0_Lb1_anti_symmetric}{Lb0\_Lb1\_anti\_symmetric}
  should be \code{true} if the contributions of the complete first
  order term to the \emph{local} element matrix are anti symmetric
  (only possible if both \code{Lb0} and \code{Lb1} are not \nil, $\bar
  b^0 = -\bar b^1$, e.g.); if the finite element spaces for rows and
  columns are the same then only the upper part of the element matrix
  for the first order term has to be computed; elements of the lower
  part can then be set using the anti symmetry; otherwise the complete
  element matrix has to be calculated;

\hyperitem{OPERATOR_INFO:Lb_degree} See the explanations for
  \hyperlink{OPERATOR_INFO:LALt_degree}{\code{LALt\_degree}} above. 

\hyperitem{OPERATOR_INFO:c}{c, c.real, c.real\_d, c.real\_dd} is a
  pointer to a function for the evaluation of $\bar c$ at quadrature
  nodes on the element; \code{c} may be a \nil pointer, if no zero
  order term has to be integrated;

  if \code{c} is not \nil, \code{c(el\_info, quad, iq, user\_data)}
  returns the value of the function $\bar c$ at
  \code{quad->lambda[iq]}; \code{quad} is the quadrature for the zero
  order term and \code{user\_data} is a pointer to user data;

\hyperitem{OPERATOR_INFO:c_type}{c\_type} see \code{LALt\_type}.
  
\hyperitem{OPERATOR_INFO:c_pw_const}{c\_pw\_const} should be
  \code{true} if the zero order term $\bar c$ is piecewise constant on
  the mesh (constant function $c$ on a non--parametric mesh, e.g.);
  thus integration of the zero order term can use pre--computed
  integrals of the basis functions on the standard element; otherwise
  integration is done by using quadrature on each element; the same
  remarks about parametric meshes as for the other \code{*\_pw\_const}
  entries hold;

\hyperitem{OPERATOR_INFO:c_degree} See the explanations for
  \hyperlink{OPERATOR_INFO:LALt_degree}{\code{LALt\_degree}} above. 

\hyperitem{OPERATOR_INFO:dirichlet_bndry}{dirichlet\_bndry} A bit
  mask flagging those components of the boundary of the triangulation
  which are subject to Dirichlet boundary conditions. See
  \secref{S:boundary}.

\hyperitem{OPERATOR_INFO:user_data}{user\_data} optional pointer to
  memory segment for ``user data'' used by \code{init\_element()},
  \code{LALt()}, \code{Lb0()}, \code{Lb1()}, and/or \code{c()} and is
  the last argument to these functions. A better name would maybe
  ``application data'', at any rate this is the channel were an
  application program can communicate data -- like the determinant of
  the transformation to the reference element -- from
  \code{init\_element()} to the operator kernels \code{LALt} \&
  friends, without using global variables. \emph{The data behind this
    pointer must be persistent for the entire life time of the
    application program. Especially, \code{user\_data} must not point
    to the stack area of some sub-routine call.}

\hyperitem{OPERATOR_INFO:fill_flag}{fill\_flag} the flag for the mesh
  traversal routine indicating which elements should be visited and
  which information should be present in the \code{EL\_INFO} structure
  for \code{init\_element()}, \code{LALt()}, \code{Lb0()},
  \code{Lb1()}, and/or \code{c()} on the visited elements.
\end{descr}

\hrulefill

Sometimes it is necessary to add contributions of boundary integrals
to the system matrix. One example are ``Robin'' boundary conditions
(see \secref{S:robin_bound}), other important examples include
capillary boundary conditions in the context of free boundary
problems, or penalty terms to penalize tangential stresses. Another
context which requires integration over the boundaries of all mesh
elements is the implementation of discontinuous Galerkin (DG) methods.
To aid these tasks there is a \code{BNDRY\_OPERATOR\_INFO} structure,
which resembles in its layout the (bulk-) \code{OPERATOR\_INFO}
structure; it is defined as follows:

\ddx{BNDRY_OPERATOR_INFO@{\code{BNDRY\_OPERATOR\_INFO}}}
\idx{assemblage tools!BNDRY_OPERATOR_INFO@{\code{BNDRY\_OPERATOR\_INFO}}}
%%
\newcommand{\BNDRYOPERATORINFO}{\hyperref[T:BNDRY_OPERATOR_INFO]{\code{BNDRY\_OPERATOR\_INFO}}\xspace}
%%
\bv\begin{lstlisting}[label=T:BNDRY_OPERATOR_INFO,name=BNDRY_OPERATOR_INFO]
typedef struct bndry_operator_info  BNDRY_OPERATOR_INFO;

struct bndry_operator_info
{
  const FE_SPACE *row_fe_space;
  const FE_SPACE *col_fe_space;

  const WALL_QUAD *quad[3];

  bool         (*init_element)(const EL_INFO *el_info, int wall,
			       const WALL_QUAD *quad[3], void *ud);

  union {
    const REAL_B *(*real)(const EL_INFO *el_info,
			  const QUAD *quad, int iq, void *apd);
    const REAL_BD *(*real_d)(const EL_INFO *el_info,
			     const QUAD *quad, int iq, void *apd);
    const REAL_BDD *(*real_dd)(const EL_INFO *el_info,
			       const QUAD *quad, int iq, void *apd);
  } LALt;
  MATENT_TYPE    LALt_type; /* MATENT_REAL, _REAL_D or _REAL_DD */
  bool           LALt_pw_const;
  bool           LALt_symmetric;
  int            LALt_degree;

  union {
    const REAL *(*real)(const EL_INFO *el_info,
			const QUAD *quad, int iq, void *apd);
    const REAL_D *(*real_d)(const EL_INFO *el_info,
			    const QUAD *quad, int iq, void *apd);
    const REAL_DD *(*real_dd)(const EL_INFO *el_info,
			      const QUAD *quad, int iq, void *apd);
  } Lb0;
  bool           Lb0_pw_const;
  union {
    const REAL *(*real)(const EL_INFO *el_info,
			const QUAD *quad, int iq, void *apd);
    const REAL_D *(*real_d)(const EL_INFO *el_info,
			    const QUAD *quad, int iq, void *apd);
    const REAL_DD *(*real_dd)(const EL_INFO *el_info,
			      const QUAD *quad, int iq, void *apd);
  } Lb1;
  bool           Lb1_pw_const;
  MATENT_TYPE    Lb_type; /* MATENT_REAL, _REAL_D or _REAL_DD */
  bool           Lb0_Lb1_anti_symmetric;
  int            Lb_degree;

  union {
    REAL (*real)(const EL_INFO *el_info,
		 const QUAD *quad, int iq, void *apd);
    const REAL *(*real_d)(const EL_INFO *el_info,
			  const QUAD *quad, int iq, void *apd);
    const REAL_D *(*real_dd)(const EL_INFO *el_info,
			     const QUAD *quad, int iq, void *apd);
  } c;
  bool           c_pw_const;
  MATENT_TYPE    c_type; /* MATENT_REAL, _REAL_D or _REAL_DD */
  int            c_degree;

  /* boundary segment(s) we belong to; if
   * BNDRY_FLAGS_IS_INTERIOR(bndry_type), then the operator is invoked
   * on all interior faces, e.g. to implement a DG-method.
   */
  BNDRY_FLAGS  bndry_type;

  bool         discontinuous; /* assemble jumps w.r.t. the neighbour */
  bool         tangential;    /* use tangential gradients */

  FLAGS        fill_flag;
  void         *user_data;
};
\end{lstlisting}\ev

Description: Because the general layout is the same as for the
bulk-\code{OPERATOR\_INFO} structure explained above we document only
the differences here. There are three additional components in the
structure:

\begin{descr}
\hyperitem{BNDRY_OPERATOR_INFO:bndry_type}{bndry\_type} This is
  bit-mask and determines on which part of the boundary the operator
  should be invoked. See also \secref{S:boundary}. If
  \code{BNDRY\_FLAGS\_IS\_INTERIOR(bndry\_type)} evaluates to
  \code{true} (i.e. if bit $0$ is set in \code{bndry\_type}, then the
  operator is invoked on all walls of the triangulation, for instance
  to implement a DG-method.
  %% 
\hyperitem{BNDRY_OPERATOR_INFO:discontinuous}{discontinuous} This is a
  boolean flag. If set to \code{true}, then the operator is treated as
  a DG-operator. This means, that it is invoked once for each wall of
  each element with the set of local basis functions on the neighbor
  element being used to define the column space (i.e. as
  ansatz-functions) and the set of local basis function on the current
  element defining the row-space (i.e. the test-functions).
  
  One instance of \code{BNDRY\_OPERATOR\_INFO} can only be used to
  either implement a jump term or a term living on a single element.
  To have both, two instances have to be defined. To this aim
  \code{fill\_matrix\_info\_ext()} accepts multiple
  \code{BNDRY\_OPERATOR\_INFO} structures. The program-code
  \code{src/Common/ellipt-dg.c} in the \code{alberta-demo}-package
  implements a very simplistic DG-method as example: jumps over
  element boundaries are penalized by zero-order term.
  %% 
\hyperitem{BNDRY_OPERATOR_INFO:tangential}{tangential} This is a
  boolean flag. If set to \code{true}, then only the tangential
  component of the gradients of the basis functions is used.
\end{descr}

\hrulefill

Information stored in \code{OPERATOR\_INFO} and
\code{BNDRY\_OPERATOR\_INFO} structures is used by the following
functions which return a pointer to a filled \code{EL\_MATRIX\_INFO}
structure; this structure can be used as an argument to the
\code{update\_matrix()} function which will then assemble the discrete
matrix corresponding to the operators defined by the
\code{OPERATOR\_INFO} and \code{BNDRY\_OPERATOR\_INFO} structures:
%%

\fdx{fill_matrix_info()@{\code{fill\_matrix\_info()}}}
\idx{assemblage tools!fill_matrix_info()@{\code{fill\_matrix\_info()}}}
\fdx{fill_matrix_info_ext()@{\code{fill\_matrix\_info\_ext()}}}
\idx{assemblage tools!fill_matrix_info_)ext()@{\code{fill\_matrix\_info\_ext()}}}
\bv\begin{lstlisting}
const EL_MATRIX_INFO *
fill_matrix_info(const OPERATOR_INFO *operator_info,
                 EL_MATRIX_INFO *matrix_info);
const EL_MATRIX_INFO *
fill_matrix_info_ext(EL_MATRIX_INFO *matrix_info,
                     const OPERATOR_INFO *operator_info,
                     const BNDRY_OPERATOR_INFO *bop_info,
                     ... /* more bndry-ops */);
\end{lstlisting}\ev
Description:
\begin{descr}
\kitem{fill\_matrix\_info(op\_info, mat\_info)}
\kitem{fill\_matrix\_info\_ext(op\_info, mat\_info, bop\_info, \dots)}
  %% 
  ~\hfill

  Return a pointer to a filled \ELMATRIXINFO structure for the
  assemblage of the system matrix for the operator defined in
  \code{op\_info} and \code{bop\_info}. The difference between the two
  functions is that the \code{\dots\_ext()}-variant (``extended'')
  allows for additional arguments describing components of the
  differential operator which have to be assembled by boundary
  integrals. Multiple such boundary-operators can be passed to
  \code{fill\_matrix\_info\_ext()}, the final boundary operator must
  be followed by a \nil pointer. So
%%
  \bv\begin{lstlisting}
fill_matrix_info_ext(mat_info, operator_info, NULL);
\end{lstlisting}\ev
  is equivalent to
  \bv\begin{lstlisting}
fill_matrix_info(operator_info, mat_info);
\end{lstlisting}\ev
  There is the artificial restriction that at most 255 different
  \BNDRYOPERATORINFO structures may be passed.

  If the argument \code{mat\_info} is a \nil pointer, a new structure
  \code{mat\_info} is allocated and filled; otherwise the structure
  \code{mat\_info} is filled; all members are newly assigned.
  
  If the underlying finite element spaces form a direct sum, then this
  is taken care of automatically, and the return
  \ELMATRIXINFO structure will assemble block-matrices,
  where each block corresponds to the pairing of one component of the
  direct sum forming the ansatz-space and one component of the direct
  sum forming the space of test functions. See also
  \secref{S:chain_impl} and \secref{S:bfcts_chains}

  The remaining part of this section is rather a description what
  happens ``back-stage'', when calling the
  \code{fill\_matrix\_info[\_ext]()} functions. The components of
  \ELMATRIXINFO are initialized like follows:

  \begin{descr}
  \kitem{row\_fe\_space, col\_fe\_space}
    \code{op\_info->row\_fe\_space} and
    \code{op\_info->col\_fe\_space} are pointers to the finite element
    spaces (and by this to the basis functions and DOFs) connected to
    the row DOFs and the column DOFs of the matrix to be assembled; if
    both pointers are \nil pointers, an error message is given, and
    the program stops; if one of these pointers is \nil, rows and
    column DOFs are connected with the same finite element space (i.e.
    \code{op\_info->row\_fe\_space = op\_info->col\_fe\_space}, or
    \code{op\_info->col\_fe\_space = op\_info->row\_fe\_space} is
    used).

  \kitem{krn\_blk\_type} Based on the matrix block-type of the zero,
    first and second order kernels
    \hyperlink{OPERATOR_INFO:LALt_type}{\code{oinfo->c\_type}},
    \hyperlink{OPERATOR_INFO:LALt_type}{\code{oinfo->Lb\_type}} and
    \hyperlink{OPERATOR_INFO:LALt_type}{\code{oinfo->LALt\_type}} and
    on the dimension of the range of the row- and column finite
    element spaces \code{krn\_blk\_type} is set to either
    \code{MATENT\_REAL}, \code{MATENT\_REAL\_D} or
    \code{MATENT\_REAL\_DD} to reflect the block-matrix structure of
    the element matrix.

  \kitem{dirichlet\_bndry} is just a copy of
    \hyperlink{OPERATOR_INFO:dirichlet_bndry}{\code{oinfo->dirichlet\_bndry}},
    see also section \code{S:boundary}.

  \kitem{factor} is initialized to \code{1.0}. Note that the structure
    returned carries the \code{const} qualifier; the clean way to
    obtain \ELMATRIXINFO structures with a modifiable \code{factor}
    component is to pass storage to \code{fill\_matrix\_info[\_ext]()}
    via the \code{matrix\_info} parameter.

  \kitem{el\_matrix\_fct} The most important member in the structure,
    namely \code{mat\_info->el\_matrix\_fct}, is adjusted to some
    general routine for the integration of the element matrix for any
    set of local basis functions; \code{fill\_matrix\_info()} tries to
    use the fastest available function for the element integration for
    the operator defined in
    \hyperref[T:OPERATOR_INFO]{\code{op\_info}}, depending on
    \hyperlink{OPERATOR_INFO:LALt_pw_const}{\code{op\_info->LALt\_pw\_const}}
    and similar hints;

    Denote by \code{row\_degree} and \code{col\_degree} the degree of
    the basis functions connected to the rows and columns. Internally,
    a three-element vector ``\code{quad}'' of quadratures rules is
    used for the element integration, if not specified by
    \code{op\_info->quad}. The quadratures are chosen according to the
    following rules: pre-computed integrals of basis functions should
    be evaluated exactly, and all terms calculated by quadrature on
    the elements should use the same quadrature formula (this is more
    efficient than to use different quadratures). To be more specific:

    \begin{itemize}
    \item If the 2nd order term has to be integrated and
      \code{op\_info->quad[2]} is not \nil, \code{quad[2] =
        op\_info->quad[2]} is used, otherwise \code{quad[2]} is a
      quadrature which is exact of degree
      \code{row\_degree+col\_degree-2+oinfo->LALt\_degree}. If the 2nd
      order term is not integrated then \code{quad[2]} is set to \nil.

    \item If the 1st order term has to be integrated and
      \code{op\_info->quad[1]} is not \nil, \code{quad[1] =
        op\_info->quad[1]} is used; otherwise: if
      \code{op\_info->Lb\_pw\_const} is zero and \code{quad[2]} is not
      \nil, \code{quad[1] = quad[2]} is used, otherwise \code{quad[1]}
      is a quadrature which is exact of degree
      \code{row\_degree+col\_degree-1+oinfo->Lb\_degree}. If the 1st
      order term is not integrated then \code{quad[1]} is set to \nil.

    \item If the zero order term has to be integrated and
      \code{op\_info->quad[0]} is not \nil, \code{quad[0] =
        op\_info->quad[0]} is used; otherwise: if
      \code{op\_info->c\_pw\_const} is zero and \code{quad[2]} is not
      \nil, \code{quad[0] = quad[2]} is used, if \code{quad[2]} is
      \nil and \code{quad[1]} is not \nil, \code{quad[0] = quad[1]} is
      used, or if both quadratures are \nil, \code{quad[0]} is a
      quadrature which is exact of degree
      \code{row\_degree+col\_degree+oinfo->c\_degree}. If the zero
      order term is not integrated then \code{quad[0]} is set to \nil.
    \end{itemize}

    \noindent
    \code{el\_matrix\_fct()} roughly works as follows:

    \begin{itemize}
    \item If \code{op\_info->init\_element} is not \nil then a call of
      \code{op\_info->init\_element(el\_info, quad,
        op\_info->user\_data)} is the first statement of
      \code{mat\_info->el\_matrix\_fct()} on each element;
      \code{el\_info} is a pointer to the \code{EL\_INFO} structure of
      the actual element, \code{quad} is the quadrature vector
      described above (now giving information about the actually used
      quadratures), and the last argument is a pointer to the
      application-data pointer
      \hyperlink{OPERATOR_INFO:user_data}{\code{oinfo->user\_data}}.
    
    \item If \code{op\_info->LALt} is not \nil, the 2nd order term is
      integrated using the quadrature \code{quad[2]}; if
      \code{op\_info->LALt\_pw\_const} is not zero, the integrals of
      the product of gradients of the basis functions on the standard
      simplex are initialized (using the quadrature \code{quad[2]} for
      the integration) and used for the computation on the elements;
      \code{op\_info->LALt()} is only called once with arguments
      \code{op\_info->LALt(el\_info, quad[2], 0,
        op\_info->user\_data)}, i.e. the matrix of the 2nd order term
      is evaluated only at the first quadrature node; otherwise the
      integrals are approximated by quadrature and
      \code{op\_info->LALt()} is called for each quadrature node of
      \code{quad[2]}; if \code{op\_info->LALt\_symmetric} is not zero,
      the symmetry of the element matrix is used, if the finite
      element spaces are the same and this term is not integrated by
      the same quadrature as the first order term.

    \item If \code{op\_info->Lb0} is not \nil, this 1st order term is
      integrated using the quadrature \code{quad[1]}; if
      \code{op\_info->Lb0\_pw\_const} is not zero, the integrals of
      the product of basis functions with gradients of basis functions
      on the standard simplex are initialized (using the quadrature
      \code{quad[1]} for the integration) and used for the computation
      on the elements; \code{op\_info->Lb0()} is only called once with
      arguments \code{op\_info->Lb0(el\_info, quad[1], 0,
        op\_info->user\_data)}, i.e. the vector of this 1st order term
      is evaluated only at the first quadrature node; otherwise the
      integrals are approximated by quadrature and
      \code{op\_info->Lb0()} is called for each quadrature node of
      \code{quad[1]};

    \item If \code{op\_info->Lb1} is not \nil, this 1st order term is
      integrated also using the quadrature \code{quad[1]}; if
      \code{op\_info->Lb1\_pw\_const} is not zero, the integrals of
      the product of gradients of basis functions with basis functions
      on the standard simplex are initialized (using the quadrature
      \code{quad[1]} for the integration) and used for the computation
      on the elements; \code{op\_info->Lb1()} is only called once with
      arguments \code{op\_info->Lb1(el\_info, quad[1], 0,
        op\_info->user\_data)}, i.e. the vector of this 1st order term
      is evaluated only at the first quadrature node; otherwise the
      integrals are approximated by quadrature and
      \code{op\_info->Lb1()} is called for each quadrature node of
      \code{quad[1]}.

    \item If both function pointers \code{op\_info->Lb0} and
      \code{op\_info->Lb1} are not \nil, the finite element spaces for
      rows and columns are the same and
      \code{Lb0\_Lb1\_anti\_symmetric} is non--zero, then the
      contributions of the 1st order term are computed using this anti
      symmetry property.

    \item If \code{op\_info->c} is not \nil, the zero order term is
      integrated using the quadrature \code{quad[0]}; if
      \code{op\_info->c\_pw\_const} is not zero, the integrals of the
      product of basis functions on the standard simplex are
      initialized (using the quadrature \code{quad[0]} for the
      integration) and used for the computation on the elements;
      \code{op\_info->c()} is only called once with arguments
      \code{op\_info->c(el\_info, quad[0], 0, op\_info->user\_data)},
      i.e. the zero order term is evaluated only at the first
      quadrature node; otherwise the integrals are approximated by
      quadrature and \code{op\_info->c()} is called for each
      quadrature node of \code{quad[0]}.
  
    \item The functions \code{LALt()}, \code{Lb0()}, \code{Lb1()}, and
      \code{c()}, can be called in an arbitrary order on the elements,
      if not \nil (this depends on the type of integration, using
      pre--computed values, using same/different quadrature for the
      second, first, and/or zero order term, e.g.) but commonly used
      data for these functions is always initialized first by
      \code{op\_info->init\_element()}, if this function pointer is
      not \nil.

    \item Using all information about the operator and quadrature, an
      ``optimal'' routine for the assemblage is chosen.  Information
      for this routine is stored at \code{mat\_info} which includes
      the pointer to user data \code{op\_info->user\_data} (the last
      argument to \code{init\_element()}, \code{LALt()}, \code{Lb0()},
      \code{Lb1()}, and/or \code{c()}).
    \end{itemize}
    
  \kitem{neigh\_el\_mat\_fcts[]} See the documentation of the
    \hyperlink{BNDRY_OPERATOR_INFO:discontinuous}{discontinuous}
    component of the \BNDRYOPERATORINFO structure.

  \kitem{fill\_flag} Finally, the flag for the mesh traversal used by
    the function \code{update\_matrix()} is set in
    \code{mat\_info->fill\_flag} to \code{op\_info->fill\_flag}; it
    indicates which elements should be visited and which information
    should be present in the \code{EL\_INFO} structure for
    \code{init\_element()}, \code{LALt()}, \code{Lb0/1()}, and/or
    \code{c()} on the visited elements.

    If the boundary bit-mask \code{op\_info->dirichlet\_bndry} has
    bits set (see also \secref{S:boundary}), then the
    \code{FILL\_BOUND} flag is added to \code{mat\_info->fill\_flag}.
  \end{descr}
\end{descr}

\begin{example}%
[Implementation of the differential operator {\boldmath$-\Delta$}]%
\label{Ex:LALt}
The following source fragment gives an example of the implementation
for the operator $-\Delta$ and the access to a \code{MATRIX\_INFO}
structure for the automatic assemblage of the system matrix for
this problem for any set of used basis functions. 

The source fragment shown here is part of the implementation for a
Poisson equation, which is the first model problem described in detail
in \secref{S:poisson-impl}. However, we will generalize the code given in 
\secref{S:poisson-impl} to include the case of parametric meshes. 
The assemblage of the discrete system including the load vector and
Dirichlet boundary values is spelled out in \secref{S:ellipt_build}.

For the Poisson equation we only have to implement a constant second
order term. For passing information about the gradient of the
barycentric coordinates (at all quadrature points) from
\code{init\_element()} to the function
\code{LALt} we define the following structure
%
\bv\begin{lstlisting}
struct app_data
{
  REAL_BD *Lambda;
  REAL    *det;
};
\end{lstlisting}\ev
%
The function \code{init\_element()} calculates the Jacobians $\Lambda$
and determinants \code{det} of the barycentric coordinates and stores these 
in the above defined
structure. In the case of a parametric mesh we fill \code{Lambda} with
the Jacobians and \code{det} with the determinants at all quadrature
points of \code{quad[2]}. For a non--parametric mesh we only fill the
zeroth entry of \code{Lambda} and \code{det}. If
\code{init\_element()} returns
\false, then \code{LALt()} is only called once for the current simplex with
\code{iq==0}, otherwise it is called for each quadrature point in 
\code{quad[2]}. Note that we need a higher order quadrature than usual 
to calculate the integral exactly for a curved parametric element.
%
\idx{init_element()@{\code{init\_element()}}!Example}
\bv\begin{lstlisting}
static bool init_element(const EL_INFO *el_info, const QUAD *quad[3], void *ud)
{
  struct app_data *data = (struct app_data *)ud;
  PARAMETRIC      *parametric = el_info->mesh->parametric;

  if (parametric && parametric->init_element(el_info, parametric)) {
    parametric->grd_lambda(el_info, quad[2], 0, NULL,
			   data->Lambda, NULL, data->det);
    return true;
  } else {
    data->det[0] = el_grd_lambda(el_info, data->Lambda[0]);
    return false;
  }
}
\end{lstlisting}\ev
%%
The function \code{LALt} now has to calculate the scaled matrix
product $|\det D F_S| \Lambda \Lambda^t$. Note that \code{LALt()} is
invoked only for the first quadrature point (\code{iq == 0}), if the
\code{OPERATOR\_INFO}-structure claims that the second-order kernel is
piece-wise constant and \code{parametric->init\_element()} returns
\code{false}, so using the index \code{iq} into the fields \code{det}
and \code{Lambda} does not access invalid data, and the assembling
linear systems remains relatively efficient, even in the context of
iso-parametric boundary approximation.
%
\idx{LALt()!parametric example}
\bv\begin{lstlisting}
const REAL_B *LALt(const EL_INFO *el_info, const QUAD *quad, 
		   int iq, void *ud)
{
  struct app_data *data = (struct app_data *)ud;
  int             i, j;
  static REAL_BB  LALt; /* mind the "static" keyword! */

  for (i = 0; i < N_VERTICES(MESH_DIM); i++) {
    LALt[i][i]  = SCP_DOW(data->Lambda[iq][i], data->Lambda[iq][i]);
    LALt[i][i] *= data->det[iq];
    for (j = i+1; j <  N_VERTICES(MESH_DIM); j++) {
      LALt[i][j]  = SCP_DOW(data->Lambda[iq][i], data->Lambda[iq][j]);
      LALt[i][j] *= data->det[iq];
      LALt[j][i]  = LALt[i][j];
    }
  }

  return (const REAL_B *)LALt;
}
\end{lstlisting}\ev
  
Before assembling the system matrix for the operator $-\Delta$, we
first have to initialize an \code{EL\_MATRIX\_INFO} structure. A
pointer to this \code{EL\_MATRIX\_INFO} structure is the second
argument to the function \code{update\_matrix()}, which finally
assembles the system matrix (compare~\secref{S:matrix_assemblage}).

For the initialization we have to fill an \code{OPERATOR\_INFO}
structure collecting all information about the differential operator.
For $-\Delta$ we only need pointers to the functions
\code{init\_element()} and \code{LALt()} described above. The
differential operator is constant and symmetric, and information about
vertex coordinates is needed for computing the Jacobian of the
barycentric coordinates. Information about Dirichlet boundary values
should be assembled into the system matrix, hence the entry
\code{operator\_info->use\_get\_bound} is set \code{true}.

The functions \code{init\_element()} and \code{LALt()} do not
depend on the finite element space which is used. This functions
can be used for any conforming finite element discretization for the
Poisson equation; all information needed about the actually used
finite elements is a pointer to this finite element space; we assume
that this pointer is stored in the variable \code{fe\_space}.
%
\bv\begin{lstlisting} 
const EL_MATRIX_INFO *matrix_info = NULL;
static struct app_data app_data; /* Must be static! */
OPERATOR_INFO  o_info = { NULL, };

if(mesh->parametric)
  quad = get_quadrature(2, 2*fe_space->bas_fcts->degree + 2);
else
  quad = get_quadrature(2, 2*fe_space->bas_fcts->degree-2);

app_data.Lambda = MEM_ALLOC(quad->n_points, REAL_BD);
app_data.det    = MEM_ALLOC(quad->n_points, REAL);
 
o_info.quad[2]        = quad;
o_info.row_fe_space   = o_info.col_fe_space = fe_space;
o_info.init_element   = init_element;
o_info.LALt.real      = LALt;
o_info.LALt_pw_const  = true;      /* pw const. assemblage is faster */
o_info.LALt_symmetric = true;      /* symmetric assemblage is faster */
o_info.user_data      = &app_data; /* application data */

/* Use, e.g., Dirichlet  boundary conditions. */
BNDRY_FLAGS_CPY(o_info.dirichlet_bndry,	dirichlet_mask);

o_info.fill_flag = CALL_LEAF_EL|FILL_COORDS;

matrix_info = fill_matrix_info(&o_info, NULL);
\end{lstlisting}\ev
%
Full information about the operator is now available via the 
\code{matrix\_info} structure and the system matrix \code{matrix} can then 
easily be assembled for the selected finite element space by
\bv\begin{lstlisting}
  clear_dof_matrix(matrix);
  update_matrix(matrix, matrix_info, NoTranspose);
\end{lstlisting}\ev
\end{example}

\subsection{Matrix assemblage for coupled second order problems}%
\label{S:matrix_assemblage_coupled}%

The corresponding mechanism for coupled vector valued problems is very
similar, except for the additional indices necessary to describe
coupling. We start by stating the form of the element matrix with the
generalized first order term and different finite element spaces (see
also \ref{book:S:Dis2ndorder}):

\begin{align*}
L_{S,\mu\nu}^{ij} &:= \int_{\Shat}\nablal
  \bar\psi^{i}(\lambda(\xhat)) \cdot \bar A^{\mu\nu}(\lambda(\xhat))\,
  \nablal \pbar^{j}(\lambda(\xhat))\,d\xhat \\
   &+ \int_{\Shat} \bar\psi^{i}(\lambda(\xhat))\; \bar
  b^{0,\mu\nu}(\lambda(\xhat)) \cdot \nablal \pbar^{j}(\lambda(\xhat))
  \,d\xhat \\
  &+ \int_{\Shat}^{} \nablal \bar\psi^{i}(\lambda(\xhat)) \cdot\; \bar
  b^{1,\mu\nu}(\lambda(\xhat))\, \pbar^{j}(\lambda(\xhat)) \,d\xhat\\
  &+ \int_{\Shat} \bar c^{\mu\nu}(\lambda(\xhat))\,
  \bar\psi^{i}(\lambda(\xhat)) \, \pbar^{j}(\lambda(\xhat))\,d\xhat,
\end{align*}
with
\begin{align*}
  \bar A^{\mu\nu}(\lambda) &:= \left(\bar
  a_{kl}^{\mu\nu}(\lambda)\right)_{k,l = 0,\ldots,d} := |\det
  DF_S(\xhat(\lambda))| \, \Lambda(x(\lambda)) \,
  A^{\mu\nu}(x(\lambda)) \, \Lambda^t(x(\lambda)),\\
  \bar b^{0,\mu\nu}(\lambda) & := \left(\bar
  b_{l}^{0,\mu\nu}(\lambda)\right)_{l = 0,\ldots,d} := |\det
  DF_S(\xhat(\lambda))| \, \Lambda(x(\lambda)) \, b^{0,\mu\nu}(x(\lambda)),\\
  \bar b^{1,\mu\nu}(\lambda) & := \left(\bar
  b_{l}^{1,\mu\nu}(\lambda)\right)_{l = 0,\ldots,d} := |\det
  DF_S(\xhat(\lambda))| \, \Lambda(x(\lambda)) \, b^{1,\mu\nu}(x(\lambda)),
  \quad\text{and}\\
  \bar c^{\mu\nu}(\lambda) & := |\det
  DF_S(\xhat(\lambda))| \, c^{\mu\nu}(x(\lambda))
\end{align*}
for $\mu,\nu=1,\dots,n$, $n=\code{DIM\_OF\_WORLD}$.

To store information about the coupled operator and finite element
spaces, we use the same \verb|OPERATOR_INFO| (see page
\pageref{T:OPERATOR_INFO}) structure as for the scalar problems, we
only have to adjust the respective \verb|MATENT_TYPE| structure
components to the correct block-matrix type. Also, the same
\verb|fill_matrix_info()| and \verb|add_element_matrix()| routines are
used for scalar and vector valued problems.

\subsection{Data structures for storing pre-computed integrals of
  basis functions}\label{S:ass_info}%
\label{S:Q11_PSI_PHI}
\label{S:Q10_PSI_PHI}
\label{S:Q01_PSI_PHI}
\label{S:Q00_PSI_PHI}

Assume a non--parametric triangulation and constant coefficient
functions $A$, $b$, and $c$. Since the Jacobian of the barycentric
coordinates is constant on $S$, the functions $\bar A_S$, $\bar b^0_S$,
$\bar b^1_S$, and $\bar c_S$ are constant on $S$ also. Now, looking at the
element matrix approximated by some quadrature $\hat Q$, we observe
\begin{equation}\label{e:psiphiformula}
\begin{split}
\hat Q\Big(\sum_{k,l = 0}^d (\bar a_{S,kl} \bar\psi^i_{,\lambda_k} \, 
    \pbar^j_{,\lambda_l})\Big)
&=
\sum_{k,l = 0}^d \bar a_{S,kl} 
\hat Q\Big(\bar\psi^i_{,\lambda_k} \, \pbar^j_{,\lambda_l}\Big),
\\
\hat Q\Big(\sum_{l = 0}^d (
\bar b^0_{S,l} \, \bar\psi^i \, \pbar^j_{,\lambda_l})\Big)
&=
\sum_{l = 0}^d \bar b^0_{S,l} \, 
\hat Q\Big(\bar\psi^i \, \pbar^j_{,\lambda_l}\Big),
\\
\hat Q\Big(\sum_{k = 0}^d (
\bar b^1_{S,k} \, \bar\psi^i_{,\lambda_k} \, \pbar^j)\Big)
&=
\sum_{k = 0}^d \bar b^1_{S,k} \, 
\hat Q\Big(\bar\psi^i_{,\lambda_k} \, \pbar^j)\Big),
\qquad\mbox{and}\\
\hat Q\Big(
 (\bar c_S \, \bar\psi^i \,\pbar^j)\Big)
&=
\bar c_S\, \hat Q\Big(\bar\psi^i \,\pbar^j\Big).
\end{split}
\end{equation}
The values of the quadrature applied to the basis functions only
depend on the basis functions and the standard element but not on the
actual simplex $S$. All information about $S$ is given by $\bar A_S$,
$\bar b^0_S$, $\bar b^1_S$, and $\bar c_S$. Thus, these quadratures
have only to be calculated once, and can then be used on each element
during the assembling.

For this we have to store for the basis functions 
$\{\bar\psi_i\}_{i=1,\dots,n}$ and $\{\bar\vphi_j\}_{j=1,\dots,m}$
the values
\begin{alignat*}{2}
\hat Q^{11}_{ij,kl} &:=
\hat Q\Big(\bar\psi^i_{,\lambda_k} \, \pbar^j_{,\lambda_l}\Big)
&\qquad &\mbox{for } 1\leq i \leq n,\; 1\leq j\leq m,\;
0 \leq k,l \leq d,
\intertext{if $A \neq 0$,}
\hat Q^{01}_{ij,l} &:=
\hat Q\Big(\bar\psi^i \, \pbar^j_{,\lambda_l}\Big)
&&\mbox{for } 1\leq i \leq n,\; 1\leq j\leq m,\;
0 \leq l \leq d,
\intertext{if $b^0 \neq 0$,}
\hat Q^{10}_{ij,k} &:=
\hat Q\Big(\bar\psi^i_{,\lambda_k} \, \pbar^j\Big)
&&\mbox{for } 1\leq i \leq n,\; 1\leq j\leq m,\;
0 \leq k \leq d
\intertext{if $b^1 \neq 0$, and}
\hat Q^{00}_{ij} &:=
\hat Q\Big(\bar\psi^i \,\pbar^j\Big)
&&\mbox{for } 1\leq i \leq n,\; 1\leq j\leq m,
\end{alignat*}
if $c \neq 0$. Many of these values are zero, especially for the 
first and second order terms (if $\bar\psi^i$ and $\pbar^j$ are the
linear nodal basis functions $\hat Q^{11}_{ij,kl} = \delta_{ij}\delta_{kl}$).
Thus, we use special data structures
for a sparse storage of the non zero values for these terms.
These are described now.

In order to ``define'' zero entries we use
\bv\begin{lstlisting}
static const REAL TOO_SMALL = 10.0 * REAL_EPSILON;
\end{lstlisting}\ev
and all computed values \code{val} with
$|\code{val}|\leq\code{TOO\_SMALL}$ are treated as zeros. As we are 
considering here integrals over the standard simplex, non-zero
integrals are usually of order one, such that the above constant is of
the order of roundoff errors for double precision.

The following data structure is used for storing values $\hat Q^{11}$
for two sets of basis functions integrated with a given quadrature.
Note that in the context of ``chained'' basis-functions (see
\secref{S:bfcts_chains} the cache-structure nevertheless hold data for
only a single component of such a multi-component chain.

\ddx{Q11_PSI_PHI@{\code{Q11\_PSI\_PHI}}}%
\idx{assemblage tools!Q11_PSI_PHI@{\code{Q11\_PSI\_PHI}}}%
\ddx{Q11_PSI_PHI_CACHE@{\code{Q11\_PSI\_PHI\_CACHE}}}%
\idx{assemblage tools!Q11_PSI_PHI_CHACHE@{\code{Q11\_PSI\_PHI\_CACHE}}}%
\bv\begin{lstlisting}[label=T:Q11_PSI_PHI]
typedef struct q11_psi_phi   Q11_PSI_PHI;

struct q11_psi_phi
{
  const BAS_FCTS    *psi;
  const BAS_FCTS    *phi;
  const QUAD        *quad;

  const Q11_PSI_PHI_CACHE *cache;

  INIT_ELEMENT_DECL;
};

typedef struct q11_psi_phi_cache
{
  int  n_psi;
  int  n_phi;

  const int  *const*n_entries;
  const REAL *const*const*values;
  const int  *const*const*k;
  const int  *const*const*l;
} Q11_PSI_PHI_CACHE;
\end{lstlisting}\ev
Description:
\begin{descr}
  %% 
  \kitem{struct q11\_psi\_phi}:\hfill
  \begin{descr}
    %% 
    \kitem{psi} Pointer to the first set of basis functions.
    %% 
    \kitem{phi} Pointer to the second set of basis functions.
    %% 
    \kitem{quad} Pointer to the quadrature which is used for the integration.
    %% 
    \kitem{cache} Pointer to the actual data in the cache. 
    %% 
    \kitem{INIT\_ELEMENT\_DECL} Optional per-element initializer. This
    entry is initialized when calling \code{get\_q11\_psi\_phi()} if
    either the underlying basis functions or the underlying quadrature
    rule has per-element initializers. See \secref{S:init_element}.
    %% 
  \end{descr}
  \kitem{struct q11\_psi\_phi\_cache}:\hfill\
  \begin{descr}
    \kitem{n\_psi} Dimension of the local space of test-functions (row
    space), equals \code{Q11\_PSI\_PHI.psi->n\_bas\_fcts}.
    %% 
    \kitem{n\_phi} Dimension of the local space of ansatz-functions
    (column space), equals \code{Q11\_PSI\_PHI.phi->n\_bas\_fcts}.
    %% 
    \kitem{n\_entries} matrix of size $\code{n\_psi}\times\code{n\_phi}$
    storing the count of non zero integrals;

    \code{n\_entries[i][j]} is the count of non zero values of $\hat
    Q^{11}_{\code{ij,kl}}$ ($0\leq\code{k,l}\leq d$) for the pair
    $(\code{psi[i]},\code{phi[j]})$,
    $\code{0}\leq\code{i}<\code{n\_psi}$,
    $\code{0}\leq\code{j}<\code{n\_phi}$.
    % 
    \kitem{values} tensor storing the non zero integrals;
    
    \code{values[i][j]} is a vector of length \code{n\_entries[i][j]}
    storing the non zero values for the pair
    $(\code{psi[i]},\code{phi[j]})$.
    %% 
    \kitem{k, l} tensor storing the indices $k$, $l$ of the non zero
    integrals;
    
    \code{k[i][j]} and \code{l[i][j]} are vectors of length
    \code{n\_entries[i][j]} storing at \code{k[i][j][r]} and
    \code{l[i][j][r]} the indices \code{k} and \code{l} of the value
    stored at \code{values[i][j][r]}, i.e.
  \end{descr}
\end{descr}
The following formulas summarize the relationship between the cache
data-structure and the formulas \eqref{e:psiphiformula} at the
beginning of this section:
\[
\code{values[i][j][r]} = 
\hat Q^{11}_{\code{ij},\code{k[i][j][r],l[i][j][r]}}
= \hat Q\Big(\bar\psi^{\code{i}}_{,\lambda_{\code{k[i][j][r]}}} \,
\pbar^{\code{j}}_{,\lambda_{\code{l[i][j][r]}}}\Big),
\]
for $\code{0}\leq\code{r}<\code{n\_entries[i][j]}$.
Using these pre--computed values we have for all elements $S$
\[
\sum_{k,l = 0}^d \bar a_{S,kl} 
\hat Q\Big(\bar\psi^i_{,\lambda_k} \, \pbar^j_{,\lambda_l}\Big)
=
\sum_{\code{r}=\code{0}}^{\code{n\_entries[i][j]-1}}
\bar a_{S,\code{k[i][j][r],l[i][j][r]}}\,\code{*}
       {\code{values[i][j][r]}}.
\]
The following function initializes a \code{Q11\_PSI\_PHI} structure:
\fdx{get_q11_psi_phi()@{\code{get\_q11\_psi\_phi())}}}%
\idx{assemblage tools!get_q11_psi_phi()@{\code{get\_q11\_psi\_phi()}}}%
\bv\begin{lstlisting}
const Q11_PSI_PHI *get_q11_psi_phi(const BAS_FCTS *psi, const BAS_FCTS *phi,
                                   const QUAD *quad);
\end{lstlisting}\ev
Description:
\begin{descr}
  \kitem{get\_q11\_psi\_phi(psi, phi, quad)} returns a pointer to a filled
  \code{Q11\_PSI\_PHI} structure.

  \code{psi} is a pointer to the first set of basis functions,
  \code{phi} is a pointer to the second set of basis functions;
  if both are \nil pointers, nothing is done and the return value
  is \nil; if one of the pointers is a \nil pointer, the structure
  is initialized using the same set of basis functions for the
  first and second set, i.e. \code{phi = psi} or \code{psi = phi}
  is used.
  
  \code{quad} is a pointer to a quadrature for the approximation
  of the integrals; if \code{quad} is \nil, then a quadrature which 
  is exact of degree \code{psi->degree+phi->degree-2} is used.

  All used \code{Q11\_PSI\_PHI} structures are stored in a linked
  list and are identified uniquely by the members \code{psi},
  \code{phi}, and \code{quad}, and for such a combination only
  one \code{Q11\_PSI\_PHI} structure is created during runtime;

  First, \code{get\_q11\_psi\_phi()}
  looks for a matching structure in the linked list; if such a 
  structure is found a pointer to this structure is returned;
  the values are not computed a second time.
  Otherwise a new structure is generated, linked to the list, and
  the values are computed using the quadrature \code{quad}; all
  values \code{val} with $|\code{val}|\leq\code{TOO\_SMALL}$ are
  treated as zeros.
\end{descr}

\begin{example}
The following example shows how to use these pre--computed values
for the integration of the 2nd order term
\[
\int_{\Shat}
 \nablal \bar\psi^{i}(\lambda(\xhat)) \cdot \bar A(\lambda(\xhat))\,
    \nablal \pbar^{j}(\lambda(\xhat))\,d\xhat
\]
for all $i,j$. We only show the body of a function for the integration
and assume that \code{LALt\_fct} returns a matrix storing $\bar A$
(compare the member \code{LALt} in the \code{OPERATOR\_INFO} structure):
\bv\begin{lstlisting}
  static Q11_PSI_PHI_CACHE *q11_psi_phi;

  if (!q11_psi_phi) {
    q11_psi_phi = get_q11_psi_phi(psi, phi, quad[2])->cache;
  }

  LALt = LALt_fct(el_info, quad, 0, user_data);
  n_entries = q11_psi_phi->n_entries;

  for (i = 0; i < q11_psi_phi->n_psi; i++)
  {
    for (j = 0; j < q11_psi_phi->n_phi; j++)
    {
      k      = q11_psi_phi->k[i][j];
      l      = q11_psi_phi->l[i][j];
      values = q11_psi_phi->values[i][j];
      for (val = m = 0; m < n_entries[i][j]; m++)
        val += values[m]*LALt[k[m]][l[m]];
      mat[i][j] += val;
    }
  }
\end{lstlisting}\ev
\end{example}

The values $\hat Q^{01}$ for the set of basis functions \code{psi}
and \code{phi} are stored in
\ddx{Q01_PSI_PHI@{\code{Q01\_PSI\_PHI}}}%
\idx{assemblage tools!Q01_PSI_PHI@{\code{Q01\_PSI\_PHI}}}%
\ddx{Q01_PSI_PHI_CACHE@{\code{Q01\_PSI\_PHI\_CACHE}}}%
\idx{assemblage tools!Q01_PSI_PHI_CACHE@{\code{Q01\_PSI\_PHI\_CACHE}}}%
\bv\begin{lstlisting}[label=T:Q01_PSI_PHI]
typedef struct q01_psi_phi   Q01_PSI_PHI;

typedef struct q01_psi_phi_cache
{
  int  n_psi;
  int  n_phi;

  const int  *const*n_entries;
  const REAL *const*const*values;
  const int  *const*const*l;
} Q01_PSI_PHI_CACHE;

struct q01_psi_phi
{
  const BAS_FCTS    *psi;
  const BAS_FCTS    *phi;
  const QUAD        *quad;

  const Q01_PSI_PHI_CACHE *cache;

  INIT_ELEMENT_DECL;
};
\end{lstlisting}\ev
Description:
\begin{descr}  
  \kitem{struct q01\_psi\_phi}:\hfill
  \begin{descr}  
    %%
    \kitem{psi} pointer to the first set of basis functions.
    %%
    \kitem{phi} pointer to the second set of basis functions.
    %% 
    \kitem{quad} pointer to the quadrature which is used for the integration.
    %% 
    \kitem{cache} Pointer to the actual data in the cache. 
    %% 
    \kitem{INIT\_ELEMENT\_DECL} Optional per-element initializer. This
    entry is initialized when calling \code{get\_q11\_psi\_phi()} if
    either the underlying basis functions or the underlying quadrature
    rule has per-element initializers. See \secref{S:init_element}.
    %% 
  \end{descr}
  \kitem{struct q01\_psi\_phi\_cache}:\hfill
  \begin{descr}
    \kitem{n\_psi} Dimension of the local space of test-functions (row
    space), equals \code{Q11\_PSI\_PHI.psi->n\_bas\_fcts}.
    %% 
    \kitem{n\_phi} Dimension of the local space of ansatz-functions
    (column space), equals \code{Q11\_PSI\_PHI.phi->n\_bas\_fcts}.
    %% 
    \kitem{n\_entries} matrix of size 
    $\code{psi->n\_bas\_fcts}\times\code{phi->n\_bas\_fcts}$
    storing the count of non zero integrals;
    
    \code{n\_entries[i][j]} is the count of non zero values
    of $\hat Q^{01}_{\code{ij,l}}$ ($0\leq\code{l}\leq d$)
    for the pair $(\code{psi[i]},\code{phi[j]})$,
    $\code{0}\leq\code{i}<\code{n\_psi}$,
    $\code{0}\leq\code{j}<\code{n\_phi}$.
    %% 
    \kitem{values} tensor storing the non zero integrals;
    
    \code{values[i][j]} is a vector of length \code{n\_entries[i][j]}
    storing the non zero values for the pair 
    $(\code{psi[i]},\code{phi[j]})$.
    %% 
    \kitem{l} tensor storing the indices $l$ of the non zero integrals;
    
    \code{l[i][j]} is a vector of length \code{n\_entries[i][j]}
    storing at \code{l[i][j][r]} the index \code{l} of the
    value stored at \code{values[i][j][r]}, i.e.
    %% 
  \end{descr}
\end{descr}
The following formulas summarize the relationship between the cache
data-structure and the formulas \eqref{e:psiphiformula} at the
beginning of this section:
\[
\code{values[i][j][r]} = 
\hat Q^{01}_{\code{ij},\code{l[i][j][r]}}
= \hat Q\Big(\bar\psi^{\code{i}} \,
\pbar^{\code{j}}_{,\lambda_{\code{l[i][j][r]}}}\Big),
\]
for $\code{0}\leq\code{r}<\code{n\_entries[i][j]}$.
Using these pre--computed values we have for all elements $S$
\[
\sum_{l = 0}^d \bar b^0_{S,l} 
\hat Q\Big(\bar\psi^i \, \pbar^j_{,\lambda_l}\Big)
=
\sum_{\code{r}=\code{0}}^{\code{n\_entries[i][j]-1}}
\bar b^0_{S,\code{l[i][j][r]}}\,\code{*} {\code{values[i][j][r]}}.
\]
The following function initializes a \code{Q01\_PSI\_PHI} structure:
\fdx{get_q01_psi_phi()@{\code{get\_q01\_psi\_phi())}}}%
\idx{assemblage tools!get_q01_psi_phi()@{\code{get\_q01\_psi\_phi()}}}%
\bv\begin{lstlisting}
const Q01_PSI_PHI *get_q01_psi_phi(const BAS_FCTS *psi, const BAS_FCTS *phi, 
                                   const QUAD *quad);
\end{lstlisting}\ev
Description:
\begin{descr}
\kitem{get\_q01\_psi\_phi(psi, phi, quad)} returns a pointer to a filled
       \code{Q01\_PSI\_PHI} structure.

      \code{psi} is a pointer to the first set of basis functions 
      \code{phi} is a pointer to the second set of basis functions;
      if both are \nil pointers, nothing is done and the return value
      is \nil; if one of the pointers is a \nil pointer, the structure
      is initialized using the same set of basis functions for the
      first and second set, i.e. \code{phi = psi} or \code{psi = phi}
      is used.
      
      \code{quad} is a pointer to a quadrature for the approximation
      of the integrals; is \code{quad} is \nil, a quadrature which 
      is exact of degree \code{psi->degree+phi->degree-1} is used.

      All used \code{Q01\_PSI\_PHI} structures are stored in a linked
      list and are identified uniquely by the members \code{psi},
      \code{phi}, and \code{quad}, and for such a combination only
      one \code{Q01\_PSI\_PHI} structure is created during runtime;

      First, \code{get\_q01\_psi\_phi()}
      looks for a matching structure in the linked list; if such a 
      structure is found a pointer to this structure is returned;
      the values are not computed a second time.
      Otherwise a new structure is generated, linked to the list, and
      the values are computed using the quadrature \code{quad}; all
      values \code{val} with $|\code{val}|\leq\code{TOO\_SMALL}$ are
      treated as zeros.
\end{descr}
The values $\hat Q^{10}$ for the set of basis functions \code{psi}
and \code{phi} are stored in
\ddx{Q10_PSI_PHI@{\code{Q10\_PSI\_PHI}}}%
\idx{assemblage tools!Q10_PSI_PHI@{\code{Q10\_PSI\_PHI}}}%
\ddx{Q10_PSI_PHI_CACHE@{\code{Q10\_PSI\_PHI\_CACHE}}}%
\idx{assemblage tools!Q10_PSI_PHI_CACHE@{\code{Q10\_PSI\_PHI\_CACHE}}}%
\bv\begin{lstlisting}[label=T:Q10_PSI_PHI]
typedef struct q10_psi_phi   Q10_PSI_PHI;

typedef struct q10_psi_phi_cache
{
  int        n_psi;
  int        n_phi;

  const int  *const*n_entries;
  const REAL *const*const*values;
  const int  *const*const*k;
} Q10_PSI_PHI_CACHE;

struct q10_psi_phi
{
  const BAS_FCTS    *psi;
  const BAS_FCTS    *phi;
  const QUAD        *quad;

  const Q10_PSI_PHI_CACHE *cache;

  INIT_ELEMENT_DECL;
};
\end{lstlisting}\ev
Description:
\begin{descr}
  \kitem{struct q10\_psi\_phi}:\hfill
  \begin{descr}
    \kitem{psi} pointer to the first set of basis functions.
    %% 
    \kitem{phi} pointer to the second set of basis functions.
    %% 
    \kitem{quad} pointer to the quadrature which is used for the integration.
    %% 
    %% 
    \kitem{cache} Pointer to the actual data in the cache. 
    %% 
    \kitem{INIT\_ELEMENT\_DECL} Optional per-element initializer. This
    entry is initialized when calling \code{get\_q11\_psi\_phi()} if
    either the underlying basis functions or the underlying quadrature
    rule has per-element initializers. See \secref{S:init_element}.
    %% 
  \end{descr}
  \kitem{struct q10\_psi\_phi\_cache}:\hfill
  \begin{descr}
    \kitem{n\_psi} Dimension of the local space of test-functions (row
    space), equals \code{Q11\_PSI\_PHI.psi->n\_bas\_fcts}.
    %% 
    \kitem{n\_phi} Dimension of the local space of ansatz-functions
    (column space), equals \code{Q11\_PSI\_PHI.phi->n\_bas\_fcts}.
    %% 
    \kitem{n\_entries} matrix of size 
    $\code{psi->n\_bas\_fcts}\times\code{phi->n\_bas\_fcts}$
    storing the count of non zero integrals;
    
    \code{n\_entries[i][j]} is the count of non zero values
    of $\hat Q^{10}_{\code{ij,k}}$ ($0\leq\code{k}\leq d$)
    for the pair $(\code{psi[i]},\code{phi[j]})$,
    $\code{0}\leq\code{i}<\code{n\_psi}$,
    $\code{0}\leq\code{j}<\code{n\_phi}$.
    %% 
    \kitem{values} tensor storing the non zero integrals;
    
    \code{values[i][j]} is a vector of length \code{n\_entries[i][j]}
    storing the non zero values for the pair 
    $(\code{psi[i]},\code{phi[j]})$.
    %% 
    \kitem{k} tensor storing the indices $k$ of the non zero integrals;
    
    \code{k[i][j]} is a vector of length \code{n\_entries[i][j]}
    storing at \code{k[i][j][r]} the index \code{k} of the
    value stored at \code{values[i][j][r]}, i.e.
    %% 
  \end{descr}
\end{descr}
The following formulas summarize the relationship between the cache
data-structure and the formulas \eqref{e:psiphiformula} at the
beginning of this section:
\[
\code{values[i][j][r]} = 
\hat Q^{10}_{\code{ij},\code{k[i][j][r]}}
= \hat Q\Big(\bar\psi^{\code{i}}_{,\lambda_{\code{k[i][j][r]}}} \,
\pbar^{\code{j}}\Big),
\]
for $\code{0}\leq\code{r}<\code{n\_entries[i][j]}$.
Using these pre--computed values we have for all elements $S$
\[
\sum_{k = 0}^d \bar b^1_{S,k} 
\hat Q\Big(\bar\psi^i_{,\lambda_k} \, \pbar^j\Big)
=
\sum_{\code{r}=\code{0}}^{\code{n\_entries[i][j]-1}}
\bar b^1_{S,\code{k[i][j][r]}}\,\code{*} {\code{values[i][j][r]}}.
\]
The following function initializes a \code{Q10\_PSI\_PHI} structure:
\fdx{get_q10_psi_phi()@{\code{get\_q10\_psi\_phi())}}}%
\idx{assemblage tools!get_q10_psi_phi()@{\code{get\_q10\_psi\_phi()}}}%
\bv\begin{lstlisting}
const Q10_PSI_PHI *get_q10_psi_phi(const BAS_FCTS *psi, const BAS_FCTS *phi, 
                                   const QUAD *quad);
\end{lstlisting}\ev
Description:
\begin{descr}
\kitem{get\_q10\_psi\_phi(psi, phi, quad)} returns a pointer to a filled
       \code{Q10\_PSI\_PHI} structure.

      \code{psi} is a pointer to the first set of basis functions 
      \code{phi} is a pointer to the second set of basis functions;
      if both are \nil pointers, nothing is done and the return value
      is \nil; if one of the pointers is a \nil pointer, the structure
      is initialized using the same set of basis functions for the
      first and second set, i.e. \code{phi = psi} or \code{psi = phi}
      is used.
      
      \code{quad} is a pointer to a quadrature for the approximation
      of the integrals; is \code{quad} is \nil, a quadrature which 
      is exact of degree \code{psi->degree+phi->degree-1} is used.

      All used \code{Q10\_PSI\_PHI} structures are stored in a linked
      list and are identified uniquely by the members \code{psi},
      \code{phi}, and \code{quad}, and for such a combination only
      one \code{Q10\_PSI\_PHI} structure is created during runtime;

      First, \code{get\_q10\_psi\_phi()}
      looks for a matching structure in the linked list; if such a 
      structure is found a pointer to this structure is returned;
      the values are not computed a second time.
      Otherwise a new structure is generated, linked to the list, and
      the values are computed using the quadrature \code{quad}; all
      values \code{val} with $|\code{val}|\leq\code{TOO\_SMALL}$ are
      treated as zeros.
\end{descr}
Finally, the values $\hat Q^{00}$ for the set of basis functions \code{psi}
and \code{phi} are stored in
\ddx{Q00_PSI_PHI@{\code{Q00\_PSI\_PHI}}}%
\idx{assemblage tools!Q00_PSI_PHI@{\code{Q00\_PSI\_PHI}}}%
\ddx{Q00_PSI_PHI_CACHE@{\code{Q00\_PSI\_PHI\_CACHE}}}%
\idx{assemblage tools!Q00_PSI_PHI_CACHE@{\code{Q00\_PSI\_PHI\_CACHE}}}%
\bv\begin{lstlisting}[label=T:Q00_PSI_PHI]
typedef struct q00_psi_phi   Q00_PSI_PHI;

typedef struct q00_psi_phi_cache
{
  int        n_psi;
  int        n_phi;

  const REAL *const*values;
} Q00_PSI_PHI_CACHE;

struct q00_psi_phi
{
  const BAS_FCTS    *psi;
  const BAS_FCTS    *phi;
  const QUAD        *quad;

  const Q00_PSI_PHI_CACHE *cache;

  INIT_ELEMENT_DECL;
};

\end{lstlisting}\ev
Description:
\begin{descr}
  \kitem{struct q00\_psi\_phi}:\hfill
  \begin{descr}
    \kitem{psi} pointer to the first set of basis functions.
    %% 
    \kitem{phi} pointer to the second set of basis functions.
    %% 
    \kitem{quad} pointer to the quadrature which is used for the integration.
    %% 
    %% 
    \kitem{cache} Pointer to the actual data in the cache. 
    %% 
    \kitem{INIT\_ELEMENT\_DECL} Optional per-element initializer. This
    entry is initialized when calling \code{get\_q11\_psi\_phi()} if
    either the underlying basis functions or the underlying quadrature
    rule has per-element initializers. See \secref{S:init_element}.
    %% 
  \end{descr}
  \kitem{struct q00\_psi\_phi\_cache}:\hfill
  \begin{descr}
    \kitem{n\_psi} Dimension of the local space of test-functions (row
    space), equals \code{Q11\_PSI\_PHI.psi->n\_bas\_fcts}.
    %% 
    \kitem{n\_phi} Dimension of the local space of ansatz-functions
    (column space), equals \code{Q11\_PSI\_PHI.phi->n\_bas\_fcts}.
    %% 
    \kitem{values} matrix storing the integrals.
  \end{descr}
\end{descr}
The following formulas summarize the relationship between the cache
data-structure and the formulas \eqref{e:psiphiformula} at the
beginning of this section:
\[
\code{values[i][j]} = 
\hat Q^{00}_{\code{ij}}
= \hat Q\Big(\bar\psi^{\code{i}} \, \pbar^{\code{j}}\Big),
\]
       for the pair $(\code{psi[i]},\code{phi[j]})$,
       $\code{0}\leq\code{i}<\code{psi->n\_bas\_fcts}$,
       $\code{0}\leq\code{j}<\code{phi->n\_bas\_fcts}$.
The following function initializes a \code{Q00\_PSI\_PHI} structure:
\fdx{get_q00_psi_phi()@{\code{get\_q00\_psi\_phi())}}}%
\idx{assemblage tools!get_q00_psi_phi()@{\code{get\_q00\_psi\_phi()}}}%
\bv\begin{lstlisting}
const Q00_PSI_PHI *get_q00_psi_phi(const BAS_FCTS *psi, const BAS_FCTS *phi, 
                                   const QUAD *quad);
\end{lstlisting}\ev
Description:
\begin{descr}
\kitem{get\_q00\_psi\_phi(psi, phi, quad)} returns a pointer to a filled
       \code{Q00\_PSI\_PHI} structure.

      \code{psi} is a pointer to the first set of basis functions 
      \code{phi} is a pointer to the second set of basis functions;
      if both are \nil pointers, nothing is done and the return value
      is \nil; if one of the pointers is a \nil pointer, the structure
      is initialized using the same set of basis functions for the
      first and second set, i.e. \code{phi = psi} or \code{psi = phi}
      is used.
      
      \code{quad} is a pointer to a quadrature for the approximation
      of the integrals; is \code{quad} is \nil, a quadrature which 
      is exact of degree \code{psi->degree+phi->degree} is used.

      All used \code{Q00\_PSI\_PHI} structures are stored in a linked
      list and are identified uniquely by the members \code{psi},
      \code{phi}, and \code{quad}, and for such a combination only
      one \code{Q00\_PSI\_PHI} structure is created during runtime;

      First, \code{get\_q00\_psi\_phi()}
      looks for a matching structure in the linked list; if such a 
      structure is found a pointer to this structure is returned;
      the values are not computed a second time.
      Otherwise a new structure is generated, linked to the list, and
      the values are computed using the quadrature \code{quad}.
\end{descr}

\subsection{Data structures and functions for updating coefficient vectors}%
\label{S:vector_update}%

\sloppypar Besides the general routines \code{update\_real\_vec()},
\code{update\_real\_d\_vec()} and \code{update\_real\_vec\_dow()},
this section presents also easy to use routines for calculation of
$L^2$ scalar products between a given function and all basis functions
of a finite element space, taken either over the interior of the mesh
or over boundary segments.

\medskip

The following structures hold full information for the assembling
of element vectors. They are used by the functions
\code{update\_real\_vec()} and \code{update\_real\_d\_vec()} described below.
\ddx{EL_VEC_INFO@{\code{EL\_VEC\_INFO}}}
\idx{assemblage tools!EL_VEC_INFO@{\code{EL\_VEC\_INFO}}}
\ddx{EL_VEC_D_INFO@{\code{EL\_VEC\_D\_INFO}}}
\idx{assemblage tools!EL_VEC_D_INFO@{\code{EL\_VEC\_D\_INFO}}}
\ddx{EL_VEC_INFO_D@{\code{EL\_VEC\_INFO\_D}}}
\idx{assemblage tools!EL_VEC_INFO_D@{\code{EL\_VEC\_INFO\_D}}}
\bv\begin{lstlisting}[name=EL_VEC_INFO,label=T:EL_VEC_INFO]
typedef struct el_vec_info    EL_VEC_INFO;
typedef struct el_vec_d_info  EL_VEC_D_INFO;
typedef struct el_vec_info_d  EL_VEC_INFO_D;

typedef const EL_REAL_VEC *
(*EL_VEC_FCT)(const EL_INFO *el_info, void *fill_info);

typedef struct el_vec_info  EL_VEC_INFO;
struct el_vec_info
{
  const FE_SPACE *fe_space;

  BNDRY_FLAGS    dirichlet_bndry;
  REAL           factor;

  EL_VEC_FCT     el_vec_fct;
  void           *fill_info;

  FLAGS          fill_flag;
};

typedef const EL_REAL_D_VEC *
(*EL_VEC_D_FCT)(const EL_INFO *el_info, void *fill_info);

typedef struct el_vec_d_info  EL_VEC_D_INFO;
struct el_vec_d_info
{
  const FE_SPACE *fe_space;

  BNDRY_FLAGS    dirichlet_bndry;
  REAL           factor;

  EL_VEC_D_FCT   el_vec_fct;
  void           *fill_info;

  FLAGS          fill_flag;
};

typedef const EL_REAL_VEC_D *
(*EL_VEC_FCT_D)(const EL_INFO *el_info, void *fill_info);

typedef struct el_vec_info_d  EL_VEC_INFO_D;
struct el_vec_info_d
{
  const FE_SPACE *fe_space;

  BNDRY_FLAGS    dirichlet_bndry;
  REAL           factor;

  EL_VEC_FCT_D   el_vec_fct;
  void           *fill_info;

  FLAGS          fill_flag;
};
\end{lstlisting}\ev
Description:
\begin{descr}
\kitem{fe\_space} the underlying finite element space
%%
\kitem{dirichlet\_bndry} a bit mask marking the boundary segments
which are subject to Dirichlet boundary conditions, see also
\secref{S:boundary}.
%%
\kitem{factor} is a multiplier for the element contributions; usually
\code{factor} is \code{1} or \code{-1}.
%%
\kitem{el\_vec\_fct} is a pointer to a function for the computation of
the element vector; \code{el\_vec\_fct(el\_info, fill\_info)} returns
a pointer to an element vector of the respective data type, see e.g.
\verb|EL_REAL_VEC| on page \pageref{T:EL_REAL_VEC}. This vector stores
the computed values for the element described by \code{el\_info}.
\code{fill\_info} is a pointer to data needed by
\code{el\_vec\_fct()}; the function has to provide memory for storing
the element vector, which can be overwritten on the next call.
%%
\kitem{fill\_info} pointer to data needed by \code{el\_vec\_fct()}; is
the second argument of this function.
%%
\kitem{fill\_flag} the flag for the mesh traversal for assembling the
vector.
\end{descr}
The following function does the update of vectors by assembling element
contributions during mesh traversal; information for computing
the element vectors is held in a \code{EL\_VEC[\_D]\_INFO} structure:
\fdx{update_real_vec()@{\code{update\_real\_vec()}}}
\idx{assemblage tools!update_real_vec()@{\code{update\_real\_vec()}}}
\fdx{update_real_d_vec()@{\code{update\_real\_d\_vec()}}}
\idx{assemblage tools!update_real_d_vec()@{\code{update\_real\_d\_vec()}}}
\fdx{update_real_vec_dow()@{\code{update\_real\_vec\_dow()}}}
\idx{assemblage tools!update_real_vec_dow()@{\code{update\_real\_vec\_dow()}}}
\bv\begin{lstlisting}
void update_real_vec(DOF_REAL_VEC *dv, const EL_VEC_INFO *vec_info);
void update_real_d_vec(DOF_REAL_D_VEC *dv, const EL_VEC_D_INFO *vec_info);
void update_real_vec_dow(DOF_REAL_VEC_D *dv, const EL_VEC_INFO_D *vec_info)
\end{lstlisting}\ev
\begin{descr}
  \kitem{update\_real[\_d]\_vec[\_dow](dv, info)} updates the vector
  \code{dr} by traversing the underlying mesh and assembling the
  element contributions into the vector; information about the
  computation of element vectors and connection of local and global
  DOFs is stored in \code{info}.

  The flags for the mesh traversal of the mesh
  \code{dv->fe\_space->mesh} are stored at \code{info->fill\_flags}
  which specifies the elements to be visited and information that
  should be present on the elements for the calculation of the element
  vectors and boundary information (if \code{info->get\_bound} is not
  \nil).

  On the elements, information about the global DOFs is accessed by
  \code{info->get\_dof} using \code{info->admin}; the boundary type of
  the DOFs is accessed by \code{info->get\_bound} if
  \code{info->get\_bound} is not a \nil pointer; then the element
  vector is computed by \code{info->el\_vec\_fct(el\_info,
    info->fill\_info)}; these contributions are finally added to
  \code{dv} multiplied by \code{info->factor} by a call of
  \code{add\_element[\_d]\_vec[\_dow]()} with all information about global
  DOFs, the element vector, and boundary types, if available;

  \code{update\_real[\_d]\_vec[\_dow]()} only adds element
  contributions; this makes several calls for the assemblage of one
  vector possible; before the first call, the vector should be set to
  zero by a call of \code{dof\_set[\_d|\_dow](0.0, dv)}.
\end{descr}

\paragraph{$L^2$- and $H^1$-scalar- products over the bulk phase}
In many applications, the load vector is just the $L^2$- or
$H^1$-scalar-product of a given function with all basis functions of
the finite element space or this scalar product is a part of the right
hand side. Such a scalar product can be directly assembled by the
following functions.

\paragraph{Prototypes}
\fdx{L2scp_fct_bas()@{\code{L2scp\_fct\_bas()}}}
\idx{assemblage tools!L2scp_fct_bas()@{\code{L2scp\_fct\_bas()}}}
\fdx{L2scp_fct_bas_d()@{\code{L2scp\_fct\_bas\_d()}}}
\idx{assemblage tools!L2scp_fct_bas_d()@{\code{L2scp\_fct\_bas\_d()}}}
\fdx{L2scp_fct_bas_dow()@{\code{L2scp\_fct\_bas\_dow()}}}
\idx{assemblage tools!L2scp_fct_bas_dow()@{\code{L2scp\_fct\_bas\_dow()}}}
\fdx{L2scp_fct_bas_loc()@{\code{L2scp\_fct\_bas\_loc()}}}
\idx{assemblage tools!L2scp_fct_bas_loc()@{\code{L2scp\_fct\_bas\_loc()}}}
\fdx{L2scp_fct_bas_loc_dow()@{\code{L2scp\_fct\_bas\_loc\_dow()}}}
\idx{assemblage tools!L2scp_fct_bas_loc_dow()@{\code{L2scp\_fct\_bas\_loc\_dow()}}}
\fdx{H1scp_fct_bas()@{\code{H1scp\_fct\_bas()}}}
\idx{assemblage tools!H1scp_fct_bas()@{\code{H1scp\_fct\_bas()}}}
\fdx{H1scp_fct_bas_dow()@{\code{H1scp\_fct\_bas\_dow()}}}
\idx{assemblage tools!H1scp_fct_bas_dow()@{\code{H1scp\_fct\_bas\_dow()}}}
%%
\bv\begin{lstlisting}
void L2scp_fct_bas(FCT_AT_X f, const QUAD *quad, DOF_REAL_VEC *fh);
void L2scp_fct_bas_d(FCT_D_AT_X f, const QUAD *, DOF_REAL_D_VEC *fh);
void L2scp_fct_bas_dow(FCT_D_AT_X f, const QUAD *quad, DOF_REAL_VEC_D *fh);

void L2scp_fct_bas_loc(DOF_REAL_VEC *fh,
                       LOC_FCT_AT_QP f, void *f_data, FLAGS fill_flag,
                       const QUAD *quad);
void L2scp_fct_bas_loc_d(DOF_REAL_D_VEC *fh,
                         LOC_FCT_D_AT_QP f, void *fd, FLAGS fill_flag,
                         const QUAD *quad);
void L2scp_fct_bas_loc_dow(DOF_REAL_VEC_D *fh,
                           LOC_FCT_D_AT_QP f, void *fd, FLAGS fill_flag,
                           const QUAD *quad);

void H1scp_fct_bas(GRD_FCT_AT_X grd_f,
                   const QUAD *quad, DOF_REAL_VEC *fh);
void H1scp_fct_bas_d(GRD_FCT_D_AT_X grd_f,
                     const QUAD *quad, DOF_REAL_D_VEC *fh);
void H1scp_fct_bas_dow(GRD_FCT_D_AT_X grd_f,
                        const QUAD *quad, DOF_REAL_VEC_D *fh);

void H1scp_fct_bas_loc(DOF_REAL_VEC *fh,
                       GRD_LOC_FCT_AT_QP grd_f, void *fd,
                       FLAGS fill_flag, const QUAD *quad);
void H1scp_fct_bas_loc_d(DOF_REAL_VEC_D *fh,
                         GRD_LOC_FCT_D_AT_QP grd_f, void *fd,
                         FLAGS fill_flag, const QUAD *quad);
void H1scp_fct_bas_loc_dow(DOF_REAL_VEC_D *fh,
                           GRD_LOC_FCT_D_AT_QP grd_f, void *fd,
                           FLAGS fill_flag, const QUAD *quad);
\end{lstlisting}\ev

\paragraph{Descriptions}
\begin{descr}
  \kitem{L2scp\_fct\_bas(f, quad, fh)}
  \kitem{L2scp\_fct\_bas\_d(f, quad, fh)}
  \kitem{L2scp\_fct\_bas\_dow(f, quad, fh)}
  %%
  Approximate the $L^2$ scalar products of a given function with all
  basis functions by numerical quadrature and add the corresponding
  values to a DOF vector

       \code{f} is a pointer for the evaluation of the given
       function in world coordinates $x$ and returns the 
       value of that function at $x$; if \code{f} is a \nil
       pointer, nothing is done;

       \code{fh} is the DOF vector where at the \code{i}--th entry
       the approximation of the $L^2$ scalar product of the given
       function with the \code{i}--th global basis function of
       \code{fh->fe\_space} is added;

       \code{quad} is the quadrature for the approximation of
       the integral on each leaf element of \code{fh->fe\_space->mesh};
       if \code{quad} is a \nil pointer, a default quadrature
       which is exact of degree \code{2*fh->fe\_space->bas\_fcts->degree-2}
       is used.

       The integrals are approximated by looping over all 
       leaf elements, computing the approximations to the element
       contributions and adding these values to the vector \code{fh}
       by  \code{add\_element\_vec()}.

       The vector \code{fh} is \emph{not} initialized with \code{0.0};
       only the new contributions are added.
       %%
\kitem{L2scp\_fct\_bas\_d(fd, quad, fhd)} approximates the $L^2$ scalar 
       products of a given vector valued function with all scalar valued
       basis functions by numerical quadrature and adds the
       corresponding values to a vector valued DOF vector;

       \code{fd} is a pointer for the evaluation of the given function
       in world coordinates $x$; \code{fd(x, fx)} returns a pointer to
       a vector storing the value at \code{x}; if \code{fx} is not
       \nil, the value is stored at \code{fx} otherwise the function
       has to provide memory for storing this vector, which can be
       overwritten on the next call; if \code{fd} is a \nil pointer,
       nothing is done;

       \code{fhd} is the DOF vector where at the \code{i}--th entry (a
       \code{REAL\_D} vector) the approximation of the $L^2$ scalar
       product of the given vector valued function with the \code{i}--th
       global (scalar valued) basis function of \code{fhd->fe\_space}
       is added;

       \code{quad} is the quadrature for the approximation of
       the integral on each leaf element of \code{fhd->fe\_space->mesh};
       if \code{quad} is a \nil pointer, a default quadrature
       which is exact of degree \code{2*fhd->fe\_space->bas\_fcts->degree-2}
       is used.

       The integrals are approximated by looping over all 
       leaf elements, computing the approximations to the element
       contributions and adding these values to the vector \code{fhd}
       by \code{add\_element\_d\_vec()}.

       The vector \code{fhd} is \emph{not} initialized with
       $(\code{0.0},\ldots,\code{0.0})$; only the new contributions are
       added.
\kitem{L2scp\_fct\_bas\_dow(fd, quad, fhd)}
\kitem{L2scp\_fct\_bas\_loc(fh, f\_at\_qp, fct\_data, fill\_flag, quad)}
\kitem{L2scp\_fct\_bas\_loc\_dow(fh, f\_at\_qp, ud, fill\_flag, quad)}
\kitem{H1scp\_fct\_bas(grd\_f, quad, fh)}
\kitem{H1scp\_fct\_bas\_dow(grd\_fd, quad, fhd)}
\end{descr}

\subsection{Boundary conditions}\label{S:boundary_conditions}
%%
The following six functions act as a front-end to the functions
explained further below, therefore we refer the reader to
\secref{S:dirichlet_bound}, \ref{S:neumann_bound} and
\ref{S:robin_bound} for a deeper discussion of the implementation of
Dirichlet, Neumann and Robin boundary conditions within \ALBERTA.

% aus l2scp.c
\fdx{boundary_conditions()@{\code{boundary\_conditions()}}}
\idx{assemblage tools!boundary_conditions()@{\code{boundary\_conditions()}}}
\fdx{boundary_conditions_loc()@{\code{boundary\_conditions\_loc()}}}
\idx{assemblage tools!boundary_conditions_loc()@{\code{boundary\_conditions\_loc()}}}
\fdx{boundary_conditions_dow()@{\code{boundary\_conditions\_dow()}}}
\idx{assemblage tools!boundary_conditions_dow()@{\code{boundary\_conditions\_dow()}}}
\fdx{boundary_conditions_loc_dow()@{\code{boundary\_conditions\_loc\_dow()}}}
\idx{assemblage tools!boundary_conditions_loc_dow()@{\code{boundary\_conditions\_loc\_dow()}}}
\bv\begin{lstlisting}
bool boundary_conditions(DOF_MATRIX *matrix, DOF_REAL_VEC *fh,
                         DOF_REAL_VEC *uh, DOF_SCHAR_VEC *bound,
                         const BNDRY_FLAGS dirichlet_segment,
                         REAL (*g)(const REAL_D x),
                         REAL (*gn)(const REAL_D x, const REAL_D normal),
                         REAL alpha_r, const WALL_QUAD *wall_quad);
bool boundary_conditions_loc(DOF_MATRIX *matrix, DOF_REAL_VEC *fh,    
                             DOF_REAL_VEC *uh, DOF_SCHAR_VEC *bound,
                             const BNDRY_FLAGS dirichlet_segment,
                             LOC_FCT_AT_QP g_at_qp, LOC_FCT_AT_QP gn_at_qp,
                             void *app_data, FLAGS fill_flags,
                             REAL alpha_r, const WALL_QUAD *wall_quad);
bool boundary_conditions_d(DOF_MATRIX *matrix, DOF_REAL_D_VEC *fh,
                           DOF_REAL_D_VEC *uh, DOF_SCHAR_VEC *bound,
                           const BNDRY_FLAGS dirichlet_segment,
                           const REAL *(*g)(const REAL_D x, REAL_D res),
                           const REAL *(*gn)(const REAL_D x,
                                             const REAL_D normal,
                                             REAL_D res),
                           REAL alpha_r, const WALL_QUAD *wall_quad);
bool boundary_conditions_loc_d(DOF_MATRIX *matrix, DOF_REAL_D_VEC *fh,
                               DOF_REAL_D_VEC *uh, DOF_SCHAR_VEC *bound,
                               const BNDRY_FLAGS dirichlet_segment,
                               LOC_FCT_D_AT_QP g_at_qp,
                               LOC_FCT_D_AT_QP gn_at_qp,
                               void *app_data, FLAGS fill_flags,
                               REAL alpha_r, const WALL_QUAD *wall_quad);
bool boundary_conditions_dow(DOF_MATRIX *matrix, DOF_REAL_VEC_D *fh,
                             DOF_REAL_VEC_D *uh, DOF_SCHAR_VEC *bound,
                             const BNDRY_FLAGS dirichlet_segment,
                             const REAL *(*g)(const REAL_D x, REAL_D res),
                             const REAL *(*gn)(const REAL_D x,
                                               const REAL_D normal,
                                               REAL_D res),
                             REAL alpha_r, const WALL_QUAD *wall_quad);
bool boundary_conditions_loc_dow(DOF_MATRIX *matrix, DOF_REAL_VEC_D *fh,
                                 DOF_REAL_VEC_D *uh, DOF_SCHAR_VEC *bound,
                                 const BNDRY_FLAGS dirichlet_segment,
                                 LOC_FCT_D_AT_QP g_at_qp,
                                 LOC_FCT_D_AT_QP gn_at_qp,
                                 void *app_data, FLAGS fill_flags,
                                 REAL alpha_r, const WALL_QUAD *wall_quad);
\end{lstlisting}\ev
Description: These ``compound'' functions implement Dirichlet, Neumann
or Robin boundary conditions, and optionally perform a mean-value
correction of the ``right hand side'' in the context of pure Neumann
boundary conditions if $\code{alpha\_r}<0$ (in order to satisfy the
conditions for the ``right hand side'' which may be violated in the
discrete context because of quadrature errors).

For the differences between the code{\dots\_loc()} and
non-\code{\dots\_loc()} versions the reader is referred to the section
dealing with \code{dirichlet\_bound\_loc()} (see
\secref{S:dirichlet_bound}). A brief discussion of the calling
convention for the various functions pointers passed to the library
functions can also be found in \secref{S:user_fcts}.
\begin{description}
\item[Parameters]~\hfill
  \begin{descr}
    \kitem{matrix} As explained in \secref{S:robin_bound}, passed on
    to \code{robin\_bound()}.
    %%
    \kitem{fh} As explained in \secref{S:dirichlet_bound},
    \secref{S:neumann_bound} and \secref{S:robin_bound}. Passed on to
    \code{dirichlet\_bound()} and \code{bndry\_L2scp\_fct\_bas()}.
    %%
    \kitem{uh} As explained in \secref{S:dirichlet_bound}. Passed on
    to \code{dirichlet\_bound()}.
    %%
    \kitem{bound} As explained in \secref{S:dirichlet_bound}. Passed on
    to \code{dirichlet\_bound()}.
    %%
    \kitem{dirichlet\_segment} As explained in
    \secref{S:dirichlet_bound}. Passed on to
    \code{dirichlet\_bound()}. The respective bit-masks passed to
    \code{bndry\_L2scp\_fct\_bas()} and \code{robin\_bound()} are
    computed as bit-wise complement of \code{dirichlet\_segment}. See
    also \secref{S:boundary}.
    %%
    \kitem{g} As explained in \secref{S:dirichlet_bound}.  Passed on
    to \code{dirichlet\_bound()}.
    %%
    \kitem{gn} As explained in \secref{S:neumann_bound},
    \secref{S:robin_bound}. Passed on to
    \code{bndry\_L2scp\_fct\_bas()}.
    %%
    \kitem{app\_data} \code{\dots\_loc()}-variants only. As explained
    in \secref{S:dirichlet_bound}, \secref{S:user_fcts}. Passed on as
    application-data pointer to the application provided function
    hooks.
    %%
    \kitem{fill\_flags} \code{\dots\_loc()}-variants only. Additional
    fill-flags needed by \code{g()} or \code{gn}.
    %%
    \kitem{alpha\_r} As explained in \secref{S:robin_bound}. Passed on
    to \code{robin\_bound()}. \code{alpha\_r} is also abused to
    request an automatic mean-value correction of the load-vector: if
    \code{alpha\_r} is negative, and Neumann boundary conditions were
    imposed on all boundary segments, then
    \code{boundary\_conditions()} will attempt such a mean-value
    correction in order to keep fulfill the compatibility condition
    for the load-vector in the discrete setting. Of course, if the
    differential operator has lower order parts, then the load-vector
    need not have mean-value $0$.

    Robin boundary conditions will only be assembled if
    \code{alpha\_r} is strictly larger than $0$.
    %%
    \kitem{wall\_quad} As explained in \secref{S:robin_bound} and
    \secref{S:neumann_bound}. Passed on to \code{robin\_bound()} and
    \code{bndry\_L2scp\_fct\_bas()}.
  \end{descr}
\item[Return Value]~\hfill

  \code{true} if any part of the boundary was subject to Dirichlet
  boundary conditions.
\end{description}

\subsubsection{Dirichlet boundary conditions}
\label{S:dirichlet_bound}

For the solution of the discrete system \mathref{book:E:LuF} on page
\pageref{book:E:LuF} derived in \secref{book:S:Dis2ndorder}, we have
to set the Dirichlet boundary values for all Dirichlet DOFs. Usually,
we take for the approximation $g_h$ of $g$ the interpolant of $g$,
i.e. $g_h = I_h g$ and we have to copy the coefficients of $g_h$ at
the Dirichlet DOFs to the start value for an iterative solver. This
way the first matrix-vector operation performed by an iterative solver
will have the effect to transform an inhomogeneous Dirichlet boundary
problem to a homogeneous one by applying the differential operator to
the boundary values and subtracting the result from the ``right hand
side''.  Whether or not it is also necessary to copy these
coefficients to the load vector depends on how the matrices act on the
coefficients:

\begin{itemize}
\item If the matrix-rows corresponding to Dirichlet-nodes
  $k_1,\,\dots,\,k_M$ have been replaced by unit-vectors $e_{k_l}$
  $(1\leq l\leq M)$, then it is also necessary to copy the Dirichlet
  nodes to the load vector (compare \mathref{book:E:Fright} on page
  \pageref{book:E:Fright}). Copying the coefficients of $g_h$ at the
  Dirichlet DOFs to the initial guess will result in an initial
  residual (and then for all subsequent residuals) which is zero at
  all Dirichlet \code{DOF}s.

  This is the case when Dirichlet bit-masks have been copied to
  \code{OPERATOR\_INFO.dirichlet\_bndry} (compare \secref{S:boundary}
  and \ref{T:OPERATOR_INFO}); the resulting \code{DOF\_MATRIX} will
  then be assembled (\secref{S:matrix_assemblage}) in this way,
  replacing any row corresponding to a Dirichlet-node by the
  corresponding unit-vector.
\item If the matrix-rows corresponding to Dirichlet-nodes have not
  been replaced by unit-vectors, then it is still possible to solve a
  Dirichlet-problem by passing a \code{DOF\_SCHAR\_VEC} to the
  matrix-vector routines (compare \secref{S:solver}, describing the
  linear solvers available in \ALBERTA).  However, in this case the
  matrix-vector subroutines will clear all Dirichlet-nodes to zero,
  see \secref{S:DOF_BLAS}.  Therefore, in this case it is necessary to
  clear the coefficients of the ``right hand side'' which correspond
  to Dirichlet-nodes. See the \exampleref{E:CLEARING_DIRICHLET_NODES}
  for simple examples how to perform this task.
\end{itemize}
Note that the matrices generated this way -- either by clearing
Dirichlet-rows or by masking out Dirichlet rows -- are not symmetric
(compare also \mathref{book:E:matrix} on page \pageref{book:E:matrix})
even if the underlying differential operator is symmetric. However,
the restriction of the matrix to the space spanned by the
non-Dirichlet \code{DOFs} \emph{is} symmetric, so any iterative solver
for symmetric matrices will work, provided one either sets the
Dirichlet-values also in the load-vector (if the matrix acts as
identity on the Dirichlet \code{DOF}s) or clears the Dirichlet-nodes
in the load-vector (if the matrix acts as zero-operator on the
Dirichlet \code{DOF}s).

The following functions will set Dirichlet boundary values
for all DOFs on the Dirichlet boundary, using an interpolation
of the boundary values $g$:
\fdx{dirichlet_bound()@{\code{dirichlet\_bound()}}}
\idx{assemblage tools!dirichlet_bound()@{\code{dirichlet\_bound()}}}
\fdx{dirichlet_bound_loc()@{\code{dirichlet\_bound\_loc()}}}
\idx{assemblage tools!dirichlet_bound_loc()@{\code{dirichlet\_bound\_loc()}}}
\fdx{dirichlet_bound_d()@{\code{dirichlet\_bound\_d()}}}
\idx{assemblage tools!dirichlet_bound_d()@{\code{dirichlet\_bound\_d()}}}
\fdx{dirichlet_bound_loc_d()@{\code{dirichlet\_bound\_loc\_d()}}}
\idx{assemblage tools!dirichlet_bound_loc_d()@{\code{dirichlet\_bound\_loc\_d()}}}
\fdx{dirichlet_bound_loc_dow()@{\code{dirichlet\_bound\_loc\_dow()}}}
\idx{assemblage tools!dirichlet_bound_loc_dow()@{\code{dirichlet\_bound\_loc\_dow()}}}
\bv\begin{lstlisting}
bool dirichlet_bound(DOF_REAL_VEC *fh, DOF_REAL_VEC *uh, 
                     DOF_SCHAR_VEC *bound,
                     const BNDRY_FLAGS dirichlet_segments,
                     REAL (*g)(const REAL_D));
bool dirichlet_bound_d(DOF_REAL_VEC_D *fh, DOF_REAL_VEC_D *uh,
                       DOF_SCHAR_VEC *bound,
                       const BNDRY_FLAGS dirichlet_segments,
                       const REAL *(*g)(const REAL_D, REAL_D));
bool dirichlet_bound_dow(DOF_REAL_VEC_D *fh, DOF_REAL_VEC_D *uh,
                         DOF_SCHAR_VEC *bound,
                         const BNDRY_FLAGS dirichlet_segments,
                         const REAL *(*g)(const REAL_D, REAL_D));
bool dirichlet_bound_loc(DOF_REAL_VEC *fh, DOF_REAL_VEC *uh,
                         DOF_SCHAR_VEC *bound,
                         const BNDRY_FLAGS dirichlet_segments,
                         LOC_FCT_AT_QP g, void *ud, FLAGS fill_flags);
bool dirichlet_bound_loc_d(DOF_REAL_VEC_D *fh, DOF_REAL_VEC_D *uh,
                           DOF_SCHAR_VEC *bound,
                           const BNDRY_FLAGS dirichlet_segments,
                           LOC_FCT_D_AT_QP g, void *ud,
                           FLAGS fill_flags);
bool dirichlet_bound_loc_dow(DOF_REAL_VEC_D *fh, DOF_REAL_VEC_D *uh,
                             DOF_SCHAR_VEC *bound,
                             const BNDRY_FLAGS dirichlet_segments,
                             LOC_FCT_D_AT_QP g, void *ud,
                             FLAGS fill_flags);
\end{lstlisting}\ev

\paragraph{Descriptions}
\begin{descr}
  \kitem{dirichlet\_bound(fh, uh, bound, dirichlet\_segments, g)} sets
  Dirichlet boundary values for all DOFs on all leaf elements of
  \code{fh->fe\_space->mesh} or \code{uh->fe\_space->mesh}; values at
  DOFs not belonging to the Dirichlet boundary are not changed by this
  function.

  Boundary values are set element-wise on the leaf elements. The
  boundary type of the walls of an element is accessed through the
  function \code{wall\_bound(el\_info, wall)}. If the corresponding
  bit is set in \code{dirichlet\_segments}, then the local
  interpolation operator of the basis functions underlying
  \code{fh/uh->fe\_space} is invoked to compute the coefficients of
  the \code{DOF}s located on that wall.

  This variant of the \code{dirichlet\_bound...()} is rather
  simplistic; the \code{dirichlet\_bound\_loc..()} pass more
  information to the function implementing the boundary values and
  also allow for manipulating the amount of information available
  while looping over the mesh.
  %%
  \begin{description}
  \item[Parameters]~\hfill
    \begin{description}
    \item[\code{fh}, \code{uh}] are vectors where Dirichlet boundary
      values should be set (usually, \code{fh} stores the load vector
      and \code{uh} an initial guess for an iterative solver); one of
      \code{fh} and \code{uh} may be a \nil pointer; if both arguments
      are \nil pointers, nothing is done, except of filling the
      \code{bound} argument, it that is non \nil; if both arguments are
      not \nil, \code{fh->fe\_space} must equal \code{uh->fe\_space}.
     \item[\code{bound}] is a vector for storing the boundary type for
       each used DOF; \code{bound} may be \nil; if it is not \nil, the
       \code{i}-th entry of the vector is filled with the boundary
       type of the \code{i}-th DOF. \code{bound->fe\_space} must be
       the same as \code{fh}'s or \code{uh}'s \code{fe\_space}.
     \item[\code{dirichlet\_segments}] Bit-mask marking those parts of
       the boundary which are subject to Dirichlet boundary
       conditions. If bit number $k>0$ is set in
       \code{dirichlet\_segments} then the part of the boundary with
       boundary classification $k$ is marked as Dirichlet boundary.
       Compare \secref{S:boundary}.
     \item[\code{REAL (*g)(const REAL\_D arg)}] is a pointer to a
       function for the evaluation of the boundary data; if \code{g}
       is a \nil pointer, all coefficients at Dirichlet DOFs are set
       to \code{0.0}. \code{arg} are the Cartesian co-ordinates where
       the value of \code{g} should be computed.
     \end{description}
   \item[Return Value] \code{true} if any part of the boundary of the
     mesh is subject to Dirichlet boundary conditions, as indicated by
     the argument \code{dirichlet\_segments}, \code{false} otherwise.
  \end{description}
  %%
  \hrulefill
  %%
  \kitem{dirichlet\_bound\_d(fh, uh, bound, dirichlet\_segments, g)}
  does the same as \code{dirichlet\_bound()}, but \code{fh} and
  \code{uh} are \code{DOF\_REAL\_D\_VEC} vectors.

  The calling convention for \code{const REAL (*g)(const REAL\_D arg,
    REAL\_D result)} is such that \code{g} must allow for
  \code{result} being a \nil-pointer. If so, a pointer to a statically
  allocated storage area must be returned, otherwise \code{result}
  must be filled with the value of \code{g} at the evaluation point
  \code{arg}, see \exampleref{E:user_fct_d} in \secref{S:user_fcts}.
  %% 
  Otherwise everything works as for \code{dirichlet\_bound()}, see
  above for the documentation.
  %%

  \hrulefill
  %%
  \kitem{dirichlet\_bound\_dow(fh, uh, bound, dirichlet\_segments, g)}
  does the same as \code{dirichlet\_bound\_d()}, but \code{fh} and
  \code{uh} are \code{DOF\_REAL\_VEC\_D} vectors, that is, \code{uh}
  and \code{fh} may belong to a finite element space which is a direct
  sum, composed of several finite element spaces (note the location of
  the \code{\_D} suffix in the data-type names
  \code{DOF\_REAL\_VEC\_D} and \code{DOF\_REAL\_D\_VEC}!). The calling
  convention for \code{const REAL (*g)(const REAL\_D arg, REAL\_D
    result)} is the same as explained above for
  \code{dirichlet\_bound\_d()}.

  \hrulefill
  %%
  \kitem{dirichlet\_bound\_loc(fh, uh, bound, dirichlet\_segments, g,
    ud, fill\_flags)}
  %%
  This function differs from its counterpart without the
  \code{\_loc}-suffix only in the calling convention for the function
  implementing the Dirichlet boundary conditions. We document only the
  differing or additional arguments here and refer the reader to the
  documentation of \code{dirichlet\_bound()} above:
  \begin{description}
  \item[\code{REAL (*g)(const EL\_INFO *el\_info, const QUAD *quad,
      int iq, void *ud)}] The function pointer to the function
    implementing the Dirichlet boundary values. In contrast to the
    corresponding function-pointer passed to \code{dirichlet\_bound()}
    this function is invoked with a co-dimension $1$ quadrature rule
    (compare the \code{interpol}-hooks in the \code{BAS\_FCTS}
    structure, \ref{T:BAS_FCTS}, and the definition of the
    \code{QUAD}-structure, \ref{T:QUAD}) and a quadrature point, and
    first and not least with a filled \code{EL\_INFO}-structure.

    This means that \code{g} has full-access to all the information
    available through the \code{EL\_INFO} element descriptor. The
    amount of data filled-in during mesh-traversal can additionally be
    controlled by setting specific fill-flags through the argument
    \code{fill\_flags}, which is passed as last argument to
    \code{dirichlet\_bound\_loc()}. The last argument \code{ud} to
    \code{g} is the same as the pointer \code{ud} passed as pre-last
    argument to \code{dirichlet\_bound\_loc()} and may be used by an
    application to reduce the amount of global variables and thus the
    potential of bugs implied by the use of global variables. The
    following simple example shows how to get hold of the Cartesian
    co-ordinates of the quadrature point, and how to use, e.g., the
    boundary classification available through the \code{EL\_INFO}
    structure:
    %%
    \begin{example}
      \label{E:loc_fct_at_qp}
    \bv\begin{lstlisting}[name=DIRICHLET_BOUND_LOC_EXAMPLE,label=C:DIRICHLET_BOUND_LOC_EXAMPLRE]
      struct g_data
      {
        BNDRY_TYPE special_wall_type; /* other stuff */
      };

      REAL g_impl(const EL_INFO *el_info, const QUAD *quad, int iq, void *ud)
      {
        struct g_data *data = (struct g_data *)ud;
        BNDRY_TYPE btype = wall_bound(el_info, quad->sub_splx);
        REAL result;
        const QUAD_EL_CACHE *qelc =
          fill_quad_el_cache(el_info, quad, FILL_EL_QUAD_WORLD);

        if (btype == data->special_wall_type) {
          ... /* do special things */
          return sin(qelc->world[iq][0];
        } else {
          ... /* do "normal" things */
          return sin(qelc->world[iq][1];
        }
      }

      ... /* 1.000.000 lines of code later ... */
      struct g_data g_data_instance = { 42 };
      dirichlet_bound_loc(fh, uh, bound, dirichlet_bits, g_impl, &g_data_instance, FILL_COORDS|FILL_MACRO_WALLS);
    \end{lstlisting}\ev
    \end{example}
    %%
  \item[\code{ud}] Application-data-pointer, forwarded as \code{ud}
    argument to the application supplied \code{g} function-pointer.
    %%
  \item[\code{fill\_flags}] Additional fill-flags for the loop over
    the mesh. The application must make sure that \code{fill\_flags}
    contains all flags corresponding to information needed by the
    function \code{g()}.
  \end{description}
  %%
  \hrulefill
  %%
  \kitem{dirichlet\_bound\_loc\_dow(}
  \kitem{~~~~~~~~~~~~fh, uh, bound, dirichlet\_segments, g, ud, fill\_flags)}
  %%
  \kitem{dirichlet\_bound\_loc(fh, uh, bound, dirichlet\_segments, g, ud, fill\_flags)}
  %%
  These two function differ from \code{dirichlet\_bound\_loc()} only
  in the calling convention for
  \bv\begin{lstlisting}
    const REAL *(*g)(REAL_D result, const EL_INFO *el_info, const QUAD *quad, int iq, void *ud).
  \end{lstlisting}\ev
  As in the example \ref{C:FCT_D_AT_X_ABI} the implementation of
  \code{g()} must allow for \code{result} being \nil, returning a
  pointer to a static storage area in this case.
\end{descr}

\begin{example}
  \label{E:CLEARING_DIRICHLET_NODES}
  This example demonstrates how to clear the Dirichlet-nodes in the
  load-vector if Dirichlet boundary conditions are implemented using a
  \code{DOF\_SCHAR\_VEC} to mask-out Dirichlet nodes. Note that this
  example applies \emph{only} if the \code{DOF\_SCHAR\_VEC} is also
  passed to the linear solvers. Otherwise Dirichlet boundary
  conditions have to be incorporated into the matrix
  \begin{description}
  \item[\emph{scalar problem}]
    \bv\begin{lstlisting}[name=CLEARING_SCALAR_DIRICHLET_NODES,label=C:CLEARING_SCALAR_DIRICHLET_NODES]
      extern REAL g(const REAL_D x); /* defined somewhere else, e.g. */
      extern DOF_REAL_VEC *uh, *fh;  /* defined somewhere else, e.g. */
      DOF_SCHAR_VEC *bound =
        get_dof_schar_vec("dirichlet mask vector", fh->fe_space);
      BNDRY_FLAGS dirichlet_bits;
      BNDRY_FLAGS_INIT(dirichlet_bits);
      BNDRY_FLAGS_SET(dirichlet_bits, 3); /* e.g. */

      ... /* other stuff */

      dirichlet_bound(NULL, uh, bound, dirichlet_bits, g);
      FOR_ALL_DOFS(fh->fe_space->admin,
                   if (bound->vec[dof] >= DIRICHLET) {
                     fh->vec[dof] = 0.0;
                   });

     ... /* other stuff */

     oem_solve_s(matrix, bound, fh, uh, ... /* other parameters */);
    \end{lstlisting}\ev
  \item[\emph{simple vector valued problem}]
    \bv\begin{lstlisting}[name=CLEARING_VECTOR_DIRICHLET_NODES,label=C:CLEARING_VECTOR_DIRICHLET_NODES]
      extern const REAL *g(const REAL_D x, REAL_D result); /* defined somewhere else, e.g. */
      extern DOF_REAL_D_VEC *uh, *fh; /* defined somewhere else, e.g. */
      extern DOF_SCHAR_VEC *bound;
      extern BNDRY_FLAGS dirichlet_bits;

      ... /* other stuff */

      dirichlet_bound_d(NULL, uh, bound, dirichlet_bits, g);
      FOR_ALL_DOFS(fh->fe_space->admin,
                   if (bound->vec[dof] >= DIRICHLET) {
                     SET_DOW(0.0, fh->vec[dof]);
                   });

     ... /* other stuff */

     oem_solve_d(matrix, bound, fh, uh, ... /* other parameters */);
    \end{lstlisting}\ev
  \item[\emph{vector valued problem, using an FE-space which is a direct sum}]
    ~\hfill

    %%
    Note the difference between a \code{DOF\_REAL\_D\_VEC} which
    contains \code{DIM\_OF\_WORLD}-sized \code{REAL\_D} vectors and a
    \code{DOF\_REAL\_VEC\_D} which contains scalars of type
    \code{REAL} if the underlying basis function are themselves
    vector-valued, or \code{REAL\_D}-vectors if the underlying basis
    functions are scalar-valued. The first code-block of the
    \code{FOREACH\_DOF\_DOW}-macro is for the case where the basis
    functions are vector-valued (and hence the coefficients are
    scalars) and the second code-block is for the case where the basis
    functions are scalar-valued (and hence the coefficients are
    vectors). The \code{FOREACH\_DOF\_DOW()} macro is further
    explained in \secref{S:chain_loops}.
    %%
    \bv\begin{lstlisting}[name=CLEARING_CHAINED_DIRICHLET_NODES,label=C:CLEARING_CHAINED_DIRICHLET_NODES]
      extern const REAL *g(const REAL_D x, REAL_D result); /* defined somewhere else, e.g. */
      extern DOF_REAL_VEC_D *uh, *fh; /* defined somewhere else, e.g. */
      extern DOF_SCHAR_VEC *bound;
      extern BNDRY_FLAGS dirichlet_bits;

      ... /* other stuff */

      dirichlet_bound_dow(NULL, uh, bound, dirichlet_bits, g);
      FOREACH_DOF_DOW(fh->fe_space,
                      if (bound->vec[dof] >= DIRICHLET) {
                        fh->vec[dof] = 0.0;
                      },
                      if (bound->vec[dof] >= DIRICHLET) {
                        SET_DOW(0.0, ((DOF_REAL_D_VEC *)fh)->vec[dof]);
                      },
                      fh = CHAIN_NEXT(fh, DOF_REAL_VEC_D);
                      bound = CHAIN_NEXT(bound, DOF_SCHAR_VEC));

     ... /* other stuff */

     oem_solve_dow(matrix, bound, fh, uh, ... /* other parameters */);
    \end{lstlisting}\ev
  \end{description}
\end{example}

\subsubsection{Neumann boundary conditions}
\label{S:neumann_bound}

For the implementation of inhomogeneous Neumann boundary conditions it
is necessary to compute $L^2$ scalar products between the
inhomogeneity and the basis functions on the Neumann boundary
segments. The following functions compute the $L^2$ scalar product
over the boundary of the domain. They return \code{true} if not all
boundary segments of the mesh belong to the segment specified by
\code{bndry\_seg}. If \code{bndry\_seg == NULL} then the scalar
product is computed over the entire boundary (i.e. over all walls
without neighbour).
%%
Besides computing the $L^2$-scalar product over boundary segments
there are also functions to compute the $L^2$-scalar-product over trace
meshes (or ``sub-meshes'', see \secref{S:submesh_implementation}).
%%
For the calling conventions for the application provided function
pointers the reader is referred to \secref{S:user_fcts}, and the
relevant part of the discussion of \code{dirichlet\_bound\_loc()} in
\secref{S:dirichlet_bound}.

All function work additive, the contributions of the per-element
integrals are added to any data already present in \code{fh}.
%%
\paragraph{Prototypes}
\fdx{bndry_L2scp_fct_bas()@{\code{bndry\_L2scp\_fct\_bas()}}}
\idx{assemblage tools!bndry_L2scp_fct_bas()@{\code{bndry\_L2scp\_fct\_bas()}}}
\fdx{bndry_L2scp_fct_bas_loc()@{\code{bndry\_L2scp\_fct\_bas\_loc()}}}
\idx{assemblage tools!bndry_L2scp_fct_bas_loc()@{\code{bndry\_L2scp\_fct\_bas\_loc()}}}
\fdx{bndry_L2scp_fct_bas_dow()@{\code{bndry\_L2scp\_fct\_bas\_dow()}}}
\idx{assemblage tools!bndry_L2scp_fct_bas_dow()@{\code{bndry\_L2scp\_fct\_bas\_dow()}}}
\fdx{bndry_L2scp_fct_bas_loc_dow()@{\code{bndry\_L2scp\_fct\_bas\_loc\_dow()}}}
\idx{assemblage tools!bndry_L2scp_fct_bas_loc_dow()@{\code{bndry\_L2scp\_fct\_bas\_loc\_dow()}}}
\fdx{trace_L2scp_fct_bas()@{\code{trace\_L2scp\_fct\_bas()}}}
\idx{evaluation tools!trace_L2scp_fct_bas()@{\code{trace\_L2scp\_fct\_bas()}}}
\fdx{trace_L2scp_fct_bas_loc()@{\code{trace\_L2scp\_fct\_bas\_loc()}}}
\idx{evaluation tools!trace_L2scp_fct_bas_loc()@{\code{trace\_L2scp\_fct\_bas\_loc()}}}
\fdx{trace_L2scp_fct_bas_dow()@{\code{trace\_L2scp\_fct\_bas\_dow()}}}
\idx{evaluation tools!trace_L2scp_fct_bas_dow()@{\code{trace\_L2scp\_fct\_bas\_dow()}}}
\fdx{trace_L2scp_fct_bas_loc_dow()@{\code{trace\_L2scp\_fct\_bas\_loc\_dow()}}}
\idx{evaluation tools!trace_L2scp_fct_bas_loc_dow()@{\code{trace\_L2scp\_fct\_bas\_loc\_dow()}}}
%%
\bv\begin{lstlisting}[label=F:BNDRY_L2SCP_PROTO]
bool bndry_L2scp_fct_bas_loc(DOF_REAL_VEC *fh,
                             LOC_FCT_AT_QP f_at_qp, void *ud, FLAGS fill_flag,
                             const BNDRY_FLAGS bndry_seg,
                             const WALL_QUAD *quad);
bool bndry_L2scp_fct_bas_loc_dow(DOF_REAL_VEC_D *fh, LOC_FCT_D_AT_QP f_at_qp,
                                 void *ud, FLAGS fill_flag,
                                 const BNDRY_FLAGS bndry_seg,
                                 const WALL_QUAD *quad);
bool bndry_L2scp_fct_bas_dow(DOF_REAL_VEC_D *fh,
                             const REAL *(*f)(const REAL_D x,
                                              const REAL_D normal,
                                              REAL_D result),
                             const BNDRY_FLAGS bndry_seg,
                             const WALL_QUAD *quad);
bool bndry_L2scp_fct_bas(DOF_REAL_VEC *fh, 
                         REAL (*f)(const REAL_D x, const REAL_D normal),
                         const BNDRY_FLAGS bndry_seg, const WALL_QUAD *quad);

void trace_L2scp_fct_bas(DOF_REAL_VEC *fh, FCT_AT_X f,
			 MESH *trace_mesh, const QUAD *quad);
void trace_L2scp_fct_bas_loc(DOF_REAL_VEC *fh,
			     LOC_FCT_AT_QP f, void *fd, FLAGS fill_flag,
			     MESH *trace_mesh,
			     const QUAD *quad);
void trace_L2scp_fct_bas_dow(DOF_REAL_VEC_D *fh, FCT_D_AT_X f,
			     MESH *trace_mesh,
			     const QUAD *quad);
void trace_L2scp_fct_bas_loc_dow(DOF_REAL_VEC_D *fh,
				 LOC_FCT_D_AT_QP f, void *fd, FLAGS fill_flag,
				 MESH *trace_mesh,
				 const QUAD *quad);
\end{lstlisting}\ev

\paragraph{Descriptions}

\begin{descr}
  \kitem{bndry\_L2scp\_fct\_bas()}~\hfill
  \begin{description}
  \item[Parameters]~\hfill
    \kitem{fh} The load-vector to add the boundary integrals to.
    %% 
    \kitem{f} Application supplied ``right hand side''.
    %% 
    \kitem{ud} Data pointer for \code{f} for the
    \code{\dots\_loc()}-variants.
    %% 
    \kitem{fill\_flags} Additional fill-flags needed by \code{f} for the 
    \code{\dots\_loc()}-variants.
    %% 
    \kitem{bndry\_seg} A bit-mask specifying the part of the boundary
    which is the domain of integration. See \secref{S:boundary}.
    %% 
    \kitem{quad} Optional application supplied quadrature rule. Maybe
    \nil, in which case a default quadrature rule is used. See
    \secref{S:WALL_QUAD} for how to obtain such a structure, function
    \code{get\_wall\_quad()}.
  \item[Return Value]~\hfill

    \code{true} if part of the boundary did
    \emph{not} belong to the domain of integration.
  \end{description}

  \hrulefill

  \kitem{trace\_L2scp\_fct\_bas()}~\hfill
  \begin{description}
  \item[Parameters]~\hfill
    \begin{descr}
      \kitem{fh} The load-vector.
      %%
      \kitem{f} The user-supplied inhomogeneity.
      %%
      \kitem{fd} The application data pointer passed on to \code{f}
      for the \code{\dots\_loc()}-variants.
      %%
      \kitem{fill\_flags} Additional fill-flags for the
      \code{\dots\_loc()}-variants.
      %%
      \kitem{trace\_mesh} The domain of integration.
      %%
      \kitem{quad} A user supplied quadrature rule. May be \nil in
      which case a default quadrature rule will be used. The
      quadrature rule must have the dimension of \code{trace\_mesh},
      naturally.
      %%
    \end{descr}
  \item[Return Value]~\hfill
  \end{description}
\end{descr}
%Description:
%\begin{descr}
%\kitem{bndry\_L2scp\_fct\_bas\_loc(fh, f\_at\_qp, ud, fill\_flag, bndry\_seg, quad)}
%\fdx{bndry_L2scp_fct_bas_loc()@{\code{bndry\_L2scp\_fct\_bas\_loc()}}}
%\idx{assemblage tools!bndry_L2scp_fct_bas_loc()@{\code{bndry\_L2scp\_fct\_bas\_loc()}}}
%%%
%\kitem{bndry\_L2scp\_fct\_bas\_loc\_dow(fh, f\_at\_qp, ud, fill\_flag, bndry\_seg, quad)}
%\fdx{bndry_L2scp_fct_bas_loc_dow()@{\code{bndry\_L2scp\_fct\_bas\_loc\_dow()}}}
%\idx{assemblage tools!bndry_L2scp_fct_bas_loc_dow()@{\code{bndry\_L2scp\_fct\_bas\_loc\_dow()}}}


%\kitem{bndry\_L2scp\_fct\_bas\_dow(fh, f, bndry\_seg, quad)}
%\fdx{bndry_L2scp_fct_bas_dow()@{\code{bndry\_L2scp\_fct\_bas\_dow()}}}
%\idx{assemblage tools!bndry_L2scp_fct_bas_dow()@{\code{bndry\_L2scp\_fct\_bas\_dow()}}}
%%%
%\kitem{bndry\_L2scp\_fct\_bas(fh, f, bndry\_seg, quad)}
%\fdx{bndry_L2scp_fct_bas()@{\code{bndry\_L2scp\_fct\_bas()}}}
%\idx{assemblage tools!bndry_L2scp_fct_bas()@{\code{bndry\_L2scp\_fct\_bas()}}}

%\kitem{trace\_L2scp\_fct\_bas(fh, f, trace\_mesh, quad)}
%\fdx{trace_L2scp_fct_bas()@{\code{trace\_L2scp\_fct\_bas()}}}
%\idx{evaluation tools!trace_L2scp_fct_bas()@{\code{trace\_L2scp\_fct\_bas()}}}

%\kitem{trace\_L2scp\_fct\_bas\_loc(fh, f, fd, fill\_flag, trace\_mesh, quad)}
%\fdx{trace_L2scp_fct_bas_loc()@{\code{trace\_L2scp\_fct\_bas\_loc()}}}
%\idx{evaluation tools!trace_L2scp_fct_bas_loc()@{\code{trace\_L2scp\_fct\_bas\_loc()}}}

%\kitem{trace\_L2scp\_fct\_bas\_dow(fh, f, trace\_mesh, quad)}
%\fdx{trace_L2scp_fct_bas_dow()@{\code{trace\_L2scp\_fct\_bas\_dow()}}}
%\idx{evaluation tools!trace_L2scp_fct_bas_dow()@{\code{trace\_L2scp\_fct\_bas\_dow()}}}

%\kitem{trace\_L2scp\_fct\_bas\_loc\_dow(fh, f, fd, fill\_flag, trace\_mesh, quad)}
%\fdx{trace_L2scp_fct_bas_loc_dow()@{\code{trace\_L2scp\_fct\_bas\_loc\_dow()}}}
%\idx{evaluation tools!trace_L2scp_fct_bas_loc_dow()@{\code{trace\_L2scp\_fct\_bas\_loc\_dow()}}}
%\end{descr}

\subsubsection{Robin boundary conditions}
\label{S:robin_bound}

\fdx{robin_bound()@{\code{robin\_bound()}}}
\idx{assemblage tools!robin_bound()@{\code{robin\_bound()}}}
\bv\begin{lstlisting}
void robin_bound(DOF_MATRIX *matrix, const BNDRY_FLAGS robin_seg,
                 REAL alpha_r, const WALL_QUAD *wall_quad,
                 REAL exponent);
\end{lstlisting}\ev
Description: Incorporate so-called ``Robin boundary'' conditions into
the matrix, i.e. a boundary condition of the form
\[
\frac{\partial u}{\partial\nu}(x) + \alpha(x)\,u(x) = g(x)\quad\text{ on }\partial\Omega.
\]
In the context of weak formulations for elliptic second-order PDEs,
this results into two boundary integrals, one has to be added to the
linear form on the ``right hand side'', and the other one is a
contribution to bilinear-form on the ``left hand side'', namely
\[
\int_{\partial\Omega} \alpha\,u\,\phi\,do
\quad\text{ and }\quad
\int_{\partial\Omega}g\,\phi\,do
\]
The contribution to the right hand side can be assembled using one of
the \code{bndry\_L2scp\_fct\_bas()}-variants, the contribution the
left hand side should be added to the system matrix.
\code{robin\_bound()} implements this for the case where $\alpha$ is
actually just a constant coefficient.
\begin{descr}
\kitem{robin\_bound(matrix, robin\_seg, alpha\_r, wall\_quad, exponent)}~\hfill
\begin{description}
\item[Parameters]~\hfill
  \begin{descr}
  \kitem{matrix} The system matrix, the contributions from the Robin
    boundary condition are added to \code{matrix}.
    %% 
  \kitem{robin\_segment} A boundary bit-mask, marking all boundary
    segments which are subject to the Robin boundary condition. The
    position of the bits set in \code{robin\_segment} correspond to
    the boundary numbers assigned to the mesh boundary in the macro
    triangulation, compare \secref{S:macro_tria} and
    \secref{S:boundary}.
    %% 
  \kitem{alpha\_r} The constant coefficient from the Robin boundary
    condition.
    %% 
  \kitem{wall\_quad} Optional. A collection of co-dimension $1$
    quadrature formulas for doing the integration. If \code{wall\_quad
      == \nil}, then \code{robin\_bound()} chooses a default
    quadrature formula, based on the polynomial degree of the
    underlying basis-functions.
    %% 
  \kitem{exponent} If \code{exponent > 0.0}, then the boundary
    integral will be weighted by the factor $h(T)^{-\code{exponent}}$,
    where $h(T)$ denotes the local mesh-width.
  \end{descr}
\end{description}
\end{descr}

\subsection{Interpolation into finite element spaces}
\label{S:I_FES}

In time dependent problems, usually the ``solve'' step in the adaptive
method for the adaptation of the initial grid is an interpolation of
initial data to the finite element space, i.e. a DOF vector is filled
with the coefficient of the interpolant. The following functions are
implemented for this task:
%%
\fdx{interpol()@{\code{interpol()}}}
\idx{assemblage tools!interpol()@{\code{interpol()}}}
\fdx{interpol_d()@{\code{interpol\_d()}}}
\idx{assemblage tools!interpol_d()@{\code{interpol\_d()}}}
\fdx{interpol_dow()@{\code{interpol\_dow()}}}
\idx{assemblage tools!interpol_dow()@{\code{interpol\_dow()}}}
%%
\fdx{interpol_loc()@{\code{interpol\_loc()}}}
\idx{assemblage tools!interpol()@{\code{interpol\_loc()}}}
\fdx{interpol_loc_d()@{\code{interpol\_loc\_d()}}}
\idx{assemblage tools!interpol_loc_d()@{\code{interpol\_loc\_d()}}}
\fdx{interpol_loc_dow()@{\code{interpol\_loc\_dow()}}}
\idx{assemblage tools!interpol_loc_dow()@{\code{interpol\_loc\_dow()}}}
%%
\bv\begin{lstlisting}
void interpol(FCT_AT_X f, DOF_REAL_VEC *fh);
void interpol_d(const REAL *(*f)(const REAL_D, REAL_D), DOF_REAL_D_VEC *fh);
void interpol_dow(FCT_D_AT_X f, DOF_REAL_VEC_D *fh);
void interpol_loc(DOF_REAL_VEC *fh,
                  LOC_FCT_AT_QP f, void *f_data, FLAGS fill_flags);
void interpol_loc_d(DOF_REAL_D_VEC *fh,
                    LOC_FCT_D_AT_QP f, void *f_data, FLAGS fill_flags);
void interpol_loc_dow(DOF_REAL_VEC_D *fh,
		      LOC_FCT_D_AT_QP f, void *f_data, FLAGS fill_flags);
\end{lstlisting}\ev
Description:
\begin{descr}
  \kitem{interpol(f, fh)} computes the coefficients of the interpolant
  of a function and stores these in a DOF vector;
  
  Interpolation is done element--wise on the leaf elements of
  \code{fh->fe\_space->mesh}; the element interpolation is done by the
  function \code{fh->fe\_space->bas\_fcts->interpol()}; the
  \code{fill\_flag} of the mesh traversal is
  \code{CALL\_LEAF\_EL|FILL\_COORDS} and the function \code{f} must
  fit to the needs of \code{fh->fe\_space->bas\_fcts->interpol()}; for
  Lagrange elements, \code{(*f)()} is evaluated for all Lagrange nodes
  on the element and has to return the value at these nodes (compare
  \secref{S:basfct_data}).

  \begin{description}
  \item[Parameters]~\hfill
    \begin{descr}
      \kitem{f} is a pointer to a function for the
      evaluation of the function to be interpolated; if \code{f} is a \nil
      pointer, all coefficients are set to \code{0.0}.
      %% 
      \kitem{fh} is a DOF vector for storing the coefficients; if
      \code{fh} is a \nil pointer, nothing is done.
    \end{descr}
  \end{description}
  %%
  \kitem{interpol\_d(fd, fhd)} computes the coefficients of the
  interpolant of a vector valued function and stores these in a DOF
  vector. This version is for the case where the underlying
  basis-functions are themselves scalars, consequently the coefficient
  vector \code{fh} has the type \code{DOF\_REAL\_D\_VEC}. Otherwise
  this function differs from the scalar counter-part only in the
  calling convention for the application supplied function \code{f},
  which is the same as for \code{dirichlet\_bound\_d()}, see also
  \exampleref{E:user_fct_d}
  %%
  \kitem{interpol\_dow(fct, uh)} same as \code{interpol\_d()}, but for
  the case where the underlying basis function are either scalar- or
  \DOW-valued and the finite-element space may optionally be a direct
  sum of finite element spaces.
  %%
  \kitem{interpol\_loc(vec, fct\_at\_qp, app\_data, fill\_flags)}
  %%
  \kitem{interpol\_loc\_d(vec, fct\_at\_qp, app\_data, fill\_flags)}
  %%
  \kitem{interpol\_loc\_dow(vec, fct\_at\_qp, app\_data, fill\_flags)}
  %%
  The \code{\dots\_loc\dots}-variants differ from the other
  \code{interpol()}-flavours only in the calling convention for the
  application supplied function and the additional \code{fill\_flags}
  argument. This has already be explained above, see also
  \exampleref{E:loc_fct_at_qp}.
\end{descr}
\idx{assemblage tools|)}

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