\section{Data structures for parametric meshes}%
\label{S:parametric_meshes}
\idx{parametric meshes}%

The current version of \ALBERTA offers support for so-called
\emph{parametric meshes} which are triangulations where some or all of
the simplices are non-linear images of the reference element.
Typically, the transformation from the reference element $\Shat$ to
the curved simplex $S$ is a polynomial, but in principle this need not
be the case. \ALBERTA has predefined polynomial parameterisations up
to polynomial degree $4$: $S = F_S(\Shat)$, $F_S\in\P_k(\Shat)$ for
$k=1,2,3,4$. The limitation $k\leq 4$ just means that piecewise
polynomial parameterisations up to the maximal degree for the Lagrange
basis functions within \ALBERTA are supported (Section
\ref{S:basfct_impl}).

The standard case for applications is the iso-parametric approximation
of curved boundaries; care has to be taken when the polynomial degree
of the parameterisation is so high that some of the Lagrange-nodes
fall into the interior of the simplex. \ALBERTA implements the
algorithm developed in \cite{Lenoir:86}. The suite of demo-programs
shipped with the \ALBERTA-package contains a program called
\code{ellipt-isoparam}, which implements the discretization of
Poisson's equation on an iso-parametric triangulation of a unit-disc.

Many other applications besides isoparametric boundary approximation
are conceivable, for example in moving finite elements, where the
positions of nodes may change with time and need to be described by a
time dependent parameterization. Stationary example programs for $1$,
$2$ and $3$ dimensional parametric meshes can again be found in the
demo-suite:
\begin{description}
\item[\code{src/Common/ellipt-sphere.c}] Poisson's equation on the
  $1$-, $2$- and $3$-dimensional unit-sphere, i.e.
  $S^k\subset\R^{k+1}$ $(1\leq k \leq 3)$.
\item[\code{src/Common/ellipt-torus.c}] Poisson's equation on the $1$-
  , $2$- and $3$-torus, i.e. $T^k\subset\R^{k+1}$ $(1\leq k \leq 3)$.
\item[\code{src/3d/ellipt-moebius.c}] Poisson's equation on an
  embedded Moebius-strip (yes, \ALBERTA can handle unorientable meshes).
\item[\code{src/4d/ellipt-klein-bottle.c}] Embedded Klein's bottle.
\item[\code{src/5d/ellipt-klein-3-bottle.c}] Embedded non-orientable
  $3$-manifold in $\R^5$, similar to a Klein's bottle, but one
  dimension higher.
\end{description}

Using parametric elements does not imply a fundamental change of data
structures within \ALBERTA. The mesh still consists of a hierarchical
collection of \code{EL} structures, however these only represent the
topological structure of the mesh. The coordinate and shape
information of all elements, standard or parametric, is stored using
an internal \code{DOF\_REAL\_D\_VEC coords} representing the global
parametrization encoded in $F_S$ for all $S$. The finite element space
containing \code{coords} is a standard Lagrange space of order $1$,
$2$, $3$ or $4$. 

A mesh may be turned into a parametric mesh with piece-wise polynomial
parameterization by calling the function
\hyperref[S:use_lagrange_parametric_fct]{\code{use\_lagrange\_parametric()}}
described below. This allocates \code{coords} and turns some or all
mesh elements into parametric simplices, depending on the options
determined by the user. The shape of the parametric simplices is
furthermore uniquely determined by the value of \code{coords} at the
Lagrange nodes. There are interface routines
\hyperref[S:get_lagrange_coords_fct]{\code{get\_lagrange\_coords()}},
\hyperref[S:copy_lagrange_coords_fct]{\code{copy\_lagrange\_coords()}}
and
\hyperref[S:get_lagrange_touched_edges_fct]{\code{get\_lagrange\_touched\_edges()}}
to give an application access to the coordinate data, see below in
Sections
\ref{S:use_lagrange_parametric_fct}-\ref{S:get_lagrange_touched_edges_fct}.

%The parametrization \code{coords} is initialized with the original mesh 
%coordinates as determined by the macro triangulation. The shape of the 
%parametric simplices is then calculated by projecting the Lagrange nodes 
%using a \code{NODE\_PROJECTION}, see \secref{S:node_projections}.

Note that on curved elements the ordinary routines to convert between
barycentric coordinates and Cartesian coordinates, or to compute their
derivatives (see \secref{S:bary_routines}), may no longer be used.
Instead, the corresponding hooks in the \code{PARAMETRIC}-structure
described below have to be called. It may be convenient in this case
to use calls to the per-element quadrature caches (see
\secref{S:fill_quad_el_cache}). An exception is the case of
affine-linear ``parametric'' meshes, or the case of affine-linear mesh
elements of only partially parametric meshes: there the standard
routines described in \secref{S:bary_routines} may still be used.

We start with a more detailed description of how to use ``standard''
piece-wise polynomial parameterizations and continue with the
description of the general interface in Section
\secref{S:PARAMETRIC_struct} further below.

%%
% Figure was created using my code "write_mesh_fig()", DK.
\begin{figure}[tbp]
\centering
\includegraphics[scale=0.75]{EPS/refine_parametric}
\caption[Parametric meshes, triangulation of a disc]{Successive
  refinements of the triangulation of a disc with \code{strategy ==
    PARAM\_STRAIGHT\_CHILDS}.  Parametric simplices are shaded in
  gray.}%
