\section{Implementation of basis functions}%
\label{S:basfct_impl}
\label{S:BAS_FCTS}

In order to construct a finite element space, we have to specify
a set of local basis functions. We follow the concept of finite
elements which are given on a single element $S$ in local coordinates:
Finite element functions on an element $S$ are defined by a finite
dimensional function space $\Pbar$ on a reference element $\Sbar$ and
the (one to one) mapping $\lambda^S: \Sbar \to S$ from the reference
element $\Sbar$ to the element $S$. In this situation the non
vanishing basis functions on an arbitrary element are given by the set
of basis functions of $\Pbar$ in local coordinates $\lambda^S$.  Also,
derivatives are given by the derivatives of basis functions on $\Pbar$
and derivatives of $\lambda^S$. 

Each local basis function on $S$ is uniquely connected to a global
degree of freedom, which can be accessed from $S$ via the DOF
administration. Together with this DOF administration and the
underlying mesh, the finite element space is given. In the following
section we describe the basic data structures for storing basis
function information.

At the moment the following finite elements are supported by \ALBERTA:
\begin{itemize}
\item standard Lagrange finite elements of order $n$, $n\in\{1,2,3,4\}$;
\item discontinuous polynomial elements of order $n$, $n\in\{0,1,2\}$; these 
are defined as arbitrary polynomials of maximal degree $n$ on each element
with no continuity restriction;
\item orthonormal discontinuous polynomial elements of order $1$ and
  $2$; these functions are normalized and orthogonal w.r.t. the
  $L^2$-scalar product on the reference element.
\end{itemize}
We present these elements in the subsequent sections.  A tselection of
more complicated basis functions is implemented in an add-on library
called \verb|libalbas|, with the focus on stable discretizations for
the Stokes-problem. This is not discussed here.

\subsection{Data structures for basis functions}%
\label{S:basfct_data}

For the handling of local basis functions, i.e. a basis of the
function space $\Pbar$ on the reference element (compare Section
\ref{book:S:FES}) we use functions of the following type. The structure
describing the set of local basis functions (\verb|BAS_FCTS|, see
page \pageref{T:BAS_FCTS}) contains arrays of such function-pointers:
%%
\idx{basis fucntions!BAS_FCT@{basis function hooks}}
\ddx{BAS_FCT@{\code{BAS\_FCT}}}%
\ddx{GRD_BAS_FCT@{\code{GRD\_BAS\_FCT}}}%
\ddx{D2_BAS_FCT@{\code{D2\_BAS\_FCT}}}%
\ddx{D3_BAS_FCT@{\code{D3\_BAS\_FCT}}}%
\ddx{D4_BAS_FCT@{\code{D4\_BAS\_FCT}}}%
\ddx{BAS_FCT_D@{\code{BAS\_FCT\_D}}}%
\ddx{GRD_BAS_FCT_D@{\code{GRD\_BAS\_FCT\_D}}}%
\ddx{D2_BAS_FCT_D@{\code{D2\_BAS\_FCT\_D}}}%
\bv
\begin{lstlisting} 
  typedef REAL
    (*BAS_FCT)(const REAL_B lambda, const BAS_FCTS *thisptr);
  typedef const REAL *
    (*GRD_BAS_FCT)(const REAL_B lambda, const BAS_FCTS *thisptr);
  typedef const REAL_B *
    (*D2_BAS_FCT)(const REAL_B lambda, const BAS_FCTS *thisptr);
  typedef const REAL_BB *
    (*D3_BAS_FCT)(const REAL_B lambda, const BAS_FCTS *thisptr);
  typedef const REAL_BBB *
    (*D4_BAS_FCT)(const REAL_B lambda, const BAS_FCTS *thisptr);

  typedef const REAL *
  (*BAS_FCT_D)(const REAL_B lambda, const BAS_FCTS *thisptr);
  typedef const REAL_B *
  (*GRD_BAS_FCT_D)(const REAL_B lambda, const BAS_FCTS *thisptr);
  typedef const REAL_BB *
  (*D2_BAS_FCT_D)(const REAL_B lambda, const BAS_FCTS *thisptr);
\end{lstlisting}\ev
Description:
\begin{descr}
%%
\kitem{BAS\_FCT} the data type for a local finite element function,
  i.e. a function $\pbar \in \Pbar$, evaluated at barycentric
  coordinates $\lambda \in \R^{d+1}$ and its return
  value $\pbar(\lambda)$ is of type \code{REAL}.
%%
\kitem{GRD\_BAS\_FCT} the data type for the gradient (with respect 
  to $\lambda$) of a local finite element function, i.e. a
  function returning a pointer to $\nablal \pbar$ for some function
  $\pbar \in \Pbar$:
  \[
  \nablal \pbar(\lambda) = \left(
    \frac{\partial\pbar(\lambda)}{\partial \lambda_0},
    \ldots,\frac{\partial\pbar(\lambda)}{\partial \lambda_d}
  \right);
  \]
  the arguments of such a function are barycentric coordinates and
  the return value is a pointer to a \code{const REAL} vector of
  length \code{N\_LAMBDA} storing $\nablal \pbar(\lambda)$; this
  vector will be overwritten during the next call of the function.
  %% 
\kitem{D2\_BAS\_FCT} the data type for the second derivatives (with
  respect to $\lambda$) of a local finite element function, i.e. a
  function returning a pointer to the matrix $D^2_\lambda \pbar$
  for some function $\pbar \in \Pbar$:
  \[
  D^2_\lambda \pbar =
  \begin{pmatrix}\ds
    \frac{\partial^2\pbar(\lambda)}{\partial \lambda_0 \partial \lambda_0} 
    & \cdots &\ds
    \frac{\partial^2\pbar(\lambda)}{\partial \lambda_0\partial\lambda_d} \\
    \vdots & & \vdots\\\ds
    \frac{\partial^2\pbar(\lambda)}{\partial \lambda_d\partial \lambda_0} 
    & \cdots &\ds
    \frac{\partial^2\pbar(\lambda)}
    {\partial\lambda_d\partial\lambda_d}
  \end{pmatrix};
  \]
  the arguments of such a function are barycentric coordinates and
  the return value is a pointer to a $\mathtt{N\_LAMBDA \times
    N\_LAMBDA}$ matrix storing $D^2_\lambda \pbar$; this matrix will be
  overwritten during the next call of the function.
  %% 
\kitem{D3\_BAS\_FCT, D4\_BAS\_FCT}
  serve primarily for debugging purposes and need not be present.
  The format is the similar to the second barycentric derivatives, the third
  derivatives are a tensor of rank $3$, and the fourth derivatives are a
  tensor of rank $4$.
  %% 
\kitem{BAS\_FCT\_D, GRD\_BAS\_FCT\_D, D2\_BAS\_FCT\_D} Basis functions
  may optionally be \REALD-valued.  In this case we factor the basis
  functions into a scalar part which is multiplied by a direction,
  with the obvious implications for the derivatives of such basis
  functions. This is further explained below where the corresponding
  components of the \code{BAS\_FCTS} structure are discussed.
\end{descr}

\mdx{PHI@{\code{PHI()}}}
\mdx{GRD_PHI@{\code{GRD\_PHI()}}}
\mdx{D2_PHI@{\code{D2\_PHI()}}}
\mdx{D3_PHI@{\code{D3\_PHI()}}}
\mdx{D4_PHI@{\code{D4\_PHI()}}}
\mdx{PHI_D@{\code{PHI\_D()}}}
\mdx{GRD_PHI_D@{\code{GRD\_PHI\_D()}}}
\mdx{D2_PHI_D@{\code{D2\_PHI\_D()}}}
\bv\begin{lstlisting}
#define PHI(bfcts, i, lambda)     (bfcts)->phi[i](lambda, bfcts)
#define GRD_PHI(bfcts, i, lambda) (bfcts)->grd_phi[i](lambda, bfcts)
#define D2_PHI(bfcts, i, lambda)  (bfcts)->D2_phi[i](lambda, bfcts)
#define D3_PHI(bfcts, i, lambda)  (bfcts)->D3_phi[i](lambda, bfcts)
#define D4_PHI(bfcts, i, lambda)  (bfcts)->D4_phi[i](lambda, bfcts)

#define PHI_D(bfcts, i, lambda)     (bfcts)->phi_d[i](lambda, bfcts)
#define GRD_PHI_D(bfcts, i, lambda) (bfcts)->grd_phi_d[i](lambda, bfcts)
#define D2_PHI_D(bfcts, i, lambda)  (bfcts)->D2_phi_d[i](lambda, bfcts) 
\end{lstlisting}\ev
Description:
\begin{descr}
\kitem{PHI(bfcts, i, lambda)} The individual basis function pointers
  expect that a pointer to the \verb|BAS_FCTS|-structure they belong
  to is passed as last argument. To decrease the potential for coding
  errors we advocate to use
  \bv\begin{lstlisting}
    value = PHI(bfcts, nr, lambda, nr);
  \end{lstlisting}\ev
  instead of using the equivalent construct
  \bv\begin{lstlisting}
    value = bfcts->phi[nr](lambda, bfcts);
  \end{lstlisting}\ev
  The other macro work analogously.
\end{descr}

For the implementation of a finite element space, we need a basis of
the function space $\Pbar$. For such a basis we need the connection of
\emph{local} and \emph{global} DOFs on each element (compare
\secref{book:S:eval_fe}), information about the interpolation of a given
function on an element, and information about
interpolation/restriction of finite element functions during
refinement/coarsening (compare \secref{book:S:inter_restrict}). Further
information includes the traces of the local finite element space on
the walls of reference element, which are needed to define
trace-meshes (AKA sub-meshes, see Section
\ref{S:tracemesh_implementation}). Also, finite element spaces may
form a direct sum, e.g. to implement stable discretisations of the
Stokes-problem (see Section \ref{S:chain_impl}); if this is the case
then the local basis-functions also reflect this fact. Finally,
basis-functions may need a per-element initializer, for example if
they depend on the geometry of the mesh-simplex.  Such information is
stored in the \code{BAS\_FCTS} data structure:
%%
\ddx{BAS_FCTS@{\code{BAS\_FCTS}}}
\idx{basis functions!BAS_FCTS@{\code{BAS\_FCTS} data type}}
%%
\bv\begin{lstlisting}[name=BAS_FCTS,label=T:BAS_FCTS]
typedef struct bas_fcts BAS_FCTS;

struct bas_fcts
{
  const char   *name;       /*  textual description */
  int          dim;         /*  dimension of the corresponding mesh. */
  int          rdim;        /*  dimension of the range, 1 or DIM_OF_WORLD */
  int          n_bas_fcts;  /*  nu_mber of basisfunctions on one el */
  int          n_bas_fcts_max;/*  max. number in presence of init_element()*/
  int          degree;      /*  maximal degree of the basis functions,
			     *	may vary on a per-element basis if
			     *	init_element() is != NULL.
			     */
  int          n_dof[N_NODE_TYPES];   /* dofs from these bas_fcts */
  int          trace_admin; /* If >= 0, then the basis function set
			     * needs a DOF_ADMIN living on a trace
			     * mesh with id TRACE_ADMIN.
			     */

  /************** link to next set of bfcts in a direct sum ***************/
  DBL_LIST_NODE  chain;

  /* A pointer to the unchained version. It simply points back to the
   * same structure if this is an unchained basis-function
   * structure.
   */
  const BAS_FCTS *unchained;

  /*************** per-element initializer (maybe NULL) *******************/

  INIT_ELEMENT_DECL;

  /*************** the basis functions themselves       *******************/
  const BAS_FCT     *phi;
  const GRD_BAS_FCT *grd_phi;
  const D2_BAS_FCT  *D2_phi;
  const D3_BAS_FCT  *D3_phi; /* Optional, implemented for Lagrange bfcts. */
  const D4_BAS_FCT  *D4_phi; /* Optional, implemented for Lagrange bfcts. */

  /* Vector valued basis functions are always factored as phi[i]() *
   * phi_d[i](). If phi_d[i]() is piece-wise constant, then
   * dir_pw_const should be true. The directions are never cached in
   * QUAD_FAST, only the scalar factor.
   */
  const BAS_FCT_D     *phi_d;
  const GRD_BAS_FCT_D *grd_phi_d;
  const D2_BAS_FCT_D  *D2_phi_d;

  bool dir_pw_const; /*Direction is p.w. constant on the reference element.*/

  /************* the trace space on the wall *****************************/
  const BAS_FCTS *trace_bas_fcts; /* The trace space */

  /* The local DOF mapping for the trace spaces,
   * < 3d:
   * [0][0][wall][slave local dof] == master local dof,
   * 3d:
   * [type > 0][orient < 0][wall][slave local dof] == master local dof.
   */
  const int      *trace_dof_map[2][2][N_WALLS_MAX];

  /* This obscure component can vary from wall to wall in the presence
   * of an INIT_ELEMENT() method. It is _always_ equal to
   * trace_bas_fcts->n_bas_fcts ... BUT ONLY after the respective
   * element initializer has been called for trace_bas_fcts on the
   * trace mesh. If an INIT_ELEMENT() method is present then it _MUST_
   * initialize trace_dof_map _AND_ n_trace_bas_fcts. Of course, in 3D
   * only the components corresponding to type and orientation of the
   * current EL_INFO object have to be taken care of by the
   * INIT_ELEMENT() method.
   */
  int            n_trace_bas_fcts[N_WALLS_MAX];

  /*************** interconnection to DOF_ADMIN and mesh   ****************/
  const EL_DOF_VEC *(*get_dof_indices)(DOF *result,
				       const EL *, const DOF_ADMIN *,
				       const BAS_FCTS *thisptr);
  const EL_BNDRY_VEC *(*get_bound)(BNDRY_FLAGS *bndry_bits,
				   const EL_INFO *eli,
				   const BAS_FCTS *thisptr);

  /*************** entries must be set for interpolation   ****************/

  void (*interpol)(EL_REAL_VEC *coeff,
		   const EL_INFO *el_info, int wall,
		   int n, const int *indices,
		   LOC_FCT_AT_QP f, void *ud,
		   const BAS_FCTS *thisptr);
  void (*interpol_d)(EL_REAL_D_VEC *coeff,
		     const EL_INFO *el_info, int wall,
		     int n, const int *indices,
		     LOC_FCT_D_AT_QP f, void *ud,
		     const BAS_FCTS *thisptr);
  void (*interpol_dow)(EL_REAL_VEC_D *coeff,
		       const EL_INFO *el_info, int wall,
		       int n, const int *indices,
		       LOC_FCT_D_AT_QP f, void *ud,
		       const BAS_FCTS *thisptr);

  /********************   optional entries  *******************************/

  const EL_INT_VEC    *(*get_int_vec)(int result[],
				      const EL *, const DOF_INT_VEC *);
  const EL_REAL_VEC   *(*get_real_vec)(REAL result[],
				       const EL *, const DOF_REAL_VEC *);
  const EL_REAL_D_VEC *(*get_real_d_vec)(REAL_D result[],
					 const EL *, const DOF_REAL_D_VEC *);
  const EL_REAL_VEC_D *(*get_real_vec_d)(REAL result[],
					 const EL *, const DOF_REAL_VEC_D *);
  const EL_UCHAR_VEC  *(*get_uchar_vec)(U_CHAR result[],
					const EL *, const DOF_UCHAR_VEC *);
  const EL_SCHAR_VEC  *(*get_schar_vec)(S_CHAR result[],
					const EL *, const DOF_SCHAR_VEC *);
  const EL_PTR_VEC    *(*get_ptr_vec)(void *result[],
				      const EL *, const DOF_PTR_VEC *);

  void  (*real_refine_inter)(DOF_REAL_VEC *, RC_LIST_EL *, int);
  void  (*real_coarse_inter)(DOF_REAL_VEC *, RC_LIST_EL *, int);
  void  (*real_coarse_restr)(DOF_REAL_VEC *, RC_LIST_EL *, int);

  void  (*real_d_refine_inter)(DOF_REAL_D_VEC *, RC_LIST_EL *, int);
  void  (*real_d_coarse_inter)(DOF_REAL_D_VEC *, RC_LIST_EL *, int);
  void  (*real_d_coarse_restr)(DOF_REAL_D_VEC *, RC_LIST_EL *, int);

  void  (*real_refine_inter_d)(DOF_REAL_VEC_D *, RC_LIST_EL *, int);
  void  (*real_coarse_inter_d)(DOF_REAL_VEC_D *, RC_LIST_EL *, int);
  void  (*real_coarse_restr_d)(DOF_REAL_VEC_D *, RC_LIST_EL *, int);
  
  void  *ext_data; /* Implementation dependent extra data */
};

/* Barycentric coordinates of Lagrange nodes. */
#define LAGRANGE_NODES(bfcts)				\
  ((const REAL_B *)(*(void **)(bfcts)->ext_data))
\end{lstlisting}\ev

The entries yield following information:
\begin{descr}
\kitem{name} string containing a textual description or \nil.
  %% 
\kitem{dim} dimension $d$ of the mesh triangulation.
  %%
\kitem{rdim} dimension of the range space. This is either $1$ or \DOW
  for vector valued basis functions like edge and face bubbles.
  %%
\kitem{n\_bas\_fcts} number of local basis functions.
  %%
\kitem{n\_bas\_fcts\_max} maximum number of local basis functions. The
  number of basis functions on a given element may vary if the
  \verb|init_element|-hook is non-\nil. In this case
  \verb|n_bas_fcts_max| gives the upper limit and can be used to
  lay-out array dimensions, e.g. In the standard case this component
  does not differ from \verb|n_bas_fcts|. See \secref{S:chain_impl}.
  %%
