\section{Data structures for the hierarchical mesh}%
\label{S:hierarchical_mesh}

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

\subsection{Dimension of the mesh}%
\label{S:mesh_dimension}
\label{S:mesh_constants}

The current version of \ALBERTA supports meshes triangulated using 
$d$-dimensional simplices where $d\in\{1,2,3\}$. These are embedded in 
$\R^n$, with $n\geq d$. For most applications we have $d=n$. However, for 
finite element methods on curves ($d=1$) or surfaces
($d=2$) embedded in $\R^n$ (like mean curvature flow \cite{Dziuk:91}),
the vertex coordinates of the simplices have $n>d$ components.
There are three principal constants which affect the storage layout of
various data-types, from \albertaH:
%%
\cdx{DIM_LIMIT@{\DIMLIM}}
\cdx{DIM_OF_WORLD@{\DOW}}
\cdx{DIM_MAX@{\DIMMAX}}
\cdx{N_LAMBDA_MAX@{\BARYMAX}}
%%
\bv\begin{lstlisting}[name=MESH_DIMENSIONS,caption={[Hard-coded dimension limits]},label=C:MESH_DIMENSIONS]
/* DIM_OF_WORLD is a compile time constant, not defined in alberta.h */
#define DIM_LIMIT    3 /* limiting mesh-dimension */
#define DIM_MAX      MIN(DIM_OF_WORLD, DIM_LIMIT)
\end{lstlisting}
\ev
%%
\begin{description}
\item[\DIMLIM] Defined to the limit for the mesh-dimension. More
  than tetrahedral meshes are not supported, so this is defined to $3$.
\item[\DOW] Defined to the dimension of the ambient space, i.e. $n$ in
  the notation used above.
\item[\DIMMAX] Defined to the maximum value of the mesh-dimension,
  given the current value of \DOW.
\end{description}
%%
\paragraph{Derived dimension dependent constants}
%%
\ALBERTA provides some expressions for the number of face-simplices of
each possible co-dimension. In \ALBERTA, the name ``face'' is reserved
for the faces of tetrahedra; to denote the co-dimension $1$
face-simplex for simplices of arbitrary dimensions the name ``wall''
is used. Besides that, there are expressions for the possible number
of neighbours, the faculty of the mesh-dimension and the number of
barycentric co-ordinates of given dimension. \albertaH defines the
following generic macros:
%%
\mdx{N_VERTICES@{\code{N\_VERTICES}}!N_VERTICES()@{\code{N\_VERTICES()}}}
\mdx{N_EDGES@{\code{N\_EDGES}}!N_EDGES()@{\code{N\_EDGES()}}}
\mdx{N_WALLS@{\code{N\_WALLS}}!N_WALLS()@{\code{N\_WALLS()}}}
\mdx{N_FACES@{\code{N\_FACES}}!N_FACES()@{\code{N\_FACES()}}}
\mdx{N_NEIGH@{\code{N\_NEIGH}}!N_NEIGH()@{\code{N\_NEIGH()}}}
\mdx{N_LAMBDA@{\code{N\_LAMBDA}}!N_LAMBDA()@{\code{N\_LAMBDA()}}}
\mdx{DIM_FAC@{\code{DIM\_FAC}}!DIM_FAC()@{\code{DIM\_FAC()}}}
\bv\begin{lstlisting}[name=SPLX_COUNTS,caption={[Macros for the enumeration of
 sub-simplices]},label=C:SPLX_COUNTS]
#define N_VERTICES(DIM) ((DIM)+1)
#define N_EDGES(DIM)    ((DIM)*((DIM)+1)/2)
#define N_WALLS(DIM)    ((DIM)+1)
#define N_FACES(DIM)    (((DIM) == 3) * N_WALLS(DIM))
#define N_NEIGH(DIM)    (((DIM) != 0) * N_WALLS(DIM))
#define N_LAMBDA(DIM)   N_VERTICES(DIM)
#define DIM_FAC(DIM)    ((DIM) < 2 ? 1 : (DIM) == 2 ? 2 : 6)
\end{lstlisting}\ev
%%
\begin{description}
\renewcommand{\itemsep}{-1ex}
\item[\code{N\_VERTICES()}] number of vertices of a simplex
\item[\code{N\_EDGES()}] number of edges of a simplex
\item[\code{N\_WALLS()}] number of co-dimension $1$ face-simplices of a simplex
\item[\code{N\_FACES()}] number of co-dimension $1$ face-simplices of a simplex of dimension $3$
\item[\code{N\_NEIGH()}] possible number of neighbour elements across walls
\item[\code{N\_LAMBDA()}] number barycentric co-ordinates
\item[\code{DIM\_FAC()}] faculty of the mesh-dimension
\end{description}
%%
From these generic macros \albertaH specializes variants with the
suffixes:
%%
\cdx{N_VERTICES@{\code{N\_VERTICES}}!N_VERTICES_0D@{\code{N\_VERTICES\_0D}}}
\mdx{N_VERTICES@{\code{N\_VERTICES}}!N_VERTICES_0D@{\code{N\_VERTICES\_0D}}}
\cdx{N_EDGES@{\code{N\_EDGES}}!N_EDGES_0D@{\code{N\_EDGES\_0D}}}
\mdx{N_EDGES@{\code{N\_EDGES}}!N_EDGES_0D@{\code{N\_EDGES\_0D}}}
\cdx{N_NEIGH@{\code{N\_NEIGH}}!N_NEIGH_0D@{\code{N\_NEIGH\_0D}}}
\mdx{N_NEIGH@{\code{N\_NEIGH}}!N_NEIGH_0D@{\code{N\_NEIGH\_0D}}}
\cdx{N_FACES@{\code{N\_FACES}}!N_FACES_0D@{\code{N\_FACES\_0D}}}
\mdx{N_FACES@{\code{N\_FACES}}!N_FACES_0D@{\code{N\_FACES\_0D}}}
\cdx{N_WALLS@{\code{N\_WALLS}}!N_WALLS_0D@{\code{N\_WALLS\_0D}}}
\mdx{N_WALLS@{\code{N\_WALLS}}!N_WALLS_0D@{\code{N\_WALLS\_0D}}}
\cdx{N_LAMBDA@{\code{N\_LAMBDA}}!N_LAMBDA_0D@{\code{N\_LAMBDA\_0D}}}
\mdx{N_LAMBDA@{\code{N\_LAMBDA}}!N_LAMBDA_0D@{\code{N\_LAMBDA\_0D}}}
\cdx{DIM_FAC@{\code{DIM\_FAC}}!DIM_FAC_0D@{\code{DIM\_FAC\_0D}}}
\mdx{DIM_FAC@{\code{DIM\_FAC}}!DIM_FAC_0D@{\code{DIM\_FAC\_0D}}}
%%
\cdx{N_VERTICES@{\code{N\_VERTICES}}!N_VERTICES_1D@{\code{N\_VERTICES\_1D}}}
\mdx{N_VERTICES@{\code{N\_VERTICES}}!N_VERTICES_1D@{\code{N\_VERTICES\_1D}}}
\cdx{N_EDGES@{\code{N\_EDGES}}!N_EDGES_1D@{\code{N\_EDGES\_1D}}}
\mdx{N_EDGES@{\code{N\_EDGES}}!N_EDGES_1D@{\code{N\_EDGES\_1D}}}
\cdx{N_NEIGH@{\code{N\_NEIGH}}!N_NEIGH_1D@{\code{N\_NEIGH\_1D}}}
\mdx{N_NEIGH@{\code{N\_NEIGH}}!N_NEIGH_1D@{\code{N\_NEIGH\_1D}}}
\cdx{N_FACES@{\code{N\_FACES}}!N_FACES_1D@{\code{N\_FACES\_1D}}}
\mdx{N_FACES@{\code{N\_FACES}}!N_FACES_1D@{\code{N\_FACES\_1D}}}
\cdx{N_WALLS@{\code{N\_WALLS}}!N_WALLS_1D@{\code{N\_WALLS\_1D}}}
\mdx{N_WALLS@{\code{N\_WALLS}}!N_WALLS_1D@{\code{N\_WALLS\_1D}}}
\cdx{N_LAMBDA@{\code{N\_LAMBDA}}!N_LAMBDA_1D@{\code{N\_LAMBDA\_1D}}}
\mdx{N_LAMBDA@{\code{N\_LAMBDA}}!N_LAMBDA_1D@{\code{N\_LAMBDA\_1D}}}
\cdx{DIM_FAC@{\code{DIM\_FAC}}!DIM_FAC_1D@{\code{DIM\_FAC\_1D}}}
\mdx{DIM_FAC@{\code{DIM\_FAC}}!DIM_FAC_1D@{\code{DIM\_FAC\_1D}}}
%%
\cdx{N_VERTICES@{\code{N\_VERTICES}}!N_VERTICES_2D@{\code{N\_VERTICES\_2D}}}
\mdx{N_VERTICES@{\code{N\_VERTICES}}!N_VERTICES_2D@{\code{N\_VERTICES\_2D}}}
\cdx{N_EDGES@{\code{N\_EDGES}}!N_EDGES_2D@{\code{N\_EDGES\_2D}}}
\mdx{N_EDGES@{\code{N\_EDGES}}!N_EDGES_2D@{\code{N\_EDGES\_2D}}}
\cdx{N_NEIGH@{\code{N\_NEIGH}}!N_NEIGH_2D@{\code{N\_NEIGH\_2D}}}
\mdx{N_NEIGH@{\code{N\_NEIGH}}!N_NEIGH_2D@{\code{N\_NEIGH\_2D}}}
\cdx{N_FACES@{\code{N\_FACES}}!N_FACES_2D@{\code{N\_FACES\_2D}}}
\mdx{N_FACES@{\code{N\_FACES}}!N_FACES_2D@{\code{N\_FACES\_2D}}}
\cdx{N_WALLS@{\code{N\_WALLS}}!N_WALLS_2D@{\code{N\_WALLS\_2D}}}
\mdx{N_WALLS@{\code{N\_WALLS}}!N_WALLS_2D@{\code{N\_WALLS\_2D}}}
\cdx{N_LAMBDA@{\code{N\_LAMBDA}}!N_LAMBDA_2D@{\code{N\_LAMBDA\_2D}}}
\mdx{N_LAMBDA@{\code{N\_LAMBDA}}!N_LAMBDA_2D@{\code{N\_LAMBDA\_2D}}}
\cdx{DIM_FAC@{\code{DIM\_FAC}}!DIM_FAC_2D@{\code{DIM\_FAC\_2D}}}
\mdx{DIM_FAC@{\code{DIM\_FAC}}!DIM_FAC_2D@{\code{DIM\_FAC\_2D}}}
%%
\cdx{N_VERTICES@{\code{N\_VERTICES}}!N_VERTICES_3D@{\code{N\_VERTICES\_3D}}}
\mdx{N_VERTICES@{\code{N\_VERTICES}}!N_VERTICES_3D@{\code{N\_VERTICES\_3D}}}
\cdx{N_EDGES@{\code{N\_EDGES}}!N_EDGES_3D@{\code{N\_EDGES\_3D}}}
\mdx{N_EDGES@{\code{N\_EDGES}}!N_EDGES_3D@{\code{N\_EDGES\_3D}}}
\cdx{N_NEIGH@{\code{N\_NEIGH}}!N_NEIGH_3D@{\code{N\_NEIGH\_3D}}}
\mdx{N_NEIGH@{\code{N\_NEIGH}}!N_NEIGH_3D@{\code{N\_NEIGH\_3D}}}
\cdx{N_FACES@{\code{N\_FACES}}!N_FACES_3D@{\code{N\_FACES\_3D}}}
\mdx{N_FACES@{\code{N\_FACES}}!N_FACES_3D@{\code{N\_FACES\_3D}}}
\cdx{N_WALLS@{\code{N\_WALLS}}!N_WALLS_3D@{\code{N\_WALLS\_3D}}}
\mdx{N_WALLS@{\code{N\_WALLS}}!N_WALLS_3D@{\code{N\_WALLS\_3D}}}
\cdx{N_LAMBDA@{\code{N\_LAMBDA}}!N_LAMBDA_3D@{\code{N\_LAMBDA\_3D}}}
\mdx{N_LAMBDA@{\code{N\_LAMBDA}}!N_LAMBDA_3D@{\code{N\_LAMBDA\_3D}}}
\cdx{DIM_FAC@{\code{DIM\_FAC}}!DIM_FAC_3D@{\code{DIM\_FAC\_3D}}}
\mdx{DIM_FAC@{\code{DIM\_FAC}}!DIM_FAC_3D@{\code{DIM\_FAC\_3D}}}
%%
\cdx{N_VERTICES@{\code{N\_VERTICES}}!N_VERTICES_MAX@{\code{N\_VERTICES\_MAX}}}
\mdx{N_VERTICES@{\code{N\_VERTICES}}!N_VERTICES_MAX@{\code{N\_VERTICES\_MAX}}}
\cdx{N_EDGES@{\code{N\_EDGES}}!N_EDGES_MAX@{\code{N\_EDGES\_MAX}}}
\mdx{N_EDGES@{\code{N\_EDGES}}!N_EDGES_MAX@{\code{N\_EDGES\_MAX}}}
\cdx{N_NEIGH@{\code{N\_NEIGH}}!N_NEIGH_MAX@{\code{N\_NEIGH\_MAX}}}
\mdx{N_NEIGH@{\code{N\_NEIGH}}!N_NEIGH_MAX@{\code{N\_NEIGH\_MAX}}}
\cdx{N_FACES@{\code{N\_FACES}}!N_FACES_MAX@{\code{N\_FACES\_MAX}}}
\mdx{N_FACES@{\code{N\_FACES}}!N_FACES_MAX@{\code{N\_FACES\_MAX}}}
\cdx{N_WALLS@{\code{N\_WALLS}}!N_WALLS_MAX@{\code{N\_WALLS\_MAX}}}
\mdx{N_WALLS@{\code{N\_WALLS}}!N_WALLS_MAX@{\code{N\_WALLS\_MAX}}}
\cdx{N_LAMBDA@{\code{N\_LAMBDA}}!N_LAMBDA_MAX@{\code{N\_LAMBDA\_MAX}}}
\mdx{N_LAMBDA@{\code{N\_LAMBDA}}!N_LAMBDA_MAX@{\code{N\_LAMBDA\_MAX}}}
\cdx{DIM_FAC@{\code{DIM\_FAC}}!DIM_FAC_MAX@{\code{DIM\_FAC\_MAX}}}
\mdx{DIM_FAC@{\code{DIM\_FAC}}!DIM_FAC_MAX@{\code{DIM\_FAC\_MAX}}}
%%
\cdx{N_VERTICES@{\code{N\_VERTICES}}!N_VERTICES_LIMIT@{\code{N\_VERTICES\_LIMIT}}}
\mdx{N_VERTICES@{\code{N\_VERTICES}}!N_VERTICES_LIMIT@{\code{N\_VERTICES\_LIMIT}}}
\cdx{N_EDGES@{\code{N\_EDGES}}!N_EDGES_LIMIT@{\code{N\_EDGES\_LIMIT}}}
\mdx{N_EDGES@{\code{N\_EDGES}}!N_EDGES_LIMIT@{\code{N\_EDGES\_LIMIT}}}
\cdx{N_NEIGH@{\code{N\_NEIGH}}!N_NEIGH_LIMIT@{\code{N\_NEIGH\_LIMIT}}}
\mdx{N_NEIGH@{\code{N\_NEIGH}}!N_NEIGH_LIMIT@{\code{N\_NEIGH\_LIMIT}}}
\cdx{N_FACES@{\code{N\_FACES}}!N_FACES_LIMIT@{\code{N\_FACES\_LIMIT}}}
\mdx{N_FACES@{\code{N\_FACES}}!N_FACES_LIMIT@{\code{N\_FACES\_LIMIT}}}
\cdx{N_WALLS@{\code{N\_WALLS}}!N_WALLS_LIMIT@{\code{N\_WALLS\_LIMIT}}}
\mdx{N_WALLS@{\code{N\_WALLS}}!N_WALLS_LIMIT@{\code{N\_WALLS\_LIMIT}}}
\cdx{N_LAMBDA@{\code{N\_LAMBDA}}!N_LAMBDA_LIMIT@{\code{N\_LAMBDA\_LIMIT}}}
\mdx{N_LAMBDA@{\code{N\_LAMBDA}}!N_LAMBDA_LIMIT@{\code{N\_LAMBDA\_LIMIT}}}
\cdx{DIM_FAC@{\code{DIM\_FAC}}!DIM_FAC_LIMIT@{\code{DIM\_FAC\_LIMIT}}}
\mdx{DIM_FAC@{\code{DIM\_FAC}}!DIM_FAC_LIMIT@{\code{DIM\_FAC\_LIMIT}}}
%%
\begin{description}
\renewcommand{\itemsep}{-1ex}
\item[\code{\_MAX}] maximum value given \DOW
\item[\code{\_LIMIT}] limiting value ever supported
\item[\code{\_0D}, \code{\_1D}, \code{\_2D}, \code{\_3D}]
special value given the mesh-dimension
\end{description}
%%
For example, the \code{N\_VERTICES} macro exists with the following variants:
\bv\begin{lstlisting}[caption={[Macros for the number of vertices]}]
#define N_VERTICES_0D       1
#define N_VERTICES_1D       2
#define N_VERTICES_2D       3
#define N_VERTICES_3D       4
#define N_VERTICES_MAX      N_VERTICES(DIM_MAX)
#define N_VERTICES_LIMIT    N_VERTICES(DIM_LIMIT)
\end{lstlisting}\ev
Finally we use the following definitions describing possible positions of 
degrees of freedom on an element:
\cdx{N_NODE_TYPES@{\code{N\_NODE\_TYPES}}}
\cdx{CENTER@{\code{CENTER}}}
\cdx{VERTEX@{\code{VERTEX}}}
\cdx{EDGE@{\code{EDGE}}}
\cdx{FACE@{\code{FACE}}}
\ddx{NODE_TYPES@{\code{NODE\_TYPES}}}
\bv\begin{lstlisting}[name=NODE_TYPES,caption={[\code{NODE\_TYPES}]},label=T:NODE_TYPES]
typedef enum node_types {
  VERTEX = 0,
  CENTER,
  EDGE,
  FACE,
  N_NODE_TYPES
} NODE_TYPES;
\end{lstlisting}\ev
The symbols refer to \DOFs located at the face-simplices with the
following meanings:
\begin{description}
  \renewcommand{\itemsep}{-1ex}
\item[\code{VERTEX}] The vertex of a simplex. In $1$d the vertices are
  treated as the ``walls'' of an element.
\item[\code{CENTER}] The interior of an element. The \DOFs of discontinuous
  basis-functions, e.g., are always treated as \code{CENTER}-\DOFs.
\item[\code{EDGE}] The edges of an element. Note that $1$d simplices do
  not have edges in \ALBERTA as long as it concerns the location of
  \DOFs. So $1$d-meshes have only \code{VERTEX} and \code{CENTER} \DOFs.
\item[\code{FACE}] The faces of an element. This is reserved for $3$d
  only. Note that the co-dimension $1$ face-simplex is denoted as
  ``wall-simplex'' within \ALBERTA.
\end{description}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{The local indexing on elements}%
\label{S:local_index}

For the handling of higher order discretizations where besides
vertices DOFs can be located at edges (in 2d and 3d), faces (in 3d),
or center, we also need a local numbering for edges, and faces. Finally, a
local numbering of neighbours for handling neighbour information is
needed, used for instance in the refinement algorithm itself and for error
estimator calculation.