\label{F:refine_parametric}
\end{figure}
%
\subsection{Piece-wise polynomial parametric meshes}%
\label{S:access_param_mesh}%
\idx{parametric meshes!accessing}

The following functions are available to access and manipulate meshes
with ``standard'' piece-wise polynomial parameterizations:
%%
\fdx{use_lagrange_parametric()@{\code{use\_lagrange\_parametric()}}}
\fdx{get_lagrange_coords()@{\code{get\_lagrange\_coords()}}}
\fdx{get_lagrange_touched_edges()@{\code{get\_lagrange\_touched\_edges()}}}
\fdx{copy_lagrange_coords()@{\code{copy\_lagrange\_coords()}}}
\bv\begin{lstlisting}
typedef enum param_strategy {
  PARAM_ALL = 0,
  PARAM_CURVED_CHILDS = 1,
  PARAM_STRAIGHT_CHILDS = 2
} PARAM_STRATEGY;
#define PARAM_PERIODIC_COORDS 0x04

void use_lagrange_parametric(MESH *mesh, int degree,
                             NODE_PROJECTION *n_proj, FLAGS flags);
DOF_REAL_D_VEC *get_lagrange_coords(MESH *mesh);
DOF_UCHAR_VEC *get_lagrange_touched_edges(MESH *mesh);
void copy_lagrange_coords(MESH *mesh, DOF_REAL_D_VEC *coords, bool to_mesh);
\end{lstlisting}\ev
%%
\begin{figure}[tbp]
\centering
\input EPS/curved_simplex
\caption[Paremetric meshes, transformation to the reference
element]{Mapping of the standard simplex under a quadratic
  transformation $F_S$ with standard numbering of the local Lagrange
  nodes. The curve $\lambda_0=\lambda_1$ is shown dashed.}