\kitem{degree} maximal polynomial degree of the basis functions; this
  entry is used by routines using numerical quadrature where no
  \code{QUAD} structure is provided; in such a case via \code{degree}
  some default numerical quadrature is chosen (see Section
  \ref{S:QUAD}); additionally, \code{degree} is used by some graphics
  routines (see \secref{S:graph_2d}).
  %%
\kitem{n\_dof} vector with the count of DOFs for this set of basis
  functions; \code{n\_dof[VERTEX,CENTER,EDGE,FACE]} is the number of
  DOFs tied to the vertices, center, edges (only 2d and 3d), faces
  (only 3d), of an element; the corresponding DOF administration of
  the finite element space uses such information.
  %%
  \kitem{trace\_admin} is used for the purpose of defining basis
  functions with DOFs attached to a
  \hyperref[S:submesh_implementation]{trace-mesh}, e.g. to define
  face-bubbles attached to part of the boundary of the mesh. In this
  case \code{trace\_admin} is the unique ID of a trace-mesh which
  carries the \hyperref[S:DOF_ADMIN]{\code{DOF\_ADMIN}} for this
  basis-function set. In this situation the basis functions ``live''
  on the bulk mesh (i.e. extend in to to bulk-phase of the element
  containing the face belonging to the trace-mesh), but the degrees of
  freedom are maintained on the trace mesh.
  %%
\kitem{chain} contains the link to the other parts if this instance
  forms part of a chain of basis functions. This is implemented as
  doubly linked list. In the standard case the list-node just points
  back to itself. Compare \secref{S:bfcts_chains} and
  \secref{S:chain_impl}.
  %%
\kitem{unchained} points to a copy of the basis function structure
  which is unaware of being part of a direct sum. In the standard case
  \code{unchained} just points back to the same basis function
  structure it is part of.
  %%
\kitem{INIT\_ELEMENT\_DECL} is used for the initialization of element
  dependent local finite element spaces; this is needed, e.g. to
  define basis functions which depend on the element geometry like
  discretizations of the $H(\text{div})$, or for some more complicated
  discretizations for the Stokes-problem which, e.g., make use of
  edge- and face-bubbles. See \secref{S:init_element}.
  %%
\kitem{phi} vector of function pointers for the evaluation of local
  basis functions in barycentric coordinates;
  \begin{center}
    \code{(*bfcts->phi[i])(lambda, bfcts)}
  \end{center}
  returns the value $\pbar^\code{i}(\lambda)$ of the \code{i}-th basis
  function at \code{lambda} for $0 \le \code{i} <
  \code{n\_bas\_fcts}$.
  We advocate the use of the \code{PHI()}-macro instead:
  \begin{center}
    \code{PHI(bfcts, i, lambda)}
  \end{center}
  %%
\kitem{grd\_phi} vector of function pointers for the evaluation of
  gradients of the basis functions in barycentric coordinates;
  \begin{center}
    \code{(*bfcts->grd\_phi[i])(lambda)}
  \end{center}
  returns a pointer to a vector of length \code{N\_LAMBDA}
  containing all first derivatives (with respect to the
  barycentric coordinates) of the \code{i}--th basis function at
  \code{lambda}, i.e.
  \code{(*grd\_phi[i])(lambda)[k]}$ =
  \pbar^\code{i}_{,\lambda_{\code{k}}} (\lambda)$
  for $0 \le \code{k} \le d$,
  $0 \le \code{i} < \code{n\_bas\_fcts}$; this vector is
  overwritten on the next call of \code{(*grd\_phi[i])()}.
  We advocate the use of the \code{GRD\_PHI()}-macro instead:
  \begin{center}
    \code{GRD\_PHI(bfcts, i, lambda)}
  \end{center}
  %%
\kitem{D2\_phi} vector of function pointers for the evaluation of
  second derivatives of the basis functions in barycentric
  coordinates;
  \begin{center}
    \code{(*bfcts->D2\_phi[i])(lambda, bfcts)}
  \end{center}
  returns a pointer to a $\code{N\_LAMBDA}\times\code{N\_LAMBDA}$
  matrix containing all second derivatives (with respect to the
  barycentric coordinates) of the \code{i}--th basis function at
  \code{lambda}, i.e.
  \code{(*D2\_phi[i])(lambda)[k][l]}$ =
  \pbar^\code{i}_{,\lambda_{\code{k}} \lambda_{\code{l}}} (\lambda)$
  $0 \le \code{k},\code{l} \le d$,
  $0 \le \code{i} < \code{n\_bas\_fcts}$; this matrix is overwritten
  on the next call of \code{(*D2\_phi[i])()}.
  We advocate the use of the \code{D2\_PHI()}-macro instead:
  \begin{center}
    \code{D2\_PHI(bfcts, i, lambda)}
  \end{center}
  %%
\kitem{D3\_PHI(), D4\_PHI()} These do similar things as the other
  hooks, however, they need not be present in a specific
  \verb|BAS_FCTS| implementation.
  %%
\kitem{PHI\_D()} the directional part of the basis functions if
  \code{rdim == }\DOW. \ALBERTA always factors vector-valued basis
  functions into a scalar factor times a directional part, so the
  actual value of the \code{i}-th basis functions has to be calculated
  as
  \begin{lstlisting}
    REAL_D value;
    AXEY_DOW(PHI(bfcts, i, lambda), PHI_D(bfcts, i, lambda), value);
  \end{lstlisting}
  Expanding all the macros and inline functions, the above is
  equivalent to
  \begin{lstlisting}
    REAL_D value;
    const REAL *vector = bfcts->phi_d[i](lambda, bfcts);
    REAL scalar = bfctgs->phi[i](lambda, bfcts);
    int k;

    for (k = 0; k < DIM_OF_WORLD; k++) {
      value][k] = scalar * vector[k];
    }
  \end{lstlisting}
  Note that \ALBERTA never caches the directional part of
  vector-valued basis functions in its \verb|QUAD_FAST| or other
  quadrature caches; it is assumed that it changes on an
  element-to-element basis. Vector-valued basis functions and the
  associated support functions are discussed in further detail below
  in Section \ref{S:vector_bfcts}.
  %%
\kitem{GRD\_PHI\_D()} The gradient of the directional part of a
  vector-valued basis-function instance. Note that \verb|GRD_PHI_D|
  may be empty if \verb|dir_pw_const == true|.
  %%
\kitem{D2\_PHI\_D()} The second derivative of the directional part of
  a vector-valued basis function instance. Note that \verb|D2_PHI_D|
  may be empty if \verb|dir_pw_const == true|.
  %%

  \hrulefill
  %%
\kitem{dir\_pw\_const} if this is set to \code{true} then the
  directional part of the vector valued \verb|BAS_FCTS|-instance (i.e.
  \code{rdim == }\DOW) is constant on each mesh element (e.g. the
  normal to a face on affine linear elements). If this is the case
  then the computation of the derivatives of the basis functions is
  drastically simplified. See below in Section \ref{S:vector_bfcts}.
  %%
\hyperitem{BAS_FCTS:trace_bas_fcts}{trace\_bas\_fcts} A pointer to a
  basis function structure describing the trace of the current
  basis-function set on the boundaries of the reference element. In
  the case of Lagrange elements, this is again a Lagrange space of the
  same degree, but one dimension lower. The trace-spaces form a chain,
  which finally is terminated by a dimensions-$0$ ``dummy'' basis
  function set.

  The trace space may, of course, be of other nature than the bulk
  basis function set. An element bubble, for instance, already has
  zero trace on the boundary. The trace-space of face-bubbles will be
  a \DOW-valued element bubble, pointing in direction of the normal on
  the lower-dimensional element and so on.

  The trace spaces play a role when boundary integrals involving basis
  functions are computed (see, for instance,
  \hyperref[T:BNDRY_OPERATOR_INFO]{\code{BNDRY\_OPERATOR\_INFO}}
  structure -- \secref{S:matrix_assemblage_scalar} -- or
  \hyperref[S:neumann_bound]{\code{bndry\_L2\_scp\_fct\_bas()}}).
  They are also used to define global traces of finite element
  function in the context of
  \hyperref[S:tracemesh_implementation]{trace meshes}, see
  \secref{S:tracemesh_implementation} and \secref{book:S:tracemeshes}
  on page \pageref{book:S:tracemeshes}.
  %% 
  \exampleref{Ex:int1_2d} shows the use of the trace-space in the
  context of an interpolation routine for linear basis functions.
  %%
\hyperitem{BAS_FCTS:trace_dof_map}{trace\_dof\_map[][][][]} The
  mapping of the local degrees of freedom from the trace-set to the
  bulk-set of local basis functions, for each co-dimension $1$
  face-simplex (``wall''). To be more concrete:
  \begin{itemize}
  \item $\code{dim} < 3$:
    \bv\begin{lstlisting}
   trace_dof_map[0][0][wall][trace_dof] == bulk_dof
\end{lstlisting}\ev
  \item $\code{dim} = 3$:
    \bv\begin{lstlisting}
   trace_dof_map[type > 0][orient < 0][wall][trace_dof] == bulk_dof
\end{lstlisting}\ev
    In this context, \code{type} and \code{orient} denote the
    respective components of the \hyperref[T:EL_INFO]{\code{EL\_INFO}}
    structure, compare also the conceptual discussion of trace-meshes
    in \secref{book:S:tracemeshes} on page
    \pageref{book:S:tracemeshes}.
  \end{itemize}
  %%
\hyperitem{BAS_FCTS:n_trace_bas_fcts}{n\_trace\_bas\_fcts[]} The
  dimension of the local trace-space, for each co-dimension $1$
  face-simplex. In the standard case, the trace-space has the same
  dimension on each wall, but in the context of
  \hyperref[S:init_element]{per-element initializers} (see
  \secref{S:init_element}) the dimension may vary from wall to wall.
  %%
\kitem{get\_dof\_indices(result, el, admin, self)}
  \idx{get_dof_indices()@{\code{get\_dof\_indices()}!entry in
      \code{BAS\_FCTS} structure}}
  pointer to a function which connects the set of local basis
  functions with its global DOFs (an implementation of the function
  $j_S$ in Section \ref{book:S:eval_fe});  \begin{description}
  \item[Parameters]~\hfill 

   \begin{descr}
    \kitem{DOF *result} Storage for the result; \code{result} must be
      the base-address of an array to \DOFs with at least
      \code{n\_bas\_fcts} elements, or \nil, in which case the result
      is returned in a statically allocated storage area which is
      overwritten on the next call to \verb|get_dof_indices()|.  On
      return \verb|result[i]| stores the global \DOF associated to the
      \code{i}--th basis function.
    \kitem{const EL *el} the current mesh-element;
    \kitem{const DOF\_ADMIN *admin} the \DOF-admin for the
      corresponding finite element space;
    \kitem{const BAS\_FCTS *self} a pointer to the current basis
      function instance; might be used if the set of local basis
      functions depends on the mesh element.
    \end{descr}
  \item[Return Value] A pointer to a \verb|const EL_DOF_VEC| element
    vector (see page \pageref{T:EL_DOF_VEC}) if
    \code{(result == NULL)} or \nil if \code{(result != NULL)}. The
    \code{vec} component of the \verb|EL_DOF_VEC| contains the data
    which otherwise would have been stored in \code{result}. The
    contents of the return value is overwritten on the next call to
    \verb|get_dof_indices()|. If the \verb|BAS_FCTGS|-instance forms
    part of a chain of basis functions (see \secref{S:bfcts_chains}),
    then only the \DOFs associated to this instance are computed. To
    get the \DOFs for all parts of the direct sum the global function
    \verb|get_dof_indices()| has to be called, see Sections
    \ref{S:fillgetelvec} and \ref{S:chain_impl}.
  \end{description}
  To reduce the potential of coding errors we advocate the use of the
  \verb|GET_DOF_INDICES()| macro:
  \begin{center}
    \verb|GET_DOF_INDICES(bfcts, result, el, admin)|
  \end{center}
  instead of calling
  \begin{center}
    \verb|bfcts->get_dof_indices(result, el, admin, bfcts)|
  \end{center}
  directly.
  %%

  \hrulefill
  %%