\idx{local numbering!edges}
\idx{refinement!local numbering!edges}
\idx{local numbering!faces}
\idx{refinement!local numbering!faces}
\idx{local numbering!neighbours}
\idx{refinement!local numbering!neighbours}
The \code{i}-th neighbour is always the element opposite the
\code{i}-th vertex. The \code{i}-th edge/face is the
edge/face opposite the \code{i}-th vertex in 2d respectively 3d; 
edges in 3d are numbered in the following way (compare
Figure~\ref{F:edge_index}):
\begin{center}
\begin{tabular}{ll}
{\bf edge 0:} between vertex \code{0} and \code{1}, &
{\bf edge 3:} between vertex \code{1} and \code{2},\\
{\bf edge 1:} between vertex \code{0} and \code{2}, &
{\bf edge 4:} between vertex \code{1} and \code{3},\\
{\bf edge 2:} between vertex \code{0} and \code{3}, &
{\bf edge 5:} between vertex \code{2} and \code{3}.
\end{tabular}
\end{center}
\begin{figure}[htbp]
\hfill\includegraphics[scale=0.5]{EPS/edge_tria}\hfill%
\includegraphics[scale=0.5]{EPS/edge_tetra}\hfill
\vspace{-2mm}
\caption{Local indices of edges/neighbours in 2d and local indices 
  of edges in 3d.}\label{F:edge_index}
\end{figure}
The data structures described in the subsequent sections are based
on this local numbering of vertices, edges, faces, and neighbours.

\subsection{BLAS-like routines for \DOW- and \BARYMAX-arrays}
%%
The term ``BLAS'' stands for ``Basic Linear Algebra Subroutines'', see
\cite{LHKK:79,DDCHH:88}.
%%
There are several vector and array data-types associated with \DOW and
\BARYMAX. The basic array types are
%%
\ddx{REAL_D@{\code{REAL\_D}}}
\ddx{REAL_B@{\code{REAL\_B}}}
\bv\begin{lstlisting}[caption={[\code{REAL\_D}, \code{REAL\_B}]},label=T:REAL_D_REAL_B]
typedef REAL          REAL_D[DIM_OF_WORLD];
typedef REAL          REAL_B[N_LAMBDA_MAX];
\end{lstlisting}\ev
%%
\begin{description}
\renewcommand{\itemsep}{0pt}
\item[\REALD] An array of the dimension of the ambient space.
\item[\REALB] An array of the size of the maximum mesh-dimension at
  given \DOW. Note that for a given mesh only the first
  \code{N\_LAMBDA(mesh->dim)} components of a \REALB-vector are
  actually used. Excess elements should be cleared to $0$.
\end{description}
%%
\mdx{INIT_BARY_?D()@{\code{INIT\_BARY\_?D()}}}
\mdx{INIT_BARY_MAX()@{\code{INIT\_BARY\_MAX()}}}
%%
To support the static initialization of \REALB-arrays there are macros
\code{INIT\_BARY\_?D()}. The definitions of these macros depend on the
values of \DIMMAX, we have the following defines in \albertaH:
\bv\begin{lstlisting}[name=INIT_BARY_?D,caption={[\code{INIT\_BARY\_?D}]},label=M:INIT_BARY]
#if DIM_MAX == 0
# define INIT_BARY_0D(a)           { 1.0 }
# define INIT_BARY_1D(a, b)        { 1.0 }
# define INIT_BARY_2D(a, b, c)     { 1.0 }
# define INIT_BARY_3D(a, b, c, d)  { 1.0 }
# define INIT_BARY_MAX(a, b, c, d) INIT_BARY_0D(a)
#elif DIM_MAX == 1
# define INIT_BARY_0D(a)           { (a), 0.0 }
# define INIT_BARY_1D(a, b)        { (a), (b) }
# define INIT_BARY_2D(a, b, c)     { (a), (b) }
# define INIT_BARY_3D(a, b, c, d)  { (a), (b) }
# define INIT_BARY_MAX(a, b, c, d) INIT_BARY_1D(a, b)
#elif DIM_MAX == 2
# define INIT_BARY_0D(a)           { (a), 0.0, 0.0 }
# define INIT_BARY_1D(a, b)        { (a), (b), 0.0 }
# define INIT_BARY_2D(a, b, c)     { (a), (b), (c) }
# define INIT_BARY_3D(a, b, c, d)  { (a), (b), (c) }
# define INIT_BARY_MAX(a, b, c, d) INIT_BARY_2D(a, b, c)
#elif DIM_MAX == 3
# define INIT_BARY_0D(a)           { (a), 0.0, 0.0, 0.0 }
# define INIT_BARY_1D(a, b)        { (a), (b), 0.0, 0.0 }
# define INIT_BARY_2D(a, b, c)     { (a), (b), (c), 0.0 }
# define INIT_BARY_3D(a, b, c, d)  { (a), (b), (c), (d) }
# define INIT_BARY_MAX(a, b, c, d) INIT_BARY_3D(a, b, c, d)
#else
# error Unsupported DIM_MAX
#endif
\end{lstlisting}\ev
%%
To have array-types for matrices like Jacobians and Hessians there is
bunch of data-types in \albertaH. The suffixes which are composed from
the two letters \code{D} and \code{B} code the ordering of the array
dimensions, e.g. a \code{REAL\_BD} is an array which's first index
ranges from $0$ to (\BARYMAX-$1$) and which's second index ranges from
$0$ to (\DOW-$1$). Currently, the following types are defined:
%%
\ddx{REAL_B@{\code{REAL\_B}}}
\ddx{REAL_BB@{\code{REAL\_BB}}}
\ddx{REAL_D@{\code{REAL\_D}}}
\ddx{REAL_DD@{\code{REAL\_DD}}}
\ddx{REAL_BD@{\code{REAL\_BD}}}
\ddx{REAL_BBD@{\code{REAL\_BBD}}}
\ddx{REAL_DDD@{\code{REAL\_DDD}}}
\ddx{REAL_BDD@{\code{REAL\_BDD}}}
\ddx{REAL_BBDD@{\code{REAL\_BBDD}}}
\ddx{REAL_DB@{\code{REAL\_DB}}}
\ddx{REAL_DBB@{\code{REAL\_DBB}}}
\ddx{REAL_BBB@{\code{REAL\_BBB}}}
\ddx{REAL_BBBB@{\code{REAL\_BBBB}}}
\ddx{REAL_DBBB@{\code{REAL\_DBBB}}}
\ddx{REAL_DBBBB@{\code{REAL\_DBBBB}}}
\ddx{REAL_BDB@{\code{REAL\_BDB}}}
\ddx{REAL_BDBB@{\code{REAL\_BDBB}}}
%%
\bv\begin{lstlisting}[name=DOWBARY_TYPES,caption={[Geometric matrix and vector types]},label=T:DOWBARY_TYPES]
typedef REAL            REAL_B[N_LAMBDA_MAX];
typedef REAL_B          REAL_BB[N_LAMBDA_MAX];
typedef REAL            REAL_D[DIM_OF_WORLD];
typedef REAL_D          REAL_DD[DIM_OF_WORLD];
typedef REAL_D          REAL_BD[N_LAMBDA_MAX];
typedef REAL_BD         REAL_BBD[N_LAMBDA_MAX];
typedef REAL_DD         REAL_DDD[DIM_OF_WORLD];
typedef REAL_DD         REAL_BDD[N_LAMBDA_MAX];
typedef REAL_BDD        REAL_BBDD[N_LAMBDA_MAX];
typedef REAL_B          REAL_DB[DIM_OF_WORLD];
typedef REAL_BB         REAL_DBB[DIM_OF_WORLD];
typedef REAL_BB         REAL_BBB[N_LAMBDA_MAX];
typedef REAL_BBB        REAL_BBBB[N_LAMBDA_MAX];
typedef REAL_BBB        REAL_DBBB[DIM_OF_WORLD];
typedef REAL_BBBB       REAL_DBBBB[DIM_OF_WORLD];
typedef REAL_DB         REAL_BDB[N_LAMBDA_MAX];
typedef REAL_DBB        REAL_BDBB[N_LAMBDA_MAX];
\end{lstlisting}\ev
%
To ease arithmetic with such vector- and matrix-types there is a
variety of inline-functions defined in \code{alberta\_inlines.h}
(\code{alberta\_inlines.h} is automatically included by \albertaH). We
describe only a selection, for the full list we refer the reader to
the header \code{alberta\_inlines.h}. Some of the following functions
are also available as matrix versions (e.g. \code{MAXEY\_DOW(a,x,y)},
\code{MSCP\_DOW(x,y)}, ...), but they aren't descripted seperately.
The prefix \code{M} means that they expect \code{REAL\_DD} matrices
instead of \code{REAL\_D} vectors. A tabular overview can be found in
Table~\ref{T:BLAS_REAL_D} and Table ~\ref{T:BLAS_REAL_DD}.
%%
\paragraph{Prototypes}
\fdx{AFFAFF_DOW()@{\code{AFFAFF\_DOW()}}}
\fdx{INVAFF_DOW()@{\code{INVAFF\_DOW()}}}
\fdx{AFFINV_DOW()@{\code{AFFINV\_DOW()}}}
\fdx{AFFINE_DOW()@{\code{AFFINE\_DOW()}}}
\fdx{COPY_DOW()@{\code{COPY\_DOW()}}}
\fdx{DIST_DOW()@{\code{DIST\_DOW()}}}
\fdx{NORM_DOW()@{\code{NORM\_DOW()}}}
\fdx{NRM2_DOW()@{\code{NRM2\_DOW()}}}
\fdx{SCAL_DOW()@{\code{SCAL\_DOW()}}}
\fdx{DST2_DOW()@{\code{DST2\_DOW()}}}
\fdx{EXPAND_DOW()@{\code{EXPAND\_DOW()}}}
\fdx{FORMAT_DOW()@{\code{FORMAT\_DOW()}}}
\fdx{SET_DOW()@{\code{SET\_DOW()}}}
\fdx{SCP_DOW()@{\code{SCP\_DOW()}}}
\fdx{GRAMSCP_DOW()@{\code{GRAMSCP\_DOW()}}}
\fdx{GRAD_DOW()@{\code{GRAD\_DOW()}}}
\fdx{GRAD_P_DOW()@{\code{GRAD\_P\_DOW()}}}
\fdx{D2_DOW()@{\code{D2\_DOW()}}}
\fdx{D2_P_DOW()@{\code{D2\_P\_DOW()}}}
\fdx{CMP_DOW()@{\code{CMP\_DOW()}}} \fdx{AX_DOW()@{\code{AX\_DOW()}}}
\fdx{AXEY_DOW()@{\code{AXEY\_DOW()}}}
\fdx{AXPY_DOW()@{\code{AXPY\_DOW()}}}
\fdx{AXPBY_DOW()@{\code{AXPBY\_DOW()}}}
\fdx{AXPBYP_DOW()@{\code{AXPBYP\_DOW()}}}
\fdx{AXPBYPCZ_DOW()@{\code{AXPBYPCZ\_DOW()}}}
\fdx{AXPBYPCZP_DOW()@{\code{AXPBYPCZP\_DOW()}}}
\fdx{MCOPY_DOW()@{\code{MCOPY\_DOW()}}}
\fdx{MDIST_DOW()@{\code{MDIST\_DOW()}}}
\fdx{MNORM_DOW()@{\code{MNORM\_DOW()}}}
\fdx{MNRM2_DOW()@{\code{MNRM2\_DOW()}}}
\fdx{MSCAL_DOW()@{\code{MSCAL\_DOW()}}}
\fdx{MDST2_DOW()@{\code{MDST2\_DOW()}}}
\fdx{MEXPAND_DOW()@{\code{MEXPAND\_DOW()}}}
\fdx{MFORMAT_DOW()@{\code{MFORMAT\_DOW()}}}
\fdx{MSET_DOW()@{\code{MSET\_DOW()}}}
\fdx{MSCP_DOW()@{\code{MSCP\_DOW()}}}
\fdx{MGRAMSCP_DOW()@{\code{MGRAMSCP\_DOW()}}}
\fdx{MGRAD_DOW()@{\code{MGRAD\_DOW()}}}
\fdx{MGRAD_P_DOW()@{\code{MGRAD\_P\_DOW()}}}
\fdx{MD2_DOW()@{\code{MD2\_DOW()}}}
\fdx{MD2_P_DOW()@{\code{MD2\_P\_DOW()}}}
\fdx{MCMP_DOW()@{\code{MCMP\_DOW()}}}
\fdx{MAX_DOW()@{\code{MAX\_DOW()}}}
\fdx{MAXEY_DOW()@{\code{MAXEY\_DOW()}}}
\fdx{MAXPY_DOW()@{\code{MAXPY\_DOW()}}}
\fdx{MAXPBY_DOW()@{\code{MAXPBY\_DOW()}}}
\fdx{MAXPBYP_DOW()@{\code{MAXPBYP\_DOW()}}}
\fdx{MAXPBYPCZ_DOW()@{\code{MAXPBYPCZ\_DOW()}}}
\fdx{MAXPBYPCZP_DOW()@{\code{MAXPBYPCZP\_DOW()}}}
\fdx{MAXTPY_DOW()@{\code{MAXTPY\_DOW()}}}
\fdx{MINVERT_DOW()@{\code{MINVERT\_DOW()}}}
\fdx{MV_DOW()@{\code{MV\_DOW()}}} \fdx{MTV_DOW()@{\code{MTV\_DOW()}}}
\fdx{MM_DOW()@{\code{MM\_DOW()}}} \fdx{MTM_DOW()@{\code{MTM\_DOW()}}}
\fdx{MMT_DOW()@{\code{MMT\_DOW()}}}
\fdx{MDET_DOW()@{\code{MDET\_DOW()}}}
\fdx{MGEMV_DOW()@{\code{MGEMV\_DOW()}}}
\fdx{MGEMTV_DOW()@{\code{MGEMTV\_DOW()}}}
\fdx{GEMV_DOW()@{\code{GEMV\_DOW()}}}
\fdx{GEMTV_DOW()@{\code{GEMTV\_DOW()}}}
\fdx{WEDGE_DOW()@{\code{WEDGE\_DOW()}}}
\fdx{DIST1_DOW()@{\code{DIST1\_DOW()}}}
\fdx{DIST8_DOW()@{\code{DIST8\_DOW()}}}
\fdx{NORM1_DOW()@{\code{NORM1\_DOW()}}}
\fdx{NORM8_DOW()@{\code{NORM8\_DOW()}}}
\fdx{NRMP_DOW()@{\code{NRMP\_DOW()}}}
\fdx{PNRMP_DOW()@{\code{PNRMP\_DOW()}}}
\fdx{SUM_DOW()@{\code{SUM\_DOW()}}}
\fdx{POW_DOW()@{\code{POW\_DOW()}}}
\fdx{SCAN_EXPAND_DOW()@{\code{SCAN\_EXPAND\_DOW()}}}
\fdx{SCAN_FORMAT_DOW()@{\code{SCAN\_FORMAT\_DOW()}}}
%%
\bv\begin{lstlisting}[name=DOW_BLAS,label=C:DOW_BLAS,caption={[BLAS-like operations for small matrices and vectors]}]
REAL SCP_DOW(const REAL_D x, const REAL_D y)
REAL GRAMSCP_DOW(const REAL_DD A, const REAL_D x, const REAL_D y)
REAL NORM_DOW(const REAL_D x)
REAL NRM2_DOW(const REAL_D x)
REAL NORM1_DOW(const REAL_D x)
REAL NORM8_DOW(const REAL_D x)
REAL NRMP_DOW(const REAL_D x, REAL p)
REAL PNRMP_DOW(const REAL_D x, REAL p)
REAL DIST_DOW(const REAL_D x, const REAL_D y)
REAL DST2_DOW(const REAL_D x, const REAL_D y)
REAL DIST1_DOW(const REAL_D x, const REAL_D y)
REAL DIST8_DOW(const REAL_D x, const REAL_D y)
REAL SUM_DOW(const REAL_D x)
REAL POW_DOW(REAL a)

REAL *SET_DOW(REAL a, REAL_D x)
REAL *COPY_DOW(const REAL_D x, REAL_D y)
REAL *SCAL_DOW(REAL a, REAL_D x)
REAL *AX_DOW(REAL a, REAL_D x)
bool CMP_DOW(REAL val, const REAL_D a)

REAL *AXEY_DOW(REAL a, const REAL_D x, REAL_D y)
REAL *AXPY_DOW(REAL a, const REAL_D x, REAL_D y)
REAL *AXPBY_DOW(REAL a, const REAL_D x,
                REAL b, const REAL_D y, REAL_D z)
REAL *AXPBYP_DOW(REAL a, const REAL_D x,
                 REAL b, const REAL_D y, REAL_D z)
REAL *AXPBYPCZ_DOW(REAL a, const REAL_D x, REAL b, const REAL_D y,
                   REAL c, const REAL_D z, REAL_D w)
REAL *AXPBYPCZP_DOW(REAL a, const REAL_D x, REAL b, const REAL_D y,
                    REAL c, const REAL_D z, REAL_D w)
REAL WEDGE_DOW(const REAL_D x, const REAL_D y)
REAL *WEDGE_DOW(const REAL_D x, const REAL_D y, REAL_D z)

EXPAND_DOW(x)
FORMAT_DOW
SCAN_EXPAND_DOW(v)
SCAN_FORMAT_DOW

REAL *GRAD_DOW(int dim, const REAL_BD Lambda, const REAL_B b_grd, REAL_D x_grd)
REAL *GRAD_P_DOW(int dim, const REAL_BD Lambda,
                 const REAL_B b_grd, REAL_D x_grd)
REAL_D *D2_DOW(int dim, const REAL_BD Lambda,
               const REAL_BB b_hesse, REAL_DD x_hesse)
REAL_D *D2_P_DOW(int dim, const REAL_BD Lambda,
                 const REAL_BB b_hesse, REAL_DD x_hesse)

REAL *MV_DOW(const REAL_DD m, const REAL_D v, REAL_D b)
REAL *MTV_DOW(const REAL_DD m, const REAL_D v, REAL_D b)
REAL *GEMV_DOW(REAL a, const REAL_DD m, const REAL_D v, REAL beta, REAL_D b)
REAL *GEMTV_DOW(REAL a, const REAL_DD m, const REAL_D v, REAL beta, REAL_D b)

REAL *AFFINE_DOW(const AFF_TRAFO *trafo, const REAL_D x, REAL_D y)
REAL *AFFINV_DOW(const AFF_TRAFO *trafo, const REAL_D x, REAL_D y)
AFF_TRAFO *AFFAFF_DOW(const AFF_TRAFO *A, const AFF_TRAFO *B, AFF_TRAFO *C)
AFF_TRAFO *INVAFF_DOW(const AFF_TRAFO *A, AFF_TRAFO *B)

REAL MSCP_DOW(const REAL_DD x, const REAL_DD y)
REAL MNORM_DOW(const REAL_DD m)
REAL MNRM2_DOW(const REAL_DD m)
REAL MDIST_DOW(const REAL_DD a, const REA_DD b)
REAL MDST2_DOW(const REAL_DD a, const REAL_DD b)

REAL_D *MSET_DOW(REAL val, REAL_DD m)
REAL_D *MCOPY_DOW(const REAL_DD x, REAL_DD y)