\label{F:curved_simplex}
\end{figure}
%%
\begin{function}{use\_lagrange\_parametric()}
\label{S:use_lagrange_parametric_fct}
%%
\fdx{use_lagrange_parametric()@{\code{use\_lagrange\_parametric()}}|(}
\idx{parametric meshes!use_lagrange_parametric()@{\code{use\_lagrange\_parametric()}}|(}

\item[Prototype] ~\hfill
%%
\bv\begin{lstlisting}
void use_lagrange_parametric(MESH *mesh, int degree,
                             NODE_PROJECTION *selective, FLAGS strategy);

typedef enum param_strategy {
  PARAM_ALL = 0, PARAM_CURVED_CHILDS = 1, PARAM_STRAIGHT_CHILDS = 2
} PARAM_STRATEGY;

#define PARAM_PERIODIC_COORDS 0x04
\end{lstlisting}\ev
\item[Synopsis] ~\hfill

\bv\begin{lstlisting}[basicstyle=\normalsize]
use_lagrange_parametric(mesh, degree, selective, strategy);
\end{lstlisting}\ev
\item[Description] ~\hfill

  Convert the given \code{mesh} into a parametric mesh. The mesh may
  already be refined. Parametric simplices will be the image of the
  reference simplex under a polynomial transformation of specified
  \hyperlink{use_lagrange_parametric:degree}{\code{degree}}. The
  maximal value of
  \hyperlink{use_lagrange_parametric:degree}{\code{degree}}
  is limited by the maximal degree of
  the Lagrange basis functions implemented in \ALBERTA (currently
  $4$).  Internally a coordinate vector \code{coords} is allocated
  within the standard Lagrange finite element space of order
  \hyperlink{use_lagrange_parametric:degree}{\code{degree}}
  Specifying \code{1} means that simplices will still
  be the images of an affine transformation, which is useful for
  special applications.

  The \code{coords} vector employs special \code{refine\_interpol} and
  \code{coarse\_restrict} entries to enable the described refinement
  of curved simplices. Concerning the coarsening of the mesh, all
  parents of parametric elements are automatically parametric elements
  themselves. The information describing the shape of children is
  passed back up to parents in a straight forward fashion.

  The function generates a filled \code{PARAMETRIC} structure and sets
  the entry \code{mesh->parametric} to point at it. Only one call of
  the function is possible per mesh. If the mesh belongs to a
  submesh-hierarchy, then \code{use\_lagrange\_parametric()} must be
  called on the top-level master mesh. The sub-meshes will then
  inherit the parametric structure from the top-level master mesh.
  Sub-meshes are discussed in Section \ref{S:submesh_implementation}.

  When \code{use\_lagrange\_parametric()} is invoked, then this will
  initiate a mesh-traversal to initialize the coordinate vector
  \code{coords} mentioned above. On all curved elements -- see the
  parameter
  \hyperlink{use_lagrange_parametric:selective}{\code{selective}} --
  the corresponding projection routine will be invoked to project the
  affine (non-curved) coordinates of the Lagrange nodes to whatever
  manifold is defined by the projection function. As described in
  \secref{S:node_projections} \ALBERTA allows for a default projection
  for the entire element, or for distinct projections attached to the
  ``walls'' of the elements.

\item[Parameters]~\hfill
  \begin{descr}
  \hyperitem{use_lagrange_parametric:mesh}{mesh} The mesh to be
    equipped with a parametric structure.
    %% 
  \hyperitem{use_lagrange_parametric:degree}{degree} The degree of the
    parameterization. Currently, the maximum degree is $4$, limited
    only by the maximum degree of the Lagrange basis functions
    implemented in \ALBERTA. \ALBERTA takes special care --
    implementing the algorithm explained in \cite{Lenoir:86} -- that
    higher degree iso-parametric boundary approximation will yield
    optimal convergence rates.
    %% 
  \hyperitem{use_lagrange_parametric:selective}{selective} Optional,
    maybe \nil. If non-\nil, then \ALBERTA only treats those elements
    as curved ones which carry exactly this
    \hyperref[T:NODE_PROJECTION]{\code{NODE\_PROJECTION}} structure.
    If \code{selective == \nil}, then all elements carrying a
    projection routine (see \secref{S:node_projections}) will be
    treated as curved elements.
    %%
  \hyperitem{use_lagrange_parametric:strategy}{strategy} 

    The parameter \code{strategy} splits in two parts: \code{(strategy
      \& PARAM\_STRATEGY\_MASK)} determines which newly created
    simplices are treated as parametric simplices during refinement of
    the mesh.  The remaining flag is \code{PARAM\_PERIODIC\_COORDS}.
    It determines whether the finite element function which holds the
    coordinate information of the parametric mesh is itself a periodic
    function.  The demo-program
    \code{demo/src/4d/ellipt-klein-bottle.c} contains an example
    application.

    The following values are defined for \code{(strategy \&
      PARAM\_STRATEGY\_MASK)}.
 \begin{descr}
    \kitem{PARAM\_ALL}
    All elements of the mesh will be treated as parametric
    elements, implying that determinants and Jacobeans will
    be calculated at all quadrature points during assembly.
    This is useful e.g. for triangulations of embedded
    curved manifolds. Please note that during refinement a
    parent element will be split along the surface defined
    by the equation $lambda_0 = lambda_1$.
 
    \kitem{PARAM\_CURVED\_CHILDS} Only those elements of the mesh
    affected by \code{n\_proj} will be treated as parametric elements.
    Simplices are split along the surface $lambda_0 = lambda_1$ during
    mesh refinement. Using \code{PARAM\_CURVED\_CHILDS} should be
    avoided for parameterisations of degree $> 2$, maybe it should not
    be used at all.

    \kitem{PARAM\_STRAIGHT\_CHILDS} Only those elements of the mesh
    affected by \code{n\_proj} will be treated as parametric elements.
    \code{PARAM\_STRAIGHT\_CHILDS} should be used for the
    approximation of curved boundaries. This keeps the number of
    curved simplices as small as possible and \ALBERTA takes care to
    position the Lagrange nodes of the parametric elements such that
    optimal approximation order can be achieved; this is not trivial,
    see \cite{Lenoir:86}.
  \end{descr}
    %%
  \end{descr}
\item[Examples] ~\hfill

  See below \exampleref{example:isoparam_unit_ball}.
\end{function}
%%
\fdx{use_lagrange_parametric()@{\code{use\_lagrange\_parametric()}}|)}
\idx{parametric meshes!use_lagrange_parametric()@{\code{use\_lagrange\_parametric()}}|)}
%%

\begin{function}{get\_lagrange\_coords()}
\label{S:get_lagrange_coords_fct}
%%
\fdx{get_lagrange_coords()@{\code{get\_lagrange\_coords()}}|(}
\idx{parametric meshes!get_lagrange_parametric()@{\code{get\_lagrange\_coords()}}|(}

\item[Prototype] ~\hfill
%%
\bv\begin{lstlisting}
DOF_REAL_D_VEC *get_lagrange_coords(MESH *mesh);
\end{lstlisting}\ev

\item[Synopsis] ~\hfill

\bv\begin{lstlisting}[basicstyle=\normalsize]
coord_dof_vec = get_lagrange_coords(mesh);
\end{lstlisting}\ev
\item[Description] ~\hfill

  Returns the internal \code{DOF\_REAL\_D\_VEC coords} used to store
  the coordinates of parametric elements. The user may change entries
  of this vector by hand, if some care is used if the parametric mesh
  was initialized with \code{strategy != PARAM\_ALL}. See below the
  description for
  \hyperref[S:get_lagrange_touched_edges_fct]{\code{get\_lagrange\_touched\_edges()}}.

  See also
  \hyperref[S:copy_lagrange_coords_fct]{\code{copy\_lagrange\_coords()}}
  for a more secure interface to the coordinate information.
\item[Parameters]~\hfill
  \begin{descr}
  \hyperitem{get_lagrange_coords:mesh}{mesh} A
    \hyperref[T:MESH]{\code{mesh}}-structure carrying a parametric
    structure previously initialized by a call to
    \hyperref[S:use_lagrange_parametric_fct]{\code{use\_lagrange\_parametric()}},
    see \secref{S:use_lagrange_parametric_fct} above.
    %%
  \end{descr}
\item[Return Value] ~\hfill
  
  A pointer to the underlying coordinate function, a
  \code{DOF\_REAL\_D\_VEC} belonging to a finite element space of the
  piece-wise polynomial degree as specified by the
  \hyperlink{use_lagrange_parametric:degree}{\code{degree}} parameter
  passed to
  \hyperref[S:use_lagrange_parametric_fct]{\code{use\_lagrange\_parametric()}}.
%%\item[Examples] ~\hfill
\end{function}
%%
\fdx{get_lagrange_coords()@{\code{get\_lagrange\_coords()}}|)}
\idx{parametric meshes!get_lagrange_parametric()@{\code{get\_lagrange\_coords()}}|)}
%%