\kitem{get\_bound(bndry\_bits, el\_info, self)}
  \idx{get_bound()@{\code{get\_bound()}!entry in \code{BAS\_FCTS}
      structure}}
  pointer to a function which fills a vector with the boundary types
  of the basis functions.
  \begin{compatibility}
    In contrast to all previous versions of \ALBERTA the boundary-type
    of a given basis function is a bit-mask, and not a mere number.
    Each bit corresponds to the number that has been assigned to a
    given boundary segment in the macro-triangulation. Basis-functions
    tied to vertex-\DOFs, e.g., may belong to different boundary
    segments in which case the bit-mask may contain more than one bit
    set. This is further discussed in Section \ref{S:boundary}.
  \end{compatibility}
  Otherwise the calling conventions are similar to the conventions for
  \verb|get_dof_indices()| (see above), there also exists a macro
  \verb|GET_BOUND(self, bndry_bits, el_info)|.  This function
  needs boundary information; thus, all routines using this function
  on the elements need the \code{FILL\_BOUND} flag during mesh
  traversal.
  \begin{description}
  \item[Parameters]~\hfill

    \begin{descr}
      \kitem{BNDRY\_BITS *bndry\_bits} storage for the reuslt;
      \kitem{const EL\_INFO *el\_info} the current elements's
        \verb|EL_INFO| descriptor;
      \kitem{const BAS\_FCTS *self} a pointer to the current basis
        function instance.
    \end{descr}
  \item[Return Value] a pointer to a statically allocated storage area
    of type \verb|const EL_BNDRY_VEC|, see page
    \pageref{T:EL_BNDRY_VEC}.  If the \verb|BAS_FCTS|-instance forms
    part of a direct sum, then only the boundary bit-masks associated
    with this instance are computed. To get the information for all
    parts of the direct sum the global function \verb|get_bound()| has
    to be called, see Sections \ref{S:fillgetelvec} and
    \ref{S:chain_impl}.
  \end{description}

  \hrulefill

  \kitem{interpol[\_d|\_dow](coeff, el\_info, wall, n, indices, f, ud, thisptr)}
  %% 
  \idx{interpol@{\code{interpol[\_d|\_dow]}! entry in \code{BAS\_FCTS} structure}}
  %% 
  When using \ALBERTA routines for the interpolation of
  \code{REAL[\_D]} valued functions the \code{interpol[\_d]} function
  pointer must be set (for example the calculation of Dirichlet
  boundary values by \code{dirichlet\_bound()} described in Section
  \ref{S:dirichlet_bound}):

  The \code{interpol}-hooks are function-pointers to functionws which
  performs the local interpolation of a \code{REAL[\_D]} valued
  function on an element. If this instance of a local basis function
  set forms part of a chain of basis functions (see
  \secref{S:bfcts_chains}, then it is possible to call the functions
  \verb|el_interpol()| -- respectively their \code{...\_d} and
  \code{...\_dow} variants -- to interpolate \code{f} onto the local
  functions space defined by the entire chain, see
  \secref{S:fillgetelvec}. The \verb|BAS_FCTS.interpol|-hook will only
  perform the interpolation for a single member of such a chain.

  \begin{description}
  \item[Parameters]~\hfill

    \begin{descr}
    \kitem{EL\_REAL\_VEC *coeff} Mandatory, storage for the result. The vector
      valued version need an \verb|EL_REAL_D_VEC| respetively an
      \verb|EL_REAL_VEC_D|.
      %% 
    \kitem{const EL\_INFO *el\_info} Element descriptor.
      %% 
    \kitem{wall} If the interpolation is to be performed over the
      boundary of the mesh, then this is the number of the wall in
      \verb|EL_INFO| to integrate over, in this case only the
      coefficients for the basis functions with non-zero trace on the
      respective wall will be computed.
      %% 
    \kitem{n} number of items in \code{indices}. If the interpolation
      is to be performed for all local \code{DOF}s, then $-1$ should
      be passed for \code{n}.
      %%
    \kitem{indices} for selective interpolation of only some of the
      local \code{DOF}s, indices may contain as many entries as
      indicated by \code{n}. The components of \code{indices} are then
      the local \code{DOF} number for which the coefficients should be
      computed. If \code{indices == NULL}, then the coefficients for
      all local basis functions will be computed.
      %%
    \kitem{REAL (*f)(const EL\_INFO *el\_info, const QUAD *quad, int
        iq, void *ud)}
      The application provided function to interpolate onto the local
      function space defined by this local basis function set. For
      simple Lagrange elements \code{QUAD} will just be a ``lumping''
      quadrature rule, with ``quadrature'' nodes on the Lagrange nodes
      of the basis function set. But interpolation may in fact require
      larger efforts, in which case \code{QUAD} may be a ``real''
      quadrature rule. If the interpolation is to be taken over a
      boundary segment, then \code{QUAD} will be a co-dimension $1$
      quadrature rule, see also \secref{S:QUAD}.
      %%

      For interpolation of \DOW-valued basis functions \code{f} must
      be a function pointer of the format
      \bv\begin{lstlisting}
        const REAL *(*f)(REAL_D result, const EL_INFO *el_info, const QUAD *quad, int iq, void *ud);
      \end{lstlisting}\ev
      %%
    \kitem{ud} Application data-pointer, forwarded to \code{f} as last
      argument.
      %%
    \kitem{thisptr} A pointer to the local basis function set.
      To reduce the potential of coding errors we advocate the use of the
      \verb|INTERPOL()| macro:
      \begin{center}
        \verb|INTERPOL(bfcts, coeff, el_info, wall, n, indices, f, ud)|
      \end{center}
      instead of calling
      \begin{center}
        \verb|bfcts->interpol(coeff, el_info, wall, n, indices, f, ud, bfcts)|
      \end{center}
      directly. Likewise for \verb|INTERPOL_D()| and \verb|INTERPOL_DOW|.
    \end{descr}
  \item[Return Value] \code{void}.
  \end{description}

  \hrulefill

  \kitem{const EL\_INT\_VEC *}
  \kitem{(*get\_int\_vec)(int res[], const EL *el, const DOF\_INT\_VEC *dv)}
  \kitem{const EL\_REAL\_VEC *}
  \kitem{(*get\_real\_vec)(REAL res[], const EL *el, const DOF\_REAL\_VEC *dv)}
  \kitem{const EL\_REAL\_D\_VEC *}
  \kitem{(*get\_real\_d\_vec)(REAL\_D res[], const EL *el, const DOF\_REAL\_D\_VEC *dv)}
  \kitem{const EL\_REAL\_VEC\_D *}
  \kitem{(*get\_real\_vec\_d)(REAL res[], const EL *el, const DOF\_REAL\_VEC\_D *dv)}
  \kitem{const EL\_UCHAR\_VEC *}
  \kitem{(*get\_uchar\_vec)(U\_CHAR res[], const EL *el, const DOF\_UCHAR\_VEC *dv)}
  \kitem{const EL\_SCHAR\_VEC *}
  \kitem{(*get\_schar\_vec)(S\_CHAR res[], const EL *el, const DOF\_SCHAR\_VEC *dv)}
  \kitem{const EL\_PTR\_VEC *}
  \kitem{(*get\_ptr\_vec)(void *res[], const EL *el, const DOF\_PTR\_VEC *dv)}

  These are pointers to functions which fills a local per-element
  coefficient vector with values of a \code{DOF\_*\_VEC} at the DOFs
  of the basis functions. The calling convention is much the same as
  for the \code{get\_dof\_indices()}- and \code{get\_bound()}-hooks.
  Also, in the context of chains of basis functions (see
  \secref{S:bfcts_chains} below) it should be noted that these
  function hooks only work on a single component of that chain.
  However, each of the hooks has global function as counter-part which
  does the job for the entire chain, see \secref{S:fillgetelvec}.
  \begin{compatibility}
    The calling convention has changed with respect to previous
    versions of \ALBERTA. In particular, there are now dedicated
    structures for storing local per-element coefficient vectors.
  \end{compatibility}
  Note that the \code{get\_real\_vec\_d}-hook accepts a \emph{scalar}
  \code{res}-argument of type \code{REAL} and returns a
  \code{EL\_REAL\_VEC\_D} because in the context of vector-valued
  basis functions the coefficient vector are scalars, see bewlow
  \secref{S:vector_bfcts}.

  A detailed description of the parameters is only given for the
  \code{get\_real\_vec()}-hook. The others work similar.
  \begin{description}
  \item[Parameters] ~\hfill

    \begin{descr}
      \kitem{res} An optional argument to store the coefficients in.
      \code{res} maybe \nil, in which case the return value if a
      pointer to \code{dv->vec\_loc}. If \code{res} is non-\nil, then
      the return value of this function is \nil.
      \begin{compatibility}
        This implies that it is safe to call
        \code{get\_real\_vec(NULL, \dots)} repeatedly with
        \emph{different} \code{DOF\_REAL\_VEC} instances, since the
        storage area for the return value is now tied to the
        argument \code{dv}, reducing the potential for coding errors.

        On the other hand, previous versions were returning a pointer
        to the argument \code{res}, if that was non-\nil. Of course,
        that can no longer work, because the return value is now a
        fully-fledged element vector, while the \code{res}-argument is
        just a flat \code{C}-array.
      \end{compatibility}
      %%
      \kitem{el} The element to compute the local coefficients for.
      %%
      \kitem{dv} The global coefficient vector to fetch the
      data-values from.
    \end{descr}
  \item[Return Value] \nil if (\code{res != \nil}), and
    \code{dv->vec\_loc} otherwise.
  \end{description}

  \hrulefill

  \kitem{void  (*real\_refine\_inter)(DOF\_REAL\_VEC *, RC\_LIST\_EL *rcl, int n)}
  %%
  \kitem{void  (*real\_coarse\_inter)(DOF\_REAL\_VEC *, RC\_LIST\_EL *rcl, int n)}
  %%
  \kitem{void  (*real\_coarse\_restr)(DOF\_REAL\_VEC *, RC\_LIST\_EL *rcl, int n)}
  %%
  \kitem{void  (*real\_d\_refine\_inter)(DOF\_REAL\_D\_VEC *, RC\_LIST\_EL *rcl, int n)}
  %%
  \kitem{void  (*real\_d\_coarse\_inter)(DOF\_REAL\_D\_VEC *, RC\_LIST\_EL *rcl, int n)}
  %%
  \kitem{void  (*real\_d\_coarse\_restr)(DOF\_REAL\_D\_VEC *, RC\_LIST\_EL *rcl, int n)}
  %%  
  \kitem{void  (*real\_refine\_inter\_d)(DOF\_REAL\_VEC\_D *, RC\_LIST\_EL *rcl, int n)}
  %%
  \kitem{void  (*real\_coarse\_inter\_d)(DOF\_REAL\_VEC\_D *, RC\_LIST\_EL *rcl, int n)}
  %%
  \kitem{void  (*real\_coarse\_restr\_d)(DOF\_REAL\_VEC\_D *, RC\_LIST\_EL *rcl, int n)}
  %%
  ~\hfill
  Since the interpolation of finite element functions during
  refinement and coarsening, as well as the restriction of functionals
  during coarsening, strongly depend on the basis functions and its
  DOFs (compare Section \ref{book:S:inter_restrict}), pointers for
  functions which perform those operations can be stored at above
  function pointers. Not all basis-function implementations may come
  with a full set of interpolation respectively restriction routines.
  
  \emph{Note also that these function-pointers are not used
    automatically; it is the responsibility of the application program
    to hook them into the global \code{DOF\_REAL\_[\_D]\_VEC[\_D]}
    coefficient vectors, only then \ALBERTA will make use of these
    function during mesh adaptation.}
  
  To give some aid in performing this job in the context of the
  more-complicated direct-sums of finite element spaces, there are
  global functions
  %%
  \code{set\_refine\_inter[\_d|\_dow]()},
  \code{set\_coarse\_inter[\_d|\_dow]()},
  \code{set\_coarse\_resrt[\_d|\_dow]()}
  %%
  which do this job for an entire chain of coefficient vectors (i.e.
  hook the corresponding function from the relevant
  \code{BAS\_FCTS}-component into the hook in the hook of the
  \code{DOF}-vector instance). See also \secref{S:chain_impl}.

  In Section \ref{S:linear-fe} and \ref{S:quadratic-fe} examples for
  the implementation of those functions are given.

  Functionally, the three different flavours have the following meaning:
  \begin{descr}
    \kitem{real[\_d]\_refine\_inter[\_d]} pointer to a function for
    interpolating a \code{REAL[\_D]} valued function during
    refinement; i.e. for interpolating the
    \code{DOF\_REAL[\_D]\_VEC[\_D]} vector \code{vec} on the
    refinement patch \code{rcl} onto the finer grid; information about
    all parents of the refinement patch is accessible in the vector
    \code{rcl} of length \code{n}.
    %%
    \kitem{real[\_d]\_coarse\_inter[\_d]} pointer to a function for
    interpolating a \code{REAL[\_D]} valued function during
    coarsening; i.e. for interpolating the \code{DOF\_REAL[\_D]\_VEC}
    vector \code{vec} on the coarsening patch \code{rcl} onto the
    coarser grid; information about all parents of the refinement
    patch is accessible in the vector \code{rcl} of length \code{n}.
    %%
    \kitem{real[\_d]\_coarse\_restr[\_d]} pointer to a function for
    restriction of \code{REAL[\_D]} valued linear functionals during
    coarsening; i.e. for restricting the \code{DOF\_REAL[\_D]\_VEC}
    vector \code{vec} on the coarsening patch \code{rcl} onto the
    coarser grid; information about all parents of the refinement
    patch is accessible in the vector \code{rcl} of length \code{n}.    
  \end{descr}
  
  \hrulefill

  \kitem{ext\_data} A \code{void}-pointer to implementation dependent
  data tied to a specific basis-functions implementation. For standard
  Lagrange basis-functions, this pointer gives access to the local
  Lagrange nodes (``local'' meaning their barycentric coordinates) via the macro
  \bv\begin{lstlisting}
    const REAL_B *nodes = LAGRANGE_NODES(bfcts);
  \end{lstlisting}\ev
  Although this is defined as a macro in \code{alberta.h}, an
  application program must not assume that this macro will not change
  in future versions of the tool-box. The ordering of data behind the
  \code{ext\_data} pointer is opaque, nothing should be assumed about it.

  \hrulefill

\end{descr}

\begin{remark}\label{R:get_int_vec}
 The access to local element vectors via the 
\code{get\_*\_vec()} routines can also be done in a standard way by using the
\code{get\_dof\_indices()} function which must be supplied; if
some of the \code{get\_*\_vec()} are pointer to \nil, \ALBERTA fills in
pointers to some standard functions using \code{get\_dof\_indices()}. But a
specialized function may be faster. An example of such a standard routine
is:
\bv\begin{lstlisting}[name=default_get_int_vec,label=C:default_get_int_vec]
  const EL_INT_VEC *
  default_get_int_vec(int *vec, const EL *el, const DOF_INT_VEC *dof_vec)
  {
    FUNCNAME("get_int_vec");
    int *rvec = vec == NULL ? dof_vec->vec_loc->vec : vec;
    const BAS_FCTS *bas_fcts = dof_vec->fe_space->bas_fcts;
    int n_bas_fcts = bas_fcts->n_bas_fcts;
    DOF index[n_bas_fcts];
    int  i;

    GET_DOF_INDICES(dof_vec->fe_space->bas_fcts,
		    el, dof_vec->fe_space->admin, index);

    for (i = 0; i < n_bas_fcts; i++) {
      rvec[i] = dof_vec->vec[index[i]];
    }

    return vec ? NULL : dof_vec->vec_loc;
  }
\end{lstlisting}\ev
A specialized implementation for linear finite elements e.g. is more
efficient:
\idx{get_int_vec()@{\code{get\_int\_vec()}}!for linear elements}
\bv\begin{lstlisting}
const EL_INT_VEC *
get_int_vec(int *ivec, const EL *el, const DOF_INT_VEC *vec)
{
  FUNCNAME("get_int_vec");
  int            i, n0;
  int            *v = vec->vec;
  int            *rvec = ivec ? ivec : vec->vec_loc->vec;
  DOF            **dof = el->dof;

  n0 = vec->fe_space->admin->n0_dof[VERTEX];

  for (i = 0; i < N_VERTICES; i++) {
    rvec[i] = v[dof[i][n0]];
  }

  return vec ? rvec : vec->vec_loc;
}
\end{lstlisting}\ev
\end{remark}

Any kind of basis functions can be implemented by filling the
above described structure for basis functions. All non-optional
entries have to be defined. Since in the functions for 
reading and writing of meshes, the basis functions are identified
by their names, all used basis functions have to be registered
before using these functions. All Lagrange
finite elements described below are already registered, with
names \code{"lagrange1\_1d"} to \code{"lagrange4\_3d"}. The discontinuous 
polynomial finite elements are registered with \code{"disc\_lagrange0\_1d"} to 
\code{"disc\_lagrange2\_3d"}. Newly defined 
basis functions must use different names.
%
\fdx{new_bas_fcts@{\code{new\_bas\_fcts}}}
\bv\begin{lstlisting}
int new_bas_fcts(const BAS_FCTS *bas_fcts);
\end{lstlisting}\ev
Description:
\begin{descr}
  \kitem{new\_bas\_fcts(bas\_fcts)} puts the new set of basis
  functions \code{bas\_fcts} to an internal list of all used basis
  functions; different sets of basis functions are identified by their
  \code{name}; thus, the member \code{name} of \code{bas\_fcts} must
  be a string with positive length holding a description; if an
  existing set of basis functions with the same name is found, the
  program stops with an error; if the entries \code{phi},
  \code{grd\_phi}, \code{get\_dof\_indices}, and \code{get\_bound} are
  not set, this also result in an error and the program stops.
\end{descr}
%
Basis functions can be accessed from that list by
\fdx{get_bas_fcts@{\code{get\_bas\_fcts}}}
\bv\begin{lstlisting}[label=F:get_bas_fcts_fct]
const BAS_FCTS *get_bas_fcts(const char *name)
\end{lstlisting}\ev
Description:
\begin{descr}
\kitem{get\_bas\_fcts(name)} looks for a set of basis functions
       with name \code{name} in the internal list of all registered
       basis functions; if such a set is found, the return value
       is a pointer to the corresponding \code{BAS\_FCTS} structure,
       otherwise the return value is \nil.
\end{descr}
%
Lagrange elements can be accessed by
a call of \code{get\_lagrange()}, see Section \ref{S:get_lagrange},
discontinuous polynomial elements by \code{get\_discontinuous\_lagrange()}, 
see Section \ref{S:get_disc_lagrange}.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Vector-valued basis functions}
\label{S:vector_bfcts}

As a new feature, the current version of this finite element toolbox
contains support for vector-valued basis functions like edge- or
face-bubbles, or Raviart-Thomas elements. The actual implementation of
those basis functions has been moved to an add-on module
\verb|libalbas|, see \secref{S:fancy_bas_fcts} below.  An example for
the implementation of face-bubbles can be found in
%%
\bv\begin{verbatim}
albertadist/add_ons/lib_albas/src/wall_bubbles.c
\end{verbatim}\ev

The current implementation assumes that it is efficient to factor
vector-valued basis functions into a scalar part which does \emph{not}
depend on the element geometry and a vector-valued part -- actually:
\DOW-valued -- which depends on the element geometry. This is
reflected by the \verb|BAS_FTCS| data-structure: the
\verb|BAS_FCTS.phi| jump-tables correspond to the scalar factor, while
the vector-valued part is stored in the jump-table
\verb|BAS_FCTS.phi_d|, and analogously for the jump-tables for the
derivatives.

\idx{Per-element initializers!vector-valued basis functions}
\idx{init_element()@{\code{init\_element()}}!vector-valued basis functions}
%%
Often vector-valued basis functions will carry a per-element
initializer, a function pointer
\verb|BAS_FCTS.init_element(el_info, self)|, which is invoked with the
current \verb|EL_INFO|-descriptor and the basis-function instance
itself. This function-hook can be used to update geometry information
like coordinates or wall-normals. In this context, the component
\verb|BAS_FCTS.fill_flags| is also of particular importance, it
contains the collection of mesh-traversal flags (see
\secref{S:traverse}) which are needed in order for the
\verb|init_element()|-hook to do its job properly. See also
\secref{S:init_element}.

If the evaluation of basis functions is computationally costly, then
it is of special importance to cache values of basis functions (and
their derivatives) at quadrature points (see \secref{S:QUAD_FAST}).
For vector-valued basis-functions, these caches are only maintained
for the scalar factor, as the vector-valued factor is assumed to vary
from element to element anyway. In order to simplify the
``recombination'' of the scalar- and the vector-factor, the library
provides three functions to perform this task (arguably this part of
the documentation would belong to \secref{S:QUAD_FAST}):

\bv\begin{lstlisting}
const REAL_D *const*get_quad_fast_phi_dow(const QUAD_FAST *cache);
const REAL_DB *const*get_quad_fast_grd_phi_dow(const QUAD_FAST *cache);
const REAL_DBB *const*get_quad_fast_D2_phi_dow(const QUAD_FAST *cache);
\end{lstlisting}\ev

These three function take a pointer to a properly initialized
\verb|QUAD_FAST|-structure as returned by \verb|get_quad_fast()| and
return arrays containing the values of the products of the scalar- and
vector-part of the basis-functions. This way an application can call
those functions, and then use the returned arrays in the same way it
used the components of the \verb|QUAD_FAST|-structure. The ordering of
indices is also the same, e.g.
\begin{example}
\bv\begin{lstlisting}
  const REAL_D *const*phi_d = get_quad_fast_phi_dow(qfast);

  for (iq = 0; iq < qfast->n_points; iq++) {
    for (b = 0; b < qfast->n_bas_fcts; b++) {
      do_something_fct(phi_d[iq][b]);
    }
  }
\end{lstlisting}\ev
\end{example}

It is also worth noting that the coefficient-vectors for
finite-element functions based on vector-valued basis-functions
contain scalars (the vector-nature of the finite element space is
induced by the vector-nature of the values of the basis functions,
thus the coefficients for the basis functions are scalars). \ALBERTA's
assemble infra-structure (see \secref{S:ass_tools}) has full support
for assembling linear systems based on vector-valued basis functions,
and to freely pair scalar- and \DOW-valued finite element spaces.
Actually, the support-functions for the assembling of the discrete
systems use up most of the compilation time during the installation of
the \ALBERTA-package.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Chains of basis function sets}
\label{S:bfcts_chains}

The current version of this finite element toolbox has support for
direct sums of finite-element spaces. Each component of such a sum is
defined by a set of local basis functions which is part of a chain of
basis functions. The chain-connectivity is implemented as a doubly
linked cyclic list, the corresponding list-link can be found in the
\verb|BAS_FCTS.chain| component. In order to chain single basis
function implementations together the support function
\verb|chain_bas_fcts()| has to be called:

\fdx{chain_bas_fcts@{\code{chain\_bas\_fcts}}}%
\bv\begin{lstlisting}
BAS_FCTS *chain_bas_fcts(const BAS_FCTS *head, BAS_FCTS *tail);
\end{lstlisting}\ev