REAL_D *MSCAL_DOW(REAL a, REAL_DD m)
REAL_D *MAX_DOW(REAL a, REAL_DD m)
bool MCMP_DOW(REAL val, const REAL_DD a)

REAL_D *MAXEY_DOW(REAL a, const REAL_DD x, REAL_DD y)
REAL_D *MAXPY_DOW(REAL a, const REAL_DD x, REAL_DD y)
REAL_D *MAXTPY_DOW(REAL a, const REAL_DD x, REAL_DD y)
REAL_D *MAXPBY_DOW(REAL a, const REAL_DD x,
                   REAL b, const REAL_DD y, REAL_DD z)
REAL_D *MAXPBYP_DOW(REAL a, const REAL_DD x,
                    REAL b, const REAL_DD y, REAL_DD z)
REAL_D *MAXPBYPCZ_DOW(REAL a, const REAL_DD x, REAL b, const REAL_DD y,
                      REAL c, const REAL_DD z, REAL_DD w)
REAL_D *MAXPBYPCZP_DOW(REAL a, const REAL_DD x, REAL b, const REAL_DD y,
                       REAL c, const REAL_DD z, REAL_DD w)

MEXPAND_DOW(m)
MFORMAT_DOW

REAL_D *MGRAD_DOW(int dim, const REAL_BD Lambda, const REAL_DB b_grd,
                  REAL_DD x_grd)
REAL_D *MGRAD_P_DOW(int dim, const REAL_BD Lambda, const REAL_DB b_grd,
                    REAL_DD x_grd)
REAL_DD *MD2_DOW(int dim, const REAL_BD Lambda, const REAL_BB *b_hesse,
                 REAL_DDD x_hesse)
REAL_DD *MD2_P_DOW(int dim, const REAL_BD Lambda, const REAL_BB *b_hesse,
                   REAL_DDD x_hesse)

REAL_D *MINVERT_DOW(const REAL_DD m, REAL_DD mi)
REAL_D *MM_DOW(const REAL_DD a, const REAL_DD b, REAL_DD c)
REAL_D *MTM_DOW(const REAL_DD a, const REAL_DD b, REAL_DD c)
REAL_D *MMT_DOW(const REAL_DD a, const REAL_DD b, REAL_DD c)
REAL MDET_DOW(const REAL_DD m)
\end{lstlisting}\ev
%%
\paragraph{Descriptions}~\hfill\\
For a more compact presentation, refer to \tableref{T:BLAS_REAL_D} and
\ref{T:BLAS_REAL_DD}.

\begin{descr}
  %% 
  \kitem{SCP\_DOW(x, y)} returns the Euclidean scalar product of the
two vectors \code{x}, \code{y}.
  %% 
  \kitem{GRAMSCP\_DOW(A, x, y)} in case \code{A} is a spd-matrix it
returns the scalar product of the two vectors \code{x}, \code{y},
defined by \code{A}: $(x,y)_A := <x, Ay> $
  %% 
  \kitem{NORM\_DOW(x)} returns the Euclidean norm of the vector
\code{x}.
  %% 
  \kitem{NRM2\_DOW(x)} returns the Euclidean scalar product of the
vector \code{x} with itself. This means it returns the square of the
Euclidean norm of the vector \code{x}.
  %% 
  \kitem{NORM1\_DOW(x)} returns the 1-norm of the vector
\code{x}. This means it returns the sum of the absolut values of the
vector entries.
  %% 
  \kitem{NORM8\_DOW(x)} returns the infinity norm or maximum norm of
the vector \code{x}.
  %% 
  \kitem{NRMP\_DOW(x, p)} returns the p-norm of the vector \code{x}.
  %% 
  \kitem{PNRMP\_DOW(x, p)} returns the p-norm to the power of p of the
vector \code{x}.
  %% 
  \kitem{DIST\_DOW(x, y)} returns the Euclidean distance between the
two vectors \code{x}, \code{y}.
  %% 
  \kitem{DST2\_DOW(x, y)} returns the square of the Euclidean distance
between two vectors \code{x}, \code{y}.
  %% 
  \kitem{DIST1\_DOW(x, y)} returns the 1-norm of the vector
\code{(x-y)}.
  %% 
  \kitem{DIST8\_DOW(x, y)} returns the infinity norm of the vector
\code{(x-y)}.
  %% 
  \kitem{SUM\_DOW(x)} returns the sum of the vector entries of the
vector \code{x}.
  %% 
  \kitem{POW\_DOW(a)} returns \code{a} to the power of
\code{DIM\_OF\_WORLD}.

  %% 
  \kitem{SET\_DOW(a, x)} set all elements of vector \code{x} to
\code{a}. Returns \code{x}
  %% 
  \kitem{COPY\_DOW(x, y)} copies all elements of vector \code{x} to
\code{y}. Returns \code{y}.
  %% 
  \kitem{SCAL\_DOW(a, x)} scales all elements of the vector \code{x}
with \code{a}. Returns \code{x}.
  %% 
  \kitem{AX\_DOW(a, x)} scales all elements of vector \code{x} with
\code{a}. Returns \code{x}.
  %% 
  \kitem{CMP\_DOW(val, a)} returns \code{true} if all elements of the
vector \code{a} have the same value \code{val}, and it returns
\code{false} if there is any element of the vector \code{a} with value
\code{!=val}.\\

  %% 
  \kitem{AXEY\_DOW(a, x, y)} scales all elements of vector \code{x}
with \code{a} and stores the resulting vector in \code{y}. Returns
\code{y}.
  %% 
  \kitem{AXPY\_DOW(a, x, y)} scales all elements of vector \code{x}
with \code{a} and add it up to the vector \code{y}.  Returns \code{y}.
  %% 
  \kitem{AXPBY\_DOW(a, x, b, y, z)} scales all elements of vector
\code{x} with \code{a}, scales all elements of the vector \code{y}
with \code{b} and add these two vectors. The result is stored in the
vector \code{z}.  Returns \code{z}.
  %% 
  \kitem{AXPBYP\_DOW(a, x, b, y, z)} scales all elements of vector
\code{x} with \code{a}, scales all elements of vector \code{y} with
\code{b} and add these two vectors up to the vector \code{z}.  Returns
\code{z}.
  %% 
  \kitem{AXPBYPCZ\_DOW(a, x, b, y, c, z, w)} scales all elements of
vector \code{x} with \code{a}, scales all elements of vector \code{y}
with \code{b}, scales all elements of vector \code{z} with \code{c}
and add these three vectors. The result is stored in the vector
\code{w}. Returns \code{w}.
  %% 
  \kitem{AXPBYPCZP\_DOW(a, x, b, y, c, z, w)} scales all elements of
vector \code{x} with \code{a}, scales all elements of vector \code{y}
with \code{b}, scales all elements of vector \code{z} with \code{c}
and add these three vectors up to the vector \code{w}. Returns
\code{w}.
  %% 
  \kitem{WEDGE\_DOW(a, b)} for \code{DIM\_OF\_WORLD==2} returns the
product \code{a[0]*b[1]-a[1]*b[0]}.
  %% 
  \kitem{WEDGE\_DOW(a, b, r)} for \code{DIM\_OF\_WORLD==3} fills
\code{r} with the vector product $a\times b\in\R^3$. Returns \code{r}.
  %% 
  \kitem{EXPAND\_DOW(x)} returns every entry of the vector \code{x}
seperated with a comma. It is used for easier print-out of
\code{REAL\_D}. An example is stated below.
  %% 
\kitem{FORMAT\_DOW} Example for \DOW\code{ == 2}:
\bv\begin{lstlisting}[label=C:FORMATEXPAND_DOW,name=FORMATEXPAND_DOW,caption={[\code{FORMAT\_DOW}, \code{EXPAND\_DOW}]}]
printf{"text" FORMAT_DOW "more text\n", EXPAND_DOW(x));
\end{lstlisting}\ev
%%
equivalent to:
%%
\bv\begin{lstlisting}[nolol,caption={}]
printf("text" "[%10.5le, %10.5le]" "more text\n", x[0], x[1]);
\end{lstlisting}\ev
%%
\kitem{SCAN\_EXPAND\_DOW(v)}
%%
\kitem{SCAN\_FORMAT\_DOW} an example will explain both (\DOW\code{ == 2}):
%%
\bv\begin{lstlisting}[caption={[\code{SCAN\_FORMAT\_DOW}, \code{SCAN\_EXPAND\_DOW}]},label=C:SCAN_FORMAT_EXPAND_DOW,name=SCAN_FORMAT_EXPAND_DOW]
printf{"text" SCAN_FORMAT_DOW "more text\n", SCAN_EXPAND_DOW(v));
\end{lstlisting}\ev
%%
equivalent to:
%%
\bv\begin{lstlisting}[nolol,caption={}]
printf("text" "%f %f" "more text\n", &v[0], &v[1]);
\end{lstlisting}\ev
%%
\kitem{GRAD\_DOW(dim, Lambda, b\_grd, x\_grd)} convert a barycentric
gradient \code{b\_grd} to a world gradient and stores the resulting
vector in \code{x\_grd}, given the gradient of the transformation to
the reference element \code{Lambda}. Whereas \code{dim} is the
dimension of the mesh. Returns \code{x\_grd}.
%%
\kitem{GRAD\_P\_DOW(dim, Lambda, b\_grd, x\_grd)} convert a
barycentric gradient \code{b\_grd} to a world gradient and add it up
to the vector \code{x\_grd}, given the gradient of the transformation
to the reference element \code{Lambda}. Whereas \code{dim} is the
dimension of the mesh. Returns \code{x\_grd}.
%%
\kitem{D2\_DOW(dim, Lambda, b\_hesse, x\_hesse)} convert a barycentric
Hesse matrix \code{b\_hesse} to a world Hesse matrix and stores the
resulting matrix in \code{x\_hesse}, given the gradient of the
transformation to the reference element \code{Lambda}. Whereas
\code{dim} is the dimension of the mesh.  Returns \code{x\_hesse}.
%%
\kitem{D2\_P\_DOW(dim, Lambda, b\_hesse, x\_hesse)} convert a
barycentric Hesse matrix \code{b\_hesse} to a world Hesse matrix and
add it up to the matrix \code{x\_hesse}, given the gradient of the
transformation to the reference element \code{Lambda}. Whereas
\code{dim} is the dimension of the mesh.  Returns \code{x\_hesse}.

%%
\kitem{MV\_DOW(m, v, b)} calculates the matrix-vector multiplication
of the matrix \code{m} and the vector \code{v}: \code{b += m*v}.
Returns \code{b}.
%%
%%
\kitem{MTV\_DOW(m, v, b)} calculates the matrix-vector multiplication
of the transpose of matrix \code{m} and the vector \code{v}: \code{b
+= $m^t$*v}. Returns \code{b}.
%%
%%
\kitem{GEMV\_DOW(a, m, v, beta, b)} returns \code{b = beta*b +
a*(m*v)}.  Where \code{a} and \code{beta} are scalar (type
\code{REAL}), \code{m} a matrix (type \code{REAL\_DD}) and \code{v}
and \code{b} are vectors (type \code{REAL\_D}).
%%
\kitem{GEMTV\_DOW(a, m, v, beta, b)} returns \code{b = beta*b +
a*($m^t$*v)}. Where \code{a} and \code{beta} are scalar (type
\code{REAL}), \code{m} a matrix (type \code{REAL\_DD}) and \code{v}
and \code{b} are vectors (type \code{REAL\_D}).

%%
\kitem{AFFINE\_DOW(trafo, x, y)} calculates the affine transformation
between the two vectors \code{x} and \code{y} and returns the vector
\code{y}. It consists of a linear transformation (matrix-vector
multiplication with the matrix \code{trafo->M}) followed by a
translation (with the translation vector \code{trafo->t}). Adequate
formular: \code{y = trafo->M * x + trafo->t}.
%%
\kitem{AFFINV\_DOW(trafo, x, y)} applies the inverse of the affine
transformation between \code{x} and \code{y}. Returns \code{y}.
%%
\kitem{AFFAFF\_DOW(A, B, C)} returns ...  \code{A}, \code{B},
\code{C}.
%%
\kitem{INVAFF\_DOW(A, B)} returns ...  \code{A}, \code{B}.
%%
\kitem{MINVERT\_DOW(m, mi)} returns the inverted matrix \code{mi} of
the matrix \code{m}.
%%
\kitem{MM\_DOW(a, b, c)} returns the matrix matrix multiplication of
\code{a} and \code{b} and stores the resulting matrix in \code{c}.
Returns \code{c}.
%%
\kitem{MTM\_DOW(a, b, c)} returns the matrix matrix multiplication of
the transposed matrix of \code{a} and \code{b} and stores the
resulting matrix in \code{c}. Returns \code{c}.
%%
\kitem{MMT\_DOW(a, b, c)} returns the matrix matrix multiplication of
\code{a} and the transposed matrix of \code{b} and stores the
resulting matrix in \code{c}. Returns \code{c}.
%%
\kitem{MDET\_DOW(m)} returns the determinant of matrix \code{m}.
\end{descr}
%%
\begin{table}[htbp]
\fdx{DIST_DOW()@{\code{DIST\_DOW()}}}
\idx{BLAS for REAL_D@{BLAS for \code{REAL\_D}}!DIST_DOW()@{\code{DIST\_DOW()}}}
\fdx{NORM_DOW()@{\code{NORM\_DOW()}}}
\idx{BLAS for REAL_D@{BLAS for \code{REAL\_D}}!NORM_DOW()@{\code{NORM\_DOW()}}}
\fdx{NRM2_DOW()@{\code{NRM2\_DOW()}}}
\idx{BLAS for REAL_D@{BLAS for \code{REAL\_D}}!NRM2_DOW()@{\code{NRM2\_DOW()}}}
\fdx{SCAL_DOW()@{\code{SCAL\_DOW()}}}
\idx{BLAS for REAL_D@{BLAS for \code{REAL\_D}}!SCAL_DOW()@{\code{SCAL\_DOW()}}}
\fdx{DST2_DOW()@{\code{DST2\_DOW()}}}
\idx{BLAS for REAL_D@{BLAS for \code{REAL\_D}}!DST2_DOW()@{\code{DST2\_DOW()}}}
\fdx{SCP_DOW()@{\code{SCP\_DOW()}}}
\idx{BLAS for REAL_D@{BLAS for \code{REAL\_D}}!SCP_DOW()@{\code{SCP\_DOW()}}}
\fdx{GRAMSCP_DOW()@{\code{GRAMSCP\_DOW()}}}
\idx{BLAS for REAL_D@{BLAS for \code{REAL\_D}}!GRAMSCP_DOW()@{\code{GRAMSCP\_DOW()}}}
\fdx{AX_DOW()@{\code{AX\_DOW()}}}
\idx{BLAS for REAL_D@{BLAS for \code{REAL\_D}}!AX_DOW()@{\code{AX\_DOW()}}}
\fdx{AXEY_DOW()@{\code{AXEY\_DOW()}}}
\idx{BLAS for REAL_D@{BLAS for \code{REAL\_D}}!AXEY_DOW()@{\code{AXEY\_DOW()}}}
\fdx{AXPY_DOW()@{\code{AXPY\_DOW()}}}
\idx{BLAS for REAL_D@{BLAS for \code{REAL\_D}}!AXPY_DOW()@{\code{AXPY\_DOW()}}}
\fdx{AXPBY_DOW()@{\code{AXPBY\_DOW()}}}
\idx{BLAS for REAL_D@{BLAS for \code{REAL\_D}}!AXPBY_DOW()@{\code{AXPBY\_DOW()}}}
\fdx{AXPBYP_DOW()@{\code{AXPBYP\_DOW()}}}
\idx{BLAS for REAL_D@{BLAS for \code{REAL\_D}}!AXPBYP_DOW()@{\code{AXPBYP\_DOW()}}}
\fdx{AXPBYPCZ_DOW()@{\code{AXPBYPCZ\_DOW()}}}
\idx{BLAS for REAL_D@{BLAS for \code{REAL\_D}}!AXPBYPCZ_DOW()@{\code{AXPBYPCZ\_DOW()}}}
\fdx{AXPBYPCZP_DOW()@{\code{AXPBYPCZP\_DOW()}}}
\idx{BLAS for REAL_D@{BLAS for \code{REAL\_D}}!AXPBYPCZP_DOW()@{\code{AXPBYPCZP\_DOW()}}}
\fdx{MDIST_DOW()@{\code{MDIST\_DOW()}}}
\idx{BLAS for REAL_D@{BLAS for \code{REAL\_D}}!MDIST_DOW()@{\code{MDIST\_DOW()}}}
\fdx{MNORM_DOW()@{\code{MNORM\_DOW()}}}
\idx{BLAS for REAL_D@{BLAS for \code{REAL\_D}}!MNORM_DOW()@{\code{MNORM\_DOW()}}}
\fdx{MNRM2_DOW()@{\code{MNRM2\_DOW()}}}
\idx{BLAS for REAL_D@{BLAS for \code{REAL\_D}}!MNRM2_DOW()@{\code{MNRM2\_DOW()}}}
\fdx{MSCAL_DOW()@{\code{MSCAL\_DOW()}}}
\idx{BLAS for REAL_D@{BLAS for \code{REAL\_D}}!MSCAL_DOW()@{\code{MSCAL\_DOW()}}}
\fdx{MDST2_DOW()@{\code{MDST2\_DOW()}}}
\idx{BLAS for REAL_D@{BLAS for \code{REAL\_D}}!MDST2_DOW()@{\code{MDST2\_DOW()}}}
\fdx{MSCP_DOW()@{\code{MSCP\_DOW()}}}
\idx{BLAS for REAL_D@{BLAS for \code{REAL\_D}}!MSCP_DOW()@{\code{MSCP\_DOW()}}}
\fdx{MGRAMSCP_DOW()@{\code{MGRAMSCP\_DOW()}}}
\idx{BLAS for REAL_D@{BLAS for \code{REAL\_D}}!MGRAMSCP_DOW()@{\code{MGRAMSCP\_DOW()}}}
\fdx{MAX_DOW()@{\code{MAX\_DOW()}}}
\idx{BLAS for REAL_D@{BLAS for \code{REAL\_D}}!MAX_DOW()@{\code{MAX\_DOW()}}}
\fdx{MAXEY_DOW()@{\code{MAXEY\_DOW()}}}
\idx{BLAS for REAL_D@{BLAS for \code{REAL\_D}}!MAXEY_DOW()@{\code{MAXEY\_DOW()}}}
\fdx{MAXPY_DOW()@{\code{MAXPY\_DOW()}}}
\idx{BLAS for REAL_D@{BLAS for \code{REAL\_D}}!MAXPY_DOW()@{\code{MAXPY\_DOW()}}}
\fdx{MAXPBY_DOW()@{\code{MAXPBY\_DOW()}}}
\idx{BLAS for REAL_D@{BLAS for \code{REAL\_D}}!MAXPBY_DOW()@{\code{MAXPBY\_DOW()}}}
\fdx{MAXPBYP_DOW()@{\code{MAXPBYP\_DOW()}}}
\idx{BLAS for REAL_D@{BLAS for \code{REAL\_D}}!MAXPBYP_DOW()@{\code{MAXPBYP\_DOW()}}}
\fdx{MAXPBYPCZ_DOW()@{\code{MAXPBYPCZ\_DOW()}}}
\idx{BLAS for REAL_D@{BLAS for \code{REAL\_D}}!MAXPBYPCZ_DOW()@{\code{MAXPBYPCZ\_DOW()}}}
\fdx{MAXPBYPCZP_DOW()@{\code{MAXPBYPCZP\_DOW()}}}
\idx{BLAS for REAL_D@{BLAS for \code{REAL\_D}}!MAXPBYPCZP_DOW()@{\code{MAXPBYPCZP\_DOW()}}}
\fdx{MAXTPY_DOW()@{\code{MAXTPY\_DOW()}}}
\idx{BLAS for REAL_D@{BLAS for \code{REAL\_D}}!MAXTPY_DOW()@{\code{MAXTPY\_DOW()}}}
\fdx{WEDGE_DOW()@{\code{WEDGE\_DOW()}}}
\idx{BLAS for REAL_D@{BLAS for \code{REAL\_D}}!WEDGE_DOW()@{\code{WEDGE\_DOW()}}}
\fdx{DIST1_DOW()@{\code{DIST1\_DOW()}}}
\idx{BLAS for REAL_D@{BLAS for \code{REAL\_D}}!DIST1_DOW()@{\code{DIST1\_DOW()}}}
\fdx{DIST8_DOW()@{\code{DIST8\_DOW()}}}
\idx{BLAS for REAL_D@{BLAS for \code{REAL\_D}}!DIST8_DOW()@{\code{DIST8\_DOW()}}}
\fdx{NORM1_DOW()@{\code{NORM1\_DOW()}}}
\idx{BLAS for REAL_D@{BLAS for \code{REAL\_D}}!NORM1_DOW()@{\code{NORM1\_DOW()}}}
\fdx{NORM8_DOW()@{\code{NORM8\_DOW()}}}
\idx{BLAS for REAL_D@{BLAS for \code{REAL\_D}}!NORM8_DOW()@{\code{NORM8\_DOW()}}}
\fdx{NRMP_DOW()@{\code{NRMP\_DOW()}}}
\idx{BLAS for REAL_D@{BLAS for \code{REAL\_D}}!NRMP_DOW()@{\code{NRMP\_DOW()}}}
\fdx{PNRMP_DOW()@{\code{PNRMP\_DOW()}}}
\idx{BLAS for REAL_D@{BLAS for \code{REAL\_D}}!PNRMP_DOW()@{\code{PNRMP\_DOW()}}}
\fdx{SUM_DOW()@{\code{SUM\_DOW()}}}
\idx{BLAS for REAL_D@{BLAS for \code{REAL\_D}}!SUM_DOW()@{\code{SUM\_DOW()}}}
\fdx{POW_DOW()@{\code{POW\_DOW()}}}
\idx{BLAS for REAL_D@{BLAS for \code{REAL\_D}}!POW_DOW()@{\code{POW\_DOW()}}}
\begin{center}{\small
\begin{tabular}{|l|l|}
\hline
\Strut\verb|REAL SCP_DOW(const REAL_D x, const REAL_D y)| &  $ (X,Y) = \sum_{i=0}^{d-1} X_i Y_i $ \\
\hline
\Strut\verb|REAL GRAMSCP_DOW(const REAL_DD A,| & $ (X,Y)_A = \sum_{i,j=0}^{d-1} X_i A_{ij} Y_j$ \\
\Strut\verb|                 const REAL_D x, const REAL_D y)| & ~ \\
\hline
\Strut\verb|REAL NORM_DOW(const REAL_D x)| & $ \|X\|_2 = \left(\sum_{i=0}^d |X_i|^2\right)^{\frac{1}{2}}$  \\
\hline
\Strut\verb|REAL NRM2_DOW(const REAL_D x)| & $  \|X\|_2^2 =  \sum_{i=0}^{d-1} |X_i|^2 $  \\
\hline
\Strut\verb|REAL NORM1_DOW(const REAL_D x)| & $  \|X\|_1 = \sum_{i=0}^{d-1} |X_i| $  \\
\hline
\Strut\verb|REAL NORM8_DOW(const REAL_D x)| & $  \|X\|_\infty = \max_{i=0}^{d-1} |X_i| $  \\
\hline
\Strut\verb|REAL NRMP_DOW(const REAL_D x, REAL p)| & $  \|X\|_p = \left( \sum_{i=0}^{d-1} |X_i|^p\right)^{\frac{1}{p}} $  \\
\hline
\Strut\verb|REAL PNRMP_DOW(const REAL_D x, REAL p)| & $  \|X\|_p^p = \sum_{i=0}^{d-1} |X_i|^p $  \\
\hline
\Strut\verb|REAL DIST_DOW(const REAL_D x, const REAL_D y)| & $ dist = \left(\sum_{i=0}^{d-1} |X_i - Y_i|^2\right)^{\frac{1}{2}} $  \\
\hline
\Strut\verb|REAL DST2_DOW(const REAL_D x, const REAL_D y)| & $ dst2 = \sum_{i=0}^{d-1} |X_i - Y_i|^2 $  \\
\hline
\Strut\verb|REAL DIST1_DOW(const REAL_D x, const REAL_D y)| & $ dist1 = \sum_{i=0}^{d-1} |X_i - Y_i| $  \\
\hline
\Strut\verb|REAL DIST8_DOW(const REAL_D x, const REAL_D y)| & $ dist8 = \max_{i=0}^{d-1} |X_i - Y_i| $  \\
\hline
\Strut\verb|REAL SUM_DOW(const REAL_D x)| & $ sum = \sum_{i=0}^{d-1} X_i $ \\
\hline
\Strut\verb|REAL POW_DOW(REAL a)| & $ pow = a^d $ \\
\hline
\Strut\verb|REAL *SCAL_DOW(REAL a, REAL_D x)| & $ X *= a$ \\
\Strut\verb|REAL *AX_DOW(REAL a, REAL_D x)| & ~ \\
\hline
\Strut\verb|REAL *AXEY_DOW(REAL a, const REAL_D x, REAL_D y)| & $ Y = aX $\\
\hline
\Strut\verb|REAL *AXPY_DOW(REAL a, const REAL_D x, REAL_D y)| & $ Y +\!= aX$ \\
\hline
\Strut\verb|REAL *AXPBY_DOW(REAL a, const REAL_D x,| & $ Z = aX + bY $ \\
\Strut\verb|                REAL b, const REAL_D y, REAL_D z)| & ~ \\
\hline
\Strut\verb|REAL *AXPBYP_DOW(REAL a, const REAL_D x,| & $ Z +\!= aX + bY $ \\
\Strut\verb|                 REAL b, const REAL_D y, REAL_D z)| & ~ \\
\hline
\Strut\verb|REAL *AXPBYPCZ_DOW(REAL a, const REAL_D x,| & $ W = aX + bY + cZ $ \\
\Strut\verb|                   REAL b, const REAL_D y, REAL c,| & ~ \\
\Strut\verb|                   const REAL_D z, REAL_D w)| & ~ \\
\hline
\Strut\verb|REAL *AXPBYPCZP_DOW(REAL a, const REAL_D x,| & $ W +\!= aX + bY + cZ $ \\
\Strut\verb|                    REAL b, const REAL_D y, REAL c,| & ~ \\
\Strut\verb|                    const REAL_D z, REAL_D w)| & ~ \\
\hline
\Strut\verb|REAL WEDGE_DOW(const REAL_D x, const REAL_D y)| & $ X[0]*Y[1]-X[1]*Y[0] $ \\
\Strut\verb|  (for DIM_OF_WORLD == 2) | & ~ \\
\hline
\Strut\verb|REAL *WEDGE_DOW(const REAL_D x, const REAL_D y, REAL_D z)| & $ Z = X\times Y $ \\
\Strut\verb|  (for DIM_OF_WORLD == 3) | & ~ \\
\hline
\end{tabular}
}\end{center}
\caption[Implemented BLAS routines for \code{REAL\_D} vectors]{Implemented BLAS routines for \code{REAL\_D} vectors ($d = $ \code{DIM\_OF\_WORLD},
  with the prefix \code{M} for \code{REAL\_DD} matrices)}