\begin{function}{copy\_lagrange\_coords()}
\label{S:copy_lagrange_coords_fct}
%%
\fdx{copy_lagrange_coords()@{\code{copy\_lagrange\_coords()}}|(}
\idx{parametric meshes!copy_lagrange_parametric()@{\code{copy\_lagrange\_coords()}}|(}

\item[Prototype] ~\hfill
%%
\bv\begin{lstlisting}
typedef enum param_copy_direction {
  COPY_FROM_MESH = false,
  COPY_TO_MESH   = true
} PARAM_COPY_DIRECTION;

void copy_lagrange_coords(MESH *mesh, DOF_REAL_D_VEC *coords, bool to_mesh);
\end{lstlisting}\ev

\item[Synopsis] ~\hfill

\bv\begin{lstlisting}[basicstyle=\normalsize]
copy_lagrange_coords(mesh, coord_copy, to_mesh);
\end{lstlisting}\ev
\item[Description] ~\hfill

  This is the recommended interface to the coordinate information for
  (Lagrange-) parametric meshes. Only the coordinate \emph{values} are
  copied; the function also makes sure that affine elements remain
  affine by using linear interpolation between the vertices of a
  simplex if that simplex has no curved edge. The state of the edges
  is determined by the \code{touched\_edges} vector returned by
  \hyperref[S:get_lagrange_touched_edges_fct]{\code{get\_lagrange\_touched\_edges()}},
  see below. \code{copy\_lagrange\_coords()} handles also a case when
  a mesh has no parametric structure, but uses \code{EL->new\_coord}
  to store coordinate information for the vertices, see
  \secref{S:node_projections}. See also
  \hyperref[S:get_lagrange_coords_fct]{\code{get\_lagrange\_coords()}}.
\item[Parameters]~\hfill
  \begin{descr}
  \hyperitem{copy_lagrange_coords:mesh}{mesh} A
    \hyperref[T:MESH]{\code{mesh}}-structure carrying a parametric
    structure previously initialized by a call to
    \hyperref[S:use_lagrange_parametric_fct]{\code{use\_lagrange\_parametric()}},
    see \secref{S:use_lagrange_parametric_fct} above.
    %%
  \hyperitem{copy_lagrange_coords:coord_copy}{coord\_copy} A
    \code{DOF\_REAL\_D\_VEC}, storage for the coordinate information.
    Note that \code{coord\_copy} is not itself installed as coordinate
    vector in the mesh, just the coordinate data is copied to and from
    \code{coord\_copy}, where the direction of the copy-operation is
    specified by the parametric \code{to\_mesh}, see below.
    %%
  \hyperitem{copy_lagrange_coords:to\_mesh}{to\_mesh} If \code{true},
    then the coordinate data is copied from \code{coord\_copy} to the
    mesh, otherwise the coordinate function of the mesh is copied to
    \code{coord\_copy}.
  \end{descr}
%\item[Return Value] ~\hfill
%%\item[Examples] ~\hfill
\end{function}
%
\fdx{copy_lagrange_coords()@{\code{copy\_lagrange\_coords()}}|)}
\idx{parametric meshes!copy_lagrange_parametric()@{\code{copy\_lagrange\_coords()}}|)}
%%

\begin{function}{get\_lagrange\_touched\_edges()}
\label{S:get_lagrange_touched_edges_fct}
%%
\fdx{get_lagrange_touched_edges()@{\code{get\_lagrange\_touched\_edges()}}|(}
\idx{parametric meshes!get_lagrange_parametric()@{\code{get\_lagrange\_touched\_edges()}}|(}

\item[Prototype] ~\hfill
%%
\bv\begin{lstlisting}
DOF_UCHAR_VEC *get_lagrange_touched_edges(MESH *mesh);
\end{lstlisting}\ev

\item[Synopsis] ~\hfill

\bv\begin{lstlisting}[basicstyle=\normalsize]
touched_edges_vec = get_lagrange_touched_edges(mesh);
\end{lstlisting}\ev
\item[Description] ~\hfill

  Returns the internally used \code{DOF\_UCHAR\_VEC touched\_edges}.
  Internally \ALBERTA maintains the ``projection-state'' of all edges.
  \code{1} means that the corresponding edge has suffered a
  projection, \code{0} means that it is still in the affine linear
  state. A simplex is treated as parametric simplex if and only if any
  of its edges has been projected. Otherwise a simplex is not curved.
  The flags vector is only used if \code{strategy != PARAM\_ALL}, this
  function will produce a warning and return \nil if \code{strategy ==
    PARAM\_ALL}.

  When changing the coordinate vector returned by
  \hyperref[S:get_lagrange_coords_fct]{\code{get\_lagrange\_coords()}}
  it falls into the responsibility of the application to also change
  the projection status of the edges.
\item[Parameters]~\hfill
  \begin{descr}
  \hyperitem{get_lagrange_touched_edges:mesh}{mesh} A
    \hyperref[T:MESH]{\code{mesh}}-structure carrying a parametric
    structure previously initialized by a call to
    \hyperref[S:use_lagrange_parametric_fct]{\code{use\_lagrange\_parametric()}},
    see \secref{S:use_lagrange_parametric_fct} above.
    %%
  \end{descr}
\item[Return Value] ~\hfill

  A pointer to a \code{DOF\_SCHAR\_VEC}, with one DOF per edge,
  indicating whether the respective edge is curved or not, with
  \code{touched\_edges->vec[dof] == true} meaning the edge is curved
  and \code{touched\_edges->vec[dof] == false} meaning the edge is \emph{not}
  curved.
%%\item[Examples] ~\hfill
\end{function}
%%
\fdx{get_lagrange_touched_edges()@{\code{get\_lagrange\_touched\_edges()}}|)}
\idx{parametric meshes!get_lagrange_parametric()@{\code{get\_lagrange\_touched\_edges()}}|)}
%%