\begin{descr}
  \kitem{chain\_bas\_fcts(head, tail)} Clone the set of
  basis-functions specified as \code{head}, if \code{tail != NULL},
  then add the copy of \code{head} as head to the list specified by
  \code{tail}.  The original instance of \code{head} is hooked into
  the copy using the \verb|BAS_FCTS.unchained| pointer; it remains a
  single, un-chained \verb|BAS_FCTS|-instance.  The chain-connectivity
  is implemented using the doubly-linked list-node \verb|BAS_FCTS.chain|.
  
  The chained basis functions will be given the name
  \begin{center}
    \verb|"HEAD_NAME#TAIL_NAME"|
  \end{center}
  
  During the construction of the name any \verb|"_Xd"|-suffixes are
  discarded in order not to make the name too complicated. At the
  end of the name-chain an appropriate \verb|_Xd|-suffix is added.
  The \verb|_Xd|-suffixes are used only for debugging, they are
  ignored everywhere else.
  
  A basis-function chain is cyclic. The trace-spaces of the given sets
  of basis functions are chained together accordingly.
  %%
  \idx{Per-element initializers!basis function chains}
  \idx{init_element()@{\code{init\_element()}}!basis function chains}
  \idx{Per-element initializers!direct sums of functions spaces}
  \idx{init_element()@{\code{init\_element()}}!direct sums of function spaces}
  %% 
  If any part of
  the chain needs a per-element initialization (see
  \secref{S:init_element}), then \code{chain\_bas\_fcts()} assigns the
  resulting chain a special \code{init\_element()} hook which follows
  the convention described in \secref{S:init_element} and which calls
  the per-element initializers of each member of the chain in turn.
  
  Note: this function does \emph{not} call \verb|new_bas_fcts()|; the
  caller has to do so after constructing the desired chain.

  More about direct sums of finite-element spaces should be found in
  \secref{S:chain_impl}. It is generally a bad idea to chain sets of
  basis functions together which are not linearly independent from
  each other \dots.

  \begin{example}
    Forming the velocity-part of the ``Mini''-element, and looping
    over the generated two-element chain, printing the name of the
    unchained instances. Note that the \verb|CHAIN_DO()|-macro rolls
    over the entire list, which is cyclic. So inside the loop
    \verb|bas_fcts| points to different components of the chain, after
    the loop \verb|bas_fcts| again points to the first instance of the
    chain.

    \emph{This is a certain potential for bugs when jumping out of the
      loop, using, e.g., a \code{break} statement. This should be
      avoided. \code{continue} statements should also be avoided
      because the trailing \code{CHAIN\_WHILE()}-part increments the
      list-pointer.}

    \bv\begin{lstlisting}
      lagrange = get_lagrange(dim, degree);
      bubble   = get_bas_fcts(dim, "Bubble");
      bas_fcts = chain_bas_fcts(lagrange, chain_bas_fcts(bubble, NULL));
      old_fcts = new_bas_fcts(bas_fcts);
      
      if (old_fcts != NULL) {
        MSG("Overriding old definition for \"%s\"\n", old_fcts->name);
      }

      MSG("New name: \"%s\"\n", bas_fcts->name);
      CHAIN_DO(bas_fcts, BAS_FCTS) {
        MSG("Rolling component name: \"%s\"\n", bas_fcts->name);
        MSG("Component name: \"%s\"\n", bas_fcts->unchained->name);
      } CHAIN_WHILE(bas_fcts, BAS_FCTS);
      MSG("Again, the name of the entire chain: \"%s\"\n", bas_fcts->name);
    \end{lstlisting}\ev
  \end{example}
\end{descr}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Lagrange finite elements}
\label{S:Lagrange_elements}

\ALBERTA provides Lagrange finite elements up to order four
which are described in the following sections. Lagrange finite
elements are given by $\Pbar = \P_p(\Sbar)$ (polynomials of degree $p
\in \N$ on $\Sbar$) and they are globally continuous.  They are
uniquely determined by the values at the associated Lagrange nodes
$\big\{x_i\big\}$.  The Lagrange basis functions $\big\{\phi_i\big\}$
satisfy
\[
\phi_i(x_j) = \delta_{ij} \qquad \mbox{for } i,j = 1,\ldots, N = \dim\ X_h.
\]
Now, consider the basis functions $\{\pbar^i\}_{i = 1}^m$ of $\Pbar$
with the associated Lagrange nodes $\{\lambda_i\}_{i = 1}^m$ given
in barycentric coordinates:
\[
\pbar^i(\lambda_j) = \delta_{ij} \qquad \mbox{for } i,j = 1,\ldots, m.
\]
Basis functions are located at the vertices, center, edges, or faces
of an element. The corresponding DOF is a vertex, center, edge,
or face DOF, respectively. The boundary type of a basis
function is the boundary type of the associated vertex (or edge or face).
Basis functions at the center are always \code{INTERIOR}.
Such boundary information is filled by the \code{get\_bound()} function
in the \code{BAS\_FCTS} structure and is straight forward.

The interpolation coefficient for a function $f$ for basis function $\pbar^i$
on element $S$ is the value of $f$ at the Lagrange node: $f(x(\lambda_i))$.
These coefficients are calculated by the \code{interpol[\_d]()} function
in the \code{BAS\_FCTS} structure. Examples for both functions
are given below for linear finite elements.

\subsubsection{Piecewise linear finite elements}\label{S:linear-fe}

Piecewise linear, continuous finite elements are uniquely defined by
their values at the vertices of the triangulation. On each element we
have \code{N\_VERTICES(dim)} basis functions which are the barycentric
coordinates of the element.  Thus, in 1d we have two, in 2d we have
three, and in 3d four basis functions for Lagrange elements of first
order; the basis functions and the corresponding Lagrange nodes in
barycentric coordinates are shown in Tables \ref{T:P1_1d},
\ref{T:P1_2d} and \ref{T:P1_3d}. The calculation of derivatives is
straight forward.
\begin{table}[htbp]
\begin{center}\begin{small}\begin{math}
\begin{array}{|rcl@{\qquad}|@{\qquad}l@{\qquad}|@{\qquad}rcl|}\hline
\multicolumn{2}{|c}{\mbox{function}} & & \mbox{position} &
\multicolumn{3}{c|}{\text{Lagrange node}}
 \\\hline\vphantom{\pbar^{\int}}
\pbar^\code{0}(\lambda) & = & \lambda_\code{0} &
\mbox{vertex \code{0}} & \lambda^\code{0} & = & (1,0)\\
\pbar^\code{1}(\lambda) & = & \lambda_\code{1} &
\mbox{vertex \code{1}} & \lambda^\code{1} & = & (0,1)\\
\hline\end{array}
\end{math}\end{small}\end{center}
\vskip -5mm
\caption{Local basis functions for linear finite elements in 1d.}
\label{T:P1_1d}
\end{table}
\begin{table}[htbp]
\begin{center}\begin{small}\begin{math}
\begin{array}{|rcl@{\qquad}|@{\qquad}l@{\qquad}|@{\qquad}rcl|}\hline
\multicolumn{2}{|l}{\mbox{function}} & & \mbox{position} &
\multicolumn{3}{l|}{\mbox{Lagrange node}}
 \\\hline\vphantom{\pbar^{\int}}
\pbar^\code{0}(\lambda) & = & \lambda_\code{0} &
\mbox{vertex \code{0}} & \lambda^\code{0} & = & (1,0,0)\\
\pbar^\code{1}(\lambda) & = & \lambda_\code{1} &
\mbox{vertex \code{1}} & \lambda^\code{1} & = & (0,1,0)\\
\pbar^\code{2}(\lambda) & = & \lambda_\code{2} &
\mbox{vertex \code{2}} & \lambda^\code{2} & = & (0,0,1)\\
\hline\end{array}
\end{math}\end{small}\end{center}
\vskip -5mm
\caption{Local basis functions for linear finite elements in 2d.}
\label{T:P1_2d}
\end{table}
\begin{table}[htbp]
\begin{center}\begin{small}\begin{math}%\displaystyle
\begin{array}{|rcl@{\qquad}|@{\qquad}l@{\qquad}|@{\qquad}rcl|}\hline
\multicolumn{2}{|l}{\mbox{function}} & & \mbox{position} &
\multicolumn{3}{l|}{\mbox{Lagrange node}}
 \\\hline\vphantom{\pbar^{\int}}
\pbar^\code{0}(\lambda) & = & \lambda_\code{0} &
\mbox{vertex \code{0}} & \lambda^\code{0} & = & (1,0,0,0)\\
\pbar^\code{1}(\lambda) & = & \lambda_\code{1} &
\mbox{vertex \code{1}} & \lambda^\code{1} & = & (0,1,0,0)\\
\pbar^\code{2}(\lambda) & = & \lambda_\code{2} &
\mbox{vertex \code{2}} & \lambda^\code{2} & = & (0,0,1,0)\\
\pbar^\code{3}(\lambda) & = & \lambda_\code{3} &
\mbox{vertex \code{3}} & \lambda^\code{3} & = & (0,0,0,1)\\
\hline\end{array}
\end{math}\end{small}\end{center}
\vskip -5mm
\caption{Local basis functions for linear finite elements in 3d.}
\label{T:P1_3d}
\end{table}
The global DOF index of the \code{i}--th basis functions on element
\code{el} is stored for linear finite elements at
\bv\begin{lstlisting}
el->dof[i][admin->n0_dof[VERTEX]]
\end{lstlisting}\ev
Setting \code{nv = admin->n0\_dof[VERTEX]} the associated DOFs are shown
in Figure~\ref{F:linear_dof}.

For linear finite elements we want to give examples
for the implementation of some routines in the corresponding 
\code{BAS\_FCTS} structure. 

\begin{example}[Accessing DOFs for piecewise linear finite elements]
The implementation of \code{get\_dof\_indices()} can be done in
as in the following 2d example code, compare Figure~\ref{F:linear_dof} and 
Remark~\ref{R:get_int_vec} with the implementation of 
the function \code{get\_int\_vec()}
for accessing a local element vector from a global \code{DOF\_INT\_VEC}
for piecewise linear finite elements.
%%
\idx{get_dof_indices()@{\code{get\_dof\_indices()}}!for linear elements}
\bv\begin{lstlisting}[name=get_dof_indices1_2d,label=get_dof_indices1_2d]
  static const EL_DOFVEC *
  get_dof_indices1_2d(DOF *vec, const EL *el, const DOF_ADMIN *admin,
                      const BAS_FCTS *thisptr)
  {
    static DEF_EL_VEC_CONST(DOF, rvec_space, N_BAS_LAG_1_2d, N_BAS_LAG_1_2d);
    DOF *rvec = vec ? vec : rvec_space->vec;
    int n0, /* node, */ ibas;
    DOF **dofptr = el->dof, dof;

    /* node = admin->mesh->node[VERTEX]; */
    n0   = admin->n0_dof[VERTEX];
    for (ibas = 0; ibas < N_BAS_LAG_1_2D; ibas++) {
      dof = dofptr[/* node+ */ibas][n0];
      body;
    }

    return vec ? NULL : rvec_space;
  }
\end{lstlisting}\ev
\begin{figure}[hbpt]
\centerline{\includegraphics[scale=1.0]{EPS/linear_dof}}
\vspace{-2mm}
\caption[DOFs and local numbering of the basis functions for linear
elements]{DOFs and local numbering of the basis functions for linear
  elements in 2d and 3d.}
\label{F:linear_dof}
\end{figure}
\end{example}

\begin{example}[Accessing the boundary types of DOFs for piecewise linear 
                finite elements in 2d]
The \code{get\_bound()} function fills the \code{bound}
vector with the boundary type of the vertices, shown here for 2d:
%%
\begin{samepage}
\idx{get_bound()@{\code{get\_bound()}}!for linear elements}
\bv\begin{lstlisting}
static const EL_BNDRY_VEC *
get_bound1_2d(BNDRY_FLAGS *vec, const EL_INFO *el_info, const BAS_FCTS *thisptr)
{
  FUNCNAME("get_bound1_2d");
  static DEF_EL_VEC_CONST(BNDRY, rvec_space, N_BAS_LAG_1_2d, N_BAS_LAG_1_2d););
  BNDRY_FLAGS *rvec = vec ? vec : rvec_space->vec;
  int            i;

  DEBUG_TEST_FLAG(FILL_BOUND, el_info);

  for (i = 0; i < N_VERTICES_2D; i++) {
    BNDRY_FLAGS_CPY(rvec[i], el_info->vertex_bound[i]);
  }

  return vec ? NULL : rvec_space;
}
\end{lstlisting}\ev
\end{samepage}
\end{example}

\begin{example}[Interpolation for piecewise linear finite elements in 2d]%
\label{Ex:int1_2d}
For the interpolation \code{interpol()} routine we have to evaluate
the given function at the vertices. Assuming a 2d mesh, the
interpolation can be implemented as follows. Note the use of a
``lumping'' quadrature rule; the application supplied function
\code{f()} has the task to return its value at the quadrature point
\code{iq}.
%%
\idx{interpol() for linear elements}
\bv\begin{lstlisting}
static void interpol1_2d(EL_REAL_VEC *el_vec,
                         const EL_INFO *el_info, int wall,
                         int no, const int *b_no,
                         REAL (*f)(const EL_INFO *el_info,
                                   const QUAD *quad, int iq,
                                   void *ud),
                         void *f_data,
                         const BAS_FCTS *thisptr)
{
  FUNCNAME("interpol1_2d");
  LAGRANGE_DATA *ld = &lag_1_2d_data;
  REAL *rvec = el_vec->vec;
  const QUAD *lq;
  const int *trace_map;
  int i;

  DEBUG_TEST_EXIT(ld->lumping_quad != NULL,
                  "called for uninitialized Lagrange basis functions\n");

  if (wall < 0) {
    lq = ld->lumping_quad;
    trace_map = NULL;
  } else {
    int type = el_info->el_type > 0;
    int orient = el_info->orientation < 0;
    lq = &ld->trace_lumping_quad[type][orient][wall];
    trace_map = thisptr->trace_dof_map[type][orient][wall];
  }

  DEBUG_TEST_EXIT(!b_no || (no >= 0 && no <= lq->n_points),
                  "not for %d points\n", no);

  el_vec->n_components = thisptr->n_bas_fcts;

  if (b_no) {
    for (i = 0; i < no; i++) {
      int ib = wall < 0 ? b_no[i] : trace_map[b_no[i]];
      rvec[ib] = f(el_info, lq, b_no[i], f_data);
    }
  } else {
    for (i = 0; i < lq->n_points; i++) {
      int ib = wall < 0 ? i : trace_map[i];
      rvec[ib] = f(el_info, lq, i, f_data);
    }
  }
}
\end{lstlisting}\ev
\end{example}

\begin{example}[Interpolation and restriction routines
                for piecewise linear finite elements in 2d]
The implementation of functions for interpolation during refinement and
restriction of linear functionals during coarsening is very simple for
linear elements; we do not have to loop over the refinement patch since
only the vertices at the refinement/coarsening edge and the new DOF at
the midpoint are involved in this process. No interpolation during
coarsening has to be done since all values at the remaining
vertices stay the same; no function has to be defined.
\idx{real_refine_inter()@{\code{real\_refine\_inter()}}!for linear elements}
\bv\begin{lstlisting}
static void real_refine_inter1_2d(DOF_REAL_VEC *drv, RC_LIST_EL *list, int n)
{
  FUNCNAME("real_refine_inter1_2d");
  EL      *el;
  REAL    *vec = nil;
  DOF     dof_new, dof0, dof1;
  int     n0;

  if (n < 1) return;
  GET_DOF_VEC(vec, drv);
  n0 = drv->fe_space->admin->n0_dof[VERTEX];
  el = list->el_info.el;
  dof0 = el->dof[0][n0];           /* 1st endpoint of refinement edge */
  dof1 = el->dof[1][n0];           /* 2nd endpoint of refinement edge */
  dof_new = el->child[0]->dof[2][n0];  /*     newest vertex is dim==2 */
  vec[dof_new] = 0.5*(vec[dof0] + vec[dof1]);

  return;
}
\end{lstlisting}\ev
\idx{real_coarse_restr()@{\code{real\_coarse\_restr()}}!for linear elements}
\bv\begin{lstlisting}
static void real_coarse_restr1_2d(DOF_REAL_VEC *drv, RC_LIST_EL *list, int n)
{
  FUNCNAME("real_coarse_restr1_2d");
  EL      *el;
  REAL    *vec = nil;
  DOF     dof_new, dof0, dof1;
  int     n0;

  if (n < 1) return;
  GET_DOF_VEC(vec, drv);
  n0 = drv->fe_space->admin->n0_dof[VERTEX];
  el = list->el_info.el;
  dof0 = el->dof[0][n0];           /* 1st endpoint of refinement edge */
  dof1 = el->dof[1][n0];           /* 2nd endpoint of refinement edge */
  dof_new = el->child[0]->dof[2][n0];    /*   newest vertex is dim==2 */
  vec[dof0] += 0.5*vec[dof_new];
  vec[dof1] += 0.5*vec[dof_new];

  return;
}
\end{lstlisting}\ev
\end{example}

\subsubsection{Piecewise quadratic finite elements}\label{S:quadratic-fe}