\label{T:BLAS_REAL_D}
\end{table}

\begin{table}[htbp]
\fdx{MV_DOW()@{\code{MV\_DOW()}}}
\fdx{MTV_DOW()@{\code{MTV\_DOW()}}}
\fdx{MM_DOW()@{\code{MM\_DOW()}}}
\fdx{MTM_DOW()@{\code{MTM\_DOW()}}}
\fdx{MMT_DOW()@{\code{MMT\_DOW()}}}
\fdx{GEMV_DOW()@{\code{GEMV\_DOW()}}}
\fdx{GEMTV_DOW()@{\code{GEMTV\_DOW()}}}
\fdx{MGEMV_DOW()@{\code{MGEMV\_DOW()}}}
\fdx{MGEMTV_DOW()@{\code{MGEMTV\_DOW()}}}
\begin{center}{\small
\begin{tabular}{|l|l|}
\hline
\Strut\verb|REAL *MV_DOW(const REAL_DD m, const REAL_D v, REAL_D b)| & $ b +\!= M*v $ \\
\hline
\Strut\verb|REAL *MTV_DOW(const REAL_DD m, const REAL_D v, REAL_D b)| & $ b +\!= M^t*v $ \\
\hline
\Strut\verb|REAL *GEMV_DOW(REAL a, const REAL_DD m,| & $ b= beta*b+ a*(M*v) $ \\
\Strut\verb|               const REAL_D v, REAL beta, REAL_D b)| & ~ \\
\hline
\Strut\verb|REAL *GEMTV_DOW(REAL a, const REAL_DD m,| & $ b= beta*b+ a*(M^t*v) $ \\
\Strut\verb|                const REAL_D v, REAL beta, REAL_D b)| & ~ \\
\hline
\end{tabular}
}\end{center}
\caption{Implemented BLAS routines for matrix-vectors multiplication.}
\label{T:BLAS_REAL_DD}
\end{table}


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

\subsection{Boundary types}%
\label{S:boundary}%
\idx{boundary types}

\begin{compatibility}
  \label{CPT:boundary_types}
  Previous versions of \ALBERTA were fixing the boundary conditions --
  Dirichlet, Neumann, others -- in the macro-data file and
  data-structures. This has changed: the new scheme is to assign only
  ``street-numbers'' to boundary segments in the macro-triangulation
  and leave the interpretation to the application program. This
  section describes some of the implications of this change. Compare
  also \compatref{CPT:boundary_street_numbers} and
  \compatref{CPT:macro_el_vertex_bound}. The reader is also referred
  to the documentation for \code{dirichlet\_bound...()}
  \secref{S:dirichlet_bound}, especially to
  \exampleref{E:CLEARING_DIRICHLET_NODES}.
\end{compatibility}

In \ALBERTA boundary conditions are first of all attached to boundary
segments -- and thus to the boundary walls of an element. Boundary
segments carry ``street-numbers'' which are defined by the
macro-triangulation. At the moment $255$ different ``boundary types''
are possible, where type $0$ is reserved for interior walls.
%%
\ddx{BNDRY_TYPE@{\code{BNDRY\_TYPE}}}
\bv\begin{lstlisting}[nolol,caption={},name=BNDRY_TYPE,label=D:BNDRY_TYPE]
typedef U_CHAR BNDRY_TYPE;
\end{lstlisting}\ev
%%
Note that this is an \emph{unsigned} value, and does not carry any
information about the nature of a boundary condition (e.g. Dirichlet
versus natural versus ...) imposed on a specific boundary segment to
``close'' a specific differential equation or system of equations.

Of course, for doing the linear algebra implied by the need to solve a
discretized PDE it is often handy to assign boundary conditions to
degrees of freedom (DOFs) of the finite element spaces. For doing so
\ALBERTA uses signed characters -- \code{S\_CHAR} -- with the
convention that positive numbers flag Dirichlet boundary conditions,
negative numbers flag natural boundary conditions and $0$ indicates an
interior node. Specifically, there are three pre-defined constants
%%
\cdx{INTERIOR@{\code{INTERIOR}}}
\cdx{DIRICHLET@{\code{DIRICHLET}}}
\cdx{NEUMANN@{\code{NEUMANN}}}
\bv\begin{lstlisting}[nolol,caption={}]
#define INTERIOR      0
#define DIRICHLET     1
#define NEUMANN      -1
\end{lstlisting}\ev
%%
and some macro which may help the to make code more readable, namely
%%
\mdx{IS_INTERIOR()@{\code{IS\_INTERIOR()}}}
\mdx{IS_DIRICHLET()@{\code{IS\_DIRICHLET()}}}
\mdx{IS_NEUMANN()@{\code{IS\_NEUMANN()}}}
\bv\begin{lstlisting}[nolol,caption={}]
#define IS_NEUMANN(bound) ((bound) <= NEUMANN)
#define IS_DIRICHLET(bound) ((bound) >= DIRICHLET)
#define IS_INTERIOR(bound) ((bound) == 0)
\end{lstlisting}\ev

There are some issues for assigning boundary conditions to
\code{DOF}s. One point is that a \code{DOF} may belong to boundary
segments with differing boundary classification, e.g. vertex
\code{DOF}s in 2d and vertex and edge \code{DOF}s in 3d. To handle
this point \ALBERTA provides a boundary bit-mask data type for such
\code{DOF}s, together with some support macros:
%%
\ddx{BNDRY_FLAGS@{\code{BNDRY\_FLAGS}}}
\bv\begin{lstlisting}[name=BNDRY_FLAGS,label=D:BNDRY_FLAGS,caption={[\code{BNDRY\_FLAGS}]}]
#define N_BNDRY_TYPES 256
typedef BITS_256          BNDRY_FLAGS;

/* Some "standard" bit-field operations, meant to hide the
 * N_BNDRY_TYPES argument.
 */
#define BNDRY_FLAGS_INIT(flags)   bitfield_zap((flags), N_BNDRY_TYPES)
#define BNDRY_FLAGS_ALL(flags)    bitfield_fill((flags), N_BNDRY_TYPES)
#define BNDRY_FLAGS_CPY(to, from) bitfield_cpy((to), (from), N_BNDRY_TYPES)
#define BNDRY_FLAGS_AND(to, from) bitfield_and((to), (from), N_BNDRY_TYPES)
#define BNDRY_FLAGS_OR(to, from)  bitfield_or((to), (from), N_BNDRY_TYPES)
#define BNDRY_FLAGS_XOR(to, from) bitfield_xor((to), (from), N_BNDRY_TYPES)
#define BNDRY_FLAGS_CMP(a, b)     bitfield_cmp((a), (b), N_BNDRY_TYPES)

/* bit 0 flags boundary segments, if not set we are in the interior */
#define BNDRY_FLAGS_IS_INTERIOR(mask) (!bitfield_tst((mask), 0))

/* Set bit 0 to mark this as a boundary bit-mask. */
#define BNDRY_FLAGS_MARK_BNDRY(flags) bitfield_set((flags), INTERIOR)

/* Return TRUE if SEGMENT has BIT set _and_ BIT != 0. */
#define BNDRY_FLAGS_IS_AT_BNDRY(segment, bit)	\
  ((bit) && bitfield_tst((segment), (bit)))

/* Set a bit in the boundary-type mask. The precise meaning of BIT:
 *
 * BIT == 0: clear the boundary mask (meaning: interior node)
 * BIT >  0: set bit BIT and also bit 0 (meaning: boundary node)
 */
#define BNDRY_FLAGS_SET(flags, bit)		\
  if ((bit) != INTERIOR) {			\
    bitfield_set((flags), INTERIOR);		\
    bitfield_set((flags), (bit));		\
  } else {					\
    BNDRY_FLAGS_INIT(flags);			\
  }

/* return TRUE if SEGMENT and MASK have non-zero overlap */
#define BNDRY_FLAGS_IS_PARTOF(segment, mask)	\
  bitfield_andp((segment), (mask), 1 /* offset */, N_BNDRY_TYPES)

/* FindFirstBoundaryBit, return INTERIOR for interior nodes, otherwise the
 * number of the first bit set in MASK.
 */
#define BNDRY_FLAGS_FFBB(mask) bitfield_ffs(mask, 1 /* offset */, N_BNDRY_TYPES)
\end{lstlisting}\ev

There is also a support function which returns for a given finite
element space on a given element the boundary classification of all
local \code{DOF}s in terms of such bit-masks, namely
\code{get\_bound()}, see \secref{S:fillgetelvec}. To collect boundary
information and interprete the information returned by
\code{get\_bound()} the function \code{dirichlet\_map()} can be used.
Omitting details like the handling of direct sums of finite element
spaces its implementation looks like follows. The effect is that
\code{bound[loc\_dof]} is set either to \code{DIRICHLET} or
\code{INTERIOR}, depending on whether the input bit-mask \code{mask}
and the boundary bit-masks of the local \code{DOF}s overlap.
%%
\bv\begin{lstlisting}[name=dirichlet_map,label=C:dirichlet_map,caption={[\code{dirichlet\_map()}]}]
void dirichlet_map(EL_SCHAR_VEC *bound,
                   const EL_BNDRY_VEC *bndry_bits,
                   const BNDRY_FLAGS mask)
{
  int loc_dof;

  for (loc_dof = 0; loc_dof < bound->n_components; loc_dof++) {
    if (BNDRY_FLAGS_IS_INTERIOR(bndry_bits->vec[loc_dof])) {
      bound->vec[loc_dof] = INTERIOR;
      continue;
    }
    if (BNDRY_FLAGS_IS_PARTOF(bndry_bits->vec[loc_dof], mask)) {
      bound->vec[loc_dof] = DIRICHLET;
    } else {
      bound->vec[loc_dof] = INTERIOR;
    }
  }
}
\end{lstlisting}\ev
%%
The use of the \code{dirichlet\_map()} function is also demonstrated
in the \code{assemble()}-function in the demo-program \code{heat.c},
see \secref{S:heat_build}.
%%

Besides the single-element mapper \code{dirichlet\_map()} there is
also support for filling an entire \code{DOF\_SCHAR\_VEC} at once with
the boundary-type interpretation for a given finite element space.
This task can be performed by the function \code{dirichlet\_bound()}
(and its variants), see \secref{S:dirichlet_bound}.

Generally, many function and structures accepts an argument (or
contain a component) of type \code{BNDRY\_FLAGS} which determines on
which part of the boundary they are acting. This concerns variants of
\code{dirichlet\_bound()} (\secref{S:dirichlet_bound}), the variants
of the support functions for Neumann or Robin boundary conditions
(\secref{S:neumann_bound}, \secref{S:robin_bound}), and the residual
error-estimator support functions (\secref{S:estimator}).
Data-structures affected are \code{DOF\_MATRIX}
(\secref{S:DOF_MATRIX}), \code{EL\_MATRIX\_INFO}
(\ref{T:EL_MATRIX_INFO}), \code{EL\_VEC\_INFO} (\ref{T:EL_VEC_INFO}),
\code{OPERATOR\_INFO} (\ref{T:OPERATOR_INFO}),
\code{BNDRY\_OPERATOR\_INFO} (\ref{T:BNDRY_OPERATOR_INFO}).

\subsection{The \code{MACRO\_EL} data structure}%
\label{S:macro_element}

We now describe the macro triangulation and data type for an element
of the macro triangulation.  The macro triangulation is stored in an
array of macro elements:
%%
\ddx{MACRO_EL@{\code{MACRO\_EL}}}
\bv\begin{lstlisting}[name=MACRO_EL,label=T:MACRO_EL,caption={[\lstname]}]
#define N_BNDRY_TYPES 256
typedef U_CHAR            BNDRY_TYPE;
typedef BITS_256          BNDRY_FLAGS;
typedef struct macro_el   MACRO_EL;