\begin{example}[Isoparametric elements for the unit ball]
\label{example:isoparam_unit_ball}
We turn again to the triangulation of the unit ball treated in Example 
\ref{Ex:unit_ball}. 
\idx{parametric meshes!isoparametric elements for the unit ball}
\bv\begin{verbatim}
static void ball_proj_func(REAL_D x,
                           const EL_INFO *el_info, const REAL_B lambda)
{
  SCAL_DOW(1.0/NORM_DOW(x), x);
}

static NODE_PROJECTION ball_proj = {ball_proj_func};

static NODE_PROJECTION *init_node_proj(MESH *mesh, MACRO_EL *mel, int c)
{
  if(c > 0 && !mel->neigh[c-1])
    return &ball_proj;
  else
    return NULL;
}

int main()
{
  MESH           *mesh;
  const BAS_FCTS *bas_fcts;
  const FE_SPACE *fe_space;
  MACRO_DATA     *data;

  ...

  data = read_macro("ball.amc");
  mesh = GET_MESH(MESH_DIM, "ALBERTA mesh", data,
		  init_node_proj, NULL /* init_wall_trafos */);
  free_macro_data(data);

  bas_fcts = get_lagrange(mesh->dim, /* degree == */ 3);
  use_lagrange_parametric(mesh, 3, NULL, PARAM_STRAIGHT_CHILDS);

  ...
}
\end{verbatim}\ev
\ALBERTA compares the node-projections of all elements with the value
of \code{\&ball\_proj}, in our example only the boundary faces will
have produce a match. Since \code{strategy ==
  PARAM\_STRAIGHT\_CHILDS}, \ALBERTA will only use parametric elements
in a narrow boundary layer, see \figref{F:refine_parametric}.
\end{example}

\subsection{The \code{PARAMETRIC} structure}%
\idx{PARAMETRIC structure@{\code{PARAMETRIC} structure}}
\label{S:PARAMETRIC_struct}

A parametric mesh is described by the structure
\ddx{PARAMETRIC@{\code{PARAMETRIC}}} \code{PARAMETRIC}. The structure
is a collection of function pointers -- ``methods'' -- which define
the parameterisation. The piecewise polynomial parameterisations
predefined in \ALBERTA work in arbitrary co-dimension, see
\secref{S:access_param_mesh} above.

\bv\begin{verbatim}
typedef struct parametric       PARAMETRIC;

struct parametric 
{
  char *name;
  bool not_all;
  bool use_reference_mesh;
  bool (*init_element)(const EL_INFO *el_info, const PARAMETRIC *parametric);