Piecewise quadratic, continuous finite elements are uniquely defined by
their values at the vertices and the edges' midpoints (center in 1d)
of the triangulation.  In 1d we have three, in 2d we have six, and in
3d we have ten basis functions for Lagrange elements of second order;
the basis functions and the corresponding Lagrange nodes in
barycentric coordinates are shown in Tables \ref{T:P2_1d},
\ref{T:P2_2d}, and \ref{T:P2_3d}.
\begin{table}[htbp]
\begin{center}\begin{small}\begin{math}%\displaystyle
\begin{array}{|rcl@{\qquad}|@{\qquad}l@{\qquad}|@{\qquad}rcl|}\hline
\multicolumn{2}{|l}{\mbox{function}} & & \mbox{position} &
\multicolumn{3}{l|}{\mbox{Lagrange node}}
 \\\hline\vphantom{\pbar^{\int}}
\pbar^\code{0}(\lambda) & = & \lambda_\code{0} (2 \lambda_\code{0} - 1) &
\mbox{vertex \code{0}} & \lambda^\code{0} & = & (1,0)\\[1mm]
\pbar^\code{1}(\lambda) & = & \lambda_\code{1} (2 \lambda_\code{1} - 1) &
\mbox{vertex \code{1}} & \lambda^\code{1} & = & (0,1)\\[1mm]
\pbar^\code{2}(\lambda) & = & 4 \lambda_\code{0}\,\lambda_\code{1} &
\mbox{center}& \lambda^\code{2} & = & (\half,\half)\\[1mm]
\hline\end{array}
\end{math}\end{small}\end{center}
\vskip -5mm
\caption{Local basis functions for quadratic finite elements in 1d.}
\label{T:P2_1d}
\end{table}
\begin{table}[htbp]
\begin{center}\begin{small}\begin{math}%\displaystyle
\begin{array}{|rcl@{\qquad}|@{\qquad}l@{\qquad}|@{\qquad}rcl|}\hline
\multicolumn{2}{|l}{\mbox{function}} & & \mbox{position} &
\multicolumn{3}{l|}{\mbox{Lagrange node}}
 \\\hline\vphantom{\pbar^{\int}}
\pbar^\code{0}(\lambda) & = & \lambda_\code{0} (2 \lambda_\code{0} - 1) &
\mbox{vertex \code{0}} & \lambda^\code{0} & = & (1,0,0)\\[1mm]
\pbar^\code{1}(\lambda) & = & \lambda_\code{1} (2 \lambda_\code{1} - 1) &
\mbox{vertex \code{1}} & \lambda^\code{1} & = & (0,1,0)\\[1mm]
\pbar^\code{2}(\lambda) & = & \lambda_\code{2} (2 \lambda_\code{2} - 1) &
\mbox{vertex \code{2}} & \lambda^\code{2} & = & (0,0,1)\\[1mm]
\pbar^\code{3}(\lambda) & = & 4 \lambda_\code{1}\,\lambda_\code{2} &
\mbox{edge \code{0}}& \lambda^\code{3} & = & (0,\half,\half)\\[1mm]
\pbar^\code{4}(\lambda) & = & 4 \lambda_\code{2}\,\lambda_\code{0} &
\mbox{edge \code{1}} & \lambda^\code{4} & = & (\half,0,\half)\\[1mm]
\pbar^\code{5}(\lambda) & = & 4 \lambda_\code{0}\,\lambda_\code{1} &
\mbox{edge \code{2}}& \lambda^\code{5} & = & (\half,\half,0)\\[1mm]
\hline\end{array}
\end{math}\end{small}\end{center}
\vskip -5mm
\caption{Local basis functions for quadratic finite elements in 2d.}
\label{T:P2_2d}
\end{table}
\begin{table}[htbp]
\begin{center}\begin{small}\begin{math}%\displaystyle
\begin{array}{|rcl@{\qquad}|@{\qquad}l@{\qquad}|@{\qquad}rcl|}\hline
\multicolumn{2}{|l}{\mbox{function}} & & \mbox{position} &
\multicolumn{3}{l|}{\mbox{Lagrange node}}
 \\\hline\vphantom{\pbar^{\int}}
\pbar^\code{0}(\lambda) & = & \lambda_\code{0} (2 \lambda_\code{0} - 1)&
\mbox{vertex \code{0}}& \lambda^\code{0} & = & (1,0,0,0)\\[1mm]
\pbar^\code{1}(\lambda) & = & \lambda_\code{1} (2 \lambda_\code{1} - 1)&
\mbox{vertex \code{1}}& \lambda^\code{1} & = & (0,1,0,0)\\[1mm]
\pbar^\code{2}(\lambda) & = & \lambda_\code{2} (2 \lambda_\code{2} - 1)&
\mbox{vertex \code{2}}& \lambda^\code{2} & = & (0,0,1,0)\\[1mm]
\pbar^\code{3}(\lambda) & = & \lambda_\code{3} (2 \lambda_\code{3} - 1)&
\mbox{vertex \code{3}}& \lambda^\code{3} & = & (0,0,0,1)\\[1mm]
\pbar^\code{4}(\lambda) & = & 4 \lambda_\code{0}\,\lambda_\code{1} &
\mbox{edge \code{0}}& \lambda^\code{4} & = & (\half,\half,0,0)\\[1mm]
\pbar^\code{5}(\lambda) & = & 4 \lambda_\code{0}\,\lambda_\code{2} &
\mbox{edge \code{1}}& \lambda^\code{5} & = & (\half,0,\half,0)\\[1mm]
\pbar^\code{6}(\lambda) & = & 4  \lambda_\code{0}\,\lambda_\code{3} &
\mbox{edge \code{2}}& \lambda^\code{6} & = & (\half,0,0,\half)\\[1mm]
\pbar^\code{7}(\lambda) & = & 4 \lambda_\code{1}\,\lambda_\code{2} &
\mbox{edge \code{3}}& \lambda^\code{7} & = & (0,\half,\half,0)\\[1mm]
\pbar^\code{8}(\lambda) & = & 4  \lambda_\code{1}\,\lambda_\code{3} &
\mbox{edge \code{4}}& \lambda^\code{8} & = & (0,\half,0,\half)\\[1mm]
\pbar^\code{9}(\lambda) & = & 4  \lambda_\code{2}\,\lambda_\code{3} &
\mbox{edge \code{5}} & \lambda^\code{9} & = & (0,0,\half,\half)\\[1mm]
\hline\end{array}
\end{math}\end{small}\end{center}
\vskip -5mm
\caption{Local basis functions for quadratic finite elements in 3d.}
\label{T:P2_3d}
\end{table}
\begin{figure}[htbp]
\centerline{\includegraphics[scale=1.0]{EPS/quadratic_dof}}
\vspace{-2mm}
\caption[DOFs and local numbering of the basis functions for quadratic
elements]{DOFs and local numbering of the basis functions for
  quadratic elements in 2d and 3d.}
\label{F:quadratic_dof}
\end{figure}

The associated DOFs for basis functions at vertices/edges are located at
the vertices/edges of the element; the entry in the vector of DOF
indices at the vertices/edges is determined by the vertex/edge offset
in the corresponding \code{admin} of the finite element space: 
the DOF index of the \code{i}--th basis functions on element \code{el} is
\bv\begin{lstlisting}
el->dof[i][admin->n0_dof[VERTEX]]
\end{lstlisting}\ev
for $\code{i = 0,\ldots,N\_VERTICES-1}$ and
\bv\begin{lstlisting}
el->dof[i][admin->n0_dof[EDGE]]
\end{lstlisting}\ev
for $\code{i = N\_VERTICES,\ldots,N\_VERTICES+N\_EDGES-1}$. Here we
used the fact, that for quadratic elements DOFs are located at the
vertices and the edges on the mesh. Thus, regardless of any other
set of  DOFs, the offset \code{mesh->node[VERTEX]} is zero and
\code{mesh->node[EDGE]} is \code{N\_VERTICES}.

Setting \code{nv = admin->n0\_dof[VERTEX]} and 
\code{ne = admin->n0\_dof[EDGE]}, the associated DOFs are shown in
Figure~\ref{F:quadratic_dof}.

\begin{example}[Accessing DOFs for piecewise quadratic finite elements]
The function \code{get\_dof\_indices()} for quadratic finite elements
can be implemented in 2d by (compare Figure~\ref{F:quadratic_dof}):
%%
\idx{get_dof_indices()@{\code{get\_dof\_indices()}}!for quadratic elements}
\bv\begin{lstlisting}
static const EL_DOF_VEC *
get_dof_indices2_2d(DOF *vec, const EL *el, const DOF_ADMIN *admin,
		    const BAS_FCTS *thisptr)
{
  static DEF_EL_VEC_CONST(DOF, rvec_space, N_BAS_LAG_2_2D, N_BAS_LAG_2_2D);
  DOF_VEC_DOF *rvec = vec ? vec : rvec_space->vec;
  int n0, ibas, inode;
  DOF **dofptr = el->dof, dof;

  n0 = (admin)->n0_dof[VERTEX];
  for (ibas = inode = 0; ibas < N_VERTICES_2D; inode++, ibas++) {
    dof = dofptr[inode][n0];
    body;
  }
  n0 = (admin)->n0_dof[EDGE];
  for (inode = 0; inode < N_EDGES_2D; inode++, ibas++) {
    dof = dofptr[N_VERTICES_2D+inode][n0];
    body;
  }

  return vec ? NULL : rvec_space;
}
\end{lstlisting}\ev
\end{example}

The boundary type of a basis functions at a vertex is the the boundary
type of the vertex, and the boundary type of a basis function at an
edge is the boundary type of the edge. The \code{i}--th interpolation
coefficient of a function $f$ on element $S$ is just
$f(x(\lambda_\code{i}))$. The implementation is similar to that for
linear finite elements and is not shown here.

The implementation of functions for interpolation during refinement
and coarsening and the restriction during coarsening becomes more
complicated and differs between the dimensions. Here we have to set
values for all elements of the refinement patch. The interpolation
during coarsening in not trivial anymore. As an example of such
implementations we present the interpolation during refinement for 2d
and 3d.

\begin{example}[Interpolation during refinement for piecewise quadratic
                finite elements in 2d]
We have to set values for the new vertex in the refinement edge,
and for the two midpoints of the bisected edge. Then we have to set
the value for the midpoint of the common edge of the two children of
the bisected triangle and we have to set the corresponding value on
the neighbor in the case that the refinement edge is not a boundary
edge:
%%
\idx{real_refine_inter()@{\code{real\_refine\_inter()}}!for quadratic elements}
\bv\begin{lstlisting}[name=real_refine_inter2_3d,label=C:real_refine_inter2_3d]
static void real_refine_inter2_3d(DOF_REAL_VEC *drv, RC_LIST_EL *list, int n)
{
  FUNCNAME("real_refine_inter2_3d");
  DOF pdof[N_BAS_LAG_2_3D];
  DOF cdof[N_BAS_LAG_2_3D];
  EL        *el;
  REAL      *v = NULL;
  DOF       cdofi;
  int       i, lr_set;
  int       node0, n0;
  const DOF_ADMIN *admin;
  const BAS_FCTS *bas_fcts;

  if (n < 1) return;
  el = list->el_info.el;

  GET_DOF_VEC(v, drv);
  if (!drv->fe_space)
  {
    ERROR("no fe_space in dof_real_vec %s\n", NAME(drv));
    return;
  }
  else if (!drv->fe_space->bas_fcts)
  {
    ERROR("no basis functions in fe_space %s\n", NAME(drv->fe_space));
    return;
  }
  GET_STRUCT(admin, drv->fe_space);
  GET_STRUCT(bas_fcts, drv->fe_space);

  get_dof_indices2_3d(pdof, el, admin, bas_fcts);

  node0 = drv->fe_space->mesh->node[EDGE];
  n0 = admin->n0_dof[EDGE];

/*--------------------------------------------------------------------------*/
/*  values on child[0]                                                      */
/*--------------------------------------------------------------------------*/

  get_dof_indices2_3d(cdof, el->child[0], admin, bas_fcts);

  v[cdof[3]] = (v[pdof[4]]);
  v[cdof[6]] = (0.375*v[pdof[0]] - 0.125*v[pdof[1]]
		+ 0.75*v[pdof[4]]);
  v[cdof[8]] = (0.125*(-v[pdof[0]] - v[pdof[1]]) + 0.25*v[pdof[4]]
		+ 0.5*(v[pdof[5]] + v[pdof[7]]));
  v[cdof[9]] = (0.125*(-v[pdof[0]] - v[pdof[1]]) + 0.25*v[pdof[4]]
		+ 0.5*(v[pdof[6]] + v[pdof[8]]));

/*--------------------------------------------------------------------------*/
/*  values on child[1]                                                      */
/*--------------------------------------------------------------------------*/

  cdofi = el->child[1]->dof[node0+2][n0];
  v[cdofi] = (-0.125*v[pdof[0]] + 0.375*v[pdof[1]]
	      + 0.75*v[pdof[4]]);

/*--------------------------------------------------------------------------*/
/*   adjust neighbour values                                                */
/*--------------------------------------------------------------------------*/

  for (i = 1; i < n; i++)
  {
    el = list[i].el_info.el;
    get_dof_indices2_3d(pdof, el, admin, bas_fcts);

    lr_set = 0;
    if (list[i].neigh[0]  &&  list[i].neigh[0]->no < i)
      lr_set = 1;

    if (list[i].neigh[1]  &&  list[i].neigh[1]->no < i)
      lr_set += 2;

    DEBUG_TEST_EXIT(lr_set, "no values set on both neighbours\n");

/*--------------------------------------------------------------------------*/
/*  values on child[0]                                                      */
/*--------------------------------------------------------------------------*/

    switch (lr_set)
    {
    case 1:
      cdofi = el->child[0]->dof[node0+4][n0];
      v[cdofi] = (0.125*(-v[pdof[0]] - v[pdof[1]]) + 0.25*v[pdof[4]]
		  + 0.5*(v[pdof[5]] + v[pdof[7]]));
      break;
    case 2:
      cdofi = el->child[0]->dof[node0+5][n0];
      v[cdofi] = (0.125*(-v[pdof[0]] - v[pdof[1]]) + 0.25*v[pdof[4]]
		  + 0.5*(v[pdof[6]] + v[pdof[8]]));
    }
  }
  return;
}
\end{lstlisting}\ev
\end{example}

\begin{example}[Interpolation during refinement for piecewise quadratic
                finite elements in 3d]
Here, we first have to set values for all DOFs that belong to the
first element of the refinement patch. Then we have to loop over the
refinement patch and set all DOFs that have not previously been set on
another patch element. In order to set values only once, by the
variable \code{lr\_set} we check, if a common DOFs with a left or
right neighbor is set by the neighbor. Such values are already set
if the neighbor is a prior element in the list. Since all values
are set on the first element for all subsequent elements there
must be DOFs which have been set by another element.
\idx{real_refine_inter@{\code{real\_refine\_inter()}}!for quadratic elements}
\bv\begin{lstlisting}
static void real_refine_inter2_3d(DOF_REAL_VEC *drv, RC_LIST_EL *list, int n)
{
  FUNCNAME("real_refine_inter2_3d");
  EL        *el;
  REAL      *v = nil;
  const DOF *cdof;
  DOF       pdof[N_BAS2_3D], cdofi;
  int       i, lr_set;
  int       node0, n0;
  const DOF       *(*get_dof_indices)(const EL *, const DOF_ADMIN *, DOF *);
  const DOF_ADMIN *admin;

  if (n < 1) return;
  el = list->el_info.el;

  GET_DOF_VEC(v, drv);
  if (!drv->fe_space)
  {
    ERROR("no fe_space in dof_real_vec %s\n", NAME(drv));
    return;
  }
  else if (!drv->fe_space->bas_fcts)
  {
    ERROR("no basis functions in fe_space %s\n", NAME(drv->fe_space));
    return;
  }
  get_dof_indices = drv->fe_space->bas_fcts->get_dof_indices;
  GET_STRUCT(admin,drv->fe_space);

  get_dof_indices(el, admin, pdof);

  node0 = drv->fe_space->mesh->node[EDGE];
  n0 = admin->n0_dof[EDGE];

/*--------------------------------------------------------------------------*/
/*  values on child[0]                                                      */
/*--------------------------------------------------------------------------*/

  cdof = get_dof_indices(el->child[0], admin, nil);

  v[cdof[3]] = (v[pdof[4]]);
  v[cdof[6]] = (0.375*v[pdof[0]] - 0.125*v[pdof[1]] 
                + 0.75*v[pdof[4]]);
  v[cdof[8]] = (0.125*(-v[pdof[0]] - v[pdof[1]]) + 0.25*v[pdof[4]]
                + 0.5*(v[pdof[5]] + v[pdof[7]]));
  v[cdof[9]] = (0.125*(-v[pdof[0]] - v[pdof[1]]) + 0.25*v[pdof[4]]
                + 0.5*(v[pdof[6]] + v[pdof[8]]));

/*--------------------------------------------------------------------------*/
/*  values on child[1]                                                      */
/*--------------------------------------------------------------------------*/
  
  cdofi = el->child[1]->dof[node0+2][n0];
  v[cdofi] = (-0.125*v[pdof[0]] + 0.375*v[pdof[1]] 
              + 0.75*v[pdof[4]]);

/*--------------------------------------------------------------------------*/
/*   adjust neighbour values                                                */
/*--------------------------------------------------------------------------*/
  
  for (i = 1; i < n; i++)
  {
    el = list[i].el_info.el;
    get_dof_indices(el, admin, pdof);

    lr_set = 0;
    if (list[i].neigh[0]  &&  list[i].neigh[0]->no < i)
      lr_set = 1;

    if (list[i].neigh[1]  &&  list[i].neigh[1]->no < i)
      lr_set += 2;

    DEBUG_TEST_EXIT(lr_set, "no values set on both neighbours\n");

/*--------------------------------------------------------------------------*/
/*  values on child[0]                                                      */
/*--------------------------------------------------------------------------*/

    switch (lr_set)
    {
    case 1:
      cdofi = el->child[0]->dof[node0+4][n0];      
      v[cdofi] = (0.125*(-v[pdof[0]] - v[pdof[1]]) + 0.25*v[pdof[4]]
                  + 0.5*(v[pdof[5]] + v[pdof[7]]));
      break;
    case 2:
      cdofi = el->child[0]->dof[node0+5][n0];      
      v[cdofi] = (0.125*(-v[pdof[0]] - v[pdof[1]]) + 0.25*v[pdof[4]]
                  + 0.5*(v[pdof[6]] + v[pdof[8]]));
    }
  }
  return;
}
\end{lstlisting}\ev
\end{example}