struct macro_el
{
  EL          *el;
  REAL_D      *coord[N_VERTICES_MAX];

  BNDRY_TYPE  wall_bound[N_WALLS_MAX];
#if DIM_MAX > 1
  BNDRY_FLAGS vertex_bound[N_VERTICES_MAX];
#endif
#if DIM_MAX > 2
  BNDRY_FLAGS edge_bound[N_EDGES_MAX];
#endif

  NODE_PROJ   *projection[N_NEIGH_MAX + 1];

  int         index;

  MACRO_EL    *neigh[N_NEIGH_MAX];
  S_CHAR      opp_vertex[N_NEIGH_MAX];
  S_CHAR      neigh_vertices[N_NEIGH_MAX][N_VERTICES(DIM_MAX-1)];
  AFF_TRAFO   *wall_trafo[N_NEIGH_MAX];
#if DIM_MAX > 1
  BNDRY_FLAGS np_vertex_bound[N_VERTICES_MAX];
#endif
#if DIM_MAX > 2
  BNDRY_FLAGS np_edge_bound[N_EDGES_MAX];
#endif

  S_CHAR      orientation;

  U_CHAR      el_type;

  struct {
    MACRO_EL    *macro_el;
    S_CHAR      opp_vertex;
  } master;
};
\end{lstlisting}\ev
% 
Description of the individual structure components:
\begin{descr}
  %% 
\kitem{el} The root of the binary tree located at this macro
  element.
  %% 
\kitem{coord} The pointer to the world coordinates of the element's
  vertices.
  %% 
\kitem{wall\_bound} The boundary classification of the respective
  wall. $0$ means this is an interior wall, any other number between
  $1$ and $255$ is a ``street number'', the boundary classification as
  read from the macro triangulation. See also
  \compatref{CPT:boundary_street_numbers}. See also
  \secref{S:boundary}.
  %% 
\kitem{vertex\_bound} Only present for $\code{DIM\_MAX}>1$. The
  boundary classification of the given vertex.
  \begin{compatibility}
    \label{CPT:macro_el_vertex_bound}
    A vertex may belong to boundary segments with differing
    classification numbers (``street numbers''). To make this
    information accessible the \code{vertex\_bound} component is now a
    bit-mask, see also \compatref{CPT:boundary_street_numbers} and
    \secref{S:boundary}. The bit-mask has $256$ slots. If bit $i$ in
    \code{vertex\_bound[v]} is set, then vertex number $v$ is located
    on the boundary segment with classification number $i$. Bit $0$
    has a special meaning: if it s \emph{not} set, then the vertex is
    an interior vertex, in order to allow for a fast check whether the
    vertex is a boundary vertex at all.

    Macros and inline functions which simplify the handling of the
    multi-bit bit-masks \code{BNDRY\_FLAGS} are described in
    \secref{S:boundary}.
  \end{compatibility}
  %% 
\kitem{edge\_bound} Only present for $\code{DIM\_MAX}>2$. The
  boundary classification of a given edge. Compare the remarks in
  \compatref{CPT:macro_el_vertex_bound} above.
  %% 
\kitem{projection} pointers for possible projection of new nodes
  during refinement. \code{projection[1]}, if set, applies to all new
  nodes. \code{projection[1+nr]} ($0\leq nr\leq\code{N\_WALLS(dim)}$)
  applies to new nodes on specific walls and overrides
  \code{projection[0]}. For details see \secref{S:node_projections}.
  \nil pointers signify no projection for the given case.
  %% 
\kitem{index} The index of this macro element.
  %% 
\kitem{neigh} \code{neigh[i]} pointer to the macro element opposite
  the \code{i}-th local vertex; it is a pointer to \nil if the
  vertex/edges/faces opposite the \code{i}-th local vertex belongs to
  the boundary.
  %% 
\kitem{opp\_vertex} \code{opp\_vertex[i]} is undefined if
  \code{neigh[i]} is a pointer to \nil; otherwise it is the local
  index of the neighbour's vertex opposite the common
  vertex/edge/face.
  %% 
\hyperitem{MACRO_EL:neigh_vertices}{neigh\_vertices} If this is a
  periodic mesh and wall number $w$ in the macro-element is part of a
  periodic boundary, then \code{neigh\_vertices[w]} is the tuple of
  local vertex numbers in the periodic neighbour the vertices of wall
  number $w$ are mapped onto. This corresponds to the combinatoric
  face-transformations specified in the macro-triangulation file
  format (see \ref{E:MACRO_FILE_TEMPLATE}) and the \code{MACRO\_DATA}
  structure (see \ref{T:MACRO_DATA}).
  %% 
\kitem{wall\_trafo} If non-\nil, then \code{wall\_trafo[w]} is the
  geometrical face-transformation which maps the current mesh onto its
  periodic neighbour across the wall number $w$.
  %% 
\hyperitem{MACRO_EL:np_vertex_bound}{np\_vertex\_bound} Non-periodic
  version of the component \code{vertex\_bound}, see above. If the
  mesh carries a periodic structure, then it is nonetheless possible
  to use a non-periodic mesh-traversal and define non-periodic finite
  element spaces.
  %% 
\hyperitem{MACRO_EL:np_edge_bound}{np\_edge\_bound} Non-periodic
  version of the structure component \code{edge\_bound}, see above.
  See also the remarks to \code{np\_vertex\_bound} above.
  %% 
\kitem{el\_type} type of the element $\mathtt{\in [0,1,2]}$ used for
  refinement and coarsening (for the definition of the element type
  see Section~\ref{book:S:refinement_algorithm}), only 3d.
  %% 
\kitem{orientation} orientation of a tetrahedron --- depending on the
  vertex numbering, this is \code{+1} or \code{-1} (only 3d). 
  %% 
\kitem{master} In the presence of trace-meshes (aka ``sub-meshes'')
  \code{master} gives the link to the macro-element of the ambient
  ``master''-mesh containing the trace-mesh this
  \code{MACRO\_EL}-structure belongs to. The current (trace)-element
  is the wall numbered \code{master.opp\_vertex} in the ambient
  \code{master.macro\_el}. See \secref{S:tracemesh_implementation}.
\end{descr}

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

The elements of the binary trees and information that should be
present for tree elements are stored in the data structure:
\ddx{EL@{\code{EL}}}
\bv\begin{lstlisting}[name=EL,label=C:EL]
typedef struct el    EL;

struct el
{
  EL         *child[2];
  DOF        **dof;
  S_CHAR     mark;
  REAL       *new_coord;

#if ALBERTA_DEBUG
  int        index;
#endif
};
\end{lstlisting}\ev
%
The members yield following information:
\begin{descr}
\kitem{child}  pointers to the two children for interior elements of the
   tree; \code{child[0]} is a pointer to \nil for leaf elements;
   \code{child[1]} is a pointer to user data on leaf elements if 
   the user is storing data on leaf elements, otherwise
   \code{child[1]} is also a pointer to \nil for leaf elements (see
   \secref{S:leaf_data_info}).
\kitem{dof}  vector of pointers to DOFs; these pointers
   may be available for the element vertices; for the edges 
  (in 2d and 3d), for the faces (in 3d), and for the barycenter; 
   they are ordered in the following way: the
   first \code{N\_VERTICES} entries correspond to the DOFs at the
   vertices; the next one are those at the edges, if present, then
   those at the faces, if present, and finally those at the
   barycenter, if present; the offsets are defined in the \code{MESH} 
   structure (see Sections \ref{S:mesh_data_structure},
   \ref{S:refinement_routines}, \ref{S:coarsening_routines}).
\kitem{mark}   marker for refinement and coarsening: if \code{mark} is 
   positive for a leaf element this element is refined \code{mark}
   times; if it is negative for a leaf element the element may be
   coarsened \code{-mark} times; (see Sections
   \ref{S:refinement_routines}, \ref{S:coarsening_routines}).
\kitem{new\_coord} if the element has a boundary edge on a curved boundary
   this is a pointer to the coordinates of the new vertex that is
   created due to the refinement of the element, otherwise it is a
   \nil pointer; thus, coordinate information can also be produced by
   the traversal routines in the case of a curved boundary.
\kitem{index}  unique global index of the element; these indices are
   not strictly ordered and may be larger than the number of elements
   in the binary tree (the list of indices may have holes after
   coarsening); the index is available only if \code{ALBERTA\_DEBUG} is
   \true.
\end{descr}

\subsection{The \code{EL\_INFO} data structure}%
\label{el_info}
\label{S:EL_INFO}

The \code{EL\_INFO} data structure has entries for all information
which is not stored on elements explicitely, but may be generated by 
the mesh traversal routines; most entries of the \code{EL\_INFO} structure
are only filled if requested (see \secref{S:traverse}).
\begin{samepage}
\ddx{EL_INFO@{\code{EL\_INFO}}}
\bv\begin{lstlisting}[name=EL_INFO,label=T:EL_INFO]
typedef struct el_info   EL_INFO;

struct el_info
{
  MESH            *mesh;
  REAL_D          coord[N_VERTICES_MAX];
  const MACRO_EL  *macro_el;
  EL              *el;
  const EL_INFO   *parent;
  FLAGS           fill_flag;
  int             level;

  S_CHAR          macro_wall[N_WALLS_MAX];

  BNDRY_TYPE      wall_bound[N_WALLS_MAX];
  BNDRY_FLAGS     vertex_bound[N_VERTICES_MAX];
#if DIM_MAX > 2
  BNDRY_FLAGS     edge_bound[N_EDGES_MAX];
#endif

  const NODE_PROJ *active_projection;

  EL              *neigh[N_NEIGH_MAX];
  S_CHAR          opp_vertex[N_NEIGH_MAX];
  REAL_D          opp_coord[N_NEIGH_MAX];

  U_CHAR          el_type;
  S_CHAR          orientation;
  
  struct {
    EL              *el;
    int             opp_vertex;
    REAL_D          opp_coord;
    U_CHAR          el_type;
    S_CHAR          orientation;
  } master, mst_neigh;

  EL_GEOM_CACHE   el_geom_cache;
};
\end{lstlisting}\ev
\end{samepage}
The members yield the following information:
\begin{descr}
  \kitem{mesh} A pointer to the current mesh, this information is
  always present.
  %% 
  \kitem{coord} \code{coord[i]} is a \code{DIM\_OF\_WORLD} vector
  storing the Cartesian coordinates of the \code{i}-th vertex.  This
  information is only present if the component \code{fill\_flag}
  contains the flag \code{FILL\_COORDS}.
  %% 
  \kitem{macro\_el} The current element belongs to the binary tree
  located at the macro element \code{macro\_el}. This information is
  always present.
  %% 
  \kitem{el} Pointer to the current element. This information is always
  present.
  %% 
  \kitem{parent} \code{el} is a child of element \code{parent}. This
  information is always present.
  \begin{compatibility}
    In previous versions \ALBERTA, \code{parent} was just a pointer of
    type \code{EL *}, now it is a pointer to the \code{EL\_INFO}
    structure of the parent element.
  \end{compatibility}
  %% 
  \kitem{fill\_flag} Actually, the bit-wise ``or'' of multiple
  fill-flags, indicating which elements are called and which
  information should be present (see \secref{S:traverse}) in the
  \code{EL\_INFO}-structure.  \emph{Note that components of the
    \code{EL\_INFO} structure which are not flagged as valid by
    \code{fill\_flag} need not be initialized and may contain random
    data.}
  %% 
  \kitem{level} level of the current element; the level is zero for
  macro elements and the level of the children is (level of the parent
  + 1); the level is always filled by the traversal routines.
  %% 
  \kitem{macro\_wall} \code{macro\_wall[nr]} contains the number of
  the wall in the ambient macro-element the wall numbered \code{nr} of
  the current element is located at, or $-1$ if that wall is an
  interior wall (with respect to the ambient macro element). This
  piece of information is only present when \code{fill\_flag} contains
  the flag \code{FILL\_MACRO\_WALLS}.
  %% 
  \kitem{wall\_bound} The boundary classification of the walls of the
  current element. See also \compatref{CPT:boundary_types}. This piece
  of information is only valid if \code{fill\_flag} contains the flag
  \code{FILL\_BOUND}. \emph{Note that is not necessary to request
  \code{FILL\_BOUND} to access the boundary classification of the
  walls of the current element; this is done more efficiently by
  requesting \code{FILL\_MACRO\_WALLS} and then calling the function
  \code{wall\_bound(el\_info, wall)}.}
  \begin{compatibility}
    In previous versions of \ALBERTA the \code{EL\_INFO} structure
    also optionally contained the boundary classification of
    ``walls'', but using the names \code{vertex\_bound} for 1d-meshes,
    \code{edge\_bound} for 2d meshes and \code{face\_bound} for 3d
    meshes. As this was extremely unhandy a new name ``wall'' was
    introduced to refer to co-dimension $1$ simplices (the name
    ``face'' was unluckily already occupied and ``defined'' to refer
    to faces of tetrahedra in 3d).
  \end{compatibility}
  %%
  \kitem{vertex\_bound} Boundary classification of the vertices.  This
  piece of information is only valid if \code{fill\_flag} contains the
  flag \code{FILL\_BOUND}.
  \begin{compatibility}
    This is now a bit-field of type \code{BNDRY\_FLAGS}. See also
    \compatref{CPT:boundary_types}.
  \end{compatibility}
  %%
  \kitem{edge\_bound} Boundary classification of the edges ($d>1$).
  This piece of information is only valid if \code{fill\_flag}
  contains the flag \code{FILL\_BOUND}.
  \begin{compatibility}
    This is now a bit-field of type \code{BNDRY\_FLAGS}. See also
    \compatref{CPT:boundary_types}.
  \end{compatibility}
  %%
  \kitem{active\_projection} If not \nil, a pointer to the projection
  function which is used to project the newly created vertex during
  refinement.
  %%
  \kitem{neigh} \code{neigh[i]} pointer to the element opposite the
  \code{i}-th local vertex; it is a pointer to \nil if the wall
  opposite the \code{i}-th local vertex belongs to the boundary. This
  piece of information is only present if \code{fill\_flag} contains
  the flag \code{FILL\_NEIGH}.
  %%
  \kitem{opp\_vertex} \code{opp\_vertex[i]} is undefined if
  \code{neigh[i]} is a pointer to \nil; otherwise it is the local
  index of the neighbour's vertex opposite the common wall. This
  piece of information is only present if \code{fill\_flag} contains
  the flag \code{FILL\_NEIGH}.
  %% 
  \kitem{opp\_coord} \code{opp\_coord[i]} coordinates of the
  \code{i-th} neighbour's vertex opposite the common wall.  This piece
  of information is only present if \code{fill\_flag} contains the
  flag \code{FILL\_OPP\_COORDS}.
  %% 
  \kitem{el\_type} The element's type (see
  Section~\ref{S:refinement_routines}); is filled automatically by the
  traversal routines (only 3d).
  %%
  \kitem{orientation} $\pm 1$: sign of the determinant of the transformation
  to the reference element\idx{reference element} with vertices
  $(0,0,0)$, $(1,1,1)$, $(1,1,0)$, $(1,0,0)$ (see 
  Figure~\ref{book:F:standard_elements}).
  %%
  \kitem{master} If the current element belongs to a co-dimension $1$
  trace-mesh (aka ``slave-mesh'', ``sub-mesh'') then this
  data-structure contains information concerning the element of the
  master-mesh the current element belongs to. This piece of
  information is only valid if \code{fill\_flag} contains the flag
  \code{FILL\_MASTER\_INFO}.
  \begin{descr}
    \kitem{el} Always filled with \code{FILL\_MASTER\_INFO}.
    %%
    \kitem{opp\_vertex} Always filled with \code{FILL\_MASTER\_INFO}.
    %%
    \kitem{opp\_coord} Only filled if \code{FILL\_COORD} is also set
    in \code{fill\_flag}
    %%
    \kitem{el\_type} Always filled with \code{FILL\_MASTER\_INFO}, if
    the master-mesh is 3d.
    %%
    \kitem{orientation} Always filled with \code{FILL\_MASTER\_INFO}, if
    the master-mesh is 3d.
  \end{descr}
  %%
  \kitem{mst\_neigh} Same information as \code{master}, but for
  neighbour across the slave element. Only filled if \code{fill\_flag}
  contains \code{FILL\_MASTER\_NEIGH}.
  %%
  \kitem{el\_geom\_cache} A storage area which is used to cache
  various geometric quantities of the current element, like the
  determinant of the transformation to the reference element, the
  normals of the walls of te current element. The data should only be
  accessed through the function \code{fill\_el\_geom\_cache(el\_info,
    fill\_flags)}, see \secref{S:fill_el_geom_cache}.
\end{descr}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Caching of geometric element quantities}
\label{S:fill_el_geom_cache}

Often it would be useful to share data like the determinant of the
transformation to the reference element or the derivative of that
transformation between pieces of program-code which are separated by
call-hierarchies, or maybe one simply does not want to blow-up the
parameter lists of application provided function hooks. To this aim
there exists a caching mechanism, called \code{EL\_GEOM\_CACHE}, which
should only be accessed and is filled by calls to
\code{fill\_el\_geom\_cache()}. The reader is also referred to the
documentation of \code{fill\_quad\_el\_cache()},
\secref{S:fill_quad_el_cache}, especially in the context of parametric
meshes of higher polynomial order. \exampleref{E:fill_quad_el_cache}
contains a simplistic example for both, \code{fill\_el\_geom\_cache()}
and \code{fill\_quad\_el\_cache()}. The element-cache structure and
the related definitions and proto-types are as follows:
%%
\idx{element geometry!fill_el_geom_cache()@{\code{fill\_el\_geom\_cache()}}}
\idx{element geometry!EL_GEOM_CACHE@{\code{EL\_GEOM\_CACHE}}}
\fdx{fill_el_geom_cache()@{\code{fill\_el\_geom\_cache()}}}
\ddx{EL_GEOM_CACHE@{\code{EL\_GEOM\_CACHE}}}
\bv\begin{lstlisting}[name=EL_GEOM_CACHE,label=T:EL_GEOM_CACHE]
typedef struct el_geom_cache EL_GEOM_CACHE;

struct el_geom_cache
{
  EL      *current_el;
  FLAGS   fill_flag;
  REAL    det;
  REAL_BD Lambda;
  int     orientation[N_WALLS_MAX][2];
  int     rel_orientation[N_WALLS_MAX];
  REAL    wall_det[N_WALLS_MAX];
  REAL_D  wall_normal[N_WALLS_MAX];
};

#define FILL_EL_DET    (1 << 0)
#define FILL_EL_LAMBDA (1 << 1)

#define FILL_EL_WALL_SHIFT(wall)       (2 + 4*(wall))
#define FILL_EL_WALL_MASK(wall)        (0x7 << FILL_EL_WALL_SHIFT(wall))

#define FILL_EL_WALL_DET(wall)             (1 << (FILL_EL_WALL_SHIFT(wall)+0))
#define FILL_EL_WALL_NORMAL(wall)          (1 << (FILL_EL_WALL_SHIFT(wall)+1))
#define FILL_EL_WALL_ORIENTATION(wall)     (1 << (FILL_EL_WALL_SHIFT(wall)+2))
#define FILL_EL_WALL_REL_ORIENTATION(wall) (1 << (FILL_EL_WALL_SHIFT(wall)+3))