  void (*coord_to_world)(const EL_INFO *info, const QUAD *quad,
                         int n, const REAL_B lambda[], REAL_D *world);
  void (*world_to_coord)(const EL_INFO *info, int n,
			 const REAL_D world[],
			 REAL_B lambda[], int *k);
  void (*det)(const EL_INFO *info, const QUAD *quad,
              int n, const REAL_B lambda[], REAL dets[]);
  void (*grd_lambda)(const EL_INFO *info,  const QUAD *quad,
                     int n, const REAL_B lambda[],
                     REAL_BD Lambda[], REAL_BDD DLambda[], REAL dets[]);
  void (*grd_world)(const EL_INFO *info,  const QUAD *quad,
		    int n, const REAL_B lambda[],
		    REAL_BD grd_Xtr[], REAL_BDB D2_Xtr[], REAL_BDBB D3_Xtr[]);
  void (*wall_normal)(const EL_INFO *el_info, int wall,
		      const QUAD *wall_quad,
		      int n, const REAL_B lambda[],
		      REAL_D nu[], REAL_DB grd_nu[], REAL_DBB D2_nu[],
		      REAL dets[]);
  void (*inherit_parametric)(MESH *slave);
  void (*unchain_parametric)(MESH *slave);
  void *data;
};
\end{verbatim}\ev
Description:
\begin{descr}
  \kitem{name} a textual description of the parametric structure,
  intended as debugging aid.

  \kitem{not\_all}, if nonzero, signifies that not all of the mesh
  elements are to be parametric (curved) simplices. This entry must
  not be changed by the application program.

  \kitem{use\_reference\_mesh}, if set, means that certain routines
  should use the reference triangulation consisting of standard
  simplices instead of the parametric mesh, see the description
  further below.  Is set to \false by default.

  \idx{Per-element initializers!PARAMETRIC@{\code{PARAMETRIC}}}
  \idx{init_element()@{\code{init\_element()}}!PARAMETRIC@{\code{PARAMETRIC}}}
  \kitem{init\_element(el\_info, parametric)} This is a per-element
  initialiser which must be called for each \code{el\_info} during a
  mesh traversal before calling any other function hook of the
  \code{PARAMETRIC} structure.  The argument \code{parametric} must
  point to the \code{PARAMETRIC} structure itself.

  A specific implementation of a parametric mesh should use the
  \code{init\_element()}-hook to perform all necessary initialisations
  needed to define the transformation from the reference element to
  the given mesh element. The return value should be \code{true} if
  the given element indeed is curved, and \code{false} if it is just
  an affine image of the reference element. In the latter case
  \code{init\_element(el\_info, \dots)} is supposed to fill
  \code{el\_info->coord} with the current element's coordinate
  information -- despite the fact that the \code{el\_info} argument
  carries the \code{const} attribute. This way the normal per-element
  functions can be used (e.g. \code{el\_det()},
  \code{el\_grd\_lambda()} etc.)  instead of the parametric
  replacements defined in the \code{PARAMETRIC} structure. This
  simplifies the program flow (and source code) for applications using
  only partially parametric meshes a lot.

  \kitem{coord\_to\_world(el\_info, quad, n, lambda, world)}
  Implements the function $F_S$ itself. Given an element
  \code{el\_info}, a vector of barycentric coordinates \code{lambda}
  of length \code{n}, this function writes the corresponding vector of
  length \code{n} of world coordinates into \code{world}.  Using this
  function on multiple sets of coordinates at once may be more
  efficient than repeatedly calling this function. If the \code{quad}
  attribute is not \nil, then \code{quad->n\_points} and
  \code{quad->lambda} instead of \code{n} and \code{lambda}.
  Additionally, a specific parametric implementation may handle the
  case \code{quad != } \nil more efficiently by using caching
  \code{QUAD\_FAST} quadratures and the like.

\kitem{world\_to\_coord()} This entry replaces the standard
  \code{world\_to\_coord()} function available for standard simplices.
  It represents the inverse $F_S^{-1}$. Currently, there is only a
  partial implementation available, which may or may not work in the
  context of iso-parametric boundary approximation.

\kitem{det(el\_info, quad, n, lambda, dets)} This function computes
  $|\det DF_S(\hat x(\lambda))|$ which is required for numerical
  integration, see Remark \ref{book:R:numerical_int}. The barycentric
  coordinates are again passed as an array \code{lambda} of length
  \code{n}. The absolute value of the determinant at each $\lambda$ is
  written into the array \code{dets}. Since this routine is mostly
  used for numerical integration the user may pass a pointer
  \code{quad} to a quadrature structure instead of \code{lambda}. The
  function will then calculate the determinants at all quadrature
  nodes of the given numerical quadrature. Additionally, a specific
  parametric implementation may handle the case \code{quad != } \nil
  more efficiently by using caching \code{QUAD\_FAST} quadratures and
  the like.  See \secref{S:quad_data} for details on using numerical
  quadrature routines and structures.

\kitem{grd\_lambda(el\_info, quad, n, lambda, Lambda, DLambda, dets)}
  This routine is similar to the entry \code{dets} above. It
  additionally fills the array \code{Lambda} with the values of the
  derivative $\Lambda_S$ of the barycentric coordinates defined in
  \secref{book:S:eval_Dfe}. Optionally, \code{grd\_lambda()} also
  computes the second derivatives of the barycentric coordinates. The
  second derivatives of the barycentric coordinates are necessary to
  compute the second derivatives of finite element functions on curved
  simplices, e.g. for the implementation of residual error estimators.
  %% 
  The arguments \code{DLambda} and \code{dets} may be \nil.

  \kitem{grd\_world(el\_info, quad, n, lambda, grd\_Xtr, D2\_Xtr,
    D3\_Xtr)} Compute the derivatives of the Cartesian coordinates with
  respect to the barycentric coordinates. The arguments \code{D2\_Xtr}
  and \code{D3\_Xtr} may be \nil, in which case the quantities are
  simply not computed. The \code{tr}-suffix stands for ``transposed'',
  meaning that actually the transposed of the Jacobians is computed.
  This way, in the affine linear case \code{grd\_Xtr} is just the
  matrix formed by the vertex coordinates as rows.

  \kitem{wall\_normal(el\_info, wall, quad, n, lambda, nu, grd\_nu,
    D2\_nu, dets)} This function hook is the parametric replacement
  for library function \code{get\_wall\_normal()}. Again,
  \code{quad->lambda} and \code{quad->n\_points} is used instead of
  \code{lambda} and \code{n} if \code{quad !=} \nil. \code{quad} must
  be a co-dimension $1$ quadrature as returned by
  \code{get\_wall\_quad()} or \code{get\_bndry\_quad()}. Either of the
  arguments \code{nu}, \code{grd\_nu}, \code{D2\_nu} or \code{dets}
  may be \nil; otherwise \code{normals} stores the outer unit normal
  field of the face opposite of the vertex with local number
  \code{wall} and \code{dets} stores the values of the surface
  element. The derivatives of the normal-field are, for instance,
  needed for vector-valued basis functions like face- or edge-bubbles
  (``wall-bubbles''). To this aim the outer normal field is extended
  into the interior of an element by setting it constant on the
  coordinate lines defined by the barycentric coordinates on the
  reference element.

  \kitem{inherit\_parametric(slave), unchain\_parametric(slave)}~\\
  \code{inherit\_parametric()} is used by \code{get\_submesh()},
  \code{unchain\_parametric()} is used by \code{unchain\_submesh()}.
  An application which defines its own \code{PARAMETRIC} structure can
  set both pointers to \nil if the sub-mesh feature is not needed.

  \kitem{data} This \code{void *} pointer is intended for the purpose
  of chaining implementation specific information to the
  \code{PARAMETRIC} structure. In a \code{C++} context the function
  hooks defined in the \code{PARAMETRIC} structure could be virtual
  methods, and implementations would just inherit the
  \code{PARAMETRIC} base-class.