\subsubsection{Piecewise cubic finite elements}\label{S:cubic-fe}

For Lagrange elements of third order we have four basis functions in
1d, ten basis functions in 2d, and 20 in 3d; the basis functions and
the corresponding Lagrange nodes in barycentric coordinates are shown
in Tables \ref{T:P3_1d}, \ref{T:P3_2d}, and \ref{T:P3_3d}.
\begin{table}[htbp]
\begin{center}\begin{small}\begin{math}
\begin{array}{|rcl@{\qquad}|@{\qquad}l@{\qquad}|@{\qquad}rcl|}\hline
\multicolumn{2}{|l}{\mbox{function}} & & \mbox{position} &
\multicolumn{3}{l|}{\mbox{Lagrange node}}
 \\\hline\vphantom{\pbar^{\int}}
%
\pbar^\code{0}(\lambda) & = & 
\half (3 \lambda_\code{0} - 1)(3 \lambda_\code{0} - 2)\lambda_\code{0} &
\mbox{vertex \code{0}} & \lambda^\code{0} & = & (1,0)\\[1mm]
%
\pbar^\code{1}(\lambda) & = &
\half (3 \lambda_\code{1} - 1)(3 \lambda_\code{1} - 2)\lambda_\code{1} &
\mbox{vertex \code{1}} & \lambda^\code{1} & = & (0,1)\\[1mm]
%
\pbar^\code{2}(\lambda) & = & 
\frac92 (3 \lambda_\code{0} - 1)\lambda_\code{0}\lambda_\code{1}&
\mbox{center}& \lambda^\code{2} & = & (\frac23,\frac13)\\[1mm]
\pbar^\code{3}(\lambda) & = & 
\frac92 (3 \lambda_\code{1} - 1)\lambda_\code{1}\lambda_\code{0}&
\mbox{center}& \lambda^\code{3} & = & (\frac13,\frac23)\\[1mm]
%
\hline\end{array}
\end{math}\end{small}\end{center}
\vskip -5mm
\caption{Local basis functions for cubic finite elements in 1d.}
\label{T:P3_1d}
\end{table}
\begin{table}[htbp]
\begin{center}\begin{small}\begin{math}
\begin{array}{|rcl@{\qquad}|@{\qquad}l@{\qquad}|@{\qquad}rcl|}\hline
\multicolumn{2}{|l}{\mbox{function}} & & \mbox{position} &
\multicolumn{3}{l|}{\mbox{Lagrange node}}
 \\\hline\vphantom{\pbar^{\int}}
%
\pbar^\code{0}(\lambda) & = & 
\half (3 \lambda_\code{0} - 1)(3 \lambda_\code{0} - 2)\lambda_\code{0} &
\mbox{vertex \code{0}} & \lambda^\code{0} & = & (1,0,0)\\[1mm]
%
\pbar^\code{1}(\lambda) & = &
\half (3 \lambda_\code{1} - 1)(3 \lambda_\code{1} - 2)\lambda_\code{1} &
\mbox{vertex \code{1}} & \lambda^\code{1} & = & (0,1,0)\\[1mm]
%
\pbar^\code{2}(\lambda) & = &
\half (3 \lambda_\code{2} - 1)(3 \lambda_\code{2} - 2)\lambda_\code{2} &
\mbox{vertex \code{2}} & \lambda^\code{2} & = & (0,0,1)\\[1mm]
%
\pbar^\code{3}(\lambda) & = & 
\frac92 (3 \lambda_\code{1} - 1)\lambda_\code{1}\lambda_\code{2}&
\mbox{edge \code{0}}& \lambda^\code{3} & = & (0,\frac23,\frac13)\\[1mm]
\pbar^\code{4}(\lambda) & = & 
\frac92 (3 \lambda_\code{2} - 1)\lambda_\code{2}\lambda_\code{1}&
\mbox{edge \code{0}}& \lambda^\code{4} & = & (0,\frac13,\frac23)\\[1mm]
%
\pbar^\code{5}(\lambda) & = & 
\frac92 (3 \lambda_\code{2} - 1)\lambda_\code{2}\lambda_\code{0}&
\mbox{edge \code{1}}& \lambda^\code{5} & = & (\frac13,0,\frac23)\\[1mm]
\pbar^\code{6}(\lambda) & = & 
\frac92 (3 \lambda_\code{0} - 1)\lambda_\code{0}\lambda_\code{2}&
\mbox{edge \code{1}}& \lambda^\code{6} & = & (\frac23,0,\frac13)\\[1mm]
%
\pbar^\code{7}(\lambda) & = & 
\frac92 (3 \lambda_\code{0} - 1)\lambda_\code{0}\lambda_\code{1}&
\mbox{edge \code{2}}& \lambda^\code{7} & = & (\frac23,\frac13,0)\\[1mm]
\pbar^\code{8}(\lambda) & = & 
\frac92 (3 \lambda_\code{1} - 1)\lambda_\code{1}\lambda_\code{0}&
\mbox{edge \code{2}}& \lambda^\code{8} & = & (\frac13,\frac23,0)\\[1mm]
%
\pbar^\code{9}(\lambda) & = & 
27 \lambda_\code{0}\lambda_\code{1}\lambda_\code{2}&
\mbox{center}& \lambda^\code{9} & = & (\frac13,\frac13,\frac13)\\[1mm]
\hline\end{array}
\end{math}\end{small}\end{center}
\vskip -5mm
\caption{Local basis functions for cubic finite elements in 2d.}
\label{T:P3_2d}
\end{table}
\begin{table}[htbp]
\begin{center}\begin{small}\begin{math}
\begin{array}{|rcl@{\qquad}|@{\qquad}l@{\qquad}|@{\qquad}rcl|}\hline
\multicolumn{2}{|l}{\mbox{function}} & & \mbox{position} &
\multicolumn{3}{l|}{\mbox{Lagrange node}}
 \\\hline\vphantom{\pbar^{\int}}
%
\pbar^\code{0}(\lambda) & = & \half (3 \lc{0} - 1)(3 \lc{0} - 2)\lc{0} 
& \mbox{vertex \code{0}} & \lambda^\code{0} & = &  (1,0,0,0)\\[1mm]
%
\pbar^\code{1}(\lambda) & = & \half (3 \lc{1} - 1)(3 \lc{1} - 2)\lc{1}  
& \mbox{vertex \code{1}} & \lambda^\code{1} & = &  (0,1,0,0)\\[1mm]
%
\pbar^\code{2}(\lambda) & = &  \half (3 \lc{2} - 1)(3 \lc{2} - 2)\lc{2} 
& \mbox{vertex \code{2}} & \lambda^\code{2} & = &  (0,0,1,0)\\[1mm]
%
\pbar^\code{3}(\lambda) & = &  \half (3 \lc{3} - 1)(3 \lc{3} - 2)\lc{3} 
& \mbox{vertex \code{3}} & \lambda^\code{3} & = &  (0,0,0,1)\\[1mm]
%
\pbar^\code{4}(\lambda) & = &  \frac92 (3 \lc{0} - 1)\lc{0}\lc{1}
& \mbox{edge \code{0}} & \lambda^\code{4} & = & (\frac23,\frac13,0,0)\\[1mm]
%
\pbar^\code{5}(\lambda) & = &  \frac92 (3 \lc{1} - 1)\lc{1}\lc{0}
& \mbox{edge \code{0}} & \lambda^\code{5} & = & (\frac13,\frac23,0,0)\\[1mm]
%
\pbar^\code{6}(\lambda) & = &  \frac92 (3 \lc{0} - 1)\lc{0}\lc{2}
& \mbox{edge \code{1}} & \lambda^\code{6} & = & (\frac23,0,\frac13,0)\\[1mm]
%
\pbar^\code{7}(\lambda) & = &  \frac92 (3 \lc{2} - 1)\lc{2}\lc{0}
& \mbox{edge \code{1}} & \lambda^\code{7} & = &  (\frac13,0,\frac23,0)\\[1mm]
%
\pbar^\code{8}(\lambda) & = &  \frac92 (3 \lc{0} - 1)\lc{0}\lc{3}
& \mbox{edge \code{2}} & \lambda^\code{8} & = &  (\frac23,0,0,\frac13)\\[1mm]
%
\pbar^\code{9}(\lambda) & = &  \frac92 (3 \lc{3} - 1)\lc{3}\lc{0}
& \mbox{edge \code{2}} & \lambda^\code{9} & = &  (\frac13,0,0,\frac23)\\[1mm]
%
\pbar^\code{10}(\lambda) & = &  \frac92 (3 \lc{1} - 1)\lc{1}\lc{2}
& \mbox{edge \code{3}} & \lambda^\code{10} & = & (0,\frac23,\frac13,0)\\[1mm]
%
\pbar^\code{11}(\lambda) & = &  \frac92 (3 \lc{2} - 1)\lc{2}\lc{1}
& \mbox{edge \code{3}} & \lambda^\code{11} & = & (0,\frac13,\frac23,0)\\[1mm]
%
\pbar^\code{12}(\lambda) & = &  \frac92 (3 \lc{1} - 1)\lc{1}\lc{3}
& \mbox{edge \code{4}} & \lambda^\code{12} & = & (0,\frac23,0,\frac13)\\[1mm]
%
\pbar^\code{13}(\lambda) & = &  \frac92 (3 \lc{3} - 1)\lc{3}\lc{1}
& \mbox{edge \code{4}} & \lambda^\code{13} & = & (0,\frac13,0,\frac23)\\[1mm]
%
\pbar^\code{14}(\lambda) & = &  \frac92 (3 \lc{2} - 1)\lc{2}\lc{3}
& \mbox{edge \code{5}} & \lambda^\code{14} & = & (0,0,\frac23,\frac13)\\[1mm]
%
\pbar^\code{15}(\lambda) & = &  \frac92 (3 \lc{3} - 1)\lc{3}\lc{2}
& \mbox{edge \code{5}} & \lambda^\code{15} & = & (0,0,\frac13,\frac23)\\[1mm]
%
\pbar^\code{16}(\lambda) & = &  27 \lc{1}\lc{2}\lc{3}
& \mbox{face \code{0}} & \lambda^\code{16} & = & (0,\frac13,\frac13,\frac13)
\\[1mm]
%
\pbar^\code{17}(\lambda) & = &  27 \lc{2}\lc{3}\lc{0}
& \mbox{face \code{1}} & \lambda^\code{17} & = & (\frac13,0,\frac13,\frac13)
\\[1mm]
%
\pbar^\code{18}(\lambda) & = &  27 \lc{3}\lc{0}\lc{1}
& \mbox{face \code{2}} & \lambda^\code{18} & = & (\frac13,\frac13,0,\frac13)
\\[1mm]
%
\pbar^\code{19}(\lambda) & = &  27 \lc{0}\lc{1}\lc{2}
& \mbox{face \code{3}} & \lambda^\code{19} & = & (\frac13,\frac13,\frac13,0)
\\[1mm]\hline\end{array}
\end{math}\end{small}\end{center}
\vskip -5mm
\caption{Local basis functions for cubic finite elements in 3d.}
\label{T:P3_3d}
\end{table}

For cubic elements we have to face a further difficulty. At each edge
two basis functions are located. The two DOFs of the \code{i}--th edge
are subsequent entries in the vector \code{el->dof[i]}. 
For two neighboring triangles the common edge has a different
orientation with respect to the local numbering of vertices
on the two triangles. In Figure~\ref{F:cubic_orient} the 3rd local 
basis function on the left and the 4th on the right triangle
built up the global basis function, e.g.; thus, both local basis function
must have access to the same global DOF.

\begin{figure}[htbp]
\centerline{\includegraphics[scale=0.5]{EPS/cubic_orient}}
\vspace{-2mm}
\caption{Cubic DOFs on a patch of two triangles.}\label{F:cubic_orient}
\end{figure}
%
In order to combine the global DOF with the local basis function in
the implementation of the \code{get\_dof\_indices()} function, we have
to give every edge a \emph{global orientation}, i.e, every edge has a
unique beginning and end point. Using the orientation of an edge we
are able to order the DOFs stored at this edge. Let for example the
common edge in Figure \ref{F:cubic_orient} be oriented from bottom to
top.  The global DOF corresponding to 3rd local DOF on the left and
the 4th local DOF on the right is then
\bv\begin{lstlisting}
el->dof[N_VERTICES+0][admin->n0_dof[EDGE]]
\end{lstlisting}\ev
and for the 4th local DOF on the left and 3rd local DOF on the right
\bv\begin{lstlisting}
el->dof[N_VERTICES+0][admin->n0_dof[EDGE]+1]
\end{lstlisting}\ev
The global orientation gives a unique access to local DOFs from global ones.

\begin{example}[Accessing DOFs for piecewise cubic finite elements]
For the implementation, we use in 2d as well as in 3d an
orientation defined by the DOF indices at the edges' vertices. The vertex
with the smaller (global) DOF index is the beginning point, the vertex
with the higher index the end point. For cubics the implementation
differs between 2d and 3d. In 2d we have one degree of freedom
at the center and in 3d one degree of freedom at each face and
none at the center. The DOFs at an edge are accessed according
to the orientation of the edge. We present the implementation in 2d:
\bv\begin{lstlisting}
#define N_BAS_LAG_3_2D (N_VERTICES_2D+2*N_EDGES_2D+1)

static const EL_DOF_VEC *
get_dof_indices3_2d(DOF *vec, const EL *el, const DOF_ADMIN *admin,
			const BAS_FCTS *thisptr)
{
  static DEF_EL_VEC_CONST(type, rvec_space, N_BAS_LAG_3_2D, N_BAS_LAG_3_2D);
  DOF *rvec = vec ? vec : rvec_space->vec;
  int n0, ibas, inode;
  DOF **dofptr = el->dof, dof;

  n0 = admin->n0_dof[VERTEX];
  for (ibas = 0; ibas < N_VERTICES_2D; ibas++) {
    dof = dofptr[ibas][n0];
    body;
  }
  n0 = admin->n0_dof[EDGE];
  for (inode = 0, ibas = N_VERTICES_2D; inode < N_EDGES_2D; inode++) {
    if (dofptr[vertex_of_edge_2d[inode][0]][0]
	< dofptr[vertex_of_edge_2d[inode][1]][0]) {
      dof = dofptr[N_VERTICES_2D+inode][n0];
      body;
      ibas++;
      dof = dofptr[N_VERTICES_2D+inode][n0+1];
      body;
      ibas++;
    } else {
      dof = dofptr[N_VERTICES_2D+inode][n0+1];
      body;
      ibas++;
      dof  = dofptr[N_VERTICES_2D+inode][n0];
      body;
      ibas++;
    }
  }
  n0 = admin->n0_dof[CENTER];
  dof = dofptr[6][n0];
  body;

  return vec ? NULL : rvec_space;
}
\end{lstlisting}\ev
\end{example}

\subsubsection{Piecewise quartic finite elements}\label{S:quartic-fe}

For Lagrange elements of fourth order we have 5 basis functions
in 1d, 15 in 2d, and 35 in 3d; the basis functions and the corresponding
Lagrange nodes in barycentric coordinates are shown in
Tables \ref{T:P4_1d}, \ref{T:P4_2d}, and \ref{T:P4_3d}.
\begin{table}[htbp]
\begin{center}\begin{small}\begin{math}
\begin{array}{|rcl@{\qquad}|@{\qquad}l@{\qquad}|@{\qquad}rcl|}\hline
\multicolumn{2}{|l}{\mbox{function}} & & \mbox{position} &
\multicolumn{3}{l|}{\mbox{Lagrange node}}
 \\\hline\vphantom{\pbar^{\int}}
%
\pbar^\code{0}(\lambda) & = & 
\frac13 (4 \lambda_\code{0} - 1)(2 \lambda_\code{0} - 1)
(4 \lambda_\code{0} - 3)\lambda_\code{0} &
\mbox{vertex \code{0}} & \lambda^\code{0} & = & (1,0)\\[1mm]
%
\pbar^\code{1}(\lambda) & = & 
\frac13 (4 \lambda_\code{1} - 1)(2 \lambda_\code{1} - 1)
(4 \lambda_\code{1} - 3)\lambda_\code{1} &
\mbox{vertex \code{1}} & \lambda^\code{1} & = & (0,1)\\[1mm]
%
\pbar^\code{2}(\lambda) & = & 
\frac{16}3 (4 \lambda_\code{0} - 1)(2 \lambda_\code{0} - 1)
\lambda_\code{0}\lambda_\code{1} &
\mbox{center} & \lambda^\code{2} & = & (\frac34,\frac14)\\[1mm]
%
\pbar^\code{3}(\lambda) & = & 
4 (4 \lambda_\code{0} - 1)(4 \lambda_\code{1} - 1)
\lambda_\code{0}\lambda_\code{1} &
\mbox{center} & \lambda^\code{3} & = & (\frac12,\frac12)\\[1mm]
%
\pbar^\code{4}(\lambda) & = & 
\frac{16}3 (4 \lambda_\code{1} - 1)(2 \lambda_\code{1} - 1)
\lambda_\code{0}\lambda_\code{1} &
\mbox{center} & \lambda^\code{4} & = & (\frac14,\frac34)\\[1mm]
%
\hline\end{array}
\end{math}\end{small}\end{center}
\vskip -5mm
\caption{Local basis functions for quartic finite elements in 1d.}
\label{T:P4_1d}
\end{table}
\begin{table}[htbp]
\begin{center}\begin{small}\begin{math}
\begin{array}{|rcl@{\qquad}|@{\qquad}l@{\qquad}|@{\qquad}rcl|}\hline
\multicolumn{2}{|l}{\mbox{function}} & & \mbox{position} &
\multicolumn{3}{l|}{\mbox{Lagrange node}}
 \\\hline\vphantom{\pbar^{\int}}