#define FILL_EL_WALL_DETS			\
  (FILL_EL_WALL_DET(0)|FILL_EL_WALL_DET(1)|	\
   FILL_EL_WALL_DET(2)|FILL_EL_WALL_DET(3))

#define FILL_EL_WALL_NORMALS				\
  (FILL_EL_WALL_NORMAL(0)|FILL_EL_WALL_NORMAL(1)|	\
   FILL_EL_WALL_NORMAL(2)|FILL_EL_WALL_NORMAL(3))

#define FILL_EL_WALL_ORIENTATIONS				\
  (FILL_EL_WALL_ORIENTATION(0)|FILL_EL_WALL_ORIENTATION(1)|	\
   FILL_EL_WALL_ORIENTATION(2)|FILL_EL_WALL_ORIENTATION(3))

#define FILL_EL_WALL_REL_ORIENTATIONS					\
  (FILL_EL_WALL_REL_ORIENTATION(0)|FILL_EL_WALL_REL_ORIENTATION(1)|	\
   FILL_EL_WALL_REL_ORIENTATION(2)|FILL_EL_WALL_REL_ORIENTATION(3))

static inline const EL_GEOM_CACHE *
fill_el_geom_cache(const EL_INFO *el_info, FLAGS fill_flag);
\end{lstlisting}\ev

The members of \code{EL\_GEOM\_CACHE} have the following meaning:
\begin{descr}
  \kitem{current\_el} For internal use only.
  %%
  \kitem{fill\_flag} A bit-mask, bit-wise or of the fill flags
  listed above (\ref{T:EL_GEOM_CACHE}).
  %%
  \kitem{det} The determinant of the transformation to the reference
  element. Filled by \code{fill\_el\_geom\_cache(..., FILL\_EL\_DET)}.
  This is the cached value of the quantity computed by
  \code{el\_det()}, see \secref{S:bary_routines}.
  %%
  \kitem{Lambda} The derivative of the barycentric coordinates w.r.t.
  the Cartesian coordinates. Filled by
  \code{fill\_el\_geom\_cache(..., FILL\_EL\_LAMBDA)}.
  This is the cached value of the quantity computed by
  \code{el\_grd\_lambda()}, see \secref{S:bary_routines}.
  %% 
  \kitem{orientation} An (absolute) orientation of the walls of the
  current element and its neighbour. \code{orientation[wall][0]} is
  the orientation of the wall of the current element,
  \code{orientation[wall][1]} is the orientation of the same wall, but
  relative to the neighbour. Filled by
  \code{fill\_el\_geom\_cache(..., FILL\_EL\_WALL\_ORIENTATION(wall))}.
  These are the cached values of two calls to
  \code{wall\_orientation()}, see \secref{S:bary_routines}.
  %%
  \kitem{rel\_orientation} \code{rel\_orientation[wall]} is the cached
  value of \code{wall\_rel\_orientation()}, see
  \secref{S:bary_routines}.  Filled by \code{fill\_el\_geom\_cache(...,
    FILL\_EL\_WALL\_REL\_ORIENTATION(wall))}.
  %%
  \kitem{wall\_det} The cached return value of
  \code{get\_wall\_normal)()}, see \secref{S:bary_routines}.  Filled by
  \code{fill\_el\_geom\_cache(..., FILL\_EL\_WALL\_DET(wall))}.
  %%
  \kitem{wall\_det} The cached value of the quantity computed by
  \code{get\_wall\_normal)()}, see \secref{S:bary_routines}.  Filled by
  \code{fill\_el\_geom\_cache(..., FILL\_EL\_WALL\_NORMAL(wall))}.
  %%
\end{descr}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{The \code{INDEX} macro}%
\label{S:index_macro}%
\idx{element indices}

It is often very helpful --- especially during program development ---
for every element to have a unique global index.
This requires an entry in the element data structure which adds to
the needed computer memory. 

On the other hand this additional amount of computer memory may be 
a disadvantage in real applications where a big number of elements
is needed, and --- after program development --- element index information
is no longer of interest.

In the debug versions of the \ALBERTA libraries (\code{ALBERTA\_DEBUG==1}) 
an element index is available. The macro
\mdx{INDEX()@{\code{INDEX()}}}
\bv\begin{lstlisting}[name=INDEX(),label=C:INDEX]
INDEX(el)
\end{lstlisting}\ev
is defined to access element indices independently of the value of
\code{ALBERTA\_DEBUG}. If no indices are available, the macro returns 
\code{-1}.

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

\subsection{Application data on leaf elements}%
\label{S:leaf_data_info}

As mentioned in \secref{book:S:hierarchical_mesh1}, it is often necessary to
provide access to special user data which is needed only on leaf elements.
Error indicators give examples for such data.

Information for leaf elements depends strongly on the application and
so it seems not to be appropriate to define a fixed data type for
storing this information. Thus, we implemented the following general
concept: The user can define his own type for data that should be
present on leaf elements. \ALBERTA only needs the size of memory that
is required to store leaf data. During refinement and coarsening
\ALBERTA automatically allocates and deallocates memory for user data
on leaf elements. Thus, after grid modifications each leaf element
possesses a memory area which is big enough to take leaf data.

To access leaf data we must have for each leaf element a pointer to
the provided memory area. This would need an additional pointer on
leaf elements. To make the element data structure as small as possible
and in order to avoid different element types for leaf and interior
elements we ``hide'' leaf data at the pointer of the second child on
leaf elements:

By definition, a leaf element is an element without children. For a
leaf element the pointers to the first \emph{and} second child are pointers
to \nil, but since we u:se a binary tree the pointer to the second
child must be \nil if the pointer to the first child is a \nil
pointer and vice versa. Thus, only testing the first child will give
correct information whether an element is a leaf element or not, and
we do not have to use the pointer of the second child for this
test. As consequence we can use the pointer of the second child as a
pointer to the allocated area for leaf data and the user can write or
read leaf data via this pointer (using casting to a self-defined data 
type defined).

The consequence is that a pointer to the second child is only a
pointer to an element if the pointer to the first child is not a
\nil pointer. Thus testing whether an element is a leaf element
or not must only be done using the pointer to the first child. If no 
leaf data is stored on the mesh then the pointer to the second child
is also a \nil pointer for leaf elements.

Finally, the user may supply routines for transforming user data
from parent to children during refinement and for transforming user data
from children to parent during coarsening. If these routines are
not supplied, information stored for the parent or the children respectively
is lost.

Leaf data storage may be initialized only once for any given
mesh. Please note that leaf data is not stored when exporting meshes
to disk (see \secref{S:file_formats}).

The following function initializes leaf data:
\fdx{init_leaf_data()@{\code{init\_leaf\_data()}}}
\bv\begin{lstlisting}[name=init_leaf_data(),label=C:init_leaf_data]
size_t init_leaf_data(MESH *mesh, size_t size, 
                      void (*refine_leaf_data)(EL *parent, EL *child[2]),
                      void (*coarsen_leaf_data)(EL *parent, EL *child[2]));
\end{lstlisting}\ev
\begin{descr}
\kitem{mesh} pointer to the mesh on which leaf data is to be stored
\kitem{size} size of memory area for storing leaf data; 
  \ALBERTA may increase the
  size of leaf data in order to guarantee an aligned memory access.
\kitem{refine\_leaf\_data} pointer to a function for transformation
  of leaf data during refinement; first,
  \code{refine\_leaf\_data(parent, child)} transforms leaf data from 
  the parent to the two children if  \code{refine\_leaf\_data} is
  not \nil; after that leaf data of the parent is destroyed.
\kitem{coarsen\_leaf\_data} pointer to a function for transformation
  of leaf data during coarsening; 
  first,
  \code{coarsen\_leaf\_data(parent, child)} transforms leaf data from 
  the two children to the parent if  \code{refine\_leaf\_data} is
  not \nil; after that leaf data the of the children is destroyed.
\end{descr}
%
The following macros for testing leaf elements and accessing leaf data are
provided:
\mdx{IS_LEAF_EL()@{\code{IS\_LEAF\_EL()}}}
\mdx{LEAF_DATA()@{\code{LEAF\_DATA()}}}
\bv\begin{lstlisting}[name={IS_LEAF_DATA(), LEAF_DATA()},label=C:IS_LEAF_DATA]
#define IS_LEAF_EL(el) (!(el)->child[0])
#define LEAF_DATA(el)  ((void *)(el)->child[1])
\end{lstlisting}\ev
The first macro \code{IS\_LEAF\_EL(el)} is true for leaf elements
and false for elements inside the binary tree; for leaf elements,
\code{LEAF\_DATA(el)} returns a pointer to leaf data hidden at
the pointer to the second child.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{The \code{RC\_LIST\_EL} data structure}%
\label{S:rc_list_el}

For refining and coarsening we need information of the elements at the
refinement and coarsening edge (compare Sections~\ref{book:S:refinement_algorithm}
and \ref{book:S:coarsening_algorithm}). Thus, we have to collect all
elements at this edge. In 1d the patch is built from the current element
only, in 2d we have at most the current element and
its neighbour across this edge, if the edge is not part of the
boundary. In 3d we have to loop around this edge to collect all the
elements. Every element at the edge has at most two neighbours sharing
the same edge.  Defining an orientation for this edge, we can define
the {\em right} and {\em left} neighbour in 3d.

For every element at the refinement/coarsening edge we have an entry
in a vector. The elements of this vector build the refinement/coarsening
patch. In 1d the vector has length \code{1}, in 2d length \code{2}, and 
in 3d length \code{mesh->max\_no\_edge\_neigh} since this is the maximal 
number of elements sharing the same edge in the mesh \code{mesh}.
\ddx{RC_LIST_EL@{\code{RC\_LIST\_EL}}}
\bv\begin{lstlisting}[name=RC_LIST_EL,label=C:RC_LIST_EL]
typedef struct rc_list_el  RC_LIST_EL;

struct rc_list_el
{
  EL_INFO      el_info;
  int          no;
  int          flag;
  RC_LIST_EL   *neigh[2];
  int          opp_vertex[2];
};
\end{lstlisting}\ev
Information that is provided for every element in this
\code{RC\_LIST\_EL} vector:
\begin{descr}
\kitem{el\_info} information for element corresponding to this 
\code{RC\_LIST\_EL} structure. This is not a pointer since \code{EL\_INFO} 
structures are often overwritten during mesh traversal.
\kitem{no} this is the \code{no}--th entry in the vector.
\kitem{flag} only used in the coarsening module: \code{flag} is \code{true}
       if the coarsening edge of the element is the coarsening edge of
       the patch, otherwise \code{flag} is \code{false}.
\kitem{neigh} \code{neigh[0/1]} neighbour of element to the right/left
       in the orientation of the edge, or a \nil pointer in the case
       of a boundary face (only 3d).
\kitem{opp\_vertex} \code{opp\_vertex[0/1]} the opposite vertex of 
       \code{neigh[0/1]} (only 3d).
\end{descr}
This \code{RC\_LIST\_EL} vector is one argument to the interpolation 
and restriction routines for DOF vectors (see Section~\ref{S:DOF_INTERPOL}).

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

All information about a triangulation is accessible via the \code{MESH}
data structure:
\ddx{MESH@{\code{MESH}}}
\bv\begin{lstlisting}[name=MESH,label=T:MESH]
struct mesh
{
  const char   *name;

  int          dim;

  int          n_vertices;
  int          n_elements;
  int          n_hier_elements;

  int          n_edges;            /* Only used for dim > 1  */
  int          n_faces;            /* Only used for dim == 3 */
  int          max_edge_neigh;     /* Only used for dim == 3 */

  bool         is_periodic;   /* true if it is possible to define periodic*/
  int          per_n_vertices;/* DOF_ADMINS on this mesh. The per_n_...   */
  int          per_n_edges;   /* entries count the number of quantities on*/
  int          per_n_faces;   /* the periodic mesh (i.e. n_faces counts   */
                              /* periodic faces twice, n_per_faces not).  */
  AFF_TRAFO    *const*wall_trafos;
  int          n_wall_trafos;

  int          n_macro_el;
  MACRO_EL     *macro_els;

  REAL_D       bbox[2]; /* bounding box for the mesh */
  REAL_D       diam;    /* bbox[1] - bbox[0] */

  PARAMETRIC   *parametric;

  DOF_ADMIN    **dof_admin;
  int          n_dof_admin;

  int          n_dof_el;            /* sum of all dofs from all admins */
  int          n_dof[N_NODE_TYPES]; /* sum of vertex/edge/... dofs from
                                     * all admins
                                     */
  int          n_node_el;   /* number of used nodes on each element */
  int          node[N_NODE_TYPES]; /* index of first vertex/edge/... node*/

  unsigned int cookie;    /* changed on each refine/coarsen. Use
                           * this to check consistency of meshes
                           * and DOF vectors when reading from
                           * files.
                           */

  void         *mem_info; /* pointer for administration; don't touch! */
};
\end{lstlisting}\ev
%%
The members yield following information:
\begin{descr}
\kitem{name} string with a textual description for the mesh, or \nil.
  Note that \code{name} will be duplicated by calling \code{strdup(3)}
  by the \code{GET\_MESH()} call.
  %% 
\kitem{dim} dimension $d$ of the mesh. May be any number from \code{0} to 
  \code{DIM\_OF\_WORLD}. Zero dimensional meshes are simply isolated vertices 
  lacking most of the features of 1d/2d/3d meshes. They were 
  implemented for completeness.
  %% 
\kitem{n\_vertices} number of vertices of the mesh.
  %% 
\kitem{n\_elements} number of leaf elements of the mesh.
  %% 
\kitem{n\_hier\_elements} number of all elements of the mesh.
  %% 
\kitem{n\_edges} number of edges of the mesh (2d and 3d).
  %% 
\kitem{n\_faces} number of faces of the mesh (3d).
  %% 
\kitem{max\_edge\_neigh} maximal number of elements that share one
  edge; used to allocate memory to store pointers to the
  neighbour at the refinement/coarsening edge (3d).
  %% 
\kitem{is\_periodic} a boolean value, set to \code{true} for
  periodic meshes, see \secref{S:periodic}.
  %% 
\kitem{per\_n\_vertices, per\_n\_edges, per\_n\_faces} the
  respective quantities, but counted taking the periodic structure
  into account, \code{n\_faces}, e.g., counts periodic faces twice,
  \code{per\_n\_faces} not.
  %% 
\kitem{wall\_trafos, n\_wall\_trafos} The geometric face
  transformation defining the periodic structure of the mesh, see
  \secref{S:periodic}.
  %% 
\kitem{n\_macro\_el} number of macro elements.
  %% 
\kitem{macro\_els} pointer to the macro element array.
  %% 
\kitem{bbox} the bounding box of the mesh.
  %% 
\kitem{diam} diameter of the mesh in the \code{DIM\_OF\_WORLD} directions.
  %% 
\kitem{parametric} is a pointer to \nil if the mesh contains no
  parametric elements; otherwise it is a pointer to a
  \code{PARAMETRIC} structure containing coefficients of the
  parameterization and related information; for more information see
  \secref{S:parametric_meshes}.
\end{descr}
The last entries are used for the administration of DOFs and are explained
in Section \ref{S:DOFs2} in detail.
\begin{descr}
\kitem{dof\_admin}  vector of \code{dof\_admin}s.
\kitem{n\_dof\_admin}  number of \code{dof\_admin}s.
\kitem{n\_node\_el} number of nodes on a single element where DOFs are
       located; needed for the (de-) allocation of the
       \code{dof}-vector on the element.
\kitem{n\_dof\_el} number of all DOFs on a single element.
\kitem{n\_dof} number of DOFs at the different positions \code{VERTEX},
       \code{EDGE}, (\code{FACE},) \code{CENTER} on an element:
\begin{descr}
\kitem{n\_dof[VERTEX]} number of DOFs at a vertex; if no DOFs are 
       associated to the barycenter, then this value is 0.
\kitem{n\_dof[CENTER]} number of DOFs at the barycenter; if no DOFs are
       associated to the barycenter, then this value is 0.
\kitem{n\_dof[EDGE]} number of DOFs at an edge; if no DOFs are
       associated to edges, then this value is 0 (2d and 3d);
\kitem{n\_dof[FACE]} number of DOFs at a face; if no DOFs are
       associated to faces, then this value is 0 (3d);
\end{descr}
\kitem{node} gives the index of the first node at vertex, edge (2d and 3d), 
       face (3d), and barycenter:
\begin{descr}
\kitem{node[VERTEX]} always has value 0; \code{dof[0],...,dof[N\_VERTICES-1]}
 are always DOFs at the vertices, if DOFs are located at vertices.
\kitem{node[CENTER]} \code{dof[node[CENTER]]} are the DOFs at the 
  barycenter, if DOFs are located at the barycenter.
\kitem{node[EDGE]} \code{dof[node[EDGE]],..., dof[node[EDGE]+N\_EDGES-1]}
 are the DOFs at the \code{N\_EDGES} edges, if DOFs are located at edges
 (2d and 3d);
\kitem{node[FACE]} \code{dof[node[FACE]],..., dof[node[FACE]+N\_FACES-1]}
 are the DOFs at the \code{N\_FACES} faces, if DOFs are located at
 faces (3d);
\end{descr}
\end{descr}
The \code{cookie} value is automatically initialized with a random
value if \code{ALBERTA\_DEBUG==0} and with a fixed number for
\code{ALBERTA\_DEBUG==1}. It is incremented on each mesh change
(refinement or coarsening). On writing meshes or finite element
coefficient vectors to disk the current cookie value is also stored. The
purpose is to provide a safety check on reading meshes and vectors; if
the cookies do not match, then \ALBERTA issues a warning message since
no guarantee can be given that coefficient vector and mesh will match.

Finally, the pointer \code{mem\_info} is used for internal memory management
and must not be changed.

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

\subsection{Initialization of meshes}%
\label{S:mesh_initialization}%
\idx{initialization of meshes}
It is possible to handle any number of meshes of any dimension $d\leq n$ in 
a given simulation. A mesh must be allocated by the following 
function or macro
\fdx{check_and_get_mesh()@{\code{check\_and\_get\_mesh()}}}
\fdx{GET_MESH()@{\code{GET\_MESH()}}}
\mdx{GET_MESH()@{\code{GET\_MESH()}}}
\bv\begin{lstlisting}[name=GET_MESH(),label=C:GET_MESH]
check_and_get_mesh(int dim, int dow, int neigh,
		   const char *version, const char *name,
		   const MACRO_DATA *macro_data,
		   NODE_PROJ *(*init_node_proj)(MESH *, MACRO_EL *, int),
		   AFF_TRAFO *(*init_wall_trafo)(MESH *, MACRO_EL *, int wall));
#define GET_MESH(dim, name, macro_data, init_node_proj, init_wall_trafo) \
  check_and_get_mesh((dim), DIM_OF_WORLD, ALBERTA_DEBUG,		\
		     ALBERTA_VERSION, (name), (macro_data),		\
		     (init_node_proj), (init_wall_trafos))