\end{descr}

\medskip

Using the flag \code{FILL\_COORDS} on a mesh traversal (see
\secref{S:traverse}) would fill the \code{EL\_INFO} structures with
coordinate information of the so-called \emph{reference mesh} based on
the original macro triangulation. The reference mesh is what is would
be used without a call to \code{use\_lagrange\_parametric()}. This
reference mesh is normally hidden from the application unless
specifically requested by setting the entry
\code{PARAMETRIC->use\_reference\_mesh} to \true. Furthermore, the
mesh traversal routines ignore the \code{FILL\_COORDS} flag unless
\code{use\_reference\_mesh} is \true. However, special applications
may profit from accessing the reference mesh. On the other hand, most
\ALBERTA routines, e.~g.~ routines to evaluate derivatives of basis
functions, will automatically use the parametric mesh structure when
present.

The function pointers \code{PARAMETRIC->coord\_to\_world},
\code{PARAMETRIC->world\_to\_coord}, \code{PARAMETRIC->det},
\code{PARAMETRIC->grd\_lambda}, \code{PARAMETRIC->wall\_normal} should
be used instead of the standard routines for standard simplicial
triangulations
\begin{itemize}
  \item \code{world\_to\_coord()}
  \item \code{coord\_to\_world()}
  \item \code{el\_det()}
  \item \code{el\_volume()}
  \item \code{el\_grd\_lambda()}
  \item \code{get\_wall\_normal()}
\end{itemize}
described in detail in \secref{S:bary_routines}. The exception are
affine elements on only partially parametric meshes: if
\code{PARAMETRIC->init\_element()} returns \code{false} then the
standard routines may be used instead of the function hooks of the
\code{PARAMETRIC} structure. The same holds when using a
``parametric'' mesh of piece-wise polynomial degree through the
\code{use\_lagrange\_parametric()} call, simply because this
implementation ``fakes'' a partially parametric mesh which is
non-curved on all elements. The use of the standard routines in the
affine-linear context can simplify application programs quite a bit.

If \code{ALBERTA\_DEBUG==1} and
\code{use\_reference\_mesh == false} then using the standard library
routines on parametric simplices will exit with an error message. This
is a safety measure to prevent accidental misuse.



\begin{example}[Use of a parametric mesh]
\label{E:parametric_mesh_traverse}
\idx{parametric meshes!use of a parametric mesh}
\idx{Per-element initializers!example for a parametric mesh}
\idx{Per-element initializers!example for extra fill-flags}
\idx{Per-element initializers!example for mesh-traversal}
\idx{init_element()@{\code{init\_element()}}!example for a parametric mesh}
\idx{init_element()@{\code{init\_element()}}!example for extra fill-flags}
\idx{init_element()@{\code{init\_element()}}!example for mesh-traversal}
%%
This example shows how to write a routine which performs a global
interpolation of a given function onto a finite element space. This is
a much simplified version of the \code{interpol()}-implementation
which can be found in \code{alberta/src/Common/eval.c} (path relative
to the top-level directory of the source distribution of \ALBERTA).
Compare also with \secref{S:I_FES}. The simplifications mostly concern
the missing support for direct sums of finite element spaces, but as
this is a scalar-only example, the restriction does not seem to be too
severe.

The function \code{interpol\_simple()} defined here takes a pointer to
an application defined function \code{REAL (*f)(const REAL\_D arg)},
and a \code{DOF\_REAL\_VEC} and loops over all mesh-elements, calling
the local interpolation routines in turn on all elements. We assume
here that the evaluation of \code{f()} is extremely costly, so we are
careful not to evaluate \code{f()} too often. There are two
helper-function, \code{inter\_fct\_loc()} and
\code{inter\_fct\_loc\_param()}, which are used as arguments to the
actual call to the \code{bfcts->interpol()} hook. Note that the code
uses the non-parametric version if either \code{mesh->parametric} is
\nil, or if \code{mesh->parametric->init\_element()} returns
\code{false}.