%
\pbar^\code{0}(\lambda) & = & 
\frac13 (4 \lambda_\code{0} - 1)(2 \lambda_\code{0} - 1)
(4 \lambda_\code{0} - 3)\lambda_\code{0} &
\mbox{vertex \code{0}} & \lambda^\code{0} & = & (1,0,0)\\[1mm]
%
\pbar^\code{1}(\lambda) & = & 
\frac13 (4 \lambda_\code{1} - 1)(2 \lambda_\code{1} - 1)
(4 \lambda_\code{1} - 3)\lambda_\code{1} &
\mbox{vertex \code{1}} & \lambda^\code{1} & = & (0,1,0)\\[1mm]
%
\pbar^\code{2}(\lambda) & = & 
\frac13 (4 \lambda_\code{2} - 1)(2 \lambda_\code{2} - 1)
(4 \lambda_\code{2} - 3)\lambda_\code{2} &
\mbox{vertex \code{2}} & \lambda^\code{2} & = & (0,0,1)\\[1mm]
%
\pbar^\code{3}(\lambda) & = & 
\frac{16}3 (4 \lambda_\code{1} - 1)(2 \lambda_\code{1} - 1)
\lambda_\code{1}\lambda_\code{2} &
\mbox{edge \code{0}} & \lambda^\code{3} & = & (0,\frac34,\frac14)\\[1mm]
%
\pbar^\code{4}(\lambda) & = & 
4 (4 \lambda_\code{1} - 1)(4 \lambda_\code{2} - 1)
\lambda_\code{1}\lambda_\code{2} &
\mbox{edge \code{0}} & \lambda^\code{4} & = & (0,\frac12,\frac12)\\[1mm]
%
\pbar^\code{5}(\lambda) & = & 
\frac{16}3 (4 \lambda_\code{2} - 1)(2 \lambda_\code{2} - 1)
\lambda_\code{1}\lambda_\code{2} &
\mbox{edge \code{0}} & \lambda^\code{5} & = & (0,\frac14,\frac34)\\[1mm]
%
\pbar^\code{6}(\lambda) & = & 
\frac{16}3 (4 \lambda_\code{2} - 1)(2 \lambda_\code{2} - 1)
\lambda_\code{0}\lambda_\code{2} &
\mbox{edge \code{1}} & \lambda^\code{6} & = & (\frac14,0,\frac34)\\[1mm]
%
\pbar^\code{7}(\lambda) & = & 
4 (4 \lambda_\code{2} - 1)(4 \lambda_\code{0} - 1)
\lambda_\code{0}\lambda_\code{2} &
\mbox{edge \code{1}} & \lambda^\code{7} & = & (\frac12,0,\frac12)\\[1mm]
%
\pbar^\code{8}(\lambda) & = & 
\frac{16}3 (4 \lambda_\code{0} - 1)(2 \lambda_\code{0} - 1)
\lambda_\code{0}\lambda_\code{2} &
\mbox{edge \code{1}} & \lambda^\code{8} & = & (\frac34,0,\frac14)\\[1mm]
%
\pbar^\code{9}(\lambda) & = & 
\frac{16}3 (4 \lambda_\code{0} - 1)(2 \lambda_\code{0} - 1)
\lambda_\code{0}\lambda_\code{1} &
\mbox{edge \code{2}} & \lambda^\code{9} & = & (\frac34,\frac14,0)\\[1mm]
%
\pbar^\code{10}(\lambda) & = & 
4 (4 \lambda_\code{0} - 1)(4 \lambda_\code{1} - 1)
\lambda_\code{0}\lambda_\code{1} &
\mbox{edge \code{2}} & \lambda^\code{10} & = & (\frac12,\frac12,0)\\[1mm]
%
\pbar^\code{11}(\lambda) & = & 
\frac{16}3 (4 \lambda_\code{1} - 1)(2 \lambda_\code{1} - 1)
\lambda_\code{0}\lambda_\code{1} &
\mbox{edge \code{2}} & \lambda^\code{11} & = & (\frac14,\frac34,0)\\[1mm]
%
\pbar^\code{12}(\lambda) & = & 
32 (4 \lambda_\code{0} - 1)\lambda_\code{0}\lambda_\code{1}\lambda_\code{2} &
\mbox{center} & \lambda^\code{12} & = & (\frac12,\frac14,\frac14)\\[1mm]
%
\pbar^\code{13}(\lambda) & = & 
32 (4 \lambda_\code{1} - 1)\lambda_\code{0}\lambda_\code{1}\lambda_\code{2} &
\mbox{center} & \lambda^\code{13} & = & (\frac14,\frac12,\frac14)\\[1mm]
%
\pbar^\code{14}(\lambda) & = & 
32 (4 \lambda_\code{2} - 1)\lambda_\code{0}\lambda_\code{1}\lambda_\code{2} &
\mbox{center} & \lambda^\code{14} & = & (\frac14,\frac14,\frac12)\\[1mm]
\hline\end{array}
\end{math}\end{small}\end{center}
\vskip -5mm
\caption{Local basis functions for quartic finite elements in 2d.}
\label{T:P4_2d}
\end{table}
\begin{table}[htbp]
\begin{center}\begin{small}\begin{math}%\displaystyle
\begin{array}{|rcl@{\qquad}|@{\qquad}l@{\qquad}|@{\qquad}rcl|}\hline
\multicolumn{2}{|l}{\mbox{function}} & & \mbox{position} &
\multicolumn{3}{l|}{\mbox{Lagrange node}}
 \\\hline\vphantom{\pbar^{\int}}
%
\pbar^\code{0}(\lambda) & = & 
\frac13 (4 \lambda_\code{0} - 1)(2 \lambda_\code{0} - 1)
(4 \lambda_\code{0} - 3)\lambda_\code{0} &
\mbox{vertex \code{0}} & \lambda^\code{0} & = & (1,0,0,0)\\[1mm]
%
\pbar^\code{1}(\lambda) & = & 
\frac13 (4 \lambda_\code{1} - 1)(2 \lambda_\code{1} - 1)
(4 \lambda_\code{1} - 3)\lambda_\code{1} &
\mbox{vertex \code{1}} & \lambda^\code{1} & = & (0,1,0,0)\\[1mm]
%
\pbar^\code{2}(\lambda) & = & 
\frac13 (4 \lambda_\code{2} - 1)(2 \lambda_\code{2} - 1)
(4 \lambda_\code{2} - 3)\lambda_\code{2} &
\mbox{vertex \code{2}} & \lambda^\code{2} & = & (0,0,1,0)\\[1mm]
%
\pbar^\code{3}(\lambda) & = & 
\frac13 (4 \lambda_\code{3} - 1)(2 \lambda_\code{3} - 1)
(4 \lambda_\code{3} - 3)\lambda_\code{3} &
\mbox{vertex \code{3}} & \lambda^\code{3} & = & (0,0,0,1)\\[1mm]
%
\pbar^\code{4}(\lambda) & = & 
\frac{16}3 (4 \lambda_\code{0} - 1)(2 \lambda_\code{0} - 1)
\lambda_\code{0}\lambda_\code{1} &
\mbox{edge \code{0}} & \lambda^\code{4} & = & (\frac34,\frac14,0,0)\\[1mm]
%
\pbar^\code{5}(\lambda) & = & 
4 (4 \lambda_\code{0} - 1)(4 \lambda_\code{1} - 1)
\lambda_\code{0}\lambda_\code{1} &
\mbox{edge \code{0}} & \lambda^\code{5} & = & (\frac12,\frac12,0,0)\\[1mm]
%
\pbar^\code{6}(\lambda) & = & 
\frac{16}3 (4 \lambda_\code{1} - 1)(2 \lambda_\code{1} - 1)
\lambda_\code{0}\lambda_\code{1} &
\mbox{edge \code{0}} & \lambda^\code{6} & = & (\frac14,\frac34,0,0)\\[1mm]
%
\pbar^\code{7}(\lambda) & = & 
\frac{16}3 (4 \lambda_\code{0} - 1)(2 \lambda_\code{0} - 1)
\lambda_\code{0}\lambda_\code{2} &
\mbox{edge \code{1}} & \lambda^\code{7} & = & (\frac34,0,\frac14,0)\\[1mm]
%
\pbar^\code{8}(\lambda) & = & 
4 (4 \lambda_\code{0} - 1)(4 \lambda_\code{2} - 1)
\lambda_\code{0}\lambda_\code{2} &
\mbox{edge \code{1}} & \lambda^\code{8} & = & (\frac12,0,\frac12,0)\\[1mm]
%
\pbar^\code{9}(\lambda) & = & 
\frac{16}3 (4 \lambda_\code{2} - 1)(2 \lambda_\code{2} - 1)
\lambda_\code{0}\lambda_\code{2} &
\mbox{edge \code{1}} & \lambda^\code{9} & = & (\frac14,0,\frac34,0)\\[1mm]
%
\pbar^\code{10}(\lambda) & = & 
\frac{16}3 (4 \lambda_\code{0} - 1)(2 \lambda_\code{0} - 1)
\lambda_\code{0}\lambda_\code{3} &
\mbox{edge \code{2}} & \lambda^\code{10} & = & (\frac34,0,0,\frac14)\\[1mm]
%
\pbar^\code{11}(\lambda) & = & 
4 (4 \lambda_\code{0} - 1)(4 \lambda_\code{3} - 1)
\lambda_\code{0}\lambda_\code{3} &
\mbox{edge \code{2}} & \lambda^\code{11} & = & (\frac12,0,0,\frac12)\\[1mm]
%
\pbar^\code{12}(\lambda) & = & 
\frac{16}3 (4 \lambda_\code{3} - 1)(2 \lambda_\code{3} - 1)
\lambda_\code{0}\lambda_\code{3} &
\mbox{edge \code{2}} & \lambda^\code{12} & = & (\frac14,0,0,\frac34)\\[1mm]
%
\pbar^\code{13}(\lambda) & = & 
\frac{16}3 (4 \lambda_\code{1} - 1)(2 \lambda_\code{1} - 1)
\lambda_\code{1}\lambda_\code{2} &
\mbox{edge \code{3}} & \lambda^\code{13} & = & (0,\frac34,\frac14,0)\\[1mm]
%
\pbar^\code{14}(\lambda) & = & 
4 (4 \lambda_\code{1} - 1)(4 \lambda_\code{2} - 1)
\lambda_\code{1}\lambda_\code{2} &
\mbox{edge \code{3}} & \lambda^\code{14} & = & (0,\frac12,\frac12,0)\\[1mm]
%
\pbar^\code{15}(\lambda) & = & 
\frac{16}3 (4 \lambda_\code{2} - 1)(2 \lambda_\code{2} - 1)
\lambda_\code{1}\lambda_\code{2} &
\mbox{edge \code{3}} & \lambda^\code{15} & = & (0,\frac14,\frac34,0)\\[1mm]
%
\pbar^\code{16}(\lambda) & = & 
\frac{16}3 (4 \lambda_\code{1} - 1)(2 \lambda_\code{1} - 1)
\lambda_\code{1}\lambda_\code{3} &
\mbox{edge \code{4}} & \lambda^\code{16} & = & (0,\frac34,0,\frac14)\\[1mm]
%
\pbar^\code{17}(\lambda) & = & 
4 (4 \lambda_\code{1} - 1)(4 \lambda_\code{3} - 1)
\lambda_\code{1}\lambda_\code{3} &
\mbox{edge \code{4}} & \lambda^\code{17} & = & (0,\frac12,0,\frac12)\\[1mm]
%
\pbar^\code{18}(\lambda) & = & 
\frac{16}3 (4 \lambda_\code{3} - 1)(2 \lambda_\code{3} - 1)
\lambda_\code{1}\lambda_\code{3} &
\mbox{edge \code{4}} & \lambda^\code{18} & = & (0,\frac14,0,\frac34)\\[1mm]
%
\pbar^\code{19}(\lambda) & = & 
\frac{16}3 (4 \lambda_\code{2} - 1)(2 \lambda_\code{2} - 1)
\lambda_\code{2}\lambda_\code{3} &
\mbox{edge \code{5}} & \lambda^\code{19} & = & (0,0,\frac34,\frac14)\\[1mm]
%
\pbar^\code{20}(\lambda) & = & 
4 (4 \lambda_\code{2} - 1)(4 \lambda_\code{3} - 1)
\lambda_\code{2}\lambda_\code{3} &
\mbox{edge \code{5}} & \lambda^\code{20} & = & (0,0,\frac12,\frac12)\\[1mm]
%
\pbar^\code{21}(\lambda) & = & 
\frac{16}3 (4 \lambda_\code{3} - 1)(2 \lambda_\code{3} - 1)
\lambda_\code{2}\lambda_\code{3} &
\mbox{edge \code{5}} & \lambda^\code{21} & = & (0,0,\frac14,\frac34)\\[1mm]
%
\pbar^\code{22}(\lambda) & = & 
32 (4 \lambda_\code{1} - 1)\lambda_\code{1}\lambda_\code{2}\lambda_\code{3} &
\mbox{face \code{0}} & \lambda^\code{22} & = & (0,\frac12,\frac14,\frac14)
\\[1mm]
%
\pbar^\code{23}(\lambda) & = & 
32 (4 \lambda_\code{2} - 1)\lambda_\code{1}\lambda_\code{2}\lambda_\code{3} &
\mbox{face \code{0}} & \lambda^\code{23} & = & (0,\frac14,\frac12,\frac14)
\\[1mm]
%
\pbar^\code{24}(\lambda) & = & 
32 (4 \lambda_\code{3} - 1)\lambda_\code{1}\lambda_\code{2}\lambda_\code{3} &
\mbox{face \code{0}} & \lambda^\code{24} & = & (0,\frac14,\frac14,\frac12)
\\[1mm]
%
\pbar^\code{25}(\lambda) & = & 
32 (4 \lambda_\code{0} - 1)\lambda_\code{0}\lambda_\code{2}\lambda_\code{3} &
\mbox{face \code{1}} & \lambda^\code{25} & = & (\frac12,0,\frac14,\frac14)
\\[1mm]
%
\pbar^\code{26}(\lambda) & = & 
32 (4 \lambda_\code{2} - 1)\lambda_\code{0}\lambda_\code{2}\lambda_\code{3} &
\mbox{face \code{1}} & \lambda^\code{26} & = & (\frac14,0,\frac12,\frac14)
\\[1mm]
%
\pbar^\code{27}(\lambda) & = & 
32 (4 \lambda_\code{3} - 1)\lambda_\code{0}\lambda_\code{2}\lambda_\code{3} &
\mbox{face \code{1}} & \lambda^\code{27} & = & (\frac14,0,\frac14,\frac12)
\\[1mm]
%
\pbar^\code{28}(\lambda) & = & 
32 (4 \lambda_\code{0} - 1)\lambda_\code{0}\lambda_\code{1}\lambda_\code{3} &
\mbox{face \code{2}} & \lambda^\code{28} & = & (\frac12,\frac14,0,\frac14)
\\[1mm]
%
\pbar^\code{29}(\lambda) & = & 
32 (4 \lambda_\code{1} - 1)\lambda_\code{0}\lambda_\code{1}\lambda_\code{3} &
\mbox{face \code{2}} & \lambda^\code{29} & = & (\frac14,\frac12,0,\frac14)
\\[1mm]
%
\pbar^\code{30}(\lambda) & = & 
32 (4 \lambda_\code{3} - 1)\lambda_\code{0}\lambda_\code{1}\lambda_\code{3} &
\mbox{face \code{2}} & \lambda^\code{30} & = & (\frac14,\frac14,0,\frac12)
\\[1mm]
%
\pbar^\code{31}(\lambda) & = & 
32 (4 \lambda_\code{0} - 1)\lambda_\code{0}\lambda_\code{1}\lambda_\code{2} &
\mbox{face \code{3}} & \lambda^\code{31} & = & (\frac12,\frac14,\frac14,0)
\\[1mm]
%
\pbar^\code{32}(\lambda) & = & 
32 (4 \lambda_\code{1} - 1)\lambda_\code{0}\lambda_\code{1}\lambda_\code{2} &
\mbox{face \code{3}} & \lambda^\code{32} & = & (\frac14,\frac12,\frac14,0)
\\[1mm]
%
\pbar^\code{33}(\lambda) & = & 
32 (4 \lambda_\code{2} - 1)\lambda_\code{0}\lambda_\code{1}\lambda_\code{2} &
\mbox{face \code{3}} & \lambda^\code{33} & = & (\frac14,\frac14,\frac12,0)
\\[1mm]
%
\pbar^\code{34}(\lambda) & = & 
256 \lambda_\code{0}\lambda_\code{1}\lambda_\code{2}\lambda_\code{3} &
\mbox{center} & \lambda^\code{34} & = & (\frac14,\frac14,\frac14,\frac14)
\\[1mm]\hline\end{array}
\end{math}\end{small}\end{center}
\vskip -5mm
\caption{Local basis functions for quartic finite elements in 3d.}
\label{T:P4_3d}
\end{table}