\end{lstlisting}\ev
\paragraph{Descriptions}
\begin{descr}
  \kitem{check\_and\_get\_mesh(dim, dow, debug, version, name, macro\_data,}
  \kitem{\hspace*{22ex}init\_node\_proj, init\_wall\_trafos)}~\\
  Return a
  pointer to a filled mesh structure; several consistency checks are
  performed. The application should not change any entry in the
  returned structure.  There is no other possibility to define new
  meshes inside \ALBERTA.  The arguments \code{dow}, \code{debug} and
  \code{version} are checked against the constants in the used
  library; if these values are identical, the mesh is allocated,
  otherwise an error message is produced and the program stops.
  %%
  \begin{description}
  \item[parameters]\hfill
    \begin{description}
    \item[\code{dim}] Desired dimension of the mesh
      ($1\leq\text{\code{dim}}\leq\min\{\text{\code{DIM\_OF\_WOLRD}},\,3\}$).
    \item[\code{dow}] Must be \code{DIM\_OF\_WORLD}.
    \item[\code{debug}] Must be \code{ALBERTA\_DEBUG}.
    \item[\code{version}] Must be \code{ALBERTA\_VERSION}.
    \item[\code{name}] A string holding a textual description of mesh
      and is duplicated at the member \code{name} of the mesh.
    \item[\code{macro\_data}] A pointer to the desired macro
      triangulation, see \secref{S:macro_tria} for details.
    \item[\code{init\_node\_proj}] Optional, may be \code{NULL}. A
      pointer to a function that will perform the initialization of
      new vertex projections, see \secref{S:node_projections}.
    \item[\code{init\_wall\_trafos}] Optional, may be \code{NULL}. A
      pointer to a function which initializes face transformations in
      the context of periodic meshes.
    \end{description}
  \end{description}
  \kitem{GET\_MESH(dim, name, macro\_data, init\_node\_proj, init\_wall\_trafos)}~\\
  %%
  Return a pointer to a filled mesh structure; this macro calls
  \code{check\_and\_get\_mesh()} and automatically supplies this
  function with the three (missing) arguments; this macro should
  always be used for creation of meshes.
\end{descr}

A mesh that is not needed any more can be freed by a call of the
function
\fdx{free_mesh()@{\code{free\_mesh()}}}
\bv\begin{lstlisting}[name=free_mesh(),label=C:free_mesh()]
void free_mesh(MESH *);
\end{lstlisting}\ev
%
Description:
\begin{descr}
\kitem{free\_mesh(mesh)} will de--allocate all memory used
  by \code{mesh} (elements, DOFs, etc.), and finally the data
  structure \code{mesh} too. Submeshes of this mesh are also freed, see 
  also \secref{S:submesh_implementation}.
\end{descr}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Projection of new nodes}%
\label{S:node_projections}%
\idx{node projection}

During refinement of simplices \ALBERTA usually places the new nodes
at the midpoint of the refinement edge. Some applications require
meshes having curved boundaries parametrized by a given continuous
function. For these it is possible to automatically project new nodes
on the boundary using this function. As the mesh is refined the curved
interface is successively better approximated.
Figure~\ref{F:node_projection} illustrates some refinements of a
triangle with one edge on the curved boundary. The projections of
refinement edge midpoints (small circles) to the curved boundary are
shown by the black dots.

\begin{figure}[htbp]
\centerline{\includegraphics[scale=0.7]{EPS/param_bound}}
\caption[Refinement at curved boundary]{Refinement at curved boundary:
  refinement edge midpoints $\circ$ are projected to the curved
  boundary $\bullet$}
\label{F:node_projection}
\end{figure}

\ALBERTA implements this in a very general way. It is possible to not only 
project nodes to boundaries, but also to arbitrary interfaces  in the interior 
of the mesh. It is even possible to project \emph{all} new nodes of the mesh 
to a given surface, making it possible to triangulate parametrized embedded 
surfaces or curves.

The following type is used to describe node projection functions:
\ddx{NODE_PROJECTION@{\code{NODE\_PROJECTION}}}
\bv\begin{lstlisting}[label=T:NODE_PROJECTION]
typedef struct node_projection  NODE_PROJECTION;
struct node_projection
{
  void (*func)(REAL_D old_coord, const EL_INFO *eli, const REAL_B lambda);
};
\end{lstlisting}\ev
The component \code{func} must overwrite the given coordinate vector 
\code{old\_coord} with the projected coordinates. As an alternative to 
world coordinates, the function may use the barycentric coordinates 
\code{lambda} describing a position on the element \code{eli}. The result 
must always be returned as world coordinates in the vector \code{old\_coord},
however.

The idea is that the user provides a callback function
\code{init\_node\_proj} during mesh initialization. This function
must decide which vertices/edges/faces (for 1d/2d/3d) of which macro
elements are to belong to the parametrized interface. All nodes
belonging to the interface are automatically projected during
refinement. \ALBERTA calls \code{init\_node\_proj} several times for
each macro element and thus builds the \code{projection} entries of the
\code{MACRO\_EL} structures, see \secref{S:macro_element}.

During the allocation of a
mesh with \code{check\_and\_get\_mesh()}, see \secref{S:mesh_initialization}, 
the user may pass the function pointer
\code{init\_node\_proj}. This function has the following form:
\idx{node projection!init_node_proj()@{\code{init\_node\_proj()}}}
\bv\begin{lstlisting} 
NODE_PROJECTION *init_node_proj(MESH *mesh, MACRO_EL *mel, int case);
\end{lstlisting}\ev
Description:
\begin{descr}
\kitem{mesh} pointer to the mesh
\kitem{mel} pointer to the macro element
\kitem{case} \ALBERTA calls \code{init\_node\_proj} once with
\code{case==0} and additionally for \code{case==1} to 
\code{case==N\_NEIGH(mesh->dim)+1} if $\mathtt{dim\in\{2,3\}}$.

If \code{init\_node\_proj} returns a \code{NODE\_PROJECTION} for 
\code{case==0}, then all new nodes will be projected. If
\code{init\_node\_proj} returns a \code{NODE\_PROJECTION} for
$\mathtt{case\in\{1,\dots,N\_NEIGH(dim)+1\}}$, $\mathtt{dim\in\{2,3\}}$, 
then all new nodes on edge/face \code{case-1} will be projected. This overrides
the \code{case==0} projection, if also set. A \nil value represents no
projection.
\end{descr}

\begin{example}[Triangulation of a unit ball]
  \label{Ex:unit_ball}
The following code demonstrates the projection of boundary nodes 
to the unit sphere in any dimension.
\idx{node projection!example of node projection}
\bv\begin{lstlisting} 
static void ball_proj_func(REAL_D vertex, const EL_INFO *eli, 
                           const REAL_B lambda)
{
  REAL norm = NORM_DOW(vertex);
  
  norm = 1.0 / MAX(1.0E-15, norm);
  SCAL_DOW(norm, vertex);

  return;
}

static NODE_PROJECTION *init_node_proj(MESH *mesh, MACRO_EL *mel, int c)
{
  static NODE_PROJECTION ball_proj = {ball_proj_func};

  if(c > 0 && !mel->neigh[c-1])
    return &ball_proj;
  else
    return nil;
}
\end{lstlisting}\ev
\end{example}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Reading and writing macro triangulations}%
\label{S:macro_tria}%
\label{S:macro_triang}%
\idx{macro triangulation!reading macro triangulations}
\idx{macro triangulation!writing macro triangulations}