The example also shows the use of another type of per-element
initializers: basis functions may also carry such a function-hook,
refer to \secref{S:init_element} for a detailed description.
%%
\bv\begin{lstlisting}
static
REAL inter_fct_loc(const EL_INFO *el_info, const QUAD *quad, int iq,
		   void *ud)
{
  FCT_AT_X fct = *(FCT_AT_X *)ud;
  REAL_D world;
  
  coord_to_world(el_info, quad->lambda[iq], world);
  
  return fct(world);
}

static
REAL inter_fct_loc_param(const EL_INFO *el_info, const QUAD *quad, int iq,
			 void *ud)
{
  const PARAMETRIC *parametric = el_info->mesh->parametric;
  FCT_AT_X fct = *(FCT_AT_X *)ud;
  REAL_D world;

  parametric->coord_to_world(el_info, NULL, 1, quad->lambda + iq, &world);

  return fct(world);
}

void interpol_simple(DOF_REAL_VEC *dv, FCT_AT_X f)
{
  /* Some abbreviations ... */
  const FE_SPACE   *fe_space = dv->fe_space;
  const BAS_FCTS   *bfcts    = fe_space->bas_fcts;
  const DOF_ADMIN  *admin    = fe_space->admin;
  MESH             *mesh     = fe_space->mesh;
  const PARAMETRIC *param    = mesh->parametric;
  EL_REAL_VEC *vec_loc;
  bool        is_param;
  FLAGS       fill_flags;
  int         indices[bfcts->n_bas_fcts_max];
  DOF         dofs[bfcts->n_bas_fcts_max];

  /* Initialize each component of vec to HUGE_VAL, misusing it as
   * flag-argument
   */
  FOR_ALL_DOFS(admin, dv->vec[dof] = HUGE_VAL);

  /* Get an element vector to store the result of the interpolation in */
  vec_loc = get_el_real_vec(bfcts);
  
  /* Basis functions may need special fill-flags */
  fill_flags = FILL_COORDS|bfcts->fill_flags;
  TRAVERSE_FIRST(mesh, -1, CALL_LEAF_EL|fill_flags) {
    int i, n_indices;
    REAL val;

    /* Basis-functions may need a per-element initialization */
    if (INIT_ELEMENT(el_info, bfcts) == INIT_EL_TAG_NULL) {
      continue;
    }

    /* Call the per-element initializer of mesh->parametric(), if needed */
    is_param = param != NULL && param->init_element(el_info, param);

    /* Determine which of the local coefficients need to be computed */
    GET_DOF_INDICES(bfcts, el_info->el, admin, dofs);
    for (i = 0, n_indices = 0; i < bfcts->n_bas_fcts; i++) {
      if ((val = dv->vec[dofs[i]]) == HUGE_VAL) {
	indices[n_indices++] = i;
      } else {
	/* "partial" interpolation may need information about the
	 * omitted DOFs nevertheless.
	 */
	vec_loc->vec[i] = val;
      }
    }

    /* Do the actual interpolation. The parametric version could be
     * handled more efficiently if n_indices == n_bas_fcts; in this
     * case we would only need a single call to
     * param->coord_to_world(). Implementing such (and other
     * optimizations) is left to the reader as an exercise).
     */
    if (n_indices == bfcts->n_bas_fcts) {
      /* Interpolation for all DOFs. The parametric version could be
       * handled more efficiently in this case: we would only need a
       * single call to param->coord_to_world(). Implementing such
       * (and other * optimizations) is left to the reader as an
       * exercise).
       */
      INTERPOL(bfcts, vec_loc, el_info, -1, -1, NULL,
	       is_param ? inter_fct_loc_param : inter_fct_loc, &f);
      /* Store the computed values in the global DOF-vector, no need
       * for the indices indirection
       */
      for (i = 0; i < bfcts->n_bas_fcts; i++) {
	dv->vec[dofs[i]] = vec_loc->vec[i];
      }
    } else {
      /* partial interpolation */

      INTERPOL(bfcts, vec_loc, el_info, -1, n_indices, indices,
	       is_param ? inter_fct_loc_param : inter_fct_loc, &f);
      /* Store the computed values in the global DOF-vector. Note that
       * BOTH, the global and the local coefficient vector, are
       * accessed indirectly over the indices array.
       */
      for (i = 0; i < n_indices; i++) {
	dv->vec[dofs[indices[i]]] = vec_loc->vec[indices[i]];
      }
    }   
  } TRAVERSE_NEXT();

  free_el_real_vec(vec_loc); /* Cleanup after ourselves */

  if (INIT_ELEMENT_NEEDED(bfcts)) {
    /* We possibly did not ran over all elementse, initialize any
     * left-over DOFs to 0.0.
     */
    FOR_ALL_DOFS(admin,
		 if (dv->vec[dof] == HUGE_VAL) {
		   dv->vec[dof] = 0.0;
		 });
  }
}
\end{lstlisting}\ev
\end{example}

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