For the implementation of \code{get\_dof\_indices()} for quartics, we
again need a global orientation of the edges on the mesh. At every
edge three DOFs are located, which can then be ordered with respect
to the orientation of the corresponding edge. In 3d, we also need
a global orientation of faces for a one to one mapping of global DOFs located
at a face to local DOFs on an element at that face. Such an orientation
can again be defined by DOF indices at the face's vertices.

\subsubsection{Access to Lagrange elements}
\label{S:get_lagrange}

The Lagrange elements described above are already implemented in \ALBERTA;
access to Lagrange elements is given by the function
\fdx{get_lagrange@{\code{get\_lagrange}}}%
\bv\begin{lstlisting}
const BAS_FCTS  *get_lagrange(int,int);
\end{lstlisting}\ev
Description:
\begin{descr}
\kitem{get\_lagrange(dim, degree)} returns a pointer to 
a filled \code{BAS\_FCTS} structure for Lagrange elements of order
\code{degree}, where $1 \leq \code{degree} \leq 4$, for dimension \code{dim};
no additional call of \code{new\_bas\_fcts()} is needed.
\end{descr}

\subsection{Discontinuous Lagrange finite elements}%
\label{S:get_disc_lagrange}

Similar to the standard Lagrange elements described above, discontinuous 
polynomial finite elements are piecewise polynomial functions. However, 
the functions are not globally continuous. This is implemented in \ALBERTA 
by assigning all DOFs to be of type \code{CENTER}, implying that they will not 
be shared among elements. At the moment these elements are available for 
polynomial degree 0 (piecewise constants), degree 1 (piecewise linears), and 
degree 2 (piecewise quadratics). Much of the implementation is similiar to 
the case of standard Lagrange elements, hence we will not provide a detailed
description.

Access to discontinuous elements is provided by the function
\fdx{get_discontinuous_lagrange@{\code{get\_discontinuous\_lagrange}}}%
\bv\begin{lstlisting}
const BAS_FCTS  *get_discontinuous_lagrange(int dim, int degree);
\end{lstlisting}\ev
Description:
\begin{descr}
\kitem{get\_discontinuous\_lagrange(dim, degree)} returns a pointer to 
a filled \code{BAS\_FCTS} structure for discontinuous polynomial elements of 
order \code{degree}, where $0 \leq \code{degree} \leq 2$ for dimension 
\code{dim}; no additional call of \code{new\_bas\_fcts()} is needed.
\end{descr}

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

\subsection{Discontinuous orthogonal finite elements}
\label{S:get_disc_ortho_poly}

In the context of discontinuous Galerkin methods it is often easier to
use basis-functions which are orthogonal w.r.t. the $L^2$ scalar
product, especially in the context of explicit time discretization
schemes where the use of orthogonal basis functions eliminates the
need for the inversion of the mass-matrix. \ALBERTA implements
discontinuous $L^2$-orthogonal basis functions of degree $1$ and $2$
in all supported mesh-dimensions.

Access to discontinuous elements is provided by the function
\fdx{get_disc_ortho_poly@{\code{get\_disc\_ortho\_poly}}}%
\bv\begin{lstlisting}
const BAS_FCTS *get_disc_ortho_poly(int dim, int degree);
\end{lstlisting}\ev
Description:
\begin{descr}
  \kitem{get\_disc\_ortho\_poly(dim, degree)} returns a pointer to a
  filled \code{BAS\_FCTS} structure for discontinuous polynomial
  elements of order \code{degree}, where $1 \leq \code{degree} \leq 2$
  for dimension \code{dim}; no additional call of
  \code{new\_bas\_fcts()} is needed.
\end{descr}

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

\subsection{Basis-function plug-in module}
\label{S:fancy_bas_fcts}

\ALBERTA also supports a rudimentary plug-in scheme, if
\code{get\_bas\_fcts(dim, name)} cannot find an instance of basis
functions requested by the parameter \code{name}, then it looks for
the environment variable \verb|ALBERTA_BAS_FCTS_LIB_XD| respectively
\verb|ALBERTA_BAS_FCTS_LIB_XD_DEBUG|, where the \code{X} has to be
replaced by the value of \DOW. The environment variable is supposed to
contain the full path to a shared library containing additional
basis-function implementations. The module must define (and export) a
function in \code{C}-calling convention named
\verb|const BAS_FCTS *bas_fcts_init(int dim, dim dow, const char *name)|,
which is called to resolve the request of the application program.  By
pointing to a basis-function plug-in module via that environment
variable it is, e.g., possible to post-process finite element data
using ``fancy'' basis function implementations with the GRAPE or
Paraview interface tools (see \secref{S:gfx}), without having to
recompile and relink the converter-tools (like
\verb|alberta2paraview_Xd|).

Additionally, the function
%%
\fdx{add_bas_fcts_plugin()@{\code{add\_bas\_fcts\_plugin()}}}
\bv\begin{lstlisting}
typedef const BAS_FCTS *
(*BAS_FCTS_INIT_FCT)(int dim, int dow, const char *name);

void add_bas_fcts_plugin(BAS_FCTS_INIT_FCT init_fct);
\end{lstlisting}\ev
%%
can be used to supply additional plug-in functions defining even more
basis-functions.
%%
During the installation of the \ALBERTA package the basis-function
add-on module in \verb|add_ons/libalbas/| is compiled and installed
under the name \bv\begin{lstlisting}
  PREFIX/lib/libalbas_Xd[_debug].DYNEXT,
\end{lstlisting}\ev
where the \code{X} again has to be replaced by the value of \DOW and
\code{DYNEXT} stands for the architecture dependent extension attached
to dynamic libraries.

\verb|libalbas| currently implements wall-bubbles (i.e. face-bubbles
and edge-bubbles), element bubbles and the corresponding trace-spaces.
There is also a very rudimentary and untested version of the lowest
order Raviart-Thomas element. Additionally, the function
\verb|stokes_pair()| implements some of the known stable mixed
discretizations for the Stokes-problem.

\fdx{stokes_pair()@{\code{stokes\_pair()}}}
\bv\begin{lstlisting}
typedef struct stokes_pair STOKES_PAIR;
struct stokes_pair 
{
  const BAS_FCTS *velocity;
  const BAS_FCTS *pressure;
  /* const BAS_FCTS *slip_stress; */
};

STOKES_PAIR stokes_pair(const char *name, unsigned dim, unsigned degree);
\end{lstlisting}\ev

The application must be linked against \verb|libalbas| and include the
header file \verb|albas.h| to use it.  \verb|stokes_pair()| can be
invoked with the symbolic names \code{"Mini"},
\code{"BernardiRaugel"}, \code{"CrouzeixRaviar"} and
\code{"TaylorHood"} and returns the requested Stokes-discretization,
which in all cases except for Taylor-Hood consists of a chain of local
basis functions where the first part of the chain contains the
Lagrange-component and the other parts the "bubbly" add-ons used to
stabilize the resulting Stokes-pair. See also 
\secref{S:bfcts_chains} and \secref{S:chain_impl}.

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

\section{Implementation of finite element spaces}
\label{S:FES-impl}

\subsection{The finite element space data structure}
\label{S:FES-data}

All information about the underlying mesh, the local basis functions,
and the DOFs are collected in the following data structure which defines
one single finite element space:
\ddx{FE_SPACE@{\code{FE\_SPACE}}}
\bv\begin{lstlisting}[name=FE_SPACE,label=T:FE_SPACE]
typedef struct fe_space         FE_SPACE;

struct fe_space
{
  const char      *name;
  const DOF_ADMIN *admin;
  const BAS_FCTS  *bas_fcts;
  MESH            *mesh;
  int             rdim;
  DBL_LIST_NODE   chain;
  const FE_SPACE  *unchained;
};
\end{lstlisting}\ev
Description:
\begin{descr}
  %% 
\kitem{name} holds a textual description of the finite element space.
  Note that \code{name} is duplicated by calling \code{strdup(3)}
  %% 
\kitem{admin} pointer to the DOF administration for the DOFs of this
  finite element space, see Section \ref{S:DOF_ADMIN}.
  %% 
\kitem{bas\_fcts} pointer to the local basis functions, see Section
  \ref{S:basfct_data}.
  %% 
\kitem{mesh} pointer to the underlying mesh, see Section
  \ref{S:mesh_data_structure}.
  %% 
\kitem{rdim} The dimension of the range of the elements of this finite
  element space, as \ALBERTA nowadays supports vector-valued basis
  functions it becomes now important whether a given finite element
  space is actually meant for scalar functions or for vector fields. See
  also \secref{S:vector_bfcts}.
  %% 
\kitem{chain} List pointer to a chain of finite element spaces which
  form a direct sum, see \secref{S:chain_impl}. Such a direct sum is
  based on a chain of local basis functions as described by
  \secref{S:bfcts_chains}.
  %% 
\kitem{unchained} If the finite element space is part of a direct sum
  of finite element spaces (and thus \code{chain} is the link to the
  other elements of this direct sum) then \code{unchained} is a copy
  of the \code{FE\_SPACE} which is unaware of this fact, i.e.
  \code{FE\_SPACE.unchained.chain} points back to itself. If the
  finite element space does not form part of a direct sum, then
  \code{unchained} simply points back to the same \verb|FE_SPACE|. See
  also \secref{S:chain_impl} and \secref{S:basfct_impl}.
\end{descr}

\noindent
Some remarks:
\begin{itemize}
\item Several finite element spaces can be handled on the same mesh.
  Different finite element spaces can use the same DOF administration,
  if they share exactly the same DOFs.
\item Using direct sums of finite element spaces which are chained
  together using the \verb|FE_SPACE.chain|-component has the effect
  that all derived structures are also chains of objects, coefficient
  vectors become chains of coefficient vectors, matrices become block
  matrices, where the blocks are chained together using chains for
  rows and columns. The same holds for the per-element vectors and
  matrices.
\item \ALBERTA provides full support for these chains in its
  infra-structure for the assembling of the discrete systems, as well
  as in the solver infra-structure and in the support functions for
  the computation of errors and error estimates.
\end{itemize}

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

\subsection{Access to finite element spaces}%
\label{S:access_fe_space}

A finite element space can only be accessed by the function
\fdx{get_fe_space()@{\code{get\_fe\_space()}}}
\bv\begin{lstlisting}
const FE_SPACE  *get_fe_space(MESH *mesh, const char *name,
                              const BAS_FCTS *bas_fcts, int rdim,
                              FLAGS adm_flags);
const FE_SPACE  *get_dof_space(MESH *mesh, const char *name,
                               const int ndof[N_NODE_TYPES],
                               FLAGS adm_flags);
\end{lstlisting}\ev
\paragraph{Descriptions}
\begin{descr}
\kitem{get\_fe\_space(mesh, name, bas\_fcts, rdim, adm\_flags)} defines a new
  finite element space on \code{mesh}; it looks for an existing
  \code{dof\_admin} defined on \code{mesh}, which manages DOFs
  uniquely defined by \code{bas\_fcts->dof\_admin->n\_dof}
  (compare Section \ref{S:DOF_ADMIN});
  if such a \code{dof\_admin} is not found, a new \code{dof\_admin} is
  created. 
  \begin{description}
  \item[Parameters]~\hfill

    \begin{descr}
    \kitem{mesh}
      A pointer to the underlying triangulation.
      %% 
    \kitem{name}
      A fancy name, used for pretty-printing and debugging.
      %% 
    \kitem{bas\_fcts} The underlying basis functions. If
      \code{bas\_fcts} is a disjoint union, or ``chain'', of local
      basis functions sets, then the resulting finite element space
      will be the direct sum of the finite element spaces defined by
      the respective components of the disjoint union of local basis
      function sets, see \secref{S:chain_impl} and
      \secref{S:basfct_impl}. The component \code{unchained} of the
      \code{FE\_SPACE} structure will -- for each component of the
      direct sum -- point to an \code{FE\_SPACE} instance which is
      ignorant of the fact that it forms part of a direct sum of
      finite element spaces.
      %% 
    \kitem{rdim} The dimension of the range of the elements of the
      finite element space. Now that \ALBERTA also support \DOW-valued
      basis functions, this parameter plays an important role; without
      attaching a ``range-dimension'' to a finite element space it
      would be hardly possible to assemble discrete systems in a
      consistent manner. Of course, if the underlying basis functions
      are vector-valued for themselves, then \code{rdim} has to equal
      \DOW as well. If the underlying local basis functions are
      scalar-valued, then \code{rdim} may be either $1$ to generate a
      finite element space consisting of scalar functions, or \DOW to
      generate a finite element space consisting of \DOW-valued
      functions. Other values for \code{rdim} are not supported.
      %% 
    \kitem{adm\_flags} Currently \verb|adm\_flags| is the bit-wise
      or of \verb|ADM_PRESERVE_COARSE_DOFS| and/or
      \verb|ADM_PERIODIC|. If the flag \verb|ADM_PRESERVE_COARSE_DOFS|
      is set, then it requests that \DOFs normally be deleted during
      refinement should be preserved instead. This is necessary for
      higher order multi-grid implementations as well as for the
      internal maintenance of submeshes, see
      \secref{S:submesh_implementation}. For a detailed description of
      which DOFs would normally be deleted, see \secref{S:DOFs2}.

      \verb|ADM_PERIODIC| requests a periodic finite element space
      where \DOFs across periodic walls are identified. See
      \secref{S:periodic}.
    \end{descr}
  \item{return value} The return value is a newly created
    \code{FE\_SPACE} structure, where \code{name} is duplicated, and
    the members \code{mesh}, \code{bas\_fcts} and \code{admin} are
    adjusted correctly.
  \end{description}
  %% 
\kitem{get\_dof\_space(mesh, name, n\_dof, adm\_flags)}
  %% 
  performs a similar task as \verb|get_fe_space()|, however, the
  resulting ``space'' is not bound to a \verb|BAS_FCTS| instance.
  Instead, the argument \verb|n_dof| determines the distribution of
  the \DOFs across the sub-simplices. The elements of \code{n\_dof}
  determine how many degrees of freedom are tied to each sub-simplex
  on each element, the mapping is defined by the \code{NODE\_TYPE}
  enumeration type, see the source code listing on page
  \pageref{T:NODE_TYPES}, i.e. \code{VERTEX == 0}, \code{CENTER == 1},
  \code{EDGE == 2}, \code{FACE == 3}. The following code-snippet
  defines an ``\code{FE\_SPACE}'' with $42$ \code{DOF}s on each face:
\bv\begin{lstlisting}
   int n_dof[N_NODE_TYPES] = {0, 0, 0, 42};
   const FE_SPACE *face_dof_space;

   face_dof_space =
     get_dof_space(mesh, "face dofs", n_dof, ADM_FLAGS_DFLT);
\end{lstlisting}\ev
\end{descr}

The selection of finite element spaces defines the DOFs that must
be present on the mesh elements. For each finite element space
there must be a corresponding DOF administration, having information
about the used DOFs. Each call of \code{get\_fe\_space()} requires the
internal adjustment of the \code{el->dof} pointer arrays which is 
potentially expensive, since information about which elements share 
a vertex/edge/face must be calculated for the current triangulation. It is 
therefore advisable (but not necessary) to allocate all finite element 
spaces before refining the mesh.

Since a mesh only gives access to the \code{DOF\_ADMIN}s
defined on it, the user has to store pointers to the
\code{FE\_SPACE} structures in some global variable; no access to
\code{FE\_SPACE}s is possible via the underlying mesh.

\begin{example}[Initializing DOFs for Stokes and Navier--Stokes]
Now, as an example we look at a possible \code{main} function.
In the example we want to define two 
finite element spaces on the mesh, for a mixed finite element
formulation of the Stokes or Navier--Stokes equations with the Taylor--Hood
element, e.g. we want to use Lagrange finite elements of
order \code{degree} (for the velocity) and $\code{degree} - 1$
(for the pressure). Pointers to the corresponding \code{FE\_SPACE} 
structures are stored in the global variables \code{u\_fe} and \code{p\_fe}.
\bv\begin{lstlisting}
static FE_SPACE *u_fe, *p_fe;
static int       degree, dim;

int main()
{
  const MESH      *mesh;
  const BAS_FCTS  *lagrange;
  MACRO_DATA      *data;

  TEST_EXIT(degree > 1)("degree must be greater than 1\n");

  ...

  data = read_macro(filename);
  mesh = GET_MESH(dim, "ALBERTA mesh", data, NULL, NULL);
  free_macro_data(data);

  lagrange = get_lagrange(mesh->dim, degree);
  u_fe = get_fe_space(
    mesh, "Velocity space", lagrange, DIM_OF_WORLD, ADM_FLAGS_DFLT);

  lagrange = get_lagrange(mesh->dim, degree-1);
  p_fe = get_fe_space(mesh, "Pressure space", lagrange, 1, ADM_FLAGS_DFLT);

  ...

  return;
}
\end{lstlisting}\ev
This will provide all DOFs for the two finite element spaces on
the elements and the corresponding \code{DOF\_ADMIN}s will have
information about the access to local DOFs for both finite element
spaces.

It is also possible to define only one or even more finite element spaces;
the use of special user defined basis functions is possible too.
These should be added to the list of all used basis functions
by a call of \code{new\_bas\_fcts()} before allocating finite element 
spaces.
\end{example}

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