Data for macro triangulations can easily be stored in an
\code{ASCII}-file (for binary macro files, see the end of this
section). For the macro triangulation file we use a similar key-data
format like for the parameter initialization (see
Section~\ref{par_files}). A line containing a \code{`:'}-character
defines a key. The key consists of all characters from the start of
line up to the \code{`:'}-char, including spaces. Everything after the
colon potentially contains data, either on the same line or on the
following lines. Data following a \code{`\#'}-character is ignored,
\code{`\#'} is the comment-character. The following template lists all
possible keys with a brief description of the data format. Luckily, an
application does not need to specify all of the key-value pairs in all
cases. A simple example is given further below, see
\exampleref{E:MACRO_UNIT_INTERVAL}, \ref{E:MACRO_UNIT_SQUARE} and
\ref{E:MACRO_UNIT_CUBE} below.
%
\paragraph{Macro-file template}
%
\idx{macro triangulation!macro triangulation file}
\bv\begin{lstlisting}[name=MACRO_FILE_TEMPLATE,label=E:MACRO_FILE_TEMPLATE]
# _This_ is a comment, introduced by a hash mark
DIM:          dim
DIM_OF_WORLD: dow

number of vertices: nv
number of elements: ne

vertex coordinates:
# Comments may be mixed with data
# _This_ line and the line above are comments
<DIM_OF_WORLD coordinates of vertex[0]>
...
<DIM_OF_WORLD coordinates of vertex[nv-1]>

element vertices:
<N_VERTICES(dim) indices of vertices of simplex[0]>
...
<N_VERTICES(dim) indices of vertices of simplex[ne-1]>

element boundaries:
<N_NEIGH(dim) boundary descriptions of simplex[0]>
...
<N_NEIGH(dim) boundary descriptions of simplex[ne-1]>

element neighbours:
<N_NEIGH(dim) neighbour indices of simplex[0]>
...
<N_NEIGH(dim) neighbour indices of of simplex[ne-1]>

element type:
<element type of simplex[0]>
...
<element type of simplex[ne-1]>

number of wall transformations: <number of generators>

wall transformations:
<data for first group generator, an affine isometry in projective notation>
...
<data for last group generator, an affine isometry in projective notation>

element wall transformations:
<N_WALLS(dim) wall-transformations for simplex[0]>
...
<N_WALLS(dim) wall-transformations for simplex[ne-1]>

number of wall vertex transformations: <number of transformations>

wall vertex transformations:
<first mapping between periodic walls, identifying vertex indices>
...
<last mapping between periodic walls, identifying vertex indices>

\end{lstlisting}\ev

\paragraph{Key-value descriptions}
Data for elements and vertices are read and stored in vectors for the
macro triangulation. Index information given in the file correspond to
this vector oriented storage of data.  Thus, index information must be
in the range \code{0,...,ne-1} for elements and \code{0,...,nv-1} for
vertices. Generally, ordering of data is of little importance except
that the \code{DIM} and \code{DIM\_OF\_WORLD} keys must come first,
and that ``natural'' dependencies must be obeyed: the number of
entities (vertices, elements, etc.) has to be specified before the
data defining those entities, and data attached to entities must be
defined after defining the entities it is attached to (e.g.
neighbourhood relations have to be defined after defining the elements
of the mesh).

\begin{description}
\item[\code{DIM}] Mandatory. The mesh dimension.
%%
\item[\code{DIM\_OF\_WORLD}] Mandatory. Dimension of the ambient
  space.  The parameter \code{DIM\_OF\_WORLD} must match the libary
  value of \code{DIM\_OF\_WORLD}. By these values it is checked
  whether the provided data matches the versions of the
  \ALBERTA-libraries in use.  \ALBERTA supports
  $\code{DIM\_OF\_WORLD}>3$ (but only meshes of dimension up to $3$).
  \ALBERTA-libraries with higher co-dimension can be selected through
  switches for the \code{configure}-script prior to compiling the
  \ALBERTA-package.
%%
\item[\code{number of vertices}] Mandatory. Number of vertex
  coordinates following the \code{vertex coordinates} keyword. The
  number of vertices must be specified prior to defining the
  coordinates themselves.
%%
\item[\code{number of elements}] Mandatory. The number of elements of
  the macro triangulation. This must be specified before defining any
  other data attached to elements, like the mesh connectivity or the
  neighbourhood relations.
%%
\item[\code{vertex coordinates}] Mandatory. The coordinates, specified
  by tuples of floating point values of dimension
  \code{DIM\_OF\_WORLD}.
%%
\item[\code{element vertices}] Mandatory. The mesh connectivity. The
  simplices are defined by their vertices, specified as offsets into
  the coordinate data defined in the \code{vertex coordinates}
  section.  Counting starts at $0$, so the first vertex has the number
  $0$. The data-line \code{0 3 4}, e.g., would define a triangle
  defined by the vertices $0$, $3$ and $4$. Note that the ordering of
  vertices defines the refinement edge, which is always located
  between the vertices with the local number $0$ and $1$.  This
  ordering of vertices (and the element type for 3d) determines the
  distribution of the refinement edges for the children.
%%
\item[\code{element boundaries}] Optional. For each element one line,
  which assigns a number between $0$ and $255$ (respectively $-128$
  and $+127$) to each co-dimension 1 sub-simplex of each element. The
  \code{element boundaries}-key may be omitted. If this is the case,
  each boundary segment is assigned a number of $1$. Note that
  interior walls have to be assigned a value of $0$.

  In the context of periodic meshes, periodic boundaries can still
  carry a classification number. This number is accessible in the
  \code{MACRO\_EL}-structure (see \ref{T:MACRO_EL}) and during
  non-periodic mesh-traversal.
  \begin{compatibility}
    \label{CPT:boundary_street_numbers}
    If the macro file contains boundary
    ``types'', then those are treated as mere ``street numbers'' by
    the current \ALBERTA version. Previous versions used positive
    numbers to indicate that a given boundary segment was subject to
    Dirichlet boundary condition and negative numbers were used to
    indicate that the respective segment carried natural boundary
    conditions.

    This was dropped because
    \begin{enumerate}
    \item the macro-triangulation should carry geometric information
      only
    \item it imposed too many restrictions, especially for the case
      were different components of systems of differential equations
      may be subject to different kind of boundary conditions on the
      same boundary segment
    \end{enumerate}
    Therefore the new scheme is now to only provide a classification
    of boundary segments by the macro triangulation. The
    interpretation of this classification is then left to the
    application program.
    
    Vertices (2d) and edges (3d) may in fact belong to boundary
    segments with different ``street numbers''. This information is
    for example accessible through the function \code{get\_bound()},
    see \secref{S:fillgetelvec}. See also \secref{S:boundary}.

    The current \ALBERTA versions prefer positive boundary-types, the
    \code{BNDRY\_TYPE} data type is in fact an \code{unsigned char} at
    the moment. Negative boundary type from ``old'' macro-data files
    are interpreted as positive numbers by the usual $2$-complement
    arithmetic.
  \end{compatibility}
%%
\item[\code{element neighbours}] Optional. Neighbourhood
  relationships. This information may be omitted from the
  macro-triangulation in which case it is computed by \ALBERTA. This
  computation is costly for large triangulations, so if neighbourhood
  information is available, it is advisable to include it in the
  macro-triangulation if the macro triangulation is a mesh with many
  simplices. If given then for each wall of each element the number of
  the neighbouring element has to be specified, or $-1$ if there is no
  such neighbour.
%%
\item[\code{element type}] Optional. This key is relevant only for 3d.
  In 3d, each element carries a ``type'' between $0$ an $2$
  (inclusive). This type influences the mesh refinement algorithm, see
  \secref{book:S:refinement_algorithm}. If the \code{element type} key is
  omitted, then \ALBERTA assigns each macro-element a type of $0$.
%%
\item[\code{number of wall transformations}] Optional. The number of
  face transformations which define a periodic structure on the mesh.
  See below under \code{wall transformations}.
%%
\item[\code{wall transformations}] Optional. For \ALBERTA, a periodic
  mesh is (part of) the fundamental domain of a crystallographic
  group. A fundamental domain of such group comes with a dedicated set
  of generators of the crystallographic group: the
  face-transformations which map the current fundamental domain to its
  neighbour across a given face (the notion ``face'' is already
  occupied within \ALBERTA, denoting co-dimension $1$ face-simplices
  in 3d, so ``wall'' denotes what ``face'' should have been used for:
  a co-dimension $1$ face-simplex, separating a simplex from its
  direct neighbour).
  
  The group-generators have to be specified in projective notation,
  acting on column vectors. For example, a simple translation by an
  amount of $2$ in $x_2$-direction in 3d would be specified as
  \bv\begin{verbatim}
 1 0 0 0
 0 1 0 2
 0 0 1 0
 0 0 0 1
\end{verbatim}\ev
  \ALBERTA assumes that the group generators are affine isometries,
  consequently, the inverses of the generators need not be specified.
  It is not necessary to format the matrices as shown above, \ALBERTA
  reads as many white-space separated numbers as it needs.
  See also \secref{S:periodic}.
%%
\item[\code{element wall transformations}] Optional. The corresponding
  information is computed by \ALBERTA when reading the macro-file if
  it is omitted. If specified, the data defines for each wall of each
  element the index into the array of wall-transformations which maps
  the current mesh to its periodic neighbour across the given wall.
  Per convention counting starts at $1$, where negative numbers denote
  the inverse of the given wall-transformation. A number of $0$
  indicates that the specific wall does not carry a face
  transformation (this applies to all interior walls as well as
  non-periodic parts of the boundary).  See also \secref{S:periodic}.
%%
\item[\code{number of wall vertex transformations}] Optional. Number
  of combinatoric face transformations. See below \code{wall vertex
    transformations}.
%%
\item[\code{wall vertex transformations}] Optional, computed on the
  fly if omitted. If specified the data following this key defines
  combinatoric face transformation by mapping boundary faces -- given
  by the global number of their vertices -- onto other boundary faces.
  For instance, to map a 2d boundary face -- an edge -- connecting
  vertex $0$ and $1$ onto the boundary edge between the vertices
  numbered $6$ and $7$ the following data would have to be specified:
  \bv\begin{verbatim}
 0 6
 1 7
\end{verbatim}\ev
  The ordering is important. Above lines implies that vertex $0$ is
  identified with vertex number $6$ and vertex number $1$ is
  identified with vertex $7$ -- or that the corresponding edges are
  identified with the orientation implied by the given ordering of the
  vertices. See also \secref{S:periodic}.
\end{description}

\begin{example}[The standard triangulation of the unit interval in $\R^1$]
  \label{E:MACRO_UNIT_INTERVAL}
\idx{macro triangulation!example of a macro triangulation in 1d}
\idx{macro triangulation!unit interval in 1d}
The easiest example is the macro triangulation for the interval $(0,1)$
in 1d. We just have one element and two vertices.
%%
\begin{minipage}{0.4\hsize}
\bv\begin{lstlisting}[name=MACRO_UNIT_INTERVAL,label=C:MACRO_UNIT_INTERVAL]
DIM: 1
DIM_OF_WORLD: 1

number of elements: 1
number of vertices: 2  

element vertices:
0 1

vertex coordinates:
 0.0 0.0  
 1.0 0.0
\end{lstlisting}\ev
\end{minipage}\hfill
\begin{minipage}{0.5\hsize}
\begin{center}
\includegraphics[scale=1.0]{EPS/macro1d}\\[2mm]
Macro triangulation of the unit interval.
\end{center}
\end{minipage}
\end{example}

\begin{example}[The standard triangulation of the unit square in $\R^2$]
  \label{E:MACRO_UNIT_SQUARE}
\idx{macro triangulation!unit square in 2d}
\idx{macro triangulation!example of a macro triangulation in 2d}
Still rather simple is the macro triangulation for the unit square
$(0,1)\times(0,1)$ in 2d. Here, we have two elements and four vertices.
The refinement edge is the diagonal for both elements.

\begin{minipage}{0.4\hsize}
\bv\begin{lstlisting}[name=MACRO_UNIT_SQUARE,label=C:MACRO_UNIT_SQUARE]
DIM: 2
DIM_OF_WORLD: 2

number of elements: 2  
number of vertices: 4  

element vertices:
2 0 1
0 2 3

vertex coordinates:
 0.0 0.0  
 1.0 0.0
 1.0 1.0
 0.0 1.0
\end{lstlisting}\ev
\end{minipage}\hfill
\begin{minipage}{0.5\hsize}
\begin{center}
\includegraphics[scale=0.5]{EPS/macro2d}\\[2mm]
Macro triangulation of the unit square.
\end{center}
\end{minipage}
\end{example}

\begin{example}[The standard triangulation of the unit cube in $\R^3$]%
\label{E:MACRO_UNIT_CUBE}
\idx{macro triangulation!example of a macro triangulation in 3d}
\idx{macro triangulation!unit cube in 3d}

More involved is already the macro triangulation for the unit cube
$(0,1)^3$ in 3d. Here, we have eight vertices and six elements, all
meeting at one diagonal; the shown specification of \code{element
vertices} prescribes this diagonal as the refinement edge for all
elements.

\begin{minipage}{0.4\hsize}
\bv\begin{lstlisting}[name=MACRO_UNIT_CUBE,label=C:MACRO_UNIT_CUBE]
DIM:          3 
DIM_OF_WORLD: 3

number of vertices: 8
number of elements: 6

vertex coordinates:
  0.0  0.0  0.0
  1.0  0.0  0.0
  0.0  0.0  1.0
  1.0  0.0  1.0
  1.0  1.0  0.0
  1.0  1.0  1.0
  0.0  1.0  0.0
  0.0  1.0  1.0

element vertices:
  0    5    4    1
  0    5    3    1
  0    5    3    2
  0    5    4    6
  0    5    7    6
  0    5    7    2
\end{lstlisting}\ev
\end{minipage}
\hfill
\begin{minipage}{0.5\hsize}
\begin{center}
\includegraphics{EPS/macro3d}\\[2mm]
Macro triangulation of the unit cube.
\end{center}
\end{minipage}
\end{example}

\begin{example}[A triangulation of three quarters of the unit disc]
  \idx{macro triangulation!example for three quarters of the unit disc}
  Here, we describe a more complex example where we are dealing with a
  curved boundary and mixed type boundary condition. Due to the curved
  boundary, we have to initialize the projection mechanism when allocating 
  a mesh as described in \secref{S:node_projections}. The actual projection is
  easy to implement, since we only have normalize the coordinates for
  nodes belonging to the curved boundary. 
  We assume that the two straight edges belong to the
  Neumann boundary, and the curved boundary is the Dirichlet boundary.
  For handling mixed boundary types we have to specify \code{element
    boundaries} in the macro triangulation file. Information about
  \code{element boundaries} is also used inside the function
  \code{init\_node\_proj}.\medskip

\begin{minipage}[b]{0.4\hsize}
\bv\begin{lstlisting}
DIM: 2 
DIM_OF_WORLD: 2

number of vertices: 5  
number of elements: 3  

vertex coordinates:
 0.0  0.0  
 1.0  0.0
 0.0  1.0
-1.0  0.0
 0.0 -1.0

element vertices:
1 2 0
2 3 0
3 4 0

element boundaries:
 0 -1  2 
 0  0  2 
-1  0  2 
\end{lstlisting}\ev
\end{minipage}
\hfill
\begin{minipage}[b]{0.5\hsize}
\begin{center}
\includegraphics[scale=0.75]{EPS/macro34}\\[2mm]
Macro triangulation of a 3/4 disc.
\end{center}
\end{minipage}
\medskip

\noindent
The function \code{init\_node\_proj()} to initialize projection of nodes
can be implemented similarly to Example \ref{Ex:unit_ball}. The projection
routine \code{ball\_proj\_func} remains the same.
\idx{node projection!example for three fourths of the unit disc}%
\bv\begin{lstlisting}
static NODE_PROJECTION *init_node_proj(MESH *mesh, MACRO_EL *mel, int c)
{
  static NODE_PROJECTION ball_proj = {ball_proj_func};

  if(c > 0 && mel->edge_bound[c-1] == 2)
    return &ball_proj;
  else
    return nil;
}
\end{lstlisting}\ev
\end{example}

\subsubsection{Reading macro triangulations from a file}
Reading data of the macro grid from these files can be done by
\fdx{read_macro()@{\code{read\_macro()}}}
\idx{macro triangulation!read_macro()@{\code{read\_macro()}}}
\bv\begin{lstlisting}
MACRO_DATA *read_macro(const char *filename);
\end{lstlisting}\ev
Description:
\begin{descr}
  \kitem{read\_macro(filename)} reads data of the macro triangulation
  from the \code{ASCII}-file \code{filename} and returns a pointer to
  a filled \code{MACRO\_DATA} structure (see \secref{S:macro_import}).
  Using index information from the file, all information concerning
  element vertices, neighbour relations can be calculated directly.
  
  During the initialization of the macro triangulation, other entries
  like \code{n\_edges}, \code{n\_faces}, and \code{max\_edge\_neigh}
  in the mesh data structure are calculated. Please note that
  projection of nodes as described in \secref{S:node_projections} is 
  only possible for new nodes arising during refinement.
\end{descr}

A binary data format allows faster import of a macro triangulation,
especially when the macro triangulation already consists of many
elements.  Macro data written previously by binary \code{write\_macro}
routines (see below) can be read in native or machine independent
binary format by the two routines
\fdx{read_macro_bin()@{\code{read\_macro\_bin()}}}
\idx{macro triangulation!read_macro_bin()@{\code{read\_macro\_bin()}}}
\fdx{read_macro_xdr()@{\code{read\_macro\_xdr()}}}
\idx{macro triangulation!read_macro_xdr()@{\code{read\_macro\_xdr()}}}
\bv\begin{lstlisting}
MACRO_DATA *read_macro_bin(const char *filename);
MACRO_DATA *read_macro_xdr(const char *filename);
\end{lstlisting}\ev
Description:
\begin{descr}
  \kitem{read\_macro\_bin(filename)} reads data of the macro
  triangulation from the native binary file
  \code{filename}; the file \code{filename} was previously generated by the
  function \code{write\_macro\_bin()}, see below.
  
  \kitem{read\_macro\_xdr(filename)} reads data of the macro
  triangulation from the machine independent binary
  file \code{filename}, the file \code{filename} was previously generated by
  the function \code{write\_macro\_xdr()}, see below.
\end{descr}

\subsubsection{Dumping macro triangulations to a file}%
\label{S:write_macro}%
\idx{macro triangulation!writing macro triangulations}

The counterpart of functions for reading macro triangulations
are functions for writing macro triangulations to file.
To be more general, it is possible to create a macro triangulation
from the triangulation given by the leaf elements of a mesh.
As mentioned above, it can be faster to use a binary format than the
textual formal for writing and reading macro triangulations with many
elements.
\fdx{write_macro()@{\code{write\_macro()}}}
\idx{macro triangulation!write_macro()@{\code{write\_macro()}}}
\fdx{write_macro_bin()@{\code{write\_macro\_bin()}}}
\idx{macro triangulation!write_macro_bin()@{\code{write\_macro\_bin()}}}
\fdx{write_macro_xdr()@{\code{write\_macro\_xdr()}}}
\idx{macro triangulation!write_macro_xdr()@{\code{write\_macro\_xdr()}}}
\bv\begin{lstlisting}
int write_macro(MESH *, const char *);
int write_macro_bin(MESH *, const char *);
int write_macro_xdr(MESH *, const char *);
\end{lstlisting}\ev
Description:
\begin{descr}
\kitem{write\_macro(mesh, name)} writes the triangulation given by 
       the leaf elements of 
       \code{mesh} as a macro triangulation to the file specified by
       \code{name} in the above described format; if the file could
       be written, the return value is \code{1}, otherwise an error
       message is produced and the return value is \code{0}.
\kitem{write\_macro\_bin(mesh, name)} writes the triangulation given by
       the leaf elements of 
       \code{mesh} as a macro triangulation to the file specified by
       \code{name} in native binary format.
\kitem{write\_macro\_xdr(mesh, name)} writes the triangulation given by
       the leaf elements of 
       \code{mesh} as a macro triangulation to the file specified by
       \code{name} in machine independent binary format.
\end{descr}
For exporting meshes including the whole hierarchy, see \secref{S:file_formats}

\subsection{Import and export of macro triangulations from/to other formats}%
\label{S:macro_import}%
\idx{macro triangulation!import macro triangulations}%
\idx{macro triangulation!export macro triangulations}

When meshes are created using a simplicial grid generation tool, then
data will usually not be in the \ALBERTA macro triangulation format
described above in Section \ref{S:macro_tria}. In order to simplify
the import of such meshes, an array-based data structure
\code{MACRO\_DATA} is provided, using flat \code{C}-arrays for storing
the data, and indirect index-arrays to bind the data to elements and
define the mesh connectivity. Such a data structure can easily be
filled by an import routine; the filled data structure can then
converted into an \ALBERTA mesh. Of course, another possibility is to
convert the data to \ALBERTA's textual macro-file format as described
in \secref{S:macro_triang} above. The \code{MACRO\_DATA} structure is
defined as
%
\begin{samepage}
\ddx{MACRO_DATA@{\code{MACRO\_DATA}}}
\bv\begin{lstlisting}[name=MACRO_DATA,label=T:MACRO_DATA]
typedef struct macro_data MACRO_DATA;

struct macro_data
{
  int dim;              /* dimension of the elements */

  int n_total_vertices;
  int n_macro_elements;

  REAL_D *coords;       /* Length will be n_total_vertices */

  int *mel_vertices;    /* mel_vertices[i*N_VERTICES(dim)+j]:
			 * global index of jth vertex of element i
			 */

  int *neigh;           /* neigh[i*N_NEIGH(dim)+j]:
			 * neighbour j of element i or -1 at boundaries
			 */
  int *opp_vertex;      /* opp_vertex[i*N_NEIGH(dim)+j]: if set (need not
			 * be) the local vertex number w.r.t. the neighbour
			 * of the vertex opposit the separating wall.
			 */
  BNDRY_TYPE *boundary; /* boundary[i*N_NEIGH(dim)+j]:
			 * boundary type of jth co-dim 1 facet of element i
			 *
			 * WARNING: In 1D the local index corresponds
			 * to vertex 1 & vice versa! (Consistent with
			 * macro_data.neigh)
			 */
  U_CHAR *el_type;      /* el_type[i]: type of element i only used in 3d! */

  /********* the remainder is only needed for periodic meshes ***********/

  int (*wall_vtx_trafos)[N_VERTICES(DIM_MAX-1)][2]; /* the wall trafos */
  /* Wall transformations are in terms of mappings between
   * vertices. i-th wall trafo: global vertex number
   * wall_vtx_trafos[i][v][0] maps to wall_vtx_trafos[i][v][1],
   * v loops through the local vertex number of the respective wall.
   */
  int n_wall_vtx_trafos;/* for periodic meshes: number of
			 * combinatorical wall trafos.
			 */
  int *el_wall_vtx_trafos;
  /* el_wall_vtx_trafos[i*N_WALLS(dim)+j] number of the wall
   * transformation of the j-th wall for the i-th element. > 0:
   * #wall_trafo+1. < 0: inverse of -(#wall_trafo+1)
   */
  AFF_TRAFO *wall_trafos; /* The group generators of the space group
			   * defining the periodic structure of the
			   * mesh.
			   */
  int n_wall_trafos;
  int *el_wall_trafos; /* N = el_wall_trafos[i*N_NEIGH(dim)+j]:
			*
			* number of the wall transformation mapping to
			* the neighbouring fundamental domain across
			* the given wall.
			*
			* If negative: inverse of generator -N-1
			* If positive:            generator +N-1
			*/
};
\end{lstlisting}\ev
\end{samepage}
The members yield following information:
\begin{descr}
  \kitem{dim} dimension of the triangulation.
  %% 
  \kitem{n\_total\_vertices} number of vertices.
  %%
  \kitem{n\_macro\_elements} number of mesh elements.
  %%
  \kitem{coords} \code{REAL\_D} array of size 
  \code{n\_total\_vertices} holding the
  point coordinates of all vertices.
  %%
  \kitem{mel\_vertices} integer array of size \code{n\_macro\_elements
    * N\_VERTICES(dim)} storing element index information;
  \code{mel\_vertices[i*N\_VERTICES[dim]+j]} is the index of the
  \code{j}th vertex of element \code{i}.
  %%
  \kitem{neigh} integer array of size
  \code{n\_macro\_elements*N\_NEIGH(dim)}, where
  \code{neigh[i*N\_NEIGH(dim)+j]} is the index of the \code{j}th
  neighbour element of element \code{i}, or \code{-1} in case of a
  boundary.
  %% 
  \kitem{boundary} \code{S\_CHAR} array of size
  \code{n\_macro\_elements*N\_NEIGH(dim)}, where 
  \code{boundary[i*N\_NEIGH(dim)+j]} is the
  boundary type of the \code{j}th vertex/edge/face of element \code{i}
  (in 1d/2d/3d). Please note that the index \code{0} corresponds to vertex $1$
  and vice versa in 1d, consistent with the numbering used for \code{neigh}.
  %%
  \kitem{el\_type} a \code{U\_CHAR} vector of size
  \code{n\_macro\_elements} holding the element type of each mesh element
  (only 3d).
  %%
  \kitem{wall\_vtx\_trafos, n\_wall\_vtx\_trafos} correspond to the
  data specified with the key \code{wall vertex transformations}, see
  \secref{S:macro_tria}. This field stores face-transformations in
  terms of mappings between vertices. For the $i$-the face
  transformation the global vertex number
  \code{wall\_vtx\_trafos[i][v][0]} maps to the global vertex number
  \code{wall\_vtx\_trafos[i][v][1]}, v loops through the \emph{local}
  vertex number of the respective wall.
  %% 
  \kitem{el\_wall\_vtx\_trafos} If
  $\code{el\_wall\_vtx\_trafos[i*N\_WALLS(dim)+j]}!=0$ then it is the
  number of the face-transformation the $j$-th wall on the for the
  $i$-th element is subject to. Negative number indicate that the
  inverse of the respective face-transformation is attached to that
  wall. Note that one has to subtract $1$ from this value before using
  it as index into \code{wall\_vtx\_trafos}, because arrays in
  \code{C} are indexed starting with $0$.
  %%
  \kitem{wall\_trafos, n\_wall\_trafos} The group generators and their
  number of the space group defining the periodic structure of the
  mesh. See \secref{S:periodic}.
  %%
  \kitem{el\_wall\_trafos} If $\code{N =
    el\_wall\_trafos[i*N\_NEIGH(dim)+j]}!=0$ then $N$ is the number of
  the face-transformation mapping the mesh to the neighboring
  fundamental domain across the given wall. If $N$ is negative, then
  the actual face-transformation is the inverse of the $N$-th
  transformation.  Note that one has to subtract $1$ from this value
  before using it as index into \code{wall\_trafos}, because arrays in
  \code{C} are indexed starting with $0$.
\end{descr}

A \code{MACRO\_DATA} structure can be allocated and freed by
\fdx{get_macro_data()@{\code{get\_macro\_data()}}}
\fdx{free_macro_data()@{\code{free\_macro\_data()}}}
\bv\begin{lstlisting}
MACRO_DATA *alloc_macro_data(int dim, int nv, int ne, FLAGS);
void free_macro_data(MACRO_DATA *);
\end{lstlisting}\ev
Description:
\begin{descr}
\kitem{alloc\_macro\_data(dim, n\_vertices, n\_elements, flags)} allocates a
  \code{dim}-dimensional
  \code{MACRO\_DATA} structure together with all arrays needed to hold
  \code{n\_vertices} vertices and \code{n\_elements} mesh
  elements. The \code{coords} and  \code{mel\_vertices} arrays are
  allocated in any case, while \code{neigh}, \code{boundary} and
  \code{el\_type} arrays are allocated only when requested as
  indicated by the corresponding flags \code{FILL\_NEIGH},
  \code{FILL\_BOUNDARY}, and \code{FILL\_EL\_TYPE} set by a
  bitwise \textsf{OR} in \code{flags}.
\kitem{free\_macro\_data(data)} frees all previously allocated storage
  for \code{MACRO\_DATA} \code{data} and all the arrays in it.
\end{descr}

Once \code{MACRO\_DATA} structure is filled, it can be
saved to file in the \ALBERTA macro triangulation format, or 
it can be directly be converted into a \code{MESH}.
\fdx{macro_data2mesh()@{\code{macro\_data2mesh()}}}
 \idx{macro triangulation!macro_data2mesh()@{\code{macro\_data2mesh()}}}
\fdx{write_macro_data()@{\code{write\_macro\_data()}}}
\idx{macro triangulation!write_macro_data()@{\code{write\_macro\_data()}}}
\fdx{write_macro_data_bin()@{\code{write\_macro\_data\_bin()}}}
\idx{macro triangulation!write_macro_data_bin()@{\code{write\_macro\_data\_bin()}}}
\fdx{write_macro_data_xdr()@{\code{write\_macro\_data\_xdr()}}}
\idx{macro triangulation!write_macro_data_xdr()@{\code{write\_macro\_data\_xdr()}}}
\bv\begin{lstlisting} 
void macro_data2mesh(MESH *mesh, const MACRO_DATA *data,
                     NODE_PROJECTION *(*n_proj)(MESH *,MACRO_EL *,int));
int write_macro_data(MACRO_DATA *, const char *);
int write_macro_data_bin(MACRO_DATA *, const char *);
int write_macro_data_xdr(MACRO_DATA *, const char *);
\end{lstlisting}\ev
Description:
\begin{descr}
\kitem{macro\_data2mesh(mesh, macro\_data, n\_proj)} converts the
  triangulation with data given in \code{macro\_data} into a
  \code{MESH} structure. It sets most entries in \code{mesh},
  allocates macro elements needed, assigns DOFs according to
  \code{mesh->n\_dof}, and calculates \code{mesh->diam}. The
  coordinates in \code{macro\_data->coords} are copied to a newly
  allocated array, thus the entire \code{MACRO\_DATA} structure can
  be freed after calling this routine. When not \code{nil}, the
  \code{n\_proj} function is used to initialize projection of new nodes.
\kitem{write\_macro\_data(macro\_data, name)} writes the macro
  triangulation with data stored in \code{macro\_data} in the \ALBERTA
  format described in Section \ref{S:macro_tria} to file \code{name}. 
  The return value is \code{0} when an error occured and
  \code{1} in case the file was written successfully.
\kitem{write\_macro\_data\_bin(macro\_data, name)} writes data of the
  macro triangulation stored in \code{macro\_data} in native
  binary format to file \code{name}; the return value is \code{0} when
  an error occured and \code{1} in case the file was written successfully.
\kitem{write\_macro\_data\_xdr(macro\_data, name)}
  writes data of the macro triangulation stored in
  \code{macro\_data} in machine independent binary format to file 
  \code{name}; the return value is \code{0} when an error occured and
  \code{1} in case the file was written successfully.
\end{descr}

It is appropriate to check whether a macro triangulation given in a
\code{MACRO\_DATA} structure allows for recursive refinement, by
testing for possible recursion cycles. An automatic correction by
choosing other refinement edges may be done, currently implemented only
in 2d.
%
\fdx{macro_test()@{\code{macro\_test()}}}
\idx{macro triangulation!macro_test()@{\code{macro\_test()}}}
\bv\begin{lstlisting}
 void macro_test(MACRO_DATA *, const char *);
\end{lstlisting}\ev
Description:
\begin{descr}
\kitem{macro\_test(macro\_data, name)} checks the triangulation
  given in \code{macro\_data} for potential cycles during recursive
  refinement. In the case that such a cycle is detected, the routine
  tries to correct this by renumbering element vertices (which is
  currently implemented only in 2d) and then writes the new, changed
  triangulation using \code{write\_macro\_data()} to a file
  \code{name}, when the second parameter is not \code{nil}.
\end{descr}

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