libalberta-dev 3.0.1-1build1 / usr / include / alberta / alberta.h

/*******************************************************************************
 * ALBERTA:  an Adaptive multi Level finite element toolbox using
 *           Bisectioning refinement and Error control by Residual
 *           Techniques for scientific Applications
 *
 * file: alberta.h
 *
 * description: public header file of the ALBERTA package
 *
 *******************************************************************************
 *
 *  authors:   Alfred Schmidt
 *             Zentrum fuer Technomathematik
 *             Fachbereich 3 Mathematik/Informatik
 *             Universitaet Bremen
 *             Bibliothekstr. 2
 *             D-28359 Bremen, Germany
 *
 *             Kunibert G. Siebert
 *             Institut fuer Mathematik
 *             Universitaet Augsburg
 *             Universitaetsstr. 14
 *             D-86159 Augsburg, Germany
 *
 *             Daniel Koester
 *             Institut fuer Mathematik
 *             Universitaet Augsburg
 *             Universitaetsstr. 14
 *             D-86159 Augsburg, Germany
 *
 *             Claus-Justus Heine
 *             Abteilung fuer Angewandte Mathematik
 *             Universitaet Freiburg
 *             Hermann-Herder-Strasse 10
 *             79104 Freiburg, Germany
 *
 *  http://www.alberta-fem.de
 *
 *  (c) by A. Schmidt and K.G. Siebert (1996-2005),
 *                          D. Koester (2002-2005),
 *                         C.-J. Heine (2002-2009).
 *
 ******************************************************************************/

/*******************************************************************************
 *  Header-File for ALBERTA utilities
 ******************************************************************************/
#ifndef _ALBERTA_H_
#define _ALBERTA_H_

#include <alberta/alberta_util.h>

#ifdef __cplusplus
extern "C" {
#elif 0
} /* some editors try to indent because of the brace above */
#endif

/* The version string used for disk-IO. Note the space after the
 * version number, the version number is expected to be of the forn
 * "X.Y " or "X.YZ". The IO version must not necessarily relate to the
 * version of the library.
 */

#define ALBERTA_MAGIC   "ALBERTA: Version "
#define ALBERTA_VERSION ALBERTA_MAGIC"2.3 "

/*******************************************************************************
 *  Definition of the space dimension and of parameters depending on the
 *  space dimension:
 *
 *  DIM_OF_WORLD:   space dimension
 *
 *  The ?D-suffix signals different simplex dimensions (formerly ==DIM).
 ******************************************************************************/

#ifndef DIM_OF_WORLD
# error DIM_OF_WORLD UNDEFINED
#endif

#ifndef ALBERTA_DEBUG
/* #warning is a GNU extension, and there is no way to get at compiler
 * switches like -Werror form here, so do NOT use it.
 */
/* # warning ALBERTA_DEBUG WAS NOT DEFINED! DEFAULTING TO 0. */
# define ALBERTA_DEBUG 0
#endif

/* meshes of at most this dimension are supported. This macro can
 * be use to lay-out static arrays.
 */
#define DIM_LIMIT 3

/* The master dimension limit, meshes of higher dimension are not
 * supported.
 */
#define DIM_MAX MIN(DIM_OF_WORLD, DIM_LIMIT)

/* Various constants for dimension dependent geometrical quantities
 *
 * A note about the terminology: caused by historic reasons the
 * terminology for simplexes is of naive nature: 1d and 2d simplexes
 * do not have faces (unluckily). Instead, 1d simplexes have vertices
 * and 2d simplexes have edges. Of course, 3d simplexes also have
 * vertices and edges, but: those are _very_ different from 1d and 2d
 * vertices and edges.
 *
 * The best solution would be to reserve the name "face" for a
 * co-dimension one face-simplex. Unluckily, there does not seem to be
 * consense about this. Therefore, the naive ALBERTA notion of "face"
 * is maintained, but supplemented by the notion of a "wall", which
 * denotes a co-dimension one (1) face-simplex. Whenever the notation
 * "wall" is used somewhere in ALBERTA then it denotes stuff
 * concerning co-dimension 1 simplexes (i.e. "faces"). Like WALL_QUAD,
 * EL_INFO->wall_bound, wall-transformation (for periodic boundaries)
 * etc.
 *
 * In my opionion this is a FIXME. "face" should be used everywhere to
 * denote a co-dimension 1 face-simplex. In German "Face" means "Seite".
 *
 * cH.
 */
#define N_VERTICES(DIM)   ((DIM)+1)
#define N_EDGES(DIM)      ((DIM)*((DIM)+1)/2)
#define N_WALLS(DIM)      ((DIM)+1) /* number of codim 1 subsimplexes */
#define N_FACES(DIM)      (((DIM) == 3) * N_WALLS(DIM))
#define N_NEIGH(DIM)      (((DIM) != 0) * N_WALLS(DIM))
#define N_LAMBDA(DIM)     N_VERTICES(DIM)
#define DIM_FAC(DIM)      ((DIM) < 2 ? 1 : (DIM) == 2 ? 2 : 6)

#define VERTEX_OF_EDGE(DIM, EDGE)		\
  ((DIM) == 1					\
   ? vertex_of_edge_1d[(EDGE)]			\
   : ((DIM == 2)				\
      ? vertex_of_edge_2d[(EDGE)]		\
      : vertex_of_edge_3d[(EDGE)]))

#define VERTEX_OF_WALL(DIM, WALL)		\
  ((DIM) == 1					\
   ? vertex_of_wall_1d[(WALL)]			\
   : ((DIM == 2)				\
      ? vertex_of_wall_2d[(WALL)]		\
      : vertex_of_wall_3d[(WALL)]))

#define N_VERTICES_0D     N_VERTICES(0)
#define N_EDGES_0D        N_EDGES(0)
#define N_FACES_0D        N_FACES(0)
#define N_NEIGH_0D        N_NEIGH(0)
#define N_WALLS_0D        N_WALLS(0)
#define N_LAMBDA_0D       N_LAMBDA(0)
#define DIM_FAC_0D        DIM_FAC(0)
#define VERTEX_OF_EDGE_0D NULL
#define VERTEX_OF_WALL_0D NULL

#define N_VERTICES_1D     N_VERTICES(1)
#define N_EDGES_1D        N_EDGES(1)
#define N_FACES_1D        N_FACES(1)
#define N_NEIGH_1D        N_NEIGH(1)
#define N_WALLS_1D        N_WALLS(1)
#define N_LAMBDA_1D       N_LAMBDA(1)
#define DIM_FAC_1D        DIM_FAC(1)
#define VERTEX_OF_EDGE_1D(E) VERTEX_OF_EDGE(1, E)
#define VERTEX_OF_WALL_1D(W) VERTEX_OF_WALL(1, W)

#define N_VERTICES_2D     N_VERTICES(2)
#define N_EDGES_2D        N_EDGES(2)
#define N_FACES_2D        N_FACES(2)
#define N_NEIGH_2D        N_NEIGH(2)
#define N_WALLS_2D        N_WALLS(2)
#define N_LAMBDA_2D       N_LAMBDA(2)
#define DIM_FAC_2D        DIM_FAC(2)
#define VERTEX_OF_EDGE_2D(E) VERTEX_OF_EDGE(2, E)
#define VERTEX_OF_WALL_2D(W) VERTEX_OF_WALL(2, W)

#define N_VERTICES_3D     N_VERTICES(3)
#define N_EDGES_3D        N_EDGES(3)
#define N_FACES_3D        N_FACES(3)
#define N_NEIGH_3D        N_NEIGH(3)
#define N_WALLS_3D        N_WALLS(3)
#define N_LAMBDA_3D       N_LAMBDA(3)
#define DIM_FAC_3D        DIM_FAC(3)
#define VERTEX_OF_EDGE_3D(E) VETEX_OF_EDGE(3, E)
#define VERTEX_OF_WALL_3D(W) VERTEX_OF_WALL(3, W)

/* The maximal number for a given DIM_OF_WORLD.
 */
#define N_LAMBDA_MAX    N_VERTICES(DIM_MAX)
#define N_VERTICES_MAX  N_VERTICES(DIM_MAX)
#define N_EDGES_MAX     N_EDGES(DIM_MAX)
#define N_FACES_MAX     N_FACES(DIM_MAX)
#define N_NEIGH_MAX     N_NEIGH(DIM_MAX)
#define N_WALLS_MAX     N_WALLS(DIM_MAX)
#define DIM_FAC_MAX     DIM_FAC(DIM_MAX)

/* The maximal number which can possibly occur within ALBERTA. These
 * macros can be used to layout static arrays. Of course, DIM_LIMIT is
 * hard-wired to 3, and that will not change.
 */
#define N_LAMBDA_LIMIT    N_VERTICES(DIM_LIMIT)
#define N_VERTICES_LIMIT  N_VERTICES(DIM_LIMIT)
#define N_EDGES_LIMIT     N_EDGES(DIM_LIMIT)
#define N_FACES_LIMIT     N_FACES(DIM_LIMIT)
#define N_NEIGH_LIMIT     N_NEIGH(DIM_LIMIT)
#define N_WALLS_LIMIT     N_WALLS(DIM_LIMIT)
#define DIM_FAC_LIMIT     DIM_FAC(DIM_LIMIT)

/* As N_LAMBDA depends on DIM_MAX, we provide some convenience
 * macros for initializing arrays etc.
 */
#if DIM_MAX == 0
# define INIT_BARY_0D(a)           { 1.0 }
# define INIT_BARY_1D(a, b)        { 1.0 }
# define INIT_BARY_2D(a, b, c)     { 1.0 }
# define INIT_BARY_3D(a, b, c, d)  { 1.0 }
# define INIT_BARY_MAX(a, b, c, d) INIT_BARY_0D(a)
#elif DIM_MAX == 1
# define INIT_BARY_0D(a)           { (a), 0.0 }
# define INIT_BARY_1D(a, b)        { (a), (b) }
# define INIT_BARY_2D(a, b, c)     { (a), (b) }
# define INIT_BARY_3D(a, b, c, d)  { (a), (b) }
# define INIT_BARY_MAX(a, b, c, d) INIT_BARY_1D(a, b)
#elif DIM_MAX == 2
# define INIT_BARY_0D(a)           { (a), 0.0, 0.0 }
# define INIT_BARY_1D(a, b)        { (a), (b), 0.0 }
# define INIT_BARY_2D(a, b, c)     { (a), (b), (c) }
# define INIT_BARY_3D(a, b, c, d)  { (a), (b), (c) }
# define INIT_BARY_MAX(a, b, c, d) INIT_BARY_2D(a, b, c)
#elif DIM_MAX == 3
# define INIT_BARY_0D(a)           { (a), 0.0, 0.0, 0.0 }
# define INIT_BARY_1D(a, b)        { (a), (b), 0.0, 0.0 }
# define INIT_BARY_2D(a, b, c)     { (a), (b), (c), 0.0 }
# define INIT_BARY_3D(a, b, c, d)  { (a), (b), (c), (d) }
# define INIT_BARY_MAX(a, b, c, d) INIT_BARY_3D(a, b, c, d)
#else
# error Unsupported DIM_MAX
#endif

/* Various matrix and vector types. Please use them as apropriate.
 *
 * A famous fortune reading (grin):
 *
 * The primary purpose of the DATA statement is to give names to
 * constants; instead of referring to pi as 3.141592653589793 at every
 * appearance, the variable PI can be given that value with a DATA
 * statement and used instead of the longer form of the constant.
 * This also simplifies modifying the program, should the value of pi
 * change.
 *
 *              -- FORTRAN manual for Xerox Computers
 *
 */
typedef REAL            REAL_B[N_LAMBDA_MAX];
typedef REAL_B          REAL_BB[N_LAMBDA_MAX];
typedef REAL            REAL_D[DIM_OF_WORLD];
typedef REAL_D          REAL_DD[DIM_OF_WORLD];
typedef REAL_D          REAL_BD[N_LAMBDA_MAX];
typedef REAL_BD         REAL_BBD[N_LAMBDA_MAX];
typedef REAL_DD         REAL_DDD[DIM_OF_WORLD];
typedef REAL_DD         REAL_BDD[N_LAMBDA_MAX];
typedef REAL_BDD        REAL_BBDD[N_LAMBDA_MAX];
typedef REAL_B          REAL_DB[DIM_OF_WORLD];
typedef REAL_BB         REAL_DBB[DIM_OF_WORLD];
typedef REAL_BB         REAL_BBB[N_LAMBDA_MAX];
typedef REAL_BBB        REAL_BBBB[N_LAMBDA_MAX];
typedef REAL_BBB        REAL_DBBB[DIM_OF_WORLD];
typedef REAL_BBBB       REAL_DBBBB[DIM_OF_WORLD];
typedef REAL_DB         REAL_BDB[N_LAMBDA_MAX];
typedef REAL_DBB        REAL_BDBB[N_LAMBDA_MAX];

/******************************************************************************/

/* The maximum number of quadrature points and local basis functions
 * for each dimension. Useful to define C99 variable size arrays.
 * Note that those fields are intentionlly are NOT labeled "const", it
 * is possible to generate quadrature formulas of arbitrary degree
 * using "get_product_quad(). The INIT_ELEMENT() frame-work provides
 * means to introduce quadratures with per-element initialization,
 * etc. An application may install new quadrature rules using
 * "register_quadrature()" or "add_quadrature()". Similar thing hold
 * for basis functions.
 */
extern int n_quad_points_max[];
extern int n_bas_fcts_max[];

/*******************************************************************************
 *  access to element index via element or element_info structure
 ******************************************************************************/

#if ALBERTA_DEBUG
#define INDEX(el)   ((el) ? (el)->index : -1)
#else
#define INDEX(el)   -1
#endif

/*******************************************************************************
 *  access to leaf data (only for leaf elements)
 ******************************************************************************/

#define IS_LEAF_EL(el) (!(el)->child[0])
#define LEAF_DATA(el)  ((void *)(el)->child[1])

/*******************************************************************************
 *  boundary types
 ******************************************************************************/

#define INTERIOR      0
#define DIRICHLET     1
#define NEUMANN      -1

#define IS_NEUMANN(bound) ((bound) <= NEUMANN)
#define IS_DIRICHLET(bound) ((bound) >= DIRICHLET)
#define IS_INTERIOR(bound) ((bound) == 0)

/*******************************************************************************
 *  node types (indices in n_dof[] vectors, e.g.)
 ******************************************************************************/

/* Be careful: in 1D we have only VERTEX and CENTER nodes (although
 * that violates the usual geometric meaning of VERTEX/EDGE/FACE:
 * looking at the 2d/3d code one really would expect 1d CENTER DOFs to
 * be 1d EDGE DOFs, but this is not the case).
 *
 * So:
 *
 * 1d: VERTEX and CENTER, EL_INFO->wall_bound refers to VERTEX boundary type
 * 2d: VERTEX and CENTER and EDGE, EL_INFO->wall_bound refers to EDGEs
 * 3d: FACE comes into play, EL_INFO->wall_bound refers to FACEs
 */
enum node_types {
  VERTEX = 0,
  CENTER,
  EDGE,
  FACE,
  N_NODE_TYPES
};

/*******************************************************************************
 *  basic types of the grid
 ******************************************************************************/

typedef signed int              DOF;
typedef enum node_types         NODE_TYPES;
#define N_BNDRY_TYPES 256
typedef U_CHAR                  BNDRY_TYPE;
typedef BITS_256                BNDRY_FLAGS;
typedef struct el               EL;
typedef struct macro_el         MACRO_EL;
typedef struct el_info          EL_INFO;
typedef struct el_geom_cache    EL_GEOM_CACHE;
typedef struct rc_list_el       RC_LIST_EL;
typedef struct mesh             MESH;

typedef struct parametric       PARAMETRIC;
typedef struct traverse_stack   TRAVERSE_STACK;

typedef struct adapt_stat       ADAPT_STAT;
typedef struct adapt_instat     ADAPT_INSTAT;

#ifndef DOF_ADMIN_DEF
typedef struct dof_admin          DOF_ADMIN;
typedef struct dof_int_vec        DOF_INT_VEC;
typedef struct dof_dof_vec        DOF_DOF_VEC;
typedef struct dof_uchar_vec      DOF_UCHAR_VEC;
typedef struct dof_schar_vec      DOF_SCHAR_VEC;
typedef struct dof_real_vec       DOF_REAL_VEC;
typedef struct dof_real_d_vec     DOF_REAL_D_VEC;
typedef struct dof_real_dd_vec    DOF_REAL_DD_VEC;
typedef struct dof_ptr_vec        DOF_PTR_VEC;
typedef struct dof_real_vec_d     DOF_REAL_VEC_D;
typedef struct dof_matrix         DOF_MATRIX;
typedef struct matrix_row         MATRIX_ROW;
typedef struct matrix_row_real    MATRIX_ROW_REAL;
typedef struct matrix_row_real_d  MATRIX_ROW_REAL_D;
typedef struct matrix_row_real_dd MATRIX_ROW_REAL_DD;
#endif

typedef struct el_matrix        EL_MATRIX;
typedef struct el_int_vec       EL_INT_VEC;
typedef struct el_dof_vec       EL_DOF_VEC;
typedef struct el_uchar_vec     EL_UCHAR_VEC;
typedef struct el_schar_vec     EL_SCHAR_VEC;
typedef struct el_bndry_vec     EL_BNDRY_VEC;
typedef struct el_ptr_vec       EL_PTR_VEC;
typedef struct el_real_vec      EL_REAL_VEC;
typedef struct el_real_dd_vec   EL_REAL_DD_VEC;
typedef struct el_real_d_vec    EL_REAL_D_VEC;
typedef struct el_real_vec_d    EL_REAL_VEC_D;

typedef struct bas_fcts         BAS_FCTS;
typedef struct fe_space         FE_SPACE;

typedef struct quadrature       QUAD;
typedef struct quadrature       QUADRATURE;
typedef struct quad_fast        QUAD_FAST;

typedef struct quad_el_cache    QUAD_EL_CACHE;

typedef struct wall_quadrature WALL_QUAD;
typedef struct wall_quad_fast  WALL_QUAD_FAST;

typedef struct macro_data       MACRO_DATA;
typedef struct node_projection  NODE_PROJ;
typedef struct node_projection  NODE_PROJECTION;
typedef struct aff_trafo        AFF_TRAFO;

typedef struct dof_comp_hook    DOF_COMP_HOOK;

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

typedef REAL (*FCT_AT_X)(const REAL_D x);
typedef const REAL *(*GRD_FCT_AT_X)(const REAL_D x, REAL_D result);
typedef const REAL_D *(*D2_FCT_AT_X)(const REAL_D x, REAL_DD result);

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


/*******************************************************************************
 * Macro for calling STRUCT->init_element() routines, per element
 * initializers for (at least) QUAD, BAS_FCTS, QUAD_FAST, WALL_QUAD,
 * WALL_QUAD_FAST
 ******************************************************************************/

/* INIT_ELEMENT(el_info, object)
 *
 * The following convention holds:
 *
 * a) This macro evaluates to INIT_EL_TAG_DFLT when no
 *    element-initializer is present.
 *
 * b) An init_element() method MUST allow a NULL pointer for the
 *    el_info argument. If called with el_info == NULL the
 *    init_element() method must restore its default behaviour. The
 *    "default case" is what the implementation defines as default;
 *    for performance reasons the default case should be the one which
 *    applies to the majority of mesh elements.
 *
 * c) The return value of the init_element() method must be
 *    INIT_EL_TAG_DFLT for the default case.
 *
 * d) The return value of the init_element() method must be
 *    INIT_EL_TAG_NULL for the NULL case, meaning, e.g., the number of
 *    basis functions is 0, or the number of quadrature points is
 *    zero. The application can assume that in the NULL case the
 *    structure does not contain any real data.
 *
 * e) In all other cases the return value is a tag which is used to
 *    efficiently cache values of intermediate computations, e.g. the
 *    values of basis functions at quadrature points. This tag should
 *    be locally unique, meaning that consecutive invocations of
 *    init_element() should return different tags for different
 *    simplexes. This can be used for optimizations: if the tag
 *    returned by an init_element() routine does not change, then the
 *    calling function may assume that the underlying object has not
 *    changed.
 */

enum {
  INIT_EL_TAG_NONE = 0, /* invalid tag */
  INIT_EL_TAG_DFLT = 1, /* default case */
  INIT_EL_TAG_NULL = 2  /* something is 0, e.g. no quad-points, basis
			 * functions are identically zero and so on.
			 */
};

typedef unsigned int INIT_EL_TAG;
typedef INIT_EL_TAG (*INIT_ELEMENT_FCT)(const EL_INFO *el_info, void *thisptr);

/* Tag context. */
typedef struct init_el_tag_ctx {
  INIT_EL_TAG  tag;
  unsigned int cnt;
} INIT_EL_TAG_CTX;

#define INIT_EL_TAG_CTX_INIT(ctx)		\
  {						\
    (ctx)->tag = INIT_EL_TAG_DFLT;		\
    (ctx)->cnt = 0;				\
  }

/* Generate a new unique tag != NULL & DFLT */
#define INIT_EL_TAG_CTX_UNIQ(ctx)			\
  {							\
    (ctx)->tag = INIT_EL_TAG_NULL + (++((ctx)->cnt));	\
    if ((ctx)->tag == INIT_EL_TAG_NONE) {		\
      (ctx)->cnt = 1;					\
      (ctx)->tag = INIT_EL_TAG_NULL + 1;		\
    }							\
  }
#define INIT_EL_TAG_CTX_NULL(ctx) (ctx)->tag = INIT_EL_TAG_NULL
#define INIT_EL_TAG_CTX_DFLT(ctx) (ctx)->tag = INIT_EL_TAG_DFLT
#define INIT_EL_TAG_CTX_TAG(ctx)  (ctx)->tag

#define INIT_ELEMENT_METHOD(obj)         (obj)->init_element
#define INIT_ELEMENT_FLAGS(obj)          (obj)->fill_flags
#define INIT_ELEMENT_DEFUN(obj, init_el, flags)			\
  {								\
    INIT_ELEMENT_METHOD(obj) = init_el;				\
    INIT_ELEMENT_FLAGS(obj) = (flags);				\
    INIT_EL_TAG_CTX_INIT(&(obj)->tag_ctx);			\
  }

#define INIT_OBJECT(object) (void)INIT_ELEMENT(NULL, object)
#define INIT_ELEMENT_DECL			\
  INIT_ELEMENT_FCT init_element;		\
  FLAGS fill_flags;				\
  INIT_EL_TAG_CTX tag_ctx

#define INIT_ELEMENT_INITIALIZER(init_el, flags)	\
  (init_el), (flags), { INIT_EL_TAG_DFLT, 0 }

#define INIT_ELEMENT(el_info, object)					\
  (INIT_ELEMENT_NEEDED(object)						\
   ? INIT_ELEMENT_METHOD(object)(el_info, (void *)object) : INIT_EL_TAG_DFLT)

#define INIT_ELEMENT_NEEDED(object) (INIT_ELEMENT_METHOD(object) != NULL)
#define INIT_ELEMENT_SETUP(el_info, object, tagvar, null_action, chg_action) \
  {									\
    INIT_EL_TAG AI_init_el_tag;						\
									\
    AI_init_el_tag = INIT_ELEMENT(el_info, object);			\
    if (AI_init_el_tag == (INIT_EL_TAG)INIT_EL_TAG_NULL) {		\
      (tagvar) = INIT_EL_TAG_NULL;					\
      null_action;							\
    } else if (AI_init_el_tag != (tagvar)) {				\
      (tagvar) = AI_init_el_tag;					\
      chg_action;							\
    }									\
  }

/* Initialization of single items of DOF-chain */
#define INIT_OBJECT_SINGLE(object) (void)INIT_ELEMENT_SINGLE(NULL, object)
#define INIT_ELEMENT_SINGLE(el_info, object)				\
  (INIT_ELEMENT_NEEDED((object)->unchained)				\
   ? INIT_ELEMENT_METHOD((object)->unchained)(el_info, (void *)object)	\
   : INIT_EL_TAG_DFLT)

/*******************************************************************************
 *  node projection descriptor:
 *
 *  a function pointer which calculates the projected location of a
 *  new vertex resulting from refinement.
 *
 ******************************************************************************/

struct node_projection
{
  void (*func)(REAL_D old_coord, const EL_INFO *eli, const REAL_B lambda);
};

/*******************************************************************************
 *  The geometric connectivity of periodic meshes is described by wall-
 *  transformations, affine isometries which map the current mesh to its
 *  periodic neighbour accross the wall of an element.
 *  The transformation operates as usual: Mx + t
 ******************************************************************************/

struct aff_trafo
{
  REAL_DD M;
  REAL_D  t;
};

/*******************************************************************************
 *  one single element (triangle) of the grid:
 *******************************************************************************
 *
 *  position of the nodes in 1d:
 *
 *  0 _____ 1   or  0 _____ 1
 *                      2
 *
 *                          child[0]     child[1]
 *  refinement:  0 _____ 1    0 ___ 1  0 ___ 1
 *                               2        2
 *
 *******************************************************************************
 *
 *  position of the nodes in 2d
 *         2                 2                   2                 2
 *        /\       or       /\        or        /\       or       /\
 *       /  \             4/  \ 3              /  \             4/  \ 3
 *      /    \            /    \              /  3 \            /  6 \
 *    0/______\1        0/______\1          0/______\1        0/______\1
 *                           5                                     5
 *
 *  refinement:           2          child[0]  0    1   child[1]
 *                       /\                   /|    |\
 *                     4/  \ 3  -->         5/ |4  3| \ 5
 *                     /  6 \               /6 |    |6 \
 *                   0/______\1           1/___|    |___\0
 *                        5                  3  2  2  4
 *
 *******************************************************************************
 *
 *  3d refinement: vertex numbering after (Baensch +) Kossaczky
 *
 *  edges:
 *  E0:  between V0, V1
 *  E1:  between V0, V2
 *  E2:  between V0, V3
 *  E3:  between V1, V2
 *  E4:  between V1, V3
 *  E5:  between V2, V3
 *
 *  Always edge 0 (between vertices 0 and 1) is bisected.
 *
 *                                        V1
 *                                      -+
 *                                 ****- ||
 *                     E0    ****--      | |
 *                     ****--           |   | E3
 *               ****--                 |    |
 *         ****--                      |      |
 *     V0 +. . . . . . . . . . . . . . | . . . |
 *         ---               (E1)     |         +  V2
 *            ---                     |         /
 *               ---                 |E4      /
 *                  ---              |      /
 *                 E2  ---          |     / E5
 *                        ---       |   /
 *                           ---   |  /
 *                              ---|/
 *                                 +
 *                                 V3
 *
 *******************************************************************************
 *  child:        pointers to the two children of the element
 *                if (child[0]==NULL) element is a leaf of the
 *                        tree
 *  dof:          vector of pointers to dof vectors :-)
 *  new_coord:    in case of curved boundary, coords of ref.edge midpoint
 *  index:        global element index (only for test purposes)
 *  mark:         element is a leaf:
 *                   mark == 0         do not refine/coarsen
 *                   mark > 0          refine (mark times)
 *                   mark < 0          may be coarsened (mark times)
 ******************************************************************************/

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

#if ALBERTA_DEBUG
  int        index;
#endif
};

/*******************************************************************************
 *   child_vertex_3d[el_type][child][i] =
 *       parent's local vertex index of new vertex i
 *       4 stands for the newly generated vertex
 *******************************************************************************
 *   child_edge_3d[el_type][child][i] =
 *       parent's local edge index of new edge i
 *       new edge 2 is half of old edge 0,
 *       new edges 4,5 are really new edges, and value is different:
 *         child_edge_3d[][][4,5] = index of same edge in other child
 *******************************************************************************
 *  vertex_of_edge_?d[edge][i], i = 1,2 are the two vertices of edge
 ******************************************************************************/
/* The numbering of the vertices of all sub-simplices is given by the
 * following fields.
 *
 * Note that we favour a lexicographical ordering of wall (face)
 * indices in favour of a cyclic ordering in 3d, while we use a cyclic
 * ordering in 2d.
 *
 * Note that all arrays contain "excess" element; this way it is
 * possible to avoid modulo calculus just to keep the array index in
 * range, i.e. it is legal to index like follows:
 *
 * w_3d = 3;
 * v_2d = 2;
 * for (w = 0; w < N_WALLS_3D; w++) {
 *   for (v = 0; v < N_VERTICES_2D; v++) {
 *     DO SOMETHING WITH vertex_of_wall_3d[w_3d + w][v_2d+v]
 *   }
 * }
 */
static const int vertex_of_edge_1d[N_EDGES_1D][2*N_VERTICES_1D-1] = {
  {0, 1, 0}
};
static const int vertex_of_wall_1d[2*N_WALLS_1D-1][N_VERTICES_0D] = {
  { 1 }, { 0 }, { 1 }
};
static const int vertex_of_edge_2d[2*N_EDGES_2D-1][2*N_VERTICES_1D-1] = {
  {1, 2, 1}, {2, 0, 2}, {0, 1, 0}, {1, 2, 1}, {2, 0, 2}
};
#define vertex_of_wall_2d vertex_of_edge_2d
static const int vertex_of_edge_3d[N_EDGES_3D][2*N_VERTICES_1D-1] = {
  {0, 1, 0}, {0, 2, 0}, {0, 3, 0}, {1, 2, 1}, {1, 3, 1}, {2, 3, 2}
};
static const int vertex_of_wall_3d[2*N_WALLS_3D-1][2*N_VERTICES_2D-1] = {
  {1, 2, 3, 1, 2}, {0, 2, 3, 0, 2}, {0, 1, 3, 0, 1}, {0, 1, 2, 0, 1},
  {1, 2, 3, 1, 2}, {0, 2, 3, 0, 2}, {0, 1, 3, 0, 1},
};

/*******************************************************************************
 *  edge_of_vertices_3d[i][j]: gives the local index of edge with vertices i, j
 ******************************************************************************/
static const int edge_of_vertices_3d[N_VERTICES_3D][N_VERTICES_3D] = {
  {-1,  0,  1,  2 },
  { 0, -1,  3,  4 },
  { 1,  3, -1,  5 },
  { 2,  4,  5, -1 }
};

/*******************************************************************************
 *  face_of_edge_3d[e][0/1]: gives the local number of the two adjacent faces
 ******************************************************************************/
static const int face_of_edge_3d[N_EDGES_3D][2] = {
  { 2, 3 }, { 1, 3 }, { 1, 2 }, { 0, 3 }, { 0, 2 }, { 0, 1 }
};

/* defined in element_Xd.c: permuted ordering of wall vertices to be
 * able to match quadrature points of neighbouring elements on the
 * separating wall.
 *
 * wall_orientation[_rel]_Xd() returns an index into those arrays.
 */
extern
const int sorted_wall_vertices_1d[N_WALLS_1D][DIM_FAC_1D][2*N_VERTICES_0D-1];
extern
const int sorted_wall_vertices_2d[N_WALLS_2D][DIM_FAC_2D][2*N_VERTICES_1D-1];
extern
const int sorted_wall_vertices_3d[N_WALLS_3D][DIM_FAC_3D][2*N_VERTICES_2D-1];

/*******************************************************************************
 * PARAMETRIC structure, entry in MESH structure
 * description of parametric meshes and elements
 ******************************************************************************/

/*
 * The values to pass as "flags" argument to
 * use_lagrange_parametric(..., flags)
 */
typedef enum param_strategy {
  PARAM_ALL = 0,            /* all elements are to be parametric */
  PARAM_CURVED_CHILDS = 1,  /* bisection along {\lambda_0 = \lambda_1} */
  PARAM_STRAIGHT_CHILDS = 2 /* bisection along straight lines/planes */
} PARAM_STRATEGY;

#define PARAM_STRATEGY_MASK				\
  (PARAM_ALL|PARAM_CURVED_CHILDS|PARAM_STRAIGHT_CHILDS)
#define PARAM_PERIODIC_COORDS 0x04 /* The coordinates themselves are
				    * periodic, normally parametric
				    * coordinates of a periodic mesh
				    * are _NOT_ periodic.
				    */

/* General hook-structure for parametric
 * meshes. use_lagrange_parametric() populates this struture for the
 * standard conforming iso-parametric geometry approximations.
 */
struct parametric
{
  const char *name; /* textual description analogous to BAS_FCTS. */

  /* true: some elements may be non-parametric */
  bool not_all;

  /* true: standard routines coord_to_world, etc. may be used to get
   * data about the reference triangulation. Set to "false" by
   * default.
   */
  bool use_reference_mesh;

  /* init_element(el_info, param) == false : non-parametric element,
   * init_element(el_info, param) == true  : parametric element.
   *
   * NOTE: If PARAMETRIC::init_element(el_info, ...) returns false,
   * then it is supposed to fill el_info->coord with the current
   * element's coordinate information despite the fact the "el_info"
   * is _CONST_. This way the normal per-element functions can be used
   * (e.g. el_det(), el_grd_lambda() etc.) instead of the parametric
   * ones. This simplifies the program flow (and source code) for
   * partially parametric meshes a _LOT_.
   */
  bool (*init_element)(const EL_INFO *el_info, const PARAMETRIC *parametric);

  void (*vertex_coords)(EL_INFO *info);

  void (*coord_to_world)(const EL_INFO *info, const QUAD *quad,
			 int n, const REAL_B lambda[], REAL_D *world);

  /* Be careful with this function, for some world coordinates there
   * are no barycentric coordinates. Works best if world is not too
   * far away from our simplex. */
  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[]);

  /* inherit_parametric is used by get_submesh(), unchain_parametric()
   * is used by unchain_submesh(). Can be left out if the sub-mesh
   * feature is not used.
   */
  void (*inherit_parametric)(MESH *slave);
  void (*unchain_parametric)(MESH *slave);

  void *data; /* private data for specific implementations */
};

/*******************************************************************************
 * EL_GEOM_CACHE structure; geometric information for non-parametric
 * meshes which is not generated during mesh-traversal, but needed in
 * several places. Data like det, Lambda, wall-normals, wall-det.
 * Access to this cache _MUST_ go through fill_el_geom_cache(). See
 * also fill_quad_el_cache(el_info), especially for parametric meshes.
 ******************************************************************************/

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

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

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

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

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

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

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

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

static inline const EL_GEOM_CACHE *
fill_el_geom_cache(const EL_INFO *el_info, FLAGS fill_flag);

/*******************************************************************************
 *  additional information to elements during hierarchy traversal
 *******************************************************************************
 *
 * mesh:           pointer to the mesh structure
 * coord:          world coordinates of the vertices. For curved
 *                 parametric meshes the corresponding information is
 *                 filled by the function hooks in the PARAMETRIC
 *                 structur.
 * macro_el:       pointer to the macro-element we belong to.
 * el:             node in the mesh-tree.
 * parent:         pointer to an EL_INFO structure describing the parent of this
 *                 element in the mesh tree.
 * fill_flag:      copy of the fill-flags used to generate the EL_INFO
 *                 structure.
 * level:          the depth in the mesh-tree, level 0 means root-level (i.e.
 *                 an element of the macro triangulation).
 *
 * macro_wall:     Mapping of the boundary facets of the element to
 *                 the boundary facets of the containing macro
 *                 element. macro_wall[w] == -1 means that the
 *                 boundary facet number "w" is located in the
 *                 interior of the containing macro-element, i.e.
 *                 EL_INFO::macro_el.
 *
 * wall_bound:     Boundary type of the co-dim 1 facets (all
 *                 dimensions). Boundary types range from 0 (interior
 *                 faces) to 255. Boundary types are just markers
 *                 without interpretation, used to group boundary
 *                 facets which share common properties. Needs
 *                 FILL_BOUND.  Boundary types can also be access via
 *                 the EL_INFO::macro_wall[] component which is filled
 *                 all the time. There is also a support function
 *                 wall_bound() for this purpose.
 * vertex_bound:   Boundary type of the vertices. This is a bit-field:
 *                 bit N is set if any of the co-dim 1 facets the
 *                 vertex belongs to has boundary-type N. Needs FILL_BOUND.
 * edge_bound:     Boundary type of the edges (only 3d). Also a bit field,
 *                 obtained in the same manner as vertex_bound. Needs
 *                 FILL_BOUND.
 *
 * active_projection: node projection function for the new vertex
 *                 which would result from a refinement of the current
 *                 element. Needs FILL_PROJECTION.
 *
 * neigh:          pointer to the adjacent elements NULL-pointer for a
 *                 part of the boundary. Needs FILL_NEIGH.
 * opp_coord:      world coordinates of opposite vertices. Needs
 *                 FILL_NEIGH|FILL_COORD.
 * opp_vertex:     local indices of opposite vertices. Needs FILL_NEIGH.
 *
 * el_type:        type of the element, 0, 1, or 2. Only meaningful in 3d.
 * orientation:    orientation of the tetrahedron relative to the macro
 *                 element. This is only set for 3d, otherwise it is fixed
 *                 at 0. For DIM == 3 this gives the orientation
 *                 w.r.t. to the standard co-ordinate frame, for higher
 *                 dimenstion "absolute" orientation makes no sense; so
 *                 "orientation" will be 1 for all macro elements, regard-
 *                 less of their actual relative orientations.
 *
 * el_geom_cache:  A cache to store derived quantities which are not
 *                 computed during mesh-traversal, but are derived
 *                 from the co-ordinate information. Access to the
 *                 cache _must_ go through fill_el_geom_cache(). This
 *                 stuff needs -- of course -- the FILL_COORDS
 *                 fill-flag.
 *
 ******************************************************************************/

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

  S_CHAR          macro_wall[N_WALLS_MAX];

  BNDRY_TYPE      wall_bound[N_WALLS_MAX];
  BNDRY_FLAGS     vertex_bound[N_VERTICES_MAX];
  BNDRY_FLAGS     edge_bound[N_EDGES_MAX];
#if DIM_MAX > 1
  BNDRY_TYPE      face_bound[MAX(1, N_FACES_MAX)];
#endif

  const NODE_PROJ *active_projection;

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

  U_CHAR          el_type;
  S_CHAR          orientation;

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

  EL_GEOM_CACHE   el_geom_cache;
};

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

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

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

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

/* Set a bit in the boundary-type mask. The precise meaning of BIT:
 *
 * BIT == 0: clear the boundary mask (meaning: interior node)
 * BIT >  0: set bit BIT and also bit 0 (meaning: boundary node)
 */
#define BNDRY_FLAGS_SET(flags, bit)		\
  if ((bit) != INTERIOR) {			\
    bitfield_set((flags), INTERIOR);		\
    bitfield_set((flags), (bit));		\
  } else {					\
    BNDRY_FLAGS_INIT(flags);			\
  }
/* return TRUE if SEGMENT and MASK have non-zero overlap */
#define BNDRY_FLAGS_IS_PARTOF(segment, mask)	\
  bitfield_andp((segment), (mask), 1 /* offset */, N_BNDRY_TYPES)
/* FindFirstBoundaryBit, return INTERIOR for interior nodes, otherwise the
 * number of the first bit set in MASK.
 */
#define BNDRY_FLAGS_FFBB(mask) bitfield_ffs(mask, 1 /* offset */, N_BNDRY_TYPES)

/*******************************************************************************
 * RC_LIST_EL structure to describe a refinement/coarsening patch.
 * el_info:        contains information about the patch element. This is not
 *                 a pointer since EL_INFO structures are often overwritten
 *                 during mesh traversal.
 * no:             index of the patch element in the patch.
 * flags:          see the RCLE_... defines below for a description.
 * neigh:          neighbours to the right/left in the orientation of the
 *                 edge, or NULL pointer for a boundary face. (dim == 3 only)
 * opp_vertex:     the opposite vertex of neigh[0/1]. (dim == 3 only)
 ******************************************************************************/

struct rc_list_el
{
  EL_INFO      el_info;
  int          no;
  FLAGS        flags;
  RC_LIST_EL   *neigh[2];
  int          opp_vertex[2];
};

/* Valid settings for RC_LIST_EL->flags. The "PERIODIC" flags can be
 * exploited in refine_interpol/coarse_restrict routines.
 */
#define RCLE_NONE               0x0      /* just nothing special */
#define RCLE_COARSE_EDGE_COMPAT (1 << 0) /* set if the coarsening edge
					  * of the patch element is
					  * the coarsening edge of the
					  * patch. Only for internal
					  * use.
					  */

/*******************************************************************************
 *  flags, which information should be present in the EL_INFO structure
 ******************************************************************************/

#define FILL_NOTHING            0x0000L
#define FILL_COORDS             0x0001L
#define FILL_BOUND              0x0002L
#define FILL_NEIGH              0x0004L
#define FILL_OPP_COORDS         0x0008L
#define FILL_ORIENTATION        0x0010L
#define FILL_PROJECTION         0x0020L
#define FILL_MACRO_WALLS        0x0040L
#define FILL_WALL_MAP           FILL_MACRO_WALLS
#define FILL_NON_PERIODIC       0x0080L
#define FILL_MASTER_INFO        0x0100L
#define FILL_MASTER_NEIGH       0x0200L

#define FILL_ANY					\
  (FILL_COORDS|FILL_BOUND|FILL_NEIGH|FILL_OPP_COORDS|	\
   FILL_ORIENTATION|FILL_PROJECTION|FILL_MACRO_WALLS|	\
   FILL_NON_PERIODIC|FILL_MASTER_INFO|FILL_MASTER_NEIGH)

/*******************************************************************************
 *  flags for mesh traversal
 ******************************************************************************/

#define CALL_EVERY_EL_PREORDER  0x010000L
#define CALL_EVERY_EL_INORDER   0x020000L
#define CALL_EVERY_EL_POSTORDER 0x040000L
#define CALL_LEAF_EL            0x080000L
#define CALL_LEAF_EL_LEVEL      0x100000L
#define CALL_EL_LEVEL           0x200000L
#define CALL_MG_LEVEL           0x400000L    /* used in multigrid methods */

#define TEST_FLAG(flags, el_info)				\
  TEST_EXIT(!((((el_info)->fill_flag)^(flags)) & (flags)),	\
	    "flag "#flags" not set\n")

#if ALBERTA_DEBUG==1
# define DEBUG_TEST_FLAG(flags, el_info)		\
  if((((el_info)->fill_flag)^(flags)) & (flags))	\
    print_error_funcname(funcName, __FILE__, __LINE__),	\
      print_error_msg_exit("flag "#flags" not set\n")
#else
# define DEBUG_TEST_FLAG(flags, el_info) do { funcName = funcName; } while (0)
#endif
/*******************************************************************************
 *  one single element of the macro triangulation:
 *******************************************************************************
 *  el:            pointer to the element data of the macro element
 *  coord:         world coordinates of the nodes on the macro element
 * wall_bound:     Boundary type of the co-dim 1 facets (all
 *                 dimensions). Boundary types range from 0 (interior
 *                 faces) to 127. Boundary types are just a markers
 *                 without interpretation, used to group boundary
 *                 facets which share common properties.
 * vertex_bound:   Boundary type of the vertices. This is a bit-field:
 *                 bit N is set if any of the co-dim 1 facets the
 *                 vertex belongs to has boundary-type N.
 *  edge_bound:    Boundary type of the edges (only 3d). Also a bit field,
 *                 obtained in the same manner as vertex_bound.
 *  projection:    possible node projection functions for all nodes [0]
 *                 or for specific edges or faces (dim > 1), which will
 *                 override entry [0].
 *  index:         unique global index of macro element
 *  neigh:         pointer to the adjacent macro elements
 *                 NULL-pointer for a part of the boundary
 *  opp_vertex:    local index of opposite vertex w.r.t. neighbour numbering
 *  neigh_vertices: local indices of common vertices of the periodic
 *                 neighbour, this component is set only for the virtual
 *                 neighbours on periodic meshes.
 *                 neigh_vertices[wall][loc_idx] is the local vertex number
 *                 on the neighbour the vertex with local number
 *                 (wall + 1 + loc_idx) % N_VERTICES(MESH_DIM) on this
 *                 element is mapped to.
 *  wall_trafo:    only for periodic meshes: the affine transformation which
 *                 maps the mesh across the corresponding wall to the
 *                 neighbour facet. The wall transformation must be
 *                 affine isometries.
 *  np_vertex_bound: boundary type of the vertices when treating a periodic
 *                 mesh as non-periodic
 *  np_edge_bound: like np_vertex_bound
 *  el_type:       type of corresponding element.
 *  orientation:   orientation of corresponding element, used in 3d.
 *
 ******************************************************************************/

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

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

  NODE_PROJ   *projection[N_NEIGH_MAX + 1];

  int         index;

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

  S_CHAR      orientation;

  U_CHAR      el_type;

  /* The chain to the master macro element if we belong to a trace-mesh */
  struct {
    MACRO_EL    *macro_el;
    S_CHAR      opp_vertex;
    BNDRY_FLAGS vertex_bound[MAX(1, N_VERTICES(DIM_MAX-1))];
    BNDRY_FLAGS np_vertex_bound[MAX(1, N_VERTICES(DIM_MAX-1))];
#if DIM_MAX > 1
    BNDRY_FLAGS edge_bound[N_EDGES(MAX(1, DIM_MAX-1))];
    BNDRY_FLAGS np_edge_bound[N_EDGES(MAX(1, DIM_MAX-1))];
#endif
  } master;
};

/* Some support functions to access boundary-facet related data only
 * stored on the macro-element level.
 */

static inline BNDRY_TYPE wall_bound(const EL_INFO *el_info, int wall)
{
  int mwall = el_info->macro_wall[wall];

  if (mwall < 0) {
    return INTERIOR;
  }

  if ((el_info->fill_flag & FILL_NON_PERIODIC)) {
    return el_info->macro_el->wall_bound[mwall];
  }

  if (el_info->macro_el->neigh_vertices[mwall][0] < 0) {
    return el_info->macro_el->wall_bound[mwall];
  } else {
    return INTERIOR;
  }
}

static inline const AFF_TRAFO *wall_trafo(const EL_INFO *el_info, int wall)
{
  int mwall;

  if ((el_info->fill_flag & FILL_NON_PERIODIC)) {
    return NULL;
  }

  mwall = el_info->macro_wall[wall];

  return mwall < 0 ? NULL : el_info->macro_el->wall_trafo[mwall];
}

static inline const NODE_PROJ *wall_proj(const EL_INFO *el_info, int wall)
{
  if (wall < 0) {
    return el_info->macro_el->projection[0];
  } else {
    int mwall = el_info->macro_wall[wall];

    return el_info->macro_el->projection[mwall+1];
  }
}

/*******************************************************************************
 * index based storage of macro triangulations
 ******************************************************************************/

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

  int n_total_vertices;
  int n_macro_elements;

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

  int *mel_vertices;    /* mel_vertices[i*N_VERTICES(dim)+j]:
			 * global index of jth vertex of element i
			 */
  int *neigh;           /* neigh[i*N_NEIGH(dim)+j]:
			 * neighbour j of element i or -1 at boundaries
			 */
  int *opp_vertex;      /* opp_vertex[i*N_NEIGH(dim)+j]: if set (need not
			 * be) the local vertex number w.r.t. the neighbour
			 * of the vertex opposit the separating wall.
			 */
  BNDRY_TYPE *boundary; /* boundary[i*N_NEIGH(dim)+j]:
			 * boundary type of jth co-dim 1 facet of element i
			 *
			 * WARNING: In 1D the local index corresponds
			 * to vertex 1 & vice versa! (Consistent with
			 * macro_data.neigh)
			 */
  U_CHAR *el_type;      /* el_type[i]: type of element i only used in 3d! */
  int (*wall_vtx_trafos)[N_VERTICES(DIM_MAX-1)][2]; /* the wall trafos */
  /* Wall transformations are in terms of mappings between
   * vertices. i-th wall trafo: global vertex number
   * wall_vtx_trafos[i][v][0] maps to wall_vtx_trafos[i][v][1], v loops
   * through the local vertex number of the respective wall.
   */
  int n_wall_vtx_trafos;/* for periodic meshes: number of
			 * combinatorical wall trafos.
			 */
  int *el_wall_vtx_trafos;
  /* el_wall_vtx_trafos[i*N_WALLS(dim)+j] number of the wall
   * transformation of the j-th wall for the i-th element. > 0:
   * #wall_trafo+1. < 0: inverse of -(#wall_trafo+1)
   */
  AFF_TRAFO *wall_trafos; /* The group generators of the space group
			   * defining the periodic structure of the
			   * mesh.
			   */
  int n_wall_trafos;
  int *el_wall_trafos; /* N = el_wall_trafos[i*N_NEIGH(dim)+j]:
			*
			* number of the wall transformation mapping to
			* the neighbouring fundamental domain across
			* the given wall.
			*
			* If negative: inverse of generator -N-1
			* If positive:            generator +N-1
			*/
#if ALBERTA_DEBUG
  char **mel_comment;    /* for debugging */
#endif
};

#ifndef DOF_ADMIN_DEF
# define DOF_ADMIN_DEF

/*******************************************************************************
 *  dof handling
 ******************************************************************************/
/* presumably the largest native integer type */
# define DOF_FREE_UNIT_TYPE long
typedef unsigned DOF_FREE_UNIT_TYPE DOF_FREE_UNIT;
# define DOF_FREE_SIZE     ((int)(8*sizeof(DOF_FREE_UNIT)))
# define DOF_UNIT_ALL_FREE (~0UL)
extern const DOF_FREE_UNIT dof_free_bit[DOF_FREE_SIZE]; /* in dof_admin.c */
# define DOF_UNUSED (-1) /* el->dof[][] == DOF_UNUSED, mark unused DOFs */

# define FOR_ALL_DOFS(admin, todo)					\
  if ((admin)->hole_count == 0) {					\
    int dof;								\
									\
    for (dof = 0; dof < (admin)->used_count; dof++) {			\
      todo;								\
    }									\
  } else {								\
    DOF_FREE_UNIT _dfu, *_dof_free = (admin)->dof_free;			\
    int _i, _ibit, dof=0;						\
    int _n= ((admin)->size_used + DOF_FREE_SIZE-1) / DOF_FREE_SIZE;	\
									\
    for (_i = 0; _i < _n; _i++) {					\
      if ((_dfu = _dof_free[_i])) {					\
	if (_dfu == DOF_UNIT_ALL_FREE) {				\
	  dof += DOF_FREE_SIZE;						\
	} else {							\
	  for (_ibit = 0;						\
	       _ibit < DOF_FREE_SIZE;					\
	       _ibit++, dof++, _dfu >>= 1) {				\
	    if ((_dfu & 1) == 0) {					\
	      todo;							\
	    }								\
	  }								\
	}								\
      } else {								\
	for (_ibit = 0; _ibit < DOF_FREE_SIZE; _ibit++, dof++) {	\
	  todo;								\
	}								\
      }									\
    }									\
  }

# define FOR_ALL_FREE_DOFS(admin, todo)					\
  if ((admin)->hole_count == 0) {					\
    int dof;								\
    for (dof = (admin)->used_count; dof < (admin)->size; dof++) {	\
      todo;								\
    }									\
  } else {								\
    DOF_FREE_UNIT _dfu, *_dof_free = (admin)->dof_free;			\
    int _i, _ibit, dof=0;						\
    int _n= ((admin)->size + DOF_FREE_SIZE-1) / DOF_FREE_SIZE;		\
									\
    for (_i = 0; _i < _n; _i++) {					\
      if ((_dfu = _dof_free[_i])) {					\
	if (_dfu == DOF_UNIT_ALL_FREE) {				\
	  for (_ibit = 0 ; _ibit < DOF_FREE_SIZE; _ibit++, dof++) {	\
	    todo;							\
	  }								\
	} else {							\
	  for (_ibit = 0;						\
	       _ibit < DOF_FREE_SIZE;					\
	       _ibit++, dof++, _dfu >>= 1) {				\
	    if ((_dfu & 1) != 0) {					\
	      todo;							\
	    }								\
	  }								\
	}								\
      } else {								\
	dof += DOF_FREE_SIZE;						\
      }									\
    }									\
  }

/* Stop if dof  >= size_used */
# define FOR_ALL_USED_FREE_DOFS(admin, todo)		\
  FOR_ALL_FREE_DOFS(admin,				\
		    if (dof >= admin->size_used) {	\
		      break;				\
		    }					\
		    todo)

/* Possible values for DOF_ADMIN->flags */
# define ADM_FLAGS_DFLT           0        /* nothing special */
# define ADM_PRESERVE_COARSE_DOFS (1 << 0) /* preserve non-leaf DOFs */
# define ADM_PERIODIC             (1 << 1) /* periodic ADMIN on a
					    * periodic mesh
					    */
#define ADM_FLAGS_MASK (ADM_PRESERVE_COARSE_DOFS | ADM_PERIODIC)

struct dof_admin
{
  MESH         *mesh;
  const char   *name;

  DOF_FREE_UNIT *dof_free;    /* flag bit vector */
  unsigned int  dof_free_size;/* flag bit vector size */
  unsigned int  first_hole;   /* index of first non-zero dof_free entry */

  FLAGS         flags;

  DOF  size;                 /* allocated size of dof_list vector */
  DOF  used_count;           /* number of used dof indices */
  DOF  hole_count;           /* number of FREED dof indices (NOT size-used)*/
  DOF  size_used;            /* > max. index of a used entry */

  int  n_dof[N_NODE_TYPES];  /* dofs from THIS dof_admin */
  int  n0_dof[N_NODE_TYPES]; /* start of THIS admin's DOFs in the mesh. */
  /****************************************************************************/
  DOF_INT_VEC     *dof_int_vec;           /* linked list of int vectors */
  DOF_DOF_VEC     *dof_dof_vec;           /* linked list of dof vectors */
  DOF_DOF_VEC     *int_dof_vec;           /* linked list of dof vectors */
  DOF_UCHAR_VEC   *dof_uchar_vec;         /* linked list of u_char vectors */
  DOF_SCHAR_VEC   *dof_schar_vec;         /* linked list of s_char vectors */
  DOF_REAL_VEC    *dof_real_vec;          /* linked list of real vectors   */
  DOF_REAL_D_VEC  *dof_real_d_vec;        /* linked list of real_d vectors */
  DOF_REAL_DD_VEC *dof_real_dd_vec;       /* linked list of real_d vectors */
  DOF_PTR_VEC     *dof_ptr_vec;           /* linked list of void * vectors */
  DOF_MATRIX      *dof_matrix;            /* linked list of matrices */

  DBL_LIST_NODE   compress_hooks;         /* linked list of custom compress
					   * handlers.
					   */

/*******************************************************************************
 * pointer for administration; don't touch!
 ******************************************************************************/

  void            *mem_info;
};

/* DOF_COMP_HOOK is a linked list rooted in
 * DOF_ADMIN->compress_hooks. The user may install arbitrary many
 * custom compress-handlers via add_dof_compress_hook(),
 * del_dof_compress_hook().
 */
struct dof_comp_hook
{
  DBL_LIST_NODE node; /* our link to the compress_hooks list */
  void (*handler)(DOF first, DOF last, const DOF *new_dof, void *app_data);
  void *application_data;
};

/*******************************************************************************
 *  dof vector structures
 *******************************************************************************
 *  next:        pointer to next structure containing vector of same type
 *  fe_space:    pointer to fe_space  structure
 *  refine_interpol: dof interpolation during refinement
 *  coarse_restrict: restriction of linear functionals evaluated on a finer
 *                   grid and stored in dof vector to the coarser grid
 *                   during coarsening
 *                or dof interpolation during coarsening
 *  size:        allocated size of vector
 *  vec[]:    vector entries (entry is used if dof index is used)
 ******************************************************************************/

#define UCHAR_name       uchar
#define uchar_VECNAME    SCHAR
#define SCHAR_name       schar
#define schar_VECNAME    SCHAR
#define INT_name         int
#define int_VECNAME      INT
#define DOF_name         dof
#define dof_VECNAME      INT
#define PTR_name         ptr
#define ptr_VECNAME      PTR
#define REAL_name        real
#define real_VECNAME     REAL
#define REAL_D_name      real_d
#define real_d_VECNAME   REAL_D
#define REAL_DD_name     real_dd
#define real_dd_VECNAME  REAL_DD
#define BNDRY_name       bndry
#define bndry_VECNAME    BNDRY

# define DECL_DOF_VEC(VECNAME, vectype)					\
  struct CPP_CONCAT3(dof_, VECNAME##_name, _vec)			\
  {									\
    DOF_##VECNAME##_VEC  *next;						\
    const FE_SPACE *fe_space;						\
									\
    const char     *name;						\
									\
    DOF            size;						\
    int            reserved;						\
									\
    vectype        *vec;						\
									\
    void  (*refine_interpol)(DOF_##VECNAME##_VEC *, RC_LIST_EL *, int n); \
    void  (*coarse_restrict)(DOF_##VECNAME##_VEC *, RC_LIST_EL *, int n); \
    void           *user_data;						\
									\
    DBL_LIST_NODE  chain;						\
    const DOF_##VECNAME##_VEC *unchained;				\
									\
    EL_##VECNAME##_VEC *vec_loc;					\
									\
    void           *mem_info; /*pointer for administration; don't touch! */ \
  };									\
  typedef vectype VECNAME##_VEC_TYPE

DECL_DOF_VEC(INT,    int);
DECL_DOF_VEC(DOF,    DOF);
DECL_DOF_VEC(UCHAR,  U_CHAR);
DECL_DOF_VEC(SCHAR,  S_CHAR);
DECL_DOF_VEC(PTR,    void *);
DECL_DOF_VEC(REAL,   REAL);
DECL_DOF_VEC(REAL_D, REAL_D);
DECL_DOF_VEC(REAL_DD, REAL_DD);

/* A finite element function/coefficient vector for DIM_OF_WORLD
 * problems. If the basis-functions are vector-valued themselves, then
 * the vector is actually REAL-valued (stride == 1), otherwise REAL_D
 * valued (stride == DIM_OF_WORLD).
 *
 * So: if stride == 1, then this is actually a DOF_REAL_VEC, if stride ==
 * DIM_OF_WORLD, then this is actually a DOF_REAL_D_VEC.
 */
struct dof_real_vec_d
{
  DOF_REAL_VEC_D   *next;
  const FE_SPACE   *fe_space;

  const char       *name;

  DOF              size;
  int              stride;    /* either 1 or DIM_OF_WORLD */

  REAL             *vec;

  void  (*refine_interpol)(DOF_REAL_VEC_D *, RC_LIST_EL *, int n);
  void  (*coarse_restrict)(DOF_REAL_VEC_D *, RC_LIST_EL *, int n);
  void             *user_data;

  DBL_LIST_NODE    chain;
  const DOF_REAL_VEC_D *unchained;

  EL_REAL_VEC_D    *vec_loc;

  void             *mem_info; /*pointer for administration; don't touch! */
};

typedef REAL REAL_VEC_D_TYPE; /* needed?? */

/*******************************************************************************
 *  sparse matrix with one row for each dof,
 *  entries are either REAL or REAL_DD
 *******************************************************************************
 *  next:        pointer to next matrix (linked list in MESH)
 *  matrix_row[]: pointers to row structures (or NULL if row index is unused)
 *  size:         currently allocated size of matrix_row[]
 ******************************************************************************/

/* "flag" values for "type" component */
typedef enum matent_type {
  MATENT_NONE    = -1,
  MATENT_REAL    =  0,
  MATENT_REAL_D  =  1,
  MATENT_REAL_DD =  2
} MATENT_TYPE;

static const size_t matent_size[4] = {
  0, sizeof(REAL), sizeof(REAL_D), sizeof(REAL_DD)
};
# define MATENT_SIZE(type) matent_size[(type)+1]

struct dof_matrix
{
  DOF_MATRIX     *next;
  const FE_SPACE *row_fe_space;
  const FE_SPACE *col_fe_space;
  const char     *name;

  MATRIX_ROW     **matrix_row;  /* lists of matrix entries   */
  DOF            size;          /* size of vector matrix_row */
  MATENT_TYPE    type;          /* type of matrix entries. */
  size_t         n_entries;     /* total number of entries in the
				 * matrix, updated by
				 * add_element_matrix(), set to 0 by
				 * clear_dof_matrix().
				 */
  bool           is_diagonal;
  union {
    DOF_REAL_VEC    *real;
    DOF_REAL_D_VEC  *real_d;
    DOF_REAL_DD_VEC *real_dd;
  } diagonal;             /* The diagonal entries, if is_diagonal == true */
  DOF_INT_VEC *diag_cols; /* The column indices of the diagonal entries   */
  union {
    DOF_REAL_VEC    *real;
    DOF_REAL_D_VEC  *real_d;
    DOF_REAL_DD_VEC *real_dd;
  } inv_diag; /* a cache for the diagonal entries, may be NULL. */

  BNDRY_FLAGS    dirichlet_bndry; /* bndry-type bit-mask for
				   * Dirichlet-boundary conditions built
				   * into the matrix
				   */

  void           (*refine_interpol)(DOF_MATRIX *, RC_LIST_EL *, int n);
  void           (*coarse_restrict)(DOF_MATRIX *, RC_LIST_EL *, int n);

  /* The list pointers for the block-matrix structure to support
   * direct sums of fe-spaces.
   */
  DBL_LIST_NODE  row_chain;
  DBL_LIST_NODE  col_chain;
  const DOF_MATRIX *unchained;

  void           *mem_info;
};

/*******************************************************************************
 *  row structure for sparse matrix, with either REAL or REAL_DD entries.
 *******************************************************************************
 *  next:        pointer to next structure containing entries of same row
 *  col[]:       column indices of entries (if >= 0; else unused)
 *  entry[]:     matrix entries
 ******************************************************************************/

# define ROW_LENGTH 9

/* The actual size of this structure is determined by the type of the
 * matrix entries. The correct length is allocated in
 * get_matrix_row().
 */
# define SIZEOF_MATRIX_ROW(type)					\
  (sizeof(MATRIX_ROW) - sizeof(REAL_DD) + ROW_LENGTH*sizeof(type))

struct matrix_row
{
  MATRIX_ROW  *next;
  MATENT_TYPE type;
  DOF         col[ROW_LENGTH];    /* column indices */
  union {
    REAL    real[1];
    REAL_D  real_d[1];
    REAL_DD real_dd[1];
  } entry;
};

struct matrix_row_real
{
  MATRIX_ROW_REAL *next;
  MATENT_TYPE     type;
  DOF             col[ROW_LENGTH];    /* column indices */
  REAL            entry[ROW_LENGTH];  /* matrix entries */
};

struct matrix_row_real_d
{
  MATRIX_ROW_REAL_D *next;
  MATENT_TYPE       type;
  DOF               col[ROW_LENGTH];    /* column indices */
  REAL_D            entry[ROW_LENGTH];  /* matrix entries */
};

struct matrix_row_real_dd
{
  MATRIX_ROW_REAL_DD *next;
  MATENT_TYPE        type;
  DOF                col[ROW_LENGTH];    /* column indices */
  REAL_DD            entry[ROW_LENGTH];  /* matrix entries */
};

# define ENTRY_USED(col)         ((col) >= 0)
# define ENTRY_NOT_USED(col)     ((col) < 0)
# define UNUSED_ENTRY    -1
# define NO_MORE_ENTRIES -2

# ifndef __CBLAS_H__
typedef enum { NoTranspose,
	       Transpose,
	       ConjugateTranspose } MatrixTranspose;
# endif

/* In C++ we would call this construct an iterator ... */
# define FOR_ALL_MAT_COLS(type, matrow, what)			\
  {								\
    MATRIX_ROW_##type *row;					\
    int col_idx;						\
    DOF col_dof;						\
								\
    for (row = (MATRIX_ROW_##type *)(matrow); row; row = row->next) {	\
      for (col_idx = 0; col_idx < ROW_LENGTH; col_idx++) {	\
	col_dof = row->col[col_idx];				\
	if (ENTRY_USED(col_dof)) {				\
	  what;							\
	} else if (col_dof == NO_MORE_ENTRIES) {		\
	  break;						\
	}							\
      }								\
      if (col_dof == NO_MORE_ENTRIES) {				\
	break;							\
      }								\
    }								\
  }

#endif  /* DOF_ADMIN_DEF */

typedef struct el_vec_head
{
  int           n_components;
  int           n_components_max;
  DBL_LIST_NODE chain;
  int           reserved;
} EL_VEC_HEAD;

#define DECL_DOF_EL_VEC(VECNAME, vectype)			\
  struct CPP_CONCAT3(el_,VECNAME##_name, _vec)			\
  {								\
    int           n_components;					\
    int           n_components_max;				\
    DBL_LIST_NODE chain;					\
    int           reserved;					\
    vectype       vec[1];					\
  };								\
  typedef vectype EL_##VECNAME##_VEC_TYPE

DECL_DOF_EL_VEC(INT,    int);
DECL_DOF_EL_VEC(DOF,    DOF);
DECL_DOF_EL_VEC(UCHAR,  U_CHAR);
DECL_DOF_EL_VEC(SCHAR,  S_CHAR);
DECL_DOF_EL_VEC(BNDRY,  BNDRY_FLAGS);
DECL_DOF_EL_VEC(PTR,    void *);
DECL_DOF_EL_VEC(REAL,   REAL);
DECL_DOF_EL_VEC(REAL_D, REAL_D);
DECL_DOF_EL_VEC(REAL_DD, REAL_DD);

struct el_real_vec_d
{
  int           n_components;
  int           n_components_max;
  DBL_LIST_NODE chain;
  int           stride; /* either 1 or DIM_OF_WORLD */
  REAL          vec[1];
};
typedef REAL EL_REAL_VEC_D_TYPE;

#undef DECL_DOF_EL_VEC

/*******************************************************************************
 *  Here comes the MESH (giving access to the whole triangulation)
 ******************************************************************************/

struct mesh
{
  const char      *name;

  int             dim;

  int             n_vertices;
  int             n_elements;
  int             n_hier_elements;

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

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

  int             n_macro_el;
  MACRO_EL        *macro_els;

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

  PARAMETRIC      *parametric;

  DOF_ADMIN       **dof_admin;
  int             n_dof_admin;

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

  unsigned int    cookie;    /* changed on each refine/coarsen. Use
			      * this to check consistency of meshes
			      * and DOF vectors when reading from
			      * files.
			      */
  int             trace_id;  /* if this is a trace-mesh (aka sub-mesh)
			      * then this is a unique id identifying
			      * it among all the other trace-meshes
			      * chained to the parent mesh.
			      */

  void            *mem_info; /* pointer for administration; don't touch! */
};

/*******************************************************************************
 * data structure for basis function representation
 ******************************************************************************/

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);

typedef void (*REF_INTER_FCT)(DOF_REAL_VEC *, RC_LIST_EL *, int);
typedef void (*REF_INTER_D_FCT)(DOF_REAL_D_VEC *, RC_LIST_EL *, int);
typedef void (*REF_INTER_FCT_D)(DOF_REAL_VEC_D *, RC_LIST_EL *, int);

#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)

#define GET_DOF_INDICES(bfcts, el, admin, dof)		\
  (bfcts)->get_dof_indices(dof, el, admin, bfcts)
#define INTERPOL(bfcts, coeff, el_info, wall,  n, indices, f, ud)	\
  (bfcts)->interpol(coeff, el_info, wall, n, indices, f, ud, bfcts)
#define INTERPOL_D(bfcts, coeff, el_info, wall, n, indices, f, ud)	\
  (bfcts)->interpol_d(coeff, el_info, wall, n, indices, f, ud, bfcts)
#define INTERPOL_DOW(bfcts, coeff, el_info, wall, n, indices, f, ud)	\
  (bfcts)->interpol_dow(coeff, el_info, wall, n, indices, f, ud, bfcts)
#define GET_BOUND(bfcts, el_info, bound)	\
  (bfcts)->get_bound(bound, el_info, bfcts)

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 (e.g. Mini-element) ************/
  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 *);
  const EL_REAL_DD_VEC *(*get_real_dd_vec)(REAL_DD result[],
					   const EL *, const DOF_REAL_DD_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))

/******************************************************************************
 * FE spaces are a triple of DOFs and BAS_FCTs on a MESH
 *
 * Fe-Spaces may be vector-valued, if either the basis functions are
 * vector-valued, or on request. rdim codes the dimension of the
 * range-space, it may be either 1 or DIM_OF_WORLD, so it is rather a
 * flag-value. Scalar fe-space should in principle not be attached to
 * DOF_REAL_D_VECs (but can be).
 *
 * Fe-spaces may be chained to constitute the Cartesian product space
 * of a couple of distinct fe-spaces, e.g. to form the Cartesian
 * product of some "trivially" vector valued fe-space with another
 * space spanned by "really" vector-valued basis functions, e.g. to
 * add edge respectively face bubbles.
 *
 *****************************************************************************/

struct fe_space
{
  const char      *name;
  const DOF_ADMIN *admin;
  const BAS_FCTS  *bas_fcts;
  MESH            *mesh;
  int             rdim;
  int             ref_cnt;
  DBL_LIST_NODE   chain;
  const FE_SPACE  *unchained;
};

/* How to check whether two FE_SPACE objects are essentially the same */
#define FE_SPACE_EQ_P(fe1, fe2)			\
  ((fe1) == (fe2) ||				\
   ((fe1)->admin    == (fe2)->admin    &&	\
    (fe1)->bas_fcts == (fe2)->bas_fcts &&	\
    (fe1)->mesh     == (fe2)->mesh     &&	\
    (fe1)->rdim     == (fe2)->rdim))

static inline bool fe_space_is_eq(const FE_SPACE *fe1, const FE_SPACE *fe2)
{
  return FE_SPACE_EQ_P(fe1, fe2);
}

/*******************************************************************************
 * data structures for numerical integration
 ******************************************************************************/

struct quadrature
{
  const char   *name;
  int          degree;

  int          dim;          /* barycentric coords have (dim+1) components */
  int          codim;        /* the co-dimension */
  int          subsplx;      /* co-dim 1: face number, co-dim 2: edge number */

  int          n_points;
  int          n_points_max; /* max. number in presence of INIT_ELEMENT() */
  const REAL_B *lambda;
  const REAL   *w;

  void         *metadata; /* for internally kept per element caches etc. */

  INIT_ELEMENT_DECL;
};

/*******************************************************************************
 * per-element quadrature cache for co-ordinates, det, Lambda, DLambda
 * uch a cache structure is associated with each quadrature and can be
 * filled by fill_quad_el_cache(el_info, quad, fill_flag). The cache
 * is incremental, repeated calls will fill in data as needed.
 ******************************************************************************/

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

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

static inline const QUAD_EL_CACHE *fill_quad_el_cache(const EL_INFO *el_info,
						      const QUAD *quad,
						      FLAGS fill);

/*******************************************************************************
 * data structure with precomputed values of basis functions at
 * quadrature nodes on the standard element
 ******************************************************************************/

#define INIT_PHI        0x01
#define INIT_GRD_PHI    0x02
#define INIT_D2_PHI     0x04
#define INIT_D3_PHI     0x08
#define INIT_D4_PHI     0x10
#define INIT_TANGENTIAL 0x80 /* derivatives are tangential, for co-dim 1
			      * quadratures.
			      */

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

  FLAGS           init_flag;

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

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

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

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

  INIT_ELEMENT_DECL;

  void *internal;
};

/*******************************************************************************
 *  data structure for adaptive methods
 ******************************************************************************/

typedef enum adaptation_strategy {
  NoStrategy = 0,
  GlobalRefinement = 1,
  GR   = GlobalRefinement,
  MinimumStrategy = 2,
  MS   = MinimumStrategy,
  EqualDistributionStrategy = 3,
  ES   = EqualDistributionStrategy,
  GuaranteedErrorReductionStrategy = 4,
  GERS = GuaranteedErrorReductionStrategy
} ADAPTATION_STRATEGY;

struct adapt_stat
{
  const char  *name;
  REAL        tolerance;
  REAL        p;                         /* power in estimator norm */
  int         max_iteration;
  int         info;

  REAL   (*estimate)(MESH *mesh, ADAPT_STAT *adapt);
  REAL   (*get_el_est)(EL *el);          /* local error estimate */
  REAL   (*get_el_estc)(EL *el);         /* local coarsening error estimate*/
  U_CHAR (*marking)(MESH *mesh, ADAPT_STAT *adapt);

  void   *est_info;                      /* estimator parameters */
  REAL   err_sum, err_max;               /* sum and max of el_est */

  void   (*build_before_refine)(MESH *mesh, U_CHAR flag);
  void   (*build_before_coarsen)(MESH *mesh, U_CHAR flag);
  void   (*build_after_coarsen)(MESH *mesh, U_CHAR flag);
  void   (*solve)(MESH *mesh);

  int    refine_bisections;
  bool   coarsen_allowed;                /* 0 : 1 */
  int    coarse_bisections;
  FLAGS  adaptation_fill_flags;          /* Fill-flags used during adaptation */

  ADAPTATION_STRATEGY strategy;          /* 1=GR, 2=MS, 3=ES, 4=GERS */
  REAL   MS_gamma, MS_gamma_c;           /* maximum strategy */
  REAL   ES_theta, ES_theta_c;           /* equidistribution strategy */
  REAL   GERS_theta_star, GERS_nu, GERS_theta_c;  /* willy's strategy */
};


struct adapt_instat
{
  const char  *name;

  ADAPT_STAT adapt_initial[1];
  ADAPT_STAT adapt_space[1];

  REAL   time;
  REAL   start_time, end_time;
  REAL   timestep;

  void   (*init_timestep)(MESH *mesh, ADAPT_INSTAT *adapt);
  void   (*set_time)(MESH *mesh, ADAPT_INSTAT *adapt);
  void   (*one_timestep)(MESH *mesh, ADAPT_INSTAT *adapt);
  REAL   (*get_time_est)(MESH *mesh, ADAPT_INSTAT *adapt);
  void   (*close_timestep)(MESH *mesh, ADAPT_INSTAT *adapt);

  int    strategy;
  int    max_iteration;

  REAL   tolerance;
  REAL   rel_initial_error;
  REAL   rel_space_error;
  REAL   rel_time_error;
  REAL   time_theta_1;
  REAL   time_theta_2;
  REAL   time_delta_1;
  REAL   time_delta_2;
  int    info;
};

#define MESH_REFINED   1
#define MESH_COARSENED 2

typedef enum norm
{
  NO_NORM = 0, /* uninitialized */
  H1_NORM = 1, /* H1-half-norm */
  L2_NORM = 2, /* L2-norm */
  L2H1_NORM = H1_NORM|L2_NORM /* full H1-norm */
} NORM;

/*******************************************************************************
 *  data structures for matrix and vector update
 ******************************************************************************/

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

typedef const EL_MATRIX *
(*EL_MATRIX_FCT)(const EL_INFO *el_info, void *fill_info);

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

  MATENT_TYPE    krn_blk_type;

  BNDRY_FLAGS    dirichlet_bndry;
  REAL           factor;

  EL_MATRIX_FCT  el_matrix_fct;
  void           *fill_info;

  const EL_MATRIX_FCT *neigh_el_mat_fcts;
  void           *neigh_fill_info;

  FLAGS          fill_flag;
};

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

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

  BNDRY_FLAGS    dirichlet_bndry;
  REAL           factor;

  EL_VEC_FCT     el_vec_fct;
  void           *fill_info;

  FLAGS          fill_flag;
};

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

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

  BNDRY_FLAGS    dirichlet_bndry;
  REAL           factor;

  EL_VEC_D_FCT   el_vec_fct;
  void           *fill_info;

  FLAGS          fill_flag;
};

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

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

  BNDRY_FLAGS    dirichlet_bndry;
  REAL           factor;

  EL_VEC_FCT_D   el_vec_fct;
  void           *fill_info;

  FLAGS          fill_flag;
};

/*******************************************************************************
 *  a matrix of quadratures to use for each block of a block-operator
 *  acting on a direct sum of FE-space.
 ******************************************************************************/
typedef struct quad_tensor
{
  const QUAD *quad;         /* the quadrature rules to use for this block */
  DBL_LIST_NODE row_chain;  /* cyclic list link for one row */
  DBL_LIST_NODE col_chain;  /* cyclic list link for one column */
  DBL_LIST_NODE dep_chain;  /* cyclic list link for the third index */
} QUAD_TENSOR;

typedef struct wall_quad_tensor
{
  const WALL_QUAD *quad;    /* the quadrature rule to use for this block */
  DBL_LIST_NODE row_chain;  /* cyclic list link for one row */
  DBL_LIST_NODE col_chain;  /* cyclic list link for on e column */
  DBL_LIST_NODE dep_chain;  /* cyclic list link for the third index */
} WALL_QUAD_TENSOR;

/*******************************************************************************
 *  data structure about the differential operator for matrix assemblage
 ******************************************************************************/

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

  const QUAD        *quad[3];
  const QUAD_TENSOR *quad_tensor[3];

  bool (*init_element)(const EL_INFO *el_info, const QUAD *quad[3], void *apd);

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

  union {
    const REAL *(*real)(const EL_INFO *el_info,
			const QUAD *quad, int iq, void *apd);
    const REAL_D *(*real_d)(const EL_INFO *el_info,
			    const QUAD *quad, int iq, void *apd);
    const REAL_DD *(*real_dd)(const EL_INFO *el_info,
			      const QUAD *quad, int iq, void *apd);
    const REAL_DDD *(*real_ddd)(const EL_INFO *el_info,
				const QUAD *quad, int iq, void *apd);
  } Lb0;
  bool           Lb0_pw_const;
  union {
    const REAL *(*real)(const EL_INFO *el_info,
			const QUAD *quad, int iq, void *apd);
    const REAL_D *(*real_d)(const EL_INFO *el_info,
			    const QUAD *quad, int iq, void *apd);
    const REAL_DD *(*real_dd)(const EL_INFO *el_info,
			      const QUAD *quad, int iq, void *apd);
    const REAL_DDD *(*real_ddd)(const EL_INFO *el_info,
				const QUAD *quad, int iq, void *apd);
  } Lb1;
  bool           Lb1_pw_const;
  MATENT_TYPE    Lb_type; /* MATENT_REAL, _REAL_D or _REAL_DD */
  bool           Lb0_Lb1_anti_symmetric;
  int            Lb_degree; /* quadrature degree for the Lb0() & Lb1() kernel */
  const EL_REAL_VEC_D *(*advection_field)(const EL_INFO *el_info, void *apd);
  const FE_SPACE *adv_fe_space;

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

  BNDRY_FLAGS    dirichlet_bndry; /* bndry-type bit-mask for
				   * Dirichlet-boundary conditions
				   * built into the matrix
				   */
  FLAGS          fill_flag;
  void           *user_data; /* application data, passed to init_element */
};

typedef struct bndry_operator_info BNDRY_OPERATOR_INFO;

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

  const WALL_QUAD        *quad[3];
  const WALL_QUAD_TENSOR *quad_tensor[3];

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

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

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

  /* The following two entries are not used yet. */
  const EL_REAL_VEC_D *(*advection_field)(const EL_INFO *el_info, void *apd);
  const FE_SPACE *adv_fe_space;

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

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

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

  FLAGS        fill_flag;
  void         *user_data;
};

/*******************************************************************************
 *
 * A structure describing the "closure" of an differential operator by
 * Dirichlet, Neumann or Robin boundary conditions.
 *
 * Note that the quadrature passed to neumann() and robin() is a
 * co-dim 1 quadrature, so quad->subsplx is the number of the boundary
 * face, and el_info->wall_bound[quad->subsplx] contains the boundary
 * classification.
 *
 ******************************************************************************/

typedef struct bndry_cond_info BNDRY_COND_INFO;
struct bndry_cond_info
{
  const FE_SPACE *fe_space;

  REAL (*dirichlet)(const EL_INFO *el_info,
		    const REAL_B lambda,
		    BNDRY_FLAGS bndry,
		    void *apd);
  BNDRY_FLAGS dirichlet_bndry;

  REAL (*neumann)(const EL_INFO *el_info, const QUAD *quad, int iq, void *apd);
  const WALL_QUAD *neumann_quad;
  BNDRY_FLAGS neumann_bndry;

  REAL (*robin)(const EL_INFO *el_info, const QUAD *quad, int iq, void *apd);
  const WALL_QUAD *robin_quad;
  REAL robin_const; /* if robin == NULL, use this constant value */
  BNDRY_FLAGS robin_bndry;

  S_CHAR (*bndry_type)(const BNDRY_FLAGS bndry_bits);

  void *user_data;
};

typedef struct bndry_cond_info_d BNDRY_COND_INFO_D;
struct bndry_cond_info_d
{
  const FE_SPACE *fe_space;

  const REAL *(*dirichlet)(REAL_D result,
			   const EL_INFO *el_info,
			   const REAL_B lambda,
			   void *apd);
  BNDRY_FLAGS dirichlet_bndry;

  const REAL *(*neumann)(REAL_D result,
			 const EL_INFO *el_info,
			 const QUAD *quad, int iq, void *apd);
  const WALL_QUAD *neumann_quad;
  BNDRY_FLAGS neumann_bound;

  const REAL *(*robin)(REAL_D result,
		       const EL_INFO *el_info,
		       const QUAD *quad, int iq, void *apd);
  const WALL_QUAD *robin_quad;
  REAL_D robin_const; /* if robin == NULL, use this constant value */
  BNDRY_FLAGS robin_bound;

  S_CHAR (*bndry_type)(const BNDRY_FLAGS bndry_bits);

  void *user_data;
};

/*******************************************************************************
 *  calculate element stiffness matrices by preevaluated integrals
 *  over the the reference element. In this notation, "PSI" means the
 *  row-space, i.e. the space of test-functions, "PHI" means the
 *  column-space, i.e. the space of ansatz functions.
 *
 *  The special tri-linear caches Q001, Q010 and Q100 are meant to
 *  support assembling tri-linear forms, e.g. for material
 *  derivatives. The caches are present, but currently there is not
 *  further support for tri-linear forms inside ALBERTA. The "1"
 *  denotes which of the three factors is differentiated.
 ******************************************************************************/

typedef struct q11_psi_phi      Q11_PSI_PHI;
typedef struct q01_psi_phi      Q01_PSI_PHI;
typedef struct q10_psi_phi      Q10_PSI_PHI;
typedef struct q00_psi_phi      Q00_PSI_PHI;
typedef struct q001_eta_psi_phi Q001_ETA_PSI_PHI;
typedef struct q010_eta_psi_phi Q010_ETA_PSI_PHI;
typedef struct q100_eta_psi_phi Q100_ETA_PSI_PHI;

/* for i = 0 ... n_psi-1
 *   for j = 0 ... n_phi-1
 *     for m = 0 ... n_entries[i][j]
 *
 *       Then we have: values[m][i][j] is the product of
 *       (\nabla\psi_i)_{k[i][j][m]} with (\nabla\phi_j)_{l[i][j][m]}
 *
 *     end for
 *   end for
 * end for
 *
 * Products yielding zero are left out.
 */
typedef struct q11_psi_phi_cache
{
  int  n_psi;
  int  n_phi;

  const int  *const*n_entries;
  const REAL *const*const*values;
  const int  *const*const*k;
  const int  *const*const*l;
} Q11_PSI_PHI_CACHE;

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

  const Q11_PSI_PHI_CACHE *cache;

  INIT_ELEMENT_DECL;
};

/* for i = 0 ... n_psi-1
 *   for j = 0 ... n_phi-1
 *     for m = 0 ... n_entries[i][j]
 *
 *       Then we have: values[m][i][j] is the product of
 *       \psi_i with (\nabla\phi_j)_{l[i][j][m]}
 *
 *     end for
 *   end for
 * end for
 *
 * Products yielding zero are left out.
 */
typedef struct q01_psi_phi_cache
{
  int  n_psi;
  int  n_phi;

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

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

  const Q01_PSI_PHI_CACHE *cache;

  INIT_ELEMENT_DECL;
};

/* for i = 0 ... n_psi-1
 *   for j = 0 ... n_phi-1
 *     for m = 0 ... n_entries[i][j]
 *
 *       Then we have: values[m][i][j] is the product of
 *       (\nabla\psi_i)_{k[i][j][m]} with \phi_j
 *
 *     end for
 *   end for
 * end for
 *
 * Products yielding zero are left out.
 */
typedef struct q10_psi_phi_cache
{
  int        n_psi;
  int        n_phi;

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

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

  const Q10_PSI_PHI_CACHE *cache;

  INIT_ELEMENT_DECL;
};

typedef struct q00_psi_phi_cache
{
  int        n_psi;
  int        n_phi;

  const REAL *const*values;
} Q00_PSI_PHI_CACHE;

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

  const Q00_PSI_PHI_CACHE *cache;

  INIT_ELEMENT_DECL;
};

/* for i = 0 ... n_psi-1
 *   for j = 0 ... n_phi-1
 *     for k = 0 ... n_eta-1
 *       for m = 0 ... n_entries[i][j][k]
 *
 *         Then we have: values[m][i][j][k] is the product of
 *         \psi_i, \phi_j with (\nabla\eta_k)_{l[i][j][k][m]}
 *       end for
 *     end for
 *   end for
 * end for
 *
 * Products yielding zero are left out. The analogue holds for the 010
 * and 100 structures, where the position of the "1" marks which
 * factor is differentitated.
 */
typedef struct q001_eta_psi_phi_cache
{
  int  n_eta;
  int  n_psi;
  int  n_phi;

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

struct q001_eta_psi_phi
{
  const BAS_FCTS    *eta;
  const BAS_FCTS    *psi;
  const BAS_FCTS    *phi;
  const QUAD        *quad;

  const Q001_ETA_PSI_PHI_CACHE *cache;

  INIT_ELEMENT_DECL;
};

typedef struct q010_eta_psi_phi_cache
{
  int  n_eta;
  int  n_psi;
  int  n_phi;

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

struct q010_eta_psi_phi
{
  const BAS_FCTS *eta;
  const BAS_FCTS *psi;
  const BAS_FCTS *phi;
  const QUAD     *quad;

  const Q010_ETA_PSI_PHI_CACHE *cache;

  INIT_ELEMENT_DECL;
};

typedef struct q100_eta_psi_phi_cache
{
  int  n_eta;
  int  n_psi;
  int  n_phi;

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

struct q100_eta_psi_phi
{
  const BAS_FCTS *eta;
  const BAS_FCTS *psi;
  const BAS_FCTS *phi;
  const QUAD     *quad;

  const Q100_ETA_PSI_PHI_CACHE *cache;

  INIT_ELEMENT_DECL;
};

/******************************************************************************/

/* Data for assembling a theta splitting scheme. */
typedef struct el_sys_info_instat EL_SYS_INFO_INSTAT;
struct el_sys_info_instat
{
  const FE_SPACE *row_fe_space;
  const FE_SPACE *col_fe_space;

  INIT_EL_TAG (*el_update_fct)(const EL_INFO *el_info,
			       REAL tau, REAL theta,
			       EL_SYS_INFO_INSTAT *thisptr);
  EL_MATRIX       *el_matrix;
  EL_REAL_VEC     *el_load;

  const EL_MATRIX   *el_stiff;
  const EL_MATRIX   *el_mass;
  const EL_REAL_VEC *u_h_loc;

  FLAGS           fill_flag;
  BNDRY_FLAGS     dirichlet_bndry; /* bndry-type bit-mask for
				    * Dirichlet-boundary conditions
				    * built into the matrix
				    */
  MATENT_TYPE     krn_blk_type;     /* MATENT_REAL for scalar problems */
};

typedef struct el_sys_info_dow_instat EL_SYS_INFO_DOW_INSTAT;
struct el_sys_info_dow_instat
{
  const FE_SPACE *row_fe_space;
  const FE_SPACE *col_fe_space;

  INIT_EL_TAG (*el_update_fct)(const EL_INFO *el_info,
			       REAL tau, REAL theta,
			       EL_SYS_INFO_DOW_INSTAT *thisptr);
  EL_MATRIX           *el_matrix;
  EL_REAL_VEC_D       *el_load;

  const EL_MATRIX     *el_stiff;
  const EL_MATRIX     *el_mass;
  const EL_REAL_VEC_D *u_h_loc;

  FLAGS           fill_flag;
  BNDRY_FLAGS     dirichlet_bndry; /* bndry-type bit-mask for
				    * Dirichlet-boundary conditions
				    * built into the matrix
				    */
  MATENT_TYPE     krn_blk_type;
};

typedef struct el_sys_info_d_instat EL_SYS_INFO_D_INSTAT;
struct el_sys_info_d_instat
{
  const FE_SPACE *row_fe_space;
  const FE_SPACE *col_fe_space;

  INIT_EL_TAG (*el_update_fct)(const EL_INFO *el_info,
			       REAL tau, REAL theta,
			       EL_SYS_INFO_D_INSTAT *thisptr);
  EL_MATRIX           *el_matrix;
  EL_REAL_D_VEC       *el_load;

  const EL_MATRIX     *el_stiff;
  const EL_MATRIX     *el_mass;
  const EL_REAL_D_VEC *u_h_loc;

  FLAGS           fill_flag;
  BNDRY_FLAGS     dirichlet_bndry; /* bndry-type bit-mask for
				    * Dirichlet-boundary conditions
				    * built into the matrix
				    */
  MATENT_TYPE     krn_blk_type;
};



/*******************************************************************************
 *  preconditioner types.
 ******************************************************************************/
typedef enum {
  PreconEnd  = -1,  /* Terminator for variable argument precon functions */
  PreconRepeat = PreconEnd,
  NoPrecon   = 0,
  DiagPrecon = 1,
  HBPrecon   = 2,
  BPXPrecon  = 3,
  SSORPrecon = 4,   /* omega == 1, n_iter = 1 */
  __SSORPrecon = 5, /* SSOR, but with variable omega and n_iter */
  ILUkPrecon = 6,   /* combinatorical ILU(k) */
  BlkDiagPrecon = 512,
  BlkSSORPrecon = 513,
} OEM_PRECON;

#define N_BLOCK_PRECON_MAX 10

/*******************************************************************************
 * A data structure which can be use to define more complex
 * preconditioners. The purpose of this structure is to avoid defining
 * functions with an endless number of arguments. This
 * "parameter-transport-structure" can be passed to
 * init_precon_from_type(), instead of calling init_oem_precon().
 *
 * The general idea is:
 *
 * type  -- one of the preconditoner types defined above.
 *
 * param -- if the preconditioner defined by "type" needs additional
 *          parameters, then the corresponding section in the "param"
 *          component has to be filled.
 *
 * Examples (using a C99-compliant C-compiler):
 *
 * type == __SSORPrecon:
 *
 * PRECON_TYPE prec = {
 *   __SSORPrecon,
 *   { .__SSOR = { 1.5, 2 } }
 * };
 *
 * type == BlkDiagPrecon, the FE-space is a direct sum of 3
 * components, e.g. the Crouzeix-Raviart-Mansfield Stokes
 * discretisation in 3d (Lagrange-2 + face-bubble + wall-bubbles):
 *
 * PRECON_TYPE = {
 *   BlkDiagPrecon,
 *   { .BlkDiag = {
 *     { __SSORPrecon, { 1.0, 1 } },
 *     { DiagPrecon },
 *     { DIagPrecon }, }
 *   }
 * }
 *
 *
 ******************************************************************************/

/* Helper struct for BlkPrecon parameters */
struct __precon_type {
  OEM_PRECON type;
  union {
    struct {
      REAL omega;
      int n_iter;
    } __SSOR;
    struct {
      int level;
    } ILUk;
  } param;
};

typedef struct precon_type
{
  OEM_PRECON type;
  union {
    struct {
      REAL omega;
      int n_iter;
    } __SSOR;
    struct {
      int level;
    } ILUk;
    struct {
      struct __precon_type precon[N_BLOCK_PRECON_MAX];
    } BlkDiag;
    struct {
      struct __precon_type precon[N_BLOCK_PRECON_MAX];
      REAL omega;
      int  n_iter;
    } BlkSSOR;
  } param;
} PRECON_TYPE;

/*******************************************************************************
 * The precon "class".
 *
 * precon_data:   Opaque data pointer.
 *
 * init_precon(): Has to be called after the operator has changed,
 *                e.g. to update the inverse of the diagonal after the
 *                matrix has changed.
 *
 * exit_precon(): Destroys the preconditioner, include the precon
 *                structure itself.
 *
 * precon():      The preconditioner itself.
 *
 ******************************************************************************/

typedef struct precon PRECON;
struct precon
{
  void    *precon_data;

  bool    (*init_precon)(void *precon_data);
  void    (*precon)(void *precon_data, int n, REAL *vec);
  void    (*exit_precon)(void *precon_data);
};

extern const PRECON *get_diag_precon(const DOF_MATRIX *A,
				     const DOF_SCHAR_VEC *bound);
extern const PRECON *get_HB_precon(const DOF_MATRIX *matrix,
				   const DOF_SCHAR_VEC *bound,
				   int info);
extern const PRECON *get_BPX_precon(const DOF_MATRIX *matrix,
				    const DOF_SCHAR_VEC *bound,
				    int info);
extern const PRECON *get_SSOR_precon(const DOF_MATRIX *A,
				     const DOF_SCHAR_VEC *bound,
				     REAL omega,
				     int n_iter);
extern const PRECON *get_ILUk_precon(const DOF_MATRIX *A,
				     const DOF_SCHAR_VEC *mask,
				     int ilu_level, int info);

/*******************************************************************************
 * abstract multigrid
 ******************************************************************************/

typedef struct multi_grid_info MULTI_GRID_INFO;

struct multi_grid_info
{
  REAL             tolerance;                     /* tol. for resid */
  REAL             exact_tolerance;               /* tol. for exact_solver  */

  int              cycle;                         /* 1=V-cycle, 2=W-cycle   */
  int              n_pre_smooth, n_in_smooth;     /* no of smoothing loops  */
  int              n_post_smooth;                 /* no of smoothing loops  */
  int              mg_levels;                     /* current no. of levels  */
  int              exact_level;                   /* level for exact_solver */
  int              max_iter;                      /* max. no of MG iter's   */
  int              info;

  int              (*init_multi_grid)(MULTI_GRID_INFO *mg_info);
  void             (*pre_smooth)(MULTI_GRID_INFO *mg_info, int level, int n);
  void             (*in_smooth)(MULTI_GRID_INFO *mg_info, int level, int n);
  void             (*post_smooth)(MULTI_GRID_INFO *mg_info, int level, int n);
  void             (*mg_restrict)(MULTI_GRID_INFO *mg_info, int level);
  void             (*mg_prolongate)(MULTI_GRID_INFO *mg_info, int level);
  void             (*exact_solver)(MULTI_GRID_INFO *mg_info, int level);
  REAL             (*mg_resid)(MULTI_GRID_INFO *mg_info, int level);
  void             (*exit_multi_grid)(MULTI_GRID_INFO *mg_info);

  void             *data;                         /* application dep. data */

};
int MG(MULTI_GRID_INFO *mg_info);

/*******************************************************************************
 * concrete multigrid
 ******************************************************************************/

typedef struct mg_s_info MG_S_INFO;
struct mg_s_info
{
  MULTI_GRID_INFO  *mg_info;                      /* abstract MG info */

  const FE_SPACE      *fe_space;
  const DOF_ADMIN     *vertex_admin;
  DOF_MATRIX          *mat;
  const DOF_REAL_VEC  *f;
  DOF_REAL_VEC        *u;
  const DOF_SCHAR_VEC *bound;

  int              smoother,     exact_solver;
  REAL             smooth_omega, exact_omega;

  int              size;                          /* current size of vectors*/
  DOF_MATRIX       **matrix;                      /* one for each level */
  REAL             **f_h;                         /* one for each level */
  REAL             **u_h;                         /* one for each level */
  REAL             **r_h;                         /* one for each level */
  int              *dofs_per_level;               /* count dofs per level */

  int              sort_size;                     /* size of sort vectors   */
  DOF              *sort_dof;                     /* dofs in order of levels*/
  DOF              *(dof_parent[2]);              /* (for linear elements)  */
  U_CHAR           *dof_level;
  S_CHAR           *sort_bound;                   /* sorted bound */

  int              sort_invers_size;              /* size of inv. sort list */
  int              *sort_dof_invers;              /* inverse sort list */
};
/*******************************************************************************
 *  sort_dof[ sorted dof ]          = unsorted dof
 *  sort_dof_invers[ unsorted dof ] = sorted dof
 ******************************************************************************/

/* file MG_s1.c DOF_sort routines *********************************************/
void MG_s_setup_levels(MG_S_INFO *mg_s_info);
void MG_s_setup_mat_b(MG_S_INFO *mg_s_info,
		      DOF_MATRIX *mat, const DOF_SCHAR_VEC *bound);
void MG_s_dof_copy_to_sparse(MG_S_INFO *mg_s_info,
			     const DOF_REAL_VEC *x, REAL *y);
void MG_s_dof_copy_from_sparse(MG_S_INFO *mg_s_info,
			       const REAL *x, DOF_REAL_VEC *y);
void MG_s_reset_mat(MG_S_INFO *mg_s_info);
void MG_s_sort_mat(MG_S_INFO *mg_s_info);
void MG_s_free_mem(MG_S_INFO *mg_s_info);

/* file MG_s2.c: DOF_sort independent routines ********************************/
void MG_s_restrict_mg_matrices(MG_S_INFO *mg_s_info);
void MG_s_restrict(MULTI_GRID_INFO *mg_info, int mg_level);
void MG_s_prolongate(MULTI_GRID_INFO *mg_info, int mg_level);
REAL MG_s_resid(MULTI_GRID_INFO *mg_info, int mg_level);
void MG_s_smoother(MULTI_GRID_INFO *mg_info, int mg_level, int n);
void MG_s_exact_solver(MULTI_GRID_INFO *mg_info, int mg_level);
void MG_s_gemv(MG_S_INFO *mg_s_info, int mg_level, MatrixTranspose transpose,
	       REAL alpha, DOF_MATRIX *a, REAL *x, REAL beta, REAL *y);


/* file MG_s.c: ***************************************************************/
int mg_s(DOF_MATRIX *matrix, DOF_REAL_VEC *u, const DOF_REAL_VEC *f,
	 const DOF_SCHAR_VEC *bound,
	 REAL tol, int max_iter, int info, char *prefix);
MG_S_INFO *mg_s_init(DOF_MATRIX *matrix, const DOF_SCHAR_VEC *bound,
		     int info, char *prefix);
int mg_s_solve(MG_S_INFO *mg_s_info,
	       DOF_REAL_VEC *u, const DOF_REAL_VEC *f, REAL tol, int max_iter);
void mg_s_exit(MG_S_INFO *mg_s_info);

/*******************************************************************************
 * Graphic output Definitions
 ******************************************************************************/

typedef void * GRAPH_WINDOW;
typedef float  GRAPH_RGBCOLOR[3];

/** flags used by graph_mesh(): ****/
#define GRAPH_MESH_BOUNDARY        1
#define GRAPH_MESH_ELEMENT_MARK    2
#define GRAPH_MESH_VERTEX_DOF      4
#define GRAPH_MESH_ELEMENT_INDEX   8

/*******************************************************************************
 *  very useful macro definitons
 ******************************************************************************/

#define GET_MESH(dim, name, macro_data, init_node_proj, init_wall_trafo) \
  check_and_get_mesh((dim), DIM_OF_WORLD, ALBERTA_DEBUG,		\
		     ALBERTA_VERSION, (name), (macro_data),		\
		     (init_node_proj), (init_wall_trafo))

#define GET_DOF_VEC(ptr, dof_vec)					\
{									\
  DEBUG_TEST_EXIT((dof_vec) && (dof_vec)->vec,				\
		  "%s == NULL\n", (dof_vec) ? NAME(dof_vec) : #dof_vec); \
  (ptr) = (dof_vec)->vec;						\
}

/*******************************************************************************
 *  defined in graphXO.c
 ******************************************************************************/
extern const GRAPH_RGBCOLOR rgb_black;
extern const GRAPH_RGBCOLOR rgb_white;
extern const GRAPH_RGBCOLOR rgb_red;
extern const GRAPH_RGBCOLOR rgb_green;
extern const GRAPH_RGBCOLOR rgb_blue;
extern const GRAPH_RGBCOLOR rgb_yellow;
extern const GRAPH_RGBCOLOR rgb_magenta;
extern const GRAPH_RGBCOLOR rgb_cyan;
extern const GRAPH_RGBCOLOR rgb_grey50;

extern const GRAPH_RGBCOLOR rgb_albert;
extern const GRAPH_RGBCOLOR rgb_alberta;


/*******************************************************************************
 *  used in wall_quad
 ******************************************************************************/

/* A collection of quadrature rules for the integration over walls
 * (3d: faces, 2d: edges) of a simplex.
 *
 * Each of the quadrature rules WALL_QUAD::quad[wall] may have its
 * own INIT_ELEMENT method.
 *
 * INIT_ELEMENT(el_info, WALL_QUAD) may or may not be called: it is
 * legal to only call INIT_ELEMENT(el_info, WALL_QUAD::quad[wall)
 * individually.
 *
 * If INIT_ELEMENT(el_info, WALL_QUAD) is called, then it has to
 * initialize all quadrature rules for all walls, so the sub-ordinate
 * initializers need not be called in this case.
 */
struct wall_quadrature
{
  const char *name;
  int  degree;
  int  dim;
  int  n_points_max;
  QUAD quad[N_WALLS_MAX];

  INIT_ELEMENT_DECL;

  void *metadata;
};

/* Convenience structure for WALL_QUAD: its is legal to call
 *
 * get_quad_fast(bas_fcts, WALL_QUAD::quad[wall], ...)
 *
 * individually, however
 *
 * get_wall_quad_fast()
 *
 * does this in a single run.
 *
 * If INIT_ELEMENT(el_info, WALL_QUAD_FAST) is called, then the
 * sub-ordinate initializers
 * INIT_ELEMENT(el_info,WALL_QUAD_FAST::quad_fast[wall]) need not be
 * called.
 */
struct wall_quad_fast
{
  const WALL_QUAD *wall_quad;
  const BAS_FCTS  *bas_fcts;

  FLAGS           init_flag;
  const QUAD_FAST *quad_fast[N_WALLS_MAX];

  INIT_ELEMENT_DECL;
};

/* initialize the meta-data for the given WALL_QUAD, no need to call
 * this if the WALL_QUAD has been aquired by get_wall_quad(), only
 * needed for externally defined extension quadrature rules.
 */
extern void register_wall_quadrature(WALL_QUAD *wall_quad);

/* Return a suitable quadrature for integrating over the given wall
 * (neigh number), but the barycentric co-ordinates of QUAD->lambda
 * are relative to the neighbour element.
 */
extern const QUAD *get_neigh_quad(const EL_INFO *el_info,
				  const WALL_QUAD *wall_quad,
				  int neigh);

/* Return a suitable QUAD_FAST structure for integrating over the
 * given wall, but relative to the neighbour element. If the returned
 * QUAD_FAST object has a per-element initializer, then it must be
 * called with an EL_INFO structure for the neighbour element.
 *
 * It is also legal to just call
 *
 * get_quad_fast(bas_fcts, get_neigh_quad(el_info, wall_quad, neigh), ...)
 *
 * but get_neigh_quad_fast() is slightly more efficient.
 */
const QUAD_FAST *get_neigh_quad_fast(const EL_INFO *el_info,
				     const WALL_QUAD_FAST *wqfast,
				     int neigh);

/*******************************************************************************
 *  functions supplied by ALBERTA
 ******************************************************************************/

/***   file coarsen.c   *******************************************************/
extern U_CHAR coarsen(MESH *mesh, FLAGS fill_flags);
extern U_CHAR global_coarsen(MESH *mesh, int no, FLAGS fill_flags);
extern int  get_max_level(MESH *mesh);

/***   file dof_admin.c   *****************************************************/
extern void free_dof_index(DOF_ADMIN *admin, int dof);
extern DOF  get_dof_index(DOF_ADMIN *admin);
extern void enlarge_dof_lists(DOF_ADMIN *admin, int minsize);

extern const DOF_ADMIN *get_vertex_admin(MESH *mesh, FLAGS flags);

extern void test_dof_matrix(DOF_MATRIX *matrix);
extern void dof_matrix_set_diagonal(DOF_MATRIX *matrix, bool diag);
extern void dof_matrix_try_diagonal(DOF_MATRIX *matrix);
extern void add_element_matrix(DOF_MATRIX *matrix,
			       REAL factor,
			       const EL_MATRIX *el_matrix,
			       MatrixTranspose transpose,
			       const EL_DOF_VEC *row_dof,
			       const EL_DOF_VEC *col_dof,
			       const EL_SCHAR_VEC *bound);
extern void add_element_vec(DOF_REAL_VEC *drv, REAL factor,
			    const EL_REAL_VEC *el_vec,
			    const EL_DOF_VEC *dof,
			    const EL_SCHAR_VEC *bound);
extern void add_element_d_vec(DOF_REAL_D_VEC *drdv, REAL factor,
			      const EL_REAL_D_VEC *el_vec,
			      const EL_DOF_VEC *dof,
			      const EL_SCHAR_VEC *bound);
extern void add_element_vec_dow(DOF_REAL_VEC_D *drdv, REAL factor,
				const EL_REAL_VEC_D *el_vec,
				const EL_DOF_VEC *dof,
				const EL_SCHAR_VEC *bound);

extern void update_matrix(DOF_MATRIX *dof_matrix,
			  const EL_MATRIX_INFO *minfo,
			  MatrixTranspose transpose);
void update_real_vec(DOF_REAL_VEC *drv, const EL_VEC_INFO *vec_info);
void update_real_d_vec(DOF_REAL_D_VEC *drdv, const EL_VEC_D_INFO *vecd_info);
void update_real_vec_dow(DOF_REAL_VEC_D *drdv, const EL_VEC_INFO_D *vecd_info);

extern void dof_compress(MESH *mesh);
extern void add_dof_compress_hook(const DOF_ADMIN *admin, DOF_COMP_HOOK *hook);
extern void del_dof_compress_hook(DOF_COMP_HOOK *hook);
extern void clear_dof_matrix(DOF_MATRIX *matrix);

extern void print_dof_matrix(const DOF_MATRIX *matrix);
extern void print_dof_real_vec(const DOF_REAL_VEC *drv);
extern void print_dof_real_d_vec(const DOF_REAL_D_VEC *drdv);
extern void print_dof_real_dd_vec(const DOF_REAL_DD_VEC *drdv);
extern void print_dof_real_vec_dow(const DOF_REAL_VEC_D *drvd);
extern void print_real_vec_maple(REAL *vector, int size, const char *vec_name);
extern void print_dof_real_vec_maple(const DOF_REAL_VEC *drv,
				     const char *vec_name);
extern void print_dof_real_d_vec_maple(const DOF_REAL_D_VEC *drdv,
				       const char *vec_name);
extern void print_dof_real_vec_dow_maple(const DOF_REAL_VEC_D *drvd,
					 const char *vec_name);
extern void print_dof_matrix_maple(const DOF_MATRIX *matrix,
				   const char *matrix_name);
extern void fprint_real_vec_maple(FILE *fp,
				  REAL *vector, int size, const char *vec_name);
extern void fprint_dof_real_vec_maple(FILE *fp,
				      const DOF_REAL_VEC *drv,
				      const char *vec_name);
extern void fprint_dof_real_d_vec_maple(FILE *fp,
					const DOF_REAL_D_VEC *drdv,
					const char *vec_name);
extern void fprint_dof_real_vec_dow_maple(FILE *fp,
					  const DOF_REAL_VEC_D *drvd,
					  const char *vec_name);
extern void fprint_dof_matrix_maple(FILE *fp,
				    const DOF_MATRIX *matrix,
				    const char *matrix_name);
extern void file_print_real_vec_maple(const char *file_name, const char *mode,
				      REAL *vector, int size,
				      const char *vec_name);
extern void file_print_dof_real_vec_maple(const char *file_name,
					  const char *mode,
					  const DOF_REAL_VEC *drv,
					  const char *vec_name);
extern void file_print_dof_real_d_vec_maple(const char *file_name,
					    const char *mode,
					    const DOF_REAL_D_VEC *drdv,
					    const char *vec_name);
extern void file_print_dof_real_vec_dow_maple(const char *file_name,
					      const char *mode,
					      const DOF_REAL_VEC_D *drvd,
					      const char *vec_name);
extern void file_print_dof_matrix_maple(const char *file_name,
					const char *mode,
					const DOF_MATRIX *matrix,
					const char *matrix_name);
extern void print_dof_ptr_vec(const DOF_PTR_VEC *dpv);
extern void print_dof_int_vec(const DOF_INT_VEC *div);
extern void print_dof_uchar_vec(const DOF_UCHAR_VEC *div);
extern void print_dof_schar_vec(const DOF_SCHAR_VEC *div);
/* BLAS 1 */
extern REAL dof_nrm2(const DOF_REAL_VEC *x);
extern REAL dof_asum(const DOF_REAL_VEC *x);
extern void dof_set(REAL alpha, DOF_REAL_VEC *x);
extern void dof_scal(REAL alpha, DOF_REAL_VEC *x);
extern REAL dof_dot(const DOF_REAL_VEC *x, const DOF_REAL_VEC *y);
extern void dof_copy(const DOF_REAL_VEC *x, DOF_REAL_VEC *y);
extern void dof_axpy(REAL alpha, const DOF_REAL_VEC *x, DOF_REAL_VEC *y);
/*  some non BLAS  */
extern void dof_xpay(REAL alpha, const DOF_REAL_VEC *x, DOF_REAL_VEC *y);
extern REAL dof_min(const DOF_REAL_VEC *x);
extern REAL dof_max(const DOF_REAL_VEC *x);
/* BLAS 2 */
extern void dof_gemv(MatrixTranspose transpose, REAL alpha,
		     const DOF_MATRIX *A, const DOF_SCHAR_VEC *mask,
		     const DOF_REAL_VEC *x,
		     REAL beta, DOF_REAL_VEC *y);
extern void dof_mv(MatrixTranspose transpose,
		   const DOF_MATRIX *A, const DOF_SCHAR_VEC *mask,
		   const DOF_REAL_VEC *x, DOF_REAL_VEC *y);

/* now the same for REAL_D */
extern void dof_axpy_d(REAL alpha, const DOF_REAL_D_VEC *x, DOF_REAL_D_VEC *y);
extern void dof_copy_d(const DOF_REAL_D_VEC *x, DOF_REAL_D_VEC *y);
extern REAL dof_dot_d(const DOF_REAL_D_VEC *x, const DOF_REAL_D_VEC *y);
extern REAL dof_nrm2_d(const DOF_REAL_D_VEC *x);
extern REAL dof_asum_d(const DOF_REAL_D_VEC *x);
extern void dof_scal_d(REAL alpha, DOF_REAL_D_VEC *x);
extern void dof_set_d(REAL alpha, DOF_REAL_D_VEC *x);
extern void dof_xpay_d(REAL alpha, const DOF_REAL_D_VEC *x, DOF_REAL_D_VEC *y);
extern REAL dof_min_d(const DOF_REAL_D_VEC *x);
extern REAL dof_max_d(const DOF_REAL_D_VEC *x);

extern void
dof_axpy_dow(REAL alpha, const DOF_REAL_VEC_D *x, DOF_REAL_VEC_D *y);
extern void dof_copy_dow(const DOF_REAL_VEC_D *x, DOF_REAL_VEC_D *y);
extern REAL dof_dot_dow(const DOF_REAL_VEC_D *x, const DOF_REAL_VEC_D *y);
extern REAL dof_nrm2_dow(const DOF_REAL_VEC_D *x);
extern REAL dof_asum_dow(const DOF_REAL_VEC_D *x);
extern void dof_scal_dow(REAL alpha, DOF_REAL_VEC_D *x);
extern void dof_set_dow(REAL alpha, DOF_REAL_VEC_D *x);
extern void
dof_xpay_dow(REAL alpha, const DOF_REAL_VEC_D *x, DOF_REAL_VEC_D *y);
extern REAL dof_min_dow(const DOF_REAL_VEC_D *x);
extern REAL dof_max_dow(const DOF_REAL_VEC_D *x);
/* BLAS 2 for REAL_D */
extern void dof_gemv_d(MatrixTranspose transpose, REAL alpha,
		       const DOF_MATRIX *A, const DOF_SCHAR_VEC *mask,
		       const DOF_REAL_D_VEC *x,
		       REAL beta, DOF_REAL_D_VEC *y);
extern void dof_gemv_rdr(MatrixTranspose transpose, REAL alpha,
			 const DOF_MATRIX *A, const DOF_SCHAR_VEC *mask,
			 const DOF_REAL_VEC *x,
			 REAL beta, DOF_REAL_D_VEC *y);
extern void dof_gemv_rrd(MatrixTranspose transpose, REAL alpha,
			 const DOF_MATRIX *A, const DOF_SCHAR_VEC *mask,
			 const DOF_REAL_D_VEC *x,
			 REAL beta, DOF_REAL_VEC *y);
extern void dof_mv_d(MatrixTranspose transpose,
		     const DOF_MATRIX *A, const DOF_SCHAR_VEC *mask,
		     const DOF_REAL_D_VEC *x, DOF_REAL_D_VEC *y);
extern void dof_mv_rdr(MatrixTranspose transpose,
		       const DOF_MATRIX *A, const DOF_SCHAR_VEC *mask,
		       const DOF_REAL_VEC *x,
		       DOF_REAL_D_VEC *y);
extern void dof_mv_rrd(MatrixTranspose transpose,
		       const DOF_MATRIX *A, const DOF_SCHAR_VEC *mask,
		       const DOF_REAL_D_VEC *x,
		       DOF_REAL_VEC *y);

extern void dof_gemv_dow(MatrixTranspose transpose, REAL alpha,
			 const DOF_MATRIX *A, const DOF_SCHAR_VEC *mask,
			 const DOF_REAL_VEC_D *x,
			 REAL beta, DOF_REAL_VEC_D *y);
extern void dof_gemv_dow_scl(MatrixTranspose transpose, REAL alpha,
			     const DOF_MATRIX *A, const DOF_SCHAR_VEC *mask,
			     const DOF_REAL_VEC *x,
			     REAL beta, DOF_REAL_VEC_D *y);
extern void dof_gemv_scl_dow(MatrixTranspose transpose, REAL alpha,
			     const DOF_MATRIX *A, const DOF_SCHAR_VEC *mask,
			     const DOF_REAL_VEC_D *x,
			     REAL beta, DOF_REAL_VEC *y);
extern void dof_mv_dow(MatrixTranspose transpose,
		       const DOF_MATRIX *A, const DOF_SCHAR_VEC *mask,
		       const DOF_REAL_VEC_D *x, DOF_REAL_VEC_D *y);
extern void dof_mv_dow_scl(MatrixTranspose transpose,
			   const DOF_MATRIX *A, const DOF_SCHAR_VEC *mask,
			   const DOF_REAL_VEC *x,
			   DOF_REAL_VEC_D *y);
extern void dof_mv_scl_dow(MatrixTranspose transpose,
			   const DOF_MATRIX *A, const DOF_SCHAR_VEC *mask,
			   const DOF_REAL_VEC_D *x,
			   DOF_REAL_VEC *y);

/* copy operation for DOF_MATRIXes */
extern void dof_matrix_copy(DOF_MATRIX *dst, const DOF_MATRIX *src);

/* low-level administration routines, use with caution */
extern void  add_dof_matrix_to_admin(DOF_MATRIX *, DOF_ADMIN *);
extern void  add_dof_int_vec_to_admin(DOF_INT_VEC *, DOF_ADMIN *);
extern void  add_int_dof_vec_to_admin(DOF_DOF_VEC *, DOF_ADMIN *);
extern void  add_dof_dof_vec_to_admin(DOF_DOF_VEC *, DOF_ADMIN *);
extern void  add_dof_uchar_vec_to_admin(DOF_UCHAR_VEC *, DOF_ADMIN *);
extern void  add_dof_schar_vec_to_admin(DOF_SCHAR_VEC *, DOF_ADMIN *);
extern void  add_dof_real_vec_to_admin(DOF_REAL_VEC *, DOF_ADMIN *);
extern void  add_dof_real_d_vec_to_admin(DOF_REAL_D_VEC *, DOF_ADMIN *);
extern void  add_dof_real_dd_vec_to_admin(DOF_REAL_DD_VEC *, DOF_ADMIN *);
extern void  add_dof_ptr_vec_to_admin(DOF_PTR_VEC *, DOF_ADMIN *);
extern void  remove_dof_matrix_from_admin(DOF_MATRIX *);
extern void  remove_dof_int_vec_from_admin(DOF_INT_VEC *);
extern void  remove_dof_dof_vec_from_admin(DOF_DOF_VEC *);
extern void  remove_int_dof_vec_from_admin(DOF_DOF_VEC *);
extern void  remove_dof_uchar_vec_from_admin(DOF_UCHAR_VEC *);
extern void  remove_dof_schar_vec_from_admin(DOF_SCHAR_VEC *);
extern void  remove_dof_real_vec_from_admin(DOF_REAL_VEC *);
extern void  remove_dof_real_vec_from_admin(DOF_REAL_VEC *);
extern void  remove_dof_real_d_vec_from_admin(DOF_REAL_D_VEC *);
extern void  remove_dof_real_dd_vec_from_admin(DOF_REAL_DD_VEC *);
extern void  remove_dof_ptr_vec_from_admin(DOF_PTR_VEC *);

/***   file wall_quad.c   *****************************************************/

extern const WALL_QUAD *wall_quad_from_quad(const QUAD *quad);
extern const WALL_QUAD *get_wall_quad(int dim, int degree);
extern const WALL_QUAD_FAST *
get_wall_quad_fast(const BAS_FCTS *, const WALL_QUAD *, FLAGS init_flag);
extern const QUAD *get_neigh_quad(const EL_INFO *el_info,
				  const WALL_QUAD *wall_quad,
				  int wall);
extern const QUAD_FAST *get_neigh_quad_fast(const EL_INFO *el_info,
					    const WALL_QUAD_FAST *bndry_qfast,
					    int wall);

/***   file macro.c   *********************************************************/
extern void macro_test(MACRO_DATA *data, const char *new_filename);

extern MACRO_DATA *read_macro(const char *name);
extern MACRO_DATA *read_macro_bin(const char *name);
extern MACRO_DATA *read_macro_xdr(const char *name);

extern bool write_macro(MESH *mesh, const char *name);
extern bool write_macro_bin(MESH *mesh, const char *name);
extern bool write_macro_xdr(MESH *mesh, const char *name);

extern bool write_macro_data(MACRO_DATA *data, const char *name);
extern bool write_macro_data_bin(MACRO_DATA *data, const char *name);
extern bool write_macro_data_xdr(MACRO_DATA *data, const char *name);

extern MACRO_DATA *alloc_macro_data(int dim, int nv, int ne);
extern void free_macro_data(MACRO_DATA *data);
extern void compute_neigh_fast(MACRO_DATA *data);
extern void default_boundary(MACRO_DATA *data, U_CHAR type, bool overwrite);

extern MACRO_DATA *mesh2macro_data(MESH *mesh);
extern void
macro_data2mesh(MESH *mesh, const MACRO_DATA *data,
		NODE_PROJECTION *(*n_proj)(MESH *,MACRO_EL *,int),
		AFF_TRAFO *(*init_wall_trafos)(MESH *, MACRO_EL *, int wall));

/***   file memory.c   ********************************************************/
extern MESH *
check_and_get_mesh(int dim, int dow, int neigh,
		   const char *version, const char *name,
		   const MACRO_DATA *macro_data,
		   NODE_PROJ *(*init_node_proj)(MESH *, MACRO_EL *, int),
		   AFF_TRAFO *(*init_wall_trafo)(MESH *, MACRO_EL *, int wall));
extern void free_dof_admin(DOF_ADMIN *admin, MESH *mesh);
extern void free_int_dof_vec(DOF_DOF_VEC *vec);
extern void free_dof_int_vec(DOF_INT_VEC *vec);
extern void free_dof_dof_vec(DOF_DOF_VEC *vec);
extern void free_dof_matrix(DOF_MATRIX *mat);
extern void free_dof_real_vec(DOF_REAL_VEC *vec);
extern void free_dof_real_d_vec(DOF_REAL_D_VEC *vec);
extern void free_dof_real_dd_vec(DOF_REAL_DD_VEC *vec);
extern void free_dof_real_vec_d(DOF_REAL_VEC_D *vec);
extern void free_dof_schar_vec(DOF_SCHAR_VEC *vec);
extern void free_dof_uchar_vec(DOF_UCHAR_VEC *vec);
extern void free_dof_ptr_vec(DOF_PTR_VEC *vec);
extern void free_fe_space(const FE_SPACE *fe_space);
extern void free_real_d(MESH *mesh, REAL *ptr);
extern void free_matrix_row(const FE_SPACE *, MATRIX_ROW *);
extern void free_real_matrix_row(const FE_SPACE *, MATRIX_ROW_REAL *);
extern void free_real_d_matrix_row(const FE_SPACE *, MATRIX_ROW_REAL_D *);
extern void free_real_dd_matrix_row(const FE_SPACE *, MATRIX_ROW_REAL_DD *);
extern void free_element(EL *el, MESH *mesh);
extern void free_rc_list(MESH *mesh, RC_LIST_EL *list); /* only for 3D */
extern void free_mesh(MESH *);
extern void free_dof(DOF *dof, MESH *mesh, int position, FLAGS flags);

extern DOF             *get_dof(MESH *mesh, int position);
extern DOF             *get_periodic_dof(MESH *mesh, int position,
					 const DOF *twin);
extern const FE_SPACE *copy_fe_space(const FE_SPACE *fe_space);
extern const FE_SPACE *clone_fe_space(const FE_SPACE *fe_space, int rdim);
extern const FE_SPACE  *get_fe_space(MESH *mesh,
				     const char *name,
				     const BAS_FCTS *bas_fcts,
				     int rdim,
				     FLAGS adm_flags);
extern const FE_SPACE  *get_dof_space(MESH *mesh, const char *name,
				      const int ndof[N_NODE_TYPES],
				      FLAGS adm_flags);
extern DOF_INT_VEC     *get_dof_int_vec(const char *name, const FE_SPACE *);
extern DOF_DOF_VEC     *get_int_dof_vec(const char *name, const FE_SPACE *);
extern DOF_DOF_VEC     *get_dof_dof_vec(const char *name, const FE_SPACE *);
extern DOF_MATRIX      *get_dof_matrix(const char *name,
				       const FE_SPACE *row_fe_space,
				       const FE_SPACE *col_fe_space);
extern DOF_REAL_VEC    *get_dof_real_vec(const char *name, const FE_SPACE *);
extern DOF_REAL_D_VEC  *get_dof_real_d_vec(const char *name, const FE_SPACE *);
extern DOF_REAL_DD_VEC *get_dof_real_dd_vec(const char *name, const FE_SPACE *);
extern DOF_REAL_VEC_D  *get_dof_real_vec_d(const char *name, const FE_SPACE *);
extern DOF_SCHAR_VEC   *get_dof_schar_vec(const char *name, const FE_SPACE *);
extern DOF_UCHAR_VEC   *get_dof_uchar_vec(const char *name, const FE_SPACE *);
extern DOF_PTR_VEC     *get_dof_ptr_vec(const char *name, const FE_SPACE *);
extern REAL            *get_real_d(MESH *mesh);
extern MATRIX_ROW      *get_matrix_row(const FE_SPACE *, MATENT_TYPE type);
extern EL              *get_element(MESH *mesh);
extern RC_LIST_EL      *get_rc_list(MESH *mesh);           /* only for 3D */
extern size_t init_leaf_data(MESH *mesh, size_t size,
			     void (*refine_leaf_data)(EL *parent,
						      EL *child[2]),
			     void (*coarsen_leaf_data)(EL *parent,
						       EL *child[2]));
extern EL_MATRIX *get_el_matrix(const FE_SPACE *row_fe_space,
				const FE_SPACE *col_fe_space,
 				MATENT_TYPE op_type);
extern void free_el_matrix(EL_MATRIX *el_mat);
extern void print_el_matrix(const EL_MATRIX *el_mat);

extern EL_INT_VEC *get_el_int_vec(const BAS_FCTS *bas_fcts);
extern void free_el_int_vec(EL_INT_VEC *el_vec);
extern EL_DOF_VEC *get_el_dof_vec(const BAS_FCTS *bas_fcts);
extern void free_el_dof_vec(EL_DOF_VEC *el_vec);
extern EL_UCHAR_VEC *get_el_uchar_vec(const BAS_FCTS *bas_fcts);
extern void free_el_uchar_vec(EL_UCHAR_VEC *el_vec);
extern EL_SCHAR_VEC *get_el_schar_vec(const BAS_FCTS *bas_fcts);
extern void free_el_schar_vec(EL_SCHAR_VEC *el_vec);
extern EL_BNDRY_VEC *get_el_bndry_vec(const BAS_FCTS *bas_fcts);
extern void free_el_bndry_vec(EL_BNDRY_VEC *el_vec);
extern EL_PTR_VEC *get_el_ptr_vec(const BAS_FCTS *bas_fcts);
extern void free_el_ptr_vec(EL_PTR_VEC *el_vec);
extern EL_REAL_VEC *get_el_real_vec(const BAS_FCTS *bas_fcts);
extern void free_el_real_vec(EL_REAL_VEC *el_vec);
extern EL_REAL_D_VEC *get_el_real_d_vec(const BAS_FCTS *bas_fcts);
extern void free_el_real_d_vec(EL_REAL_D_VEC *el_vec);
extern EL_REAL_DD_VEC *get_el_real_dd_vec(const BAS_FCTS *bas_fcts);
extern void free_el_real_dd_vec(EL_REAL_DD_VEC *el_vec);
extern EL_REAL_VEC_D *get_el_real_vec_d(const BAS_FCTS *bas_fcts);
extern void free_el_real_vec_d(EL_REAL_VEC_D *el_vec);

/***   file submesh.c  ********************************************************/
extern MESH *get_submesh(MESH *master, const char *name,
			 bool (*binding_method)(MESH *master, MACRO_EL *el,
						int wall, void *data),
			 void *data);
extern MESH *get_bndry_submesh(MESH *master, const char *name);
extern MESH *read_bndry_submesh_xdr(MESH *master, const char *slave_filename);
extern MESH *get_bndry_submesh_by_type(MESH *master, const char *name,
				       BNDRY_TYPE type);
extern MESH *read_bndry_submesh_by_type_xdr(MESH *master,
					    const char *slave_filename,
					    BNDRY_TYPE type);

extern MESH *get_bndry_submesh_by_segment(MESH *master,
					  const char *name,
					  const BNDRY_FLAGS segment);
extern MESH *read_bndry_submesh_by_segment(MESH *master,
					   const char *slave_filename,
					   const BNDRY_FLAGS segment);
extern MESH *read_bndry_submesh_by_segment_xdr(MESH *master,
					       const char *slave_filename,
					       const BNDRY_FLAGS segment);

extern
MESH *lookup_submesh_by_binding(MESH *master, 
				bool (*binding_method)(MESH *master,
						       MACRO_EL *el, int wall,
						       void *data),
				void *data);
extern
MESH *lookup_bndry_submesh_by_type(MESH *master, BNDRY_TYPE type);
extern
MESH *lookup_bndry_submesh_by_segment(MESH *master, const BNDRY_FLAGS segment);
extern
MESH *lookup_bndry_submesh(MESH *master);

extern void unchain_submesh(MESH *slave);

extern void bind_submesh(MESH *master,
			 MESH *slave,
			 bool (*binding_method)(MESH *master, MACRO_EL *el,
						int wall, void *data),
			 void *data);

extern MESH *read_submesh(MESH *master,
			  const char *slave_filename,
			  bool (*binding_method)(MESH *master, MACRO_EL *el,
						 int wall, void *data),
			  void *data);

extern MESH *read_submesh_xdr(MESH *master,
			      const char *slave_filename,
			      bool (*binding_method)(MESH *master, MACRO_EL *el,
						     int wall, void *data),
			      void *data);
extern MESH *lookup_submesh_by_id(MESH *mesh, int id);
extern MESH *lookup_submesh_by_name(MESH *mesh, const char *name);

extern void get_slave_dof_mapping(const FE_SPACE *m_fe_space,
				  DOF_INT_VEC *s_map);
extern MESH *get_master(MESH *slave);
extern const EL_DOF_VEC *get_master_dof_indices(EL_DOF_VEC *result,
						const EL_INFO *s_el_info,
						const FE_SPACE *m_fe_space);
extern const EL_BNDRY_VEC *get_master_bound(EL_BNDRY_VEC *result,
					    const EL_INFO *s_el_info,
					    const BAS_FCTS *m_bas_fcts);
extern void
fill_master_el_info(EL_INFO *mst_el_info,
		    const EL_INFO *slv_el_info, FLAGS fill_flags);
extern const EL *get_slave_el(const EL *el, int wall, MESH *trace_mesh);
extern void fill_slave_el_info(EL_INFO *slv_el_info,
			       const EL_INFO *el_info, int wall,
			       MESH *trace_mesh);

void trace_to_bulk_coords_2d(REAL_B result,
			     const REAL_B lambda,
			     const EL_INFO *el_info);
void trace_to_bulk_coords_1d(REAL_B result,
			     const REAL_B lambda,
			     const EL_INFO *el_info);
void trace_to_bulk_coords_0d(REAL_B result,
			     const REAL_B lambda,
			     const EL_INFO *el_info);

void bulk_to_trace_coords_2d(REAL_B result,
			     const REAL_B lambda,
			     const EL_INFO *el_info);
void bulk_to_trace_coords_1d(REAL_B result,
			     const REAL_B lambda,
			     const EL_INFO *el_info);
void bulk_to_trace_coords_0d(REAL_B result,
			     const REAL_B lambda,
			     const EL_INFO *el_info);

static inline
void trace_to_bulk_coords(REAL_B result,
			  const REAL_B lambda,
			  const EL_INFO *el_info)
{
  FUNCNAME("trace_to_bulk_coords");
  switch (el_info->mesh->dim) {
  case 2: trace_to_bulk_coords_2d(result, lambda, el_info); break;
  case 1: trace_to_bulk_coords_1d(result, lambda, el_info); break;
  case 0: trace_to_bulk_coords_0d(result, lambda, el_info); break;
  default:
    ERROR_EXIT("Illegal dimension: %d\n", el_info->mesh->dim);
    break;
  }
}

static inline
void trace_to_bulk_coords_dim(int dim,
			      REAL_B result,
			      const REAL_B lambda,
			      const EL_INFO *el_info)
{
  FUNCNAME("trace_to_bulk_coords_dim");
  switch (dim) {
  case 2: trace_to_bulk_coords_2d(result, lambda, el_info); break;
  case 1: trace_to_bulk_coords_1d(result, lambda, el_info); break;
  case 0: trace_to_bulk_coords_0d(result, lambda, el_info); break;
  default:
    ERROR_EXIT("Illegal dimension: %d\n", el_info->mesh->dim);
    break;
  }
}

static inline
void bulk_to_trace_coords(REAL_B result,
			  const REAL_B lambda,
			  const EL_INFO *el_info)
{
  FUNCNAME("bulk_to_trace_coords");
  switch (el_info->mesh->dim) {
  case 2: bulk_to_trace_coords_2d(result, lambda, el_info); break;
  case 1: bulk_to_trace_coords_1d(result, lambda, el_info); break;
  case 0: bulk_to_trace_coords_0d(result, lambda, el_info); break;
  default:
    ERROR_EXIT("Illegal dimension: %d\n", el_info->mesh->dim);
    break;
  }
}

static inline
void bulk_to_trace_coords_dim(int dim,
			      REAL_B result,
			      const REAL_B lambda,
			      const EL_INFO *el_info)
{
  FUNCNAME("bulk_to_trace_coords_dim");
  switch (dim) {
  case 2: bulk_to_trace_coords_2d(result, lambda, el_info); break;
  case 1: bulk_to_trace_coords_1d(result, lambda, el_info); break;
  case 0: bulk_to_trace_coords_0d(result, lambda, el_info); break;
  default:
    ERROR_EXIT("Illegal dimension: %d\n", el_info->mesh->dim);
    break;
  }
}

#define TRACE_DOF_VEC_PROTO(TYPE, typename)				\
  extern void trace_##typename##_vec(DOF_##TYPE##_VEC *svec,		\
				     const DOF_##TYPE##_VEC *mvec)
TRACE_DOF_VEC_PROTO(REAL, dof_real);
TRACE_DOF_VEC_PROTO(REAL_D, dof_real_d);
TRACE_DOF_VEC_PROTO(INT, dof_int);
TRACE_DOF_VEC_PROTO(DOF, dof_dof);
TRACE_DOF_VEC_PROTO(DOF, int_dof);
TRACE_DOF_VEC_PROTO(UCHAR, dof_uchar);
TRACE_DOF_VEC_PROTO(SCHAR, dof_schar);
TRACE_DOF_VEC_PROTO(PTR, dof_ptr);
extern
void trace_dof_real_vec_d(DOF_REAL_VEC_D *svec, const DOF_REAL_VEC_D *mvec);

extern void update_master_matrix(DOF_MATRIX *m_dof_matrix,
				 const EL_MATRIX_INFO *s_minfo,
				 MatrixTranspose transpose);
extern void update_master_real_vec(DOF_REAL_VEC *m_drv,
				   const EL_VEC_INFO *s_vec_info);
extern void update_master_real_d_vec(DOF_REAL_D_VEC *m_drdv,
				     const EL_VEC_D_INFO *s_vecd_info);

/***   file level.c    ********************************************************/
extern REAL level_element_det_2d(const REAL_D coord[]);
extern void level_coord_to_world_2d(const REAL_D coord[],
				    const REAL_B lambda,
				    REAL_D world);
extern void level_coord_to_el_coord_2d(const REAL_B v_lambda[],
				       const REAL_B lambda,
				       REAL_B el_lambda);
extern REAL level_element_det_3d(const REAL_D coord[]);
extern void level_coord_to_world_3d(const REAL_D coord[],
				    const REAL_B lambda,
				    REAL_D world);
extern void level_coord_to_el_coord_3d(const REAL_B v_lambda[],
				       const REAL_B lambda,
				       REAL_B el_lambda);

extern int find_level(MESH *mesh, FLAGS fill_flag, const DOF_REAL_VEC *Level,
		      REAL value,
		      int (*init)(const EL_INFO *el_info,
				  REAL v[],
				  int N, int wall, const REAL_B lambda[]),
		      void (*cal)(const EL_INFO *el_info,
				  REAL v[],
				  int i,
				  int wall, const REAL_B lambda[],
				  const REAL_D coord[]));
extern void set_element_mark(MESH *mesh, FLAGS fill_flag, S_CHAR mark);

/***   file numint.c   ********************************************************/
const QUAD *get_quadrature(int dim, int degree);
void register_quadrature(QUAD *quad);
bool new_quadrature(const QUAD *quad);
const QUAD *get_product_quad(const QUAD *quad);
const QUAD *get_lumping_quadrature(int dim);
static inline const QUAD_EL_CACHE *fill_quad_el_cache(const EL_INFO *el_info,
						      const QUAD *quad,
						      FLAGS need);
void print_quadrature(const QUAD *quad);
void check_quadrature(const QUAD *quad);

REAL  integrate_std_simp(const QUAD *quad, REAL (*f)(const REAL_B lambda));

static inline
const REAL   *f_at_qp(REAL quad_vec[],
		      const QUAD *quad, REAL (*f)(const REAL_B lambda));
static inline
const REAL_D *f_d_at_qp(REAL_D quad_vec[],
			const QUAD *quad,
			const REAL *(*f)(const REAL_B lambda));
static inline
const REAL_D *grd_f_at_qp(REAL_D quad_vec[],
			  const QUAD *quad,
			  const REAL *(*grd_f)(const REAL_B lambda));
static inline
const REAL_DD *grd_f_d_at_qp(REAL_DD quad_vec[],
			     const QUAD *quad,
			     const REAL_D *(*grd_f)(const REAL_B lambda));

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

static inline
const REAL *f_loc_at_qp(REAL quad_vec[],
			const EL_INFO *el_info, const QUAD *quad,
			LOC_FCT_AT_QP f_at_qp, void *ud);
static inline
const REAL_D *f_loc_d_at_qp(REAL_D quad_vec[],
			    const EL_INFO *el_info, const QUAD *quad,
			    LOC_FCT_D_AT_QP f_at_qp, void *ud);
static inline
const REAL_D *grd_f_loc_at_qp(REAL_D quad_vec[],
			      const EL_INFO *el_info, const QUAD *quad,
			      const REAL_BD Lambda,
			      GRD_LOC_FCT_AT_QP grd_f_at_qp, void *ud);
static inline
const REAL_DD *grd_f_loc_d_at_qp(REAL_DD quad_vec[],
				 const EL_INFO *el_info, const QUAD *quad,
				 const REAL_BD Lambda,
				 GRD_LOC_FCT_D_AT_QP grd_f_at_qp, void *ud);
static inline
const REAL_D *
param_grd_f_loc_at_qp(REAL_D quad_vec[],
		      const EL_INFO *el_info,
		      const QUAD *quad, const REAL_BD Lambda[],
		      GRD_LOC_FCT_AT_QP grd_f_at_qp, void *ud);
static inline
const REAL_DD *
param_grd_f_loc_d_at_qp(REAL_DD quad_vec[],
			const EL_INFO *el_info,
			const QUAD *quad, const REAL_BD Lambda[],
			GRD_LOC_FCT_D_AT_QP grd_f_at_qp, void *ud);

const QUAD_FAST *get_quad_fast(const BAS_FCTS *, const QUAD *, FLAGS init_flag);

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);

/***   file refine.c   ********************************************************/
extern U_CHAR refine(MESH *mesh, FLAGS fill_flags);
extern U_CHAR global_refine(MESH *mesh, int mark, FLAGS fill_flags);

/******************************************************************************/

/***   file adapt.c   *********************************************************/
extern void adapt_method_stat(MESH *mesh, ADAPT_STAT *adapt);
extern void adapt_method_instat(MESH *mesh, ADAPT_INSTAT *adapt);
extern int marking(MESH *mesh, ADAPT_STAT *adapt);
extern ADAPT_INSTAT *get_adapt_instat(int dim, const char *name,
				      const char *prefix,
				      int info, ADAPT_INSTAT *adapt_instat);
extern ADAPT_STAT *get_adapt_stat(int dim, const char *name,
				  const char *prefix,
				  int info, ADAPT_STAT *adapt_stat);

/***   file quad_cache.c ******************************************************/
extern const Q00_PSI_PHI *get_q00_psi_phi(const BAS_FCTS *psi,
					  const BAS_FCTS *phi,
					  const QUAD *quad);
extern const Q01_PSI_PHI *get_q01_psi_phi(const BAS_FCTS *psi,
					  const BAS_FCTS *phi,
					  const QUAD *quad);
extern const Q10_PSI_PHI *get_q10_psi_phi(const BAS_FCTS *psi,
					  const BAS_FCTS *phi,
					  const QUAD *quad);
extern const Q11_PSI_PHI *get_q11_psi_phi(const BAS_FCTS *psi,
					  const BAS_FCTS *phi,
					  const QUAD *quad);

extern const Q001_ETA_PSI_PHI *get_q001_eta_psi_phi(const BAS_FCTS *eta,
						    const BAS_FCTS *psi,
						    const BAS_FCTS *phi,
						    const QUAD *quad);
extern const Q010_ETA_PSI_PHI *get_q010_eta_psi_phi(const BAS_FCTS *eta,
						    const BAS_FCTS *psi,
						    const BAS_FCTS *phi,
						    const QUAD *quad);
extern const Q100_ETA_PSI_PHI *get_q100_eta_psi_phi(const BAS_FCTS *eta,
						    const BAS_FCTS *psi,
						    const BAS_FCTS *phi,
						    const QUAD *quad);

/***   file assemble.c   ******************************************************/
extern const EL_MATRIX_INFO *fill_matrix_info(const OPERATOR_INFO *,
					      EL_MATRIX_INFO *res);
extern const EL_MATRIX_INFO *fill_matrix_info_ext(EL_MATRIX_INFO *res,
						  const OPERATOR_INFO *,
						  const BNDRY_OPERATOR_INFO *,
						  ...);
extern const QUAD_TENSOR *get_quad_matrix(const FE_SPACE *row_fe_space,
					  const FE_SPACE *col_fe_space,
					  int krn_degree,
					  int n_derivatives);
extern const QUAD_TENSOR *get_quad_tensor(const FE_SPACE *row_fe_space,
					  const FE_SPACE *col_fe_space,
					  const FE_SPACE *depth_fe_space,
					  int krn_degree,
					  int n_derivatives);
extern void free_quad_tensor(const QUAD_TENSOR *qtensor);

/***   file assemble-instat.c   ***********************************************/

extern EL_SYS_INFO_INSTAT *
fill_sys_info_instat(const OPERATOR_INFO *stiff_info,
		     const OPERATOR_INFO *mass_info,
		     const DOF_REAL_VEC *u_h);
extern EL_SYS_INFO_DOW_INSTAT *
fill_sys_info_instat_dow(const OPERATOR_INFO *stiff_info,
			 const OPERATOR_INFO *mass_info,
			 const DOF_REAL_VEC_D *u_h);
static inline EL_SYS_INFO_D_INSTAT *
fill_sys_info_instat_d(const OPERATOR_INFO *stiff_info,
		       const OPERATOR_INFO *mass_info,
		       const DOF_REAL_D_VEC *u_h)
{
  return (EL_SYS_INFO_D_INSTAT *)
    fill_sys_info_instat_dow(
      stiff_info, mass_info, (const DOF_REAL_VEC_D *)u_h);
}
extern void free_sys_info_instat(EL_SYS_INFO_INSTAT *elsii);
extern void free_sys_info_d_instat(EL_SYS_INFO_D_INSTAT *elsii);
extern void update_system_instat(DOF_MATRIX *dof_matrix,
				 DOF_REAL_VEC *f_h,
				 REAL tau,
				 REAL theta,
				 EL_SYS_INFO_INSTAT *elsii);
extern void update_system_instat_dow(DOF_MATRIX *dof_matrix,
				     DOF_REAL_VEC_D *f_h,
				     REAL tau,
				     REAL theta,
				     EL_SYS_INFO_DOW_INSTAT *elsii);
static inline
void update_system_instat_d(DOF_MATRIX *dof_matrix,
			    DOF_REAL_D_VEC *f_h,
			    REAL tau,
			    REAL theta,
			    EL_SYS_INFO_D_INSTAT *elsii)
{
  update_system_instat_dow(dof_matrix, (DOF_REAL_VEC_D *)f_h, tau, theta,
			   (EL_SYS_INFO_DOW_INSTAT *)elsii);
}

/***   file bas_fct.c   *******************************************************/
extern const BAS_FCTS *get_bas_fcts(int dim, const char *name);
extern const BAS_FCTS *get_discontinuous_lagrange(int dim, int degree);
extern const BAS_FCTS *get_lagrange(int dim, int degree);
extern const BAS_FCTS *get_disc_ortho_poly(int dim, int degree);
extern const BAS_FCTS *new_bas_fcts(const BAS_FCTS * bas_fcts);
extern BAS_FCTS *chain_bas_fcts(const BAS_FCTS *head, BAS_FCTS *tail);

typedef const BAS_FCTS *
(*BAS_FCTS_INIT_FCT)(int dim, int dow, const char *name);

extern void add_bas_fcts_plugin(BAS_FCTS_INIT_FCT init_fct);

extern const EL_INT_VEC *
default_get_int_vec(int rvec[], const EL *el, const DOF_INT_VEC *vec);
extern const EL_REAL_VEC *
default_get_real_vec(REAL rvec[], const EL *el, const DOF_REAL_VEC *vec);
extern const EL_REAL_D_VEC *
default_get_real_d_vec(REAL_D rvec[],
		       const EL *el, const DOF_REAL_D_VEC *vec);
extern const EL_REAL_DD_VEC *
default_get_real_dd_vec(REAL_DD rvec[],
			const EL *el, const DOF_REAL_DD_VEC *vec);
extern const EL_REAL_VEC_D *
default_get_real_vec_d(REAL rvec[],
		       const EL *el, const DOF_REAL_VEC_D *vec);
extern const EL_UCHAR_VEC *
default_get_uchar_vec(U_CHAR rvec[],
		      const EL *el, const DOF_UCHAR_VEC *vec);
extern const EL_SCHAR_VEC *
default_get_schar_vec(S_CHAR rvec[],
		      const EL *el, const DOF_SCHAR_VEC *vec);
extern const EL_PTR_VEC *
default_get_ptr_vec(void *rvec[],
		    const EL *el, const DOF_PTR_VEC *vec);

/***   file check.c   *********************************************************/
extern void check_mesh(MESH *mesh);

/***   file element.c   *******************************************************/

/* These routines are partially available as _?d-versions to avoid looking
 * up the dimension. This should be a small efficiency bonus.
 */
extern void
fill_neigh_el_info(EL_INFO *neigh_info,
		   const EL_INFO *el_info, int wall, int rel_perm);
static inline const EL_GEOM_CACHE *
fill_el_geom_cache(const EL_INFO *el_info, FLAGS fill_flag);

#if 0
/* implemented as inline-functions further below */
extern int wall_orientation(int dim, const EL *el, int wall);
extern int world_to_coord(const EL_INFO *el_info, const REAL *,
			  REAL_B);
extern const REAL *coord_to_world(const EL_INFO *, const REAL_B, REAL_D);
extern REAL el_det(const EL_INFO *el_info);
extern REAL el_volume(const EL_INFO *el_info);
extern REAL el_grd_lambda(const EL_INFO *el_info,
			  REAL_BD grd_lam);
#endif /* inlines further below */

/* Dimension dependent routines, 0d, just dummies in most cases. */
extern int wall_orientation_0d(const EL *el, int wall);
extern int wall_rel_orientation_0d(const EL *el, const EL *neigh,
				   int wall, int ov);
extern int world_to_coord_0d(const EL_INFO *el_info, const REAL *, REAL_B);
extern const REAL *coord_to_world_0d(const EL_INFO *, const REAL_B, REAL_D);
extern REAL el_det_0d(const EL_INFO *);
extern REAL el_volume_0d(const EL_INFO *el_info);
extern REAL el_grd_lambda_0d(const EL_INFO *el_info,
			     REAL_BD grd_lam);
extern REAL get_wall_normal_0d(const EL_INFO *el_info, int wall, REAL *normal);

/* Dimension dependent routines, 1d */
extern int wall_orientation_1d(const EL *el, int wall);
extern int wall_rel_orientation_1d(const EL *el, const EL *neigh,
				   int wall, int ov);
extern int world_to_coord_1d(const EL_INFO *el_info, const REAL *,
			     REAL_B);
extern const REAL *coord_to_world_1d(const EL_INFO *, const REAL_B, REAL_D);
extern REAL el_det_1d(const EL_INFO *);
extern REAL el_volume_1d(const EL_INFO *el_info);
extern REAL el_grd_lambda_1d(const EL_INFO *,
			     REAL_BD grd_lam);
extern REAL get_wall_normal_1d(const EL_INFO *el_info, int wall, REAL *normal);

#if DIM_MAX > 1
/* Dimension dependent routines, 2d */
extern int wall_orientation_2d(const EL *el, int wall);
extern int wall_rel_orientation_2d(const EL *el, const EL *neigh,
				   int wall, int ov);
extern int world_to_coord_2d(const EL_INFO *el_info, const REAL *,
			     REAL_B);
extern const REAL *coord_to_world_2d(const EL_INFO *, const REAL_B, REAL_D);
extern REAL el_det_2d(const EL_INFO *);
extern REAL el_volume_2d(const EL_INFO *el_info);
extern REAL el_grd_lambda_2d(const EL_INFO *,
			     REAL_BD grd_lam);
extern REAL get_wall_normal_2d(const EL_INFO *el_info, int wall, REAL *normal);

#if DIM_MAX > 2
/* Dimension dependent routines, 3d */
extern int wall_orientation_3d(const EL *el, int wall);
extern int wall_rel_orientation_3d(const EL *el, const EL *neigh,
				   int wall, int oppv);
extern int world_to_coord_3d(const EL_INFO *el_info, const REAL *, REAL_B);
extern const REAL *coord_to_world_3d(const EL_INFO *, const REAL_B, REAL_D);
extern REAL el_det_3d(const EL_INFO *);
extern REAL el_volume_3d(const EL_INFO *el_info);
extern REAL el_grd_lambda_3d(const EL_INFO *,
			     REAL_BD grd_lam);
extern REAL get_wall_normal_3d(const EL_INFO *el_info, int wall, REAL *normal);
#endif
#endif

/* Below we provide wrapper functions which distinguish the dimension
 * dependent routines by the co-dimension rather than by the dimension
 * of the underlying mesh. We start by defining a preprocessor macro
 * which spares us some typing and especially typos.
 *
 * In addition, we provide wrapper functions which decide by looking
 * at el_info->mesh->dim what to do.
 *
 */
#if DIM_OF_WORLD == 1

# define ALBERTA_CODIM_WRAPPER(dim, ret, name, suf, argtypes, argnames)	\
  static inline ret name##suf argtypes					\
  {									\
    FUNCNAME(#name);							\
									\
    switch (dim) {							\
    case 0: return name##_0d argnames;					\
    case 1: return name##_1d argnames;					\
    default:								\
      ERROR_EXIT("Illegal dim!\n");					\
      return (ret)0; /* shut-off a compiler warning */			\
    }									\
  }									\
  struct _AI_semicolon_dummy

# define ALBERTA_CODIM_ALIAS(ret, name, argtypes, argnames)		\
  static inline ret name##_0cd argtypes { return name##_1d argnames; }	\
  static inline ret name##_1cd argtypes { return name##_0d argnames; }	\
  struct _AI_semicolon_dummy

/* Variants which start at DOW == 2 and thus are empty here */
# define ALBERTA_CODIM_ALIAS_2(ret, name, argtypes, argnames)	\
  struct _AI_semicolon_dummy
# define ALBERTA_VOID_CODIM_ALIAS_2(name, argtypes, argnames)	\
  struct _AI_semicolon_dummy

#elif DIM_OF_WORLD == 2

# define ALBERTA_CODIM_WRAPPER(dim, ret, name, suf, argtypes, argnames)	\
  static inline ret name##suf argtypes					\
  {									\
    FUNCNAME(#name);							\
									\
    switch (dim) {							\
    case 0: return name##_0d argnames;					\
    case 1: return name##_1d argnames;					\
    case 2: return name##_2d argnames;					\
    default:								\
      ERROR_EXIT("Illegal dim!\n");					\
      return (ret)0L; /* shut-off a compiler warning */			\
    }									\
  }									\
  struct _AI_semicolon_dummy

# define ALBERTA_CODIM_ALIAS(ret, name, argtypes, argnames)		\
  static inline ret name##_0cd argtypes { return name##_2d argnames; }	\
  static inline ret name##_1cd argtypes { return name##_1d argnames; }	\
  static inline ret name##_2cd argtypes { return name##_0d argnames; }	\
  struct _AI_semicolon_dummy

/* Variants which start at DOW == 2 */
# define ALBERTA_CODIM_ALIAS_2(ret, name, argtypes, argnames)		\
  static inline ret name##_0cd argtypes { return name##_2d argnames; }	\
  struct _AI_semicolon_dummy
# define ALBERTA_VOID_CODIM_ALIAS_2(name, argtypes, argnames)		\
  static inline void name##_0cd argtypes { name##_2d argnames; }	\
  struct _AI_semicolon_dummy

#elif DIM_OF_WORLD == 3

# define ALBERTA_CODIM_WRAPPER(dim, ret, name, suf, argtypes, argnames)	\
  static inline ret name##suf argtypes					\
  {									\
    FUNCNAME(#name);							\
									\
    switch (dim) {							\
    case 0: return name##_0d argnames;					\
    case 1: return name##_1d argnames;					\
    case 2: return name##_2d argnames;					\
    case 3: return name##_3d argnames;					\
    default:								\
      ERROR_EXIT("Illegal dim!\n");					\
      return (ret)0L; /* shut-off a compiler warning */			\
    }									\
  }									\
  struct _AI_semicolon_dummy

# define ALBERTA_CODIM_ALIAS(ret, name, argtypes, argnames)		\
  static inline ret name##_0cd argtypes { return name##_3d argnames; }	\
  static inline ret name##_1cd argtypes { return name##_2d argnames; }	\
  static inline ret name##_2cd argtypes { return name##_1d argnames; }	\
  static inline ret name##_3cd argtypes { return name##_0d argnames; }	\
  struct _AI_semicolon_dummy

/* Variants which start at DOW == 2 */
# define ALBERTA_CODIM_ALIAS_2(ret, name, argtypes, argnames)		\
  static inline ret name##_0cd argtypes { return name##_3d argnames; }	\
  static inline ret name##_1cd argtypes { return name##_2d argnames; }	\
  struct _AI_semicolon_dummy
# define ALBERTA_VOID_CODIM_ALIAS_2(name, argtypes, argnames)		\
  static inline void name##_0cd argtypes { name##_3d argnames; }	\
  static inline void name##_1cd argtypes { name##_2d argnames; }	\
  struct _AI_semicolon_dummy

#elif DIM_OF_WORLD > 3

# define ALBERTA_CODIM_WRAPPER(dim, ret, name, suf, argtypes, argnames)	\
  static inline ret name##suf argtypes					\
  {									\
    FUNCNAME(#name);							\
									\
    switch (dim) {							\
    case 0: return name##_0d argnames;					\
    case 1: return name##_1d argnames;					\
    case 2: return name##_2d argnames;					\
    case 3: return name##_3d argnames;					\
    default:								\
      ERROR_EXIT("Illegal dim!\n");					\
      return (ret)0L; /* shut-off a compiler warning */			\
    }									\
  }									\
  struct _AI_semicolon_dummy

# define ALBERTA_CODIM_ALIAS(ret, name, argtypes, argnames)	\
  struct _AI_semicolon_dummy
# define ALBERTA_CODIM_ALIAS_2(ret, name, argtypes, argnames)	\
  struct _AI_semicolon_dummy
# define ALBERTA_VOID_CODIM_ALIAS_2(name, argtypes, argnames)	\
  struct _AI_semicolon_dummy

#endif

/* ..._Xcd() alias definitions */
ALBERTA_CODIM_ALIAS(int, world_to_coord,
		    (const EL_INFO *el_info,
		     const REAL *xy,
		     REAL_B lambda),
		    (el_info, xy, lambda));
ALBERTA_CODIM_ALIAS(const REAL *, coord_to_world,
		    (const EL_INFO *el_info, const REAL_B l, REAL_D w),
		    (el_info, l, w));
ALBERTA_CODIM_ALIAS(REAL, el_volume, (const EL_INFO *el_info), (el_info));
ALBERTA_CODIM_ALIAS(REAL, el_det, (const EL_INFO *el_info), (el_info));
ALBERTA_CODIM_ALIAS(REAL, el_grd_lambda,
		    (const EL_INFO *el_info,
		     REAL_BD grd_lam),
		    (el_info, grd_lam));
ALBERTA_CODIM_ALIAS(REAL, get_wall_normal,
		    (const EL_INFO *el_info, int i0, REAL *normal),
		    (el_info, i0, normal));
ALBERTA_CODIM_ALIAS(int, wall_orientation,
		    (const EL *el, int wall),
		    (el, wall));
ALBERTA_CODIM_ALIAS(int, wall_rel_orientation,
		    (const EL *el, const EL *neigh, int wall, int oppv),
		    (el, neigh, wall, oppv));
static const int sorted_wall_vertices_0d[1][1][1] = {{{ 0 }}}; /* dummy */
ALBERTA_CODIM_ALIAS(const int *, sorted_wall_vertices,
		    (int wall, int permno),
		    [wall][permno]);
static const int vertex_of_wall_0d[1][1] = {{ 0 }}; /* dummy */
ALBERTA_CODIM_ALIAS(const int *, vertex_of_wall,
		    (int wall),
		    [wall]);
static const int vertex_of_edge_0d[1][1] = {{ 0 }}; /* dummy */
ALBERTA_CODIM_ALIAS(const int *, vertex_of_edge,
		    (int edge),
		    [edge]);

/* Wrappers which look at el_info->mesh->dim */
ALBERTA_CODIM_WRAPPER(el_info->mesh->dim,
		      int, world_to_coord, /**/,
		      (const EL_INFO *el_info, const REAL *x, REAL_B lambda),
		      (el_info, x, lambda));
ALBERTA_CODIM_WRAPPER(el_info->mesh->dim,
		      const REAL *, coord_to_world, /**/,
		      (const EL_INFO *el_info, const REAL_B lambda, REAL_D x),
		      (el_info, lambda, x));
ALBERTA_CODIM_WRAPPER(el_info->mesh->dim,
		      REAL, el_volume, /**/,
		      (const EL_INFO *el_info), (el_info));
ALBERTA_CODIM_WRAPPER(el_info->mesh->dim,
		      REAL, el_det, /**/,
		      (const EL_INFO *el_info), (el_info));
ALBERTA_CODIM_WRAPPER(el_info->mesh->dim,
		      REAL, el_grd_lambda, /**/,
		      (const EL_INFO *el_info,
		       REAL_BD grd_lam),
		      (el_info, grd_lam));
ALBERTA_CODIM_WRAPPER(el_info->mesh->dim,
		      REAL, get_wall_normal, /**/,
		      (const EL_INFO *el_info, int wall, REAL *normal),
		      (el_info, wall, normal));

/* Wrappers with addtional "dim" as argument */
ALBERTA_CODIM_WRAPPER(dim, int, wall_orientation, /**/,
		      (int dim, const EL *el, int wall),
		      (el, wall));
ALBERTA_CODIM_WRAPPER(dim, int, wall_rel_orientation, /**/,
		      (int dim,
		       const EL *el, const EL *neigh, int wall, int oppv),
		      (el, neigh, wall, oppv));
ALBERTA_CODIM_WRAPPER(dim, const int *, sorted_wall_vertices, /**/,
		      (int dim, int wall, int permno), [wall][permno]);
ALBERTA_CODIM_WRAPPER(dim, const int *, vertex_of_wall, /**/,
		      (int dim, int wall), [wall]);
ALBERTA_CODIM_WRAPPER(dim, const int *, vertex_of_edge, /**/,
		      (int dim, int edge), [edge]);
ALBERTA_CODIM_WRAPPER(dim, int, world_to_coord, _dim,
		      (int dim,
		       const EL_INFO *el_info, const REAL *x, REAL_B lambda),
		      (el_info, x, lambda));
ALBERTA_CODIM_WRAPPER(dim, const REAL *, coord_to_world, _dim,
		      (int dim,
		       const EL_INFO *el_info, const REAL_B lambda, REAL_D x),
		      (el_info, lambda, x));
ALBERTA_CODIM_WRAPPER(dim, REAL, el_volume, _dim,
		      (int dim, const EL_INFO *el_info), (el_info));
ALBERTA_CODIM_WRAPPER(dim, REAL, el_det, _dim,
		      (int dim, const EL_INFO *el_info), (el_info));
ALBERTA_CODIM_WRAPPER(dim, REAL, el_grd_lambda, _dim,
		      (int dim, const EL_INFO *el_info, REAL_BD grd_lam),
		      (el_info, grd_lam));
ALBERTA_CODIM_WRAPPER(dim, REAL, get_wall_normal, _dim,
		      (int dim, const EL_INFO *el_info, int wall, REAL *normal),
		      (el_info, wall, normal));

/* Some special wrapper functions, used for some stuff defined in
 * level.c
 */
ALBERTA_CODIM_ALIAS_2(REAL, level_element_det, (const REAL_D coord[]), (coord));
ALBERTA_VOID_CODIM_ALIAS_2(level_coord_to_world,
			   (const REAL_D coord[],
			    const REAL_B lambda,
			    REAL_D world),
			   (coord, lambda, world));
ALBERTA_VOID_CODIM_ALIAS_2(level_coord_to_el_coord,
			   (const REAL_B v_lambda[],
			    const REAL_B lambda,
			    REAL_B el_lambda),
			   (v_lambda, lambda, el_lambda));

/***   file estimator{_dowb}.c   **********************************************/

/* The values accepted by the f_flags arguments of ellipt_est() &
 * friends.
 */
#define INIT_UH     1
#define INIT_GRD_UH 2

/*
 {
 const void *est_handle;
 REAL est_el;
 const EL_GEOM_CACHE *elgc;

 est_handle = estimator_init(...);
 TRAVERSE_FIRST(mesh, -1, fill_flag) {
 est_el = element_est(el_info, est_handle);
 #if needed
 uh = element_est_uh(est_handle);
 grd_uh = element_est_grd_uh(est_handle);
 el_gc = fill_el_geom_cache(el_info, FILL_EL_DET);
 est_el += el_gc->det * additional_stuff(el_cache,...);
 #endif
 element_est_finish(est_el, est_handle);
 } TRAVERSE_NEXT();

 estimate = estimator_finish(..., est_handle)
 }
*/

extern const void *ellipt_est_init(const DOF_REAL_VEC *uh,
				   ADAPT_STAT *adapt,
				   REAL *(*rw_est)(EL *),
				   REAL *(*rw_estc)(EL *),
				   const QUAD *quad,
				   const WALL_QUAD *wall_quad,
				   NORM norm,
				   REAL C[3],
				   const REAL_DD A,
				   const BNDRY_FLAGS dirichlet_bndry,
				   REAL (*f)(const EL_INFO *el_info,
					     const QUAD *quad,
					     int qp,
					     REAL uh_qp,
					     const REAL_D grd_uh_gp),
				   FLAGS f_flags,
				   REAL (*gn)(const EL_INFO *el_info,
					      const QUAD *quad,
					      int qp,
					      REAL uh_qp,
					      const REAL_D normal),
				   FLAGS gn_flags);
extern const void *heat_est_init(const DOF_REAL_VEC *uh,
				 const DOF_REAL_VEC *uh_old,
				 ADAPT_INSTAT *adapt,
				 REAL *(*rw_est)(EL *),
				 REAL *(*rw_estc)(EL *),
				 const QUAD *quad,
				 const WALL_QUAD *wall_quad,
				 REAL C[4],
				 const REAL_DD A,
				 const BNDRY_FLAGS dirichlet_bndry,
				 REAL (*f)(const EL_INFO *el_info,
					   const QUAD *quad,
					   int qp,
					   REAL uh_qp,
					   const REAL_D grd_uh_gp,
					   REAL time),
				 FLAGS f_flags,
				 REAL (*gn)(const EL_INFO *el_info,
					    const QUAD *quad,
					    int qp,
					    REAL uh_qp,
					    const REAL_D normal,
					    REAL time),
				 FLAGS gn_flags);

extern REAL element_est(const EL_INFO *el_info, const void *est_handle);
extern void element_est_finish(const EL_INFO *el_info,
			       REAL est_el, const void *est_handle);
extern REAL ellipt_est_finish(ADAPT_STAT *adapt, const void *est_handle);
extern REAL heat_est_finish(ADAPT_INSTAT *adapt, const void *est_handle);
extern const REAL *element_est_uh(const void *est_handle);
extern const REAL_D *element_est_grd_uh(const void *est_handle);
extern const EL_REAL_VEC *element_est_uh_loc(const void *est_handle);
extern const EL_REAL_VEC *element_est_uh_old_loc(const void *est_handle);

extern REAL ellipt_est(const DOF_REAL_VEC *uh,
		       ADAPT_STAT *adapt,
		       REAL *(*rw_est)(EL *),
		       REAL *(*rw_estc)(EL *),
		       int quad_degree,
		       NORM norm,
		       REAL C[3],
		       const REAL_DD A,
		       const BNDRY_FLAGS dirichlet_bndry,
		       REAL (*f)(const EL_INFO *el_info,
				 const QUAD *quad,
				 int qp,
				 REAL uh_qp,
				 const REAL_D grd_uh_gp),
		       FLAGS f_flags,
		       REAL (*gn)(const EL_INFO *el_info,
				  const QUAD *quad,
				  int qp,
				  REAL uh_qp,
				  const REAL_D normal),
		       FLAGS gn_flags);
extern REAL heat_est(const DOF_REAL_VEC *uh,
		     const DOF_REAL_VEC *uh_old,
		     ADAPT_INSTAT *adapt,
		     REAL *(*rw_est)(EL *),
		     REAL *(*rw_estc)(EL *),
		     int quad_degree,
		     REAL C[4],
		     const REAL_DD A,
		     const BNDRY_FLAGS dirichlet_bndry,
		     REAL (*f)(const EL_INFO *el_info,
			       const QUAD *quad,
			       int qp,
			       REAL uh_qp,
			       const REAL_D grd_uh_gp,
			       REAL time),
		     FLAGS f_flags,
		     REAL (*gn)(const EL_INFO *el_info,
				const QUAD *quad,
				int qp,
				REAL uh_qp,
				const REAL_D normal,
				REAL time),
		     FLAGS gn_flags);

/***   file estimator_dowb.c **************************************************/
extern const void *ellipt_est_dow_init(const DOF_REAL_VEC_D *uh,
				       ADAPT_STAT *adapt,
				       REAL *(*rw_est)(EL *),
				       REAL *(*rw_estc)(EL *),
				       const QUAD *quad,
				       const WALL_QUAD *wall_quad,
				       NORM norm,
				       REAL C[3],
				       const void *A,
				       MATENT_TYPE A_type,
				       MATENT_TYPE A_blocktype,
				       bool sym_grad,
				       const BNDRY_FLAGS dirichlet_bndry,
				       const REAL *(*f)(REAL_D result,
							const EL_INFO *el_info,
							const QUAD *quad,
							int qp,
							const REAL_D uh_qp,
							const REAL_DD grd_uh_gp),
				       FLAGS f_flags,
				       const REAL *(*gn)(REAL_D result,
							 const EL_INFO *el_info,
							 const QUAD *quad,
							 int qp,
							 const REAL_D uh_qp,
							 const REAL_D normal),
				       FLAGS gn_flags);
extern const void *heat_est_dow_init(const DOF_REAL_VEC_D *uh,
				     const DOF_REAL_VEC_D *uh_old,
				     ADAPT_INSTAT *adapt,
				     REAL *(*rw_est)(EL *),
				     REAL *(*rw_estc)(EL *),
				     const QUAD *quad,
				     const WALL_QUAD *wall_quad,
				     REAL C[4],
				     const void *A,
				     MATENT_TYPE A_type,
				     MATENT_TYPE A_blocktype,
				     bool sym_grad,
				     const BNDRY_FLAGS dirichlet_bndry,
				     const REAL *(*f)(REAL_D result,
						      const EL_INFO *el_info,
						      const QUAD *quad,
						      int qp,
						      const REAL_D uh_qp,
						      const REAL_DD grd_uh_gp,
						      REAL time),
				     FLAGS f_flags,
				     const REAL *(*gn)(REAL_D result,
						       const EL_INFO *el_info,
						       const QUAD *quad,
						       int qp,
						       const REAL_D uh_qp,
						       const REAL_D normal,
						       REAL time),
				     FLAGS gn_flags);

extern REAL element_est_dow(const EL_INFO *el_info, const void *est_handle);
extern void element_est_dow_finish(const EL_INFO *el_info,
				   REAL est_el, const void *est_handle);
extern REAL ellipt_est_dow_finish(ADAPT_STAT *adapt, const void *est_handle);
extern REAL heat_est_dow_finish(ADAPT_INSTAT *adapt, const void *est_handle);
extern const REAL_D *element_est_uh_dow(const void *est_handle);
extern const REAL_DD *element_est_grd_uh_dow(const void *est_handle);
extern const EL_REAL_VEC_D *element_est_uh_loc_dow(const void *est_handle);
extern const EL_REAL_VEC_D *element_est_uh_old_loc_dow(const void *est_handle);

extern REAL ellipt_est_dow(const DOF_REAL_VEC_D *uh,
			   ADAPT_STAT *adapt,
			   REAL *(*rw_est)(EL *),
			   REAL *(*rw_estc)(EL *),
			   int quad_degree,
			   NORM norm,
			   REAL C[3],
			   const void *A,
			   MATENT_TYPE A_type,
			   MATENT_TYPE A_blocktype,
			   bool sym_grad,
			   const BNDRY_FLAGS dirichlet_bndry,
			   const REAL *(*f)(REAL_D result,
					    const EL_INFO *el_info,
					    const QUAD *quad,
					    int qp,
					    const REAL_D uh_qp,
					    const REAL_DD grd_uh_gp),
			   FLAGS f_flags,
			   const REAL *(*gn)(REAL_D result,
					     const EL_INFO *el_info,
					     const QUAD *quad,
					     int qp,
					     const REAL_D uh_qp,
					     const REAL_D normal),
			   FLAGS gn_flags);
extern REAL heat_est_dow(const DOF_REAL_VEC_D *uh,
			 const DOF_REAL_VEC_D *uh_old,
			 ADAPT_INSTAT *adapt,
			 REAL *(*rw_est)(EL *),
			 REAL *(*rw_estc)(EL *),
			 int quad_degree,
			 REAL C[4],
			 const void *A,
			 MATENT_TYPE A_type,
			 MATENT_TYPE A_blocktype,
			 bool sym_grad,
			 const BNDRY_FLAGS dirichlet_bndry,
			 const REAL *(*f)(REAL_D result,
					  const EL_INFO *el_info,
					  const QUAD *quad,
					  int qp,
					  const REAL_D uh_qp,
					  const REAL_DD grd_uh_gp,
					  REAL time),
			 FLAGS f_flags,
			 const REAL *(*gn)(REAL_D result,
					   const EL_INFO *el_info,
					   const QUAD *quad,
					   int qp,
					   const REAL_D uh_qp,
					   const REAL_D normal,
					   REAL time),
			 FLAGS gn_flags);

/* DOF_REAL_D_VEC versions */
static inline
const void *ellipt_est_d_init(const DOF_REAL_D_VEC *uh,
			      ADAPT_STAT *adapt,
			      REAL *(*rw_est)(EL *),
			      REAL *(*rw_estc)(EL *),
			      const QUAD *quad,
			      const WALL_QUAD *wall_quad,
			      NORM norm,
			      REAL C[3],
			      const void *A,
			      MATENT_TYPE A_type,
			      MATENT_TYPE A_blocktype,
			      bool sym_grad,
			      const BNDRY_FLAGS dirichlet_bndry,
			      const REAL *(*f)(REAL_D result,
					       const EL_INFO *el_info,
					       const QUAD *quad,
					       int qp,
					       const REAL_D uh_qp,
					       const REAL_DD grd_uh_gp),
			      FLAGS f_flags,
			      const REAL *(*gn)(REAL_D result,
						const EL_INFO *el_info,
						const QUAD *quad,
						int qp,
						const REAL_D uh_qp,
						const REAL_D normal),
			      FLAGS gn_flags)
{
  return ellipt_est_dow_init(
    (const DOF_REAL_VEC_D *)uh,
    adapt, rw_est, rw_estc, quad, wall_quad, norm, C, A, A_type, A_blocktype,
    sym_grad, dirichlet_bndry, f, f_flags, gn, gn_flags);
}

static inline
REAL ellipt_est_d(const DOF_REAL_D_VEC *uh,
		  ADAPT_STAT *adapt,
		  REAL *(*rw_est)(EL *),
		  REAL *(*rw_estc)(EL *),
		  int quad_degree,
		  NORM norm,
		  REAL C[3],
		  const void *A,
		  MATENT_TYPE A_type,
		  MATENT_TYPE A_blocktype,
		  bool sym_grad,
		  const BNDRY_FLAGS dirichlet_bndry,
		  const REAL *(*f)(REAL_D result,
				   const EL_INFO *el_info,
				   const QUAD *quad,
				   int qp,
				   const REAL_D uh_qp,
				   const REAL_DD grd_uh_gp),
		  FLAGS f_flags,
		  const REAL *(*gn)(REAL_D result,
				    const EL_INFO *el_info,
				    const QUAD *quad,
				    int qp,
				    const REAL_D uh_qp,
				    const REAL_D normal),
		  FLAGS gn_flags)
{
  return ellipt_est_dow(
    (const DOF_REAL_VEC_D *)uh,
    adapt, rw_est, rw_estc, quad_degree, norm, C, A, A_type, A_blocktype,
    sym_grad, dirichlet_bndry, f, f_flags, gn, gn_flags);
}

static inline
const void *heat_est_d_init(const DOF_REAL_D_VEC *uh,
			    const DOF_REAL_D_VEC *uh_old,
			    ADAPT_INSTAT *adapt,
			    REAL *(*rw_est)(EL *),
			    REAL *(*rw_estc)(EL *),
			    const QUAD *quad,
			    const WALL_QUAD *wall_quad,
			    REAL C[4],
			    const void *A,
			    MATENT_TYPE A_type,
			    MATENT_TYPE A_blocktype,
			    bool sym_grad,
			    const BNDRY_FLAGS dirichlet_bndry,
			    const REAL *(*f)(REAL_D result,
					     const EL_INFO *el_info,
					     const QUAD *quad,
					     int qp,
					     const REAL_D uh_qp,
					     const REAL_DD grd_uh_gp,
					     REAL time),
			    FLAGS f_flags,
			    const REAL *(*gn)(REAL_D result,
					      const EL_INFO *el_info,
					      const QUAD *quad,
					      int qp,
					      const REAL_D uh_qp,
					      const REAL_D normal,
					      REAL time),
			    FLAGS gn_flags)
{
  return heat_est_dow_init(
    (const DOF_REAL_VEC_D *)uh, (const DOF_REAL_VEC_D *)uh_old,
    adapt, rw_est, rw_estc, quad, wall_quad, C, A, A_type, A_blocktype,
    sym_grad, dirichlet_bndry, f, f_flags, gn, gn_flags);
}

static inline
REAL heat_est_d(const DOF_REAL_D_VEC *uh,
		const DOF_REAL_D_VEC *uh_old,
		ADAPT_INSTAT *adapt,
		REAL *(*rw_est)(EL *),
		REAL *(*rw_estc)(EL *),
		int quad_degree,
		REAL C[4],
		const void *A,
		MATENT_TYPE A_type,
		MATENT_TYPE A_blocktype,
		bool sym_grad,
		const BNDRY_FLAGS dirichlet_bndry,
		const REAL *(*f)(REAL_D result,
				 const EL_INFO *el_info,
				 const QUAD *quad,
				 int qp,
				 const REAL_D uh_qp,
				 const REAL_DD grd_uh_gp,
				 REAL time),
		FLAGS f_flags,
		const REAL *(*gn)(REAL_D result,
				  const EL_INFO *el_info,
				  const QUAD *quad,
				  int qp,
				  const REAL_D uh_qp,
				  const REAL_D normal,
				  REAL time),
		FLAGS gn_flags)
{
  return heat_est_dow(
    (const DOF_REAL_VEC_D *)uh, (const DOF_REAL_VEC_D *)uh_old,
    adapt, rw_est, rw_estc, quad_degree, C, A, A_type, A_blocktype,
    sym_grad, dirichlet_bndry, f, f_flags, gn, gn_flags);
}

static inline
REAL element_est_d(const EL_INFO *el_info, const void *est_handle)
{
  return element_est_dow(el_info, est_handle);
}

static inline
void element_est_d_finish(const EL_INFO *el_info,
			  REAL est_el, const void *est_handle)
{
  element_est_dow_finish(el_info, est_el, est_handle);
}

static inline
REAL ellipt_est_d_finish(ADAPT_STAT *adapt, const void *est_handle)
{
  return  ellipt_est_dow_finish(adapt, est_handle);
}

static inline
REAL heat_est_d_finish(ADAPT_INSTAT *adapt, const void *est_handle)
{
  return heat_est_dow_finish(adapt, est_handle);
}

static inline
const REAL_D *element_est_uh_d(const void *est_handle)
{
  return element_est_uh_dow(est_handle);
}

static inline
const REAL_DD *element_est_grd_uh_d(const void *est_handle) {
  return element_est_grd_uh_dow(est_handle);
}

static inline
const EL_REAL_D_VEC *element_est_uh_loc_d(const void *est_handle)
{
  return (const EL_REAL_D_VEC *)element_est_uh_loc_dow(est_handle);
}

static inline
const EL_REAL_D_VEC *element_est_uh_old_loc_d(const void *est_handle)
{
  return (const EL_REAL_D_VEC *)element_est_uh_old_loc_dow(est_handle);
}

/***   file error.c   *********************************************************/
REAL max_err_at_qp(FCT_AT_X u, const DOF_REAL_VEC *, const QUAD *);
REAL max_err_dow_at_qp(FCT_D_AT_X u, const DOF_REAL_VEC_D *, const QUAD *);
REAL max_err_at_vert(FCT_AT_X f, const DOF_REAL_VEC *);
REAL max_err_dow_at_vert(FCT_D_AT_X u, const DOF_REAL_VEC_D *uh);

REAL L2_err(FCT_AT_X u, const DOF_REAL_VEC *uh,
	    const QUAD *quad,
	    bool rel_err, bool mean_value_adjust,
	    REAL *(*rw_err_el)(EL *el), REAL *max_l2_err2);
REAL L2_err_dow(FCT_D_AT_X u, const DOF_REAL_VEC_D *uh,
		const QUAD *quad,
		bool rel_err, bool mean_value_adjust,
		REAL *(*rw_err_el)(EL *el), REAL *max_l2_err2);
REAL L2_err_weighted(FCT_AT_X weight, FCT_AT_X u, const DOF_REAL_VEC *uh,
		     const QUAD *quad,
		     bool rel_err, bool mean_value_adjust,
		     REAL *(*rw_err_el)(EL *el), REAL *max_l2_err2);
REAL L2_err_dow_weighted(FCT_AT_X wieght, FCT_D_AT_X u,
			 const DOF_REAL_VEC_D *uh,
			 const QUAD *quad,
			 bool rel_err, bool mean_value_adjust,
			 REAL *(*rw_err_el)(EL *el), REAL *max_l2_err2);

REAL H1_err(GRD_FCT_AT_X grd_u, const DOF_REAL_VEC *u_h,
	    const QUAD *quad, bool rel_err,
	    REAL *(*rw_err_el)(EL *el), REAL *max_h1_err2);
REAL H1_err_dow(GRD_FCT_D_AT_X grd_u, const DOF_REAL_VEC_D *uh,
		const QUAD *quad, bool rel_err,
		REAL *(*rw_err_esl)(EL *el), REAL *max_h1_err2);
REAL deform_err(GRD_FCT_D_AT_X grd_u, const DOF_REAL_VEC_D *uh,
		const QUAD *quad,
		bool rel_err, REAL *(*rw_err_el)(EL *), REAL *max_el_err2);
REAL H1_err_weighted(FCT_AT_X weight, GRD_FCT_AT_X grd_u,
		     const DOF_REAL_VEC *u_h,
		     const QUAD *quad, bool rel_err,
		     REAL *(*rw_err_el)(EL *el), REAL *max_h1_err2);
REAL H1_err_dow_weighted(FCT_AT_X weight, GRD_FCT_D_AT_X grd_u,
			 const DOF_REAL_VEC_D *uh,
			 const QUAD *quad, bool rel_err,
			 REAL *(*rw_err_esl)(EL *el), REAL *max_h1_err2);
REAL deform_err_weighted(FCT_AT_X weight, GRD_FCT_D_AT_X grd_u,
			 const DOF_REAL_VEC_D *uh,
			 const QUAD *quad,
			 bool rel_err, REAL *(*rw_err_el)(EL *),
			 REAL *max_el_err2);

REAL max_err_at_qp_loc(LOC_FCT_AT_QP u_at_qp, void *ud, FLAGS fill_flag,
		       const DOF_REAL_VEC *uh, const QUAD *quad);
REAL max_err_dow_at_qp_loc(LOC_FCT_D_AT_QP u_at_qp, void *ud, FLAGS fill_flag,
			   const DOF_REAL_VEC_D *uh, const QUAD *quad);
REAL max_err_at_vert_loc(LOC_FCT_AT_QP u_at_qp, void *ud, FLAGS fill_flag,
			 const DOF_REAL_VEC *uh);
REAL max_err_dow_at_vert_loc(LOC_FCT_D_AT_QP u_at_qp, void *ud, FLAGS fill_flag,
			     const DOF_REAL_VEC_D *uh);

REAL L2_err_loc(LOC_FCT_AT_QP u_at_qp, void *ud, FLAGS fill_flag,
		const DOF_REAL_VEC *uh,
		const QUAD *quad,
		bool rel_err, bool mean_value_adjust,
		REAL *(*rw_err_el)(EL *el), REAL *max_l2_err2);
REAL H1_err_loc(GRD_LOC_FCT_AT_QP grd_u_at_qp, void *ud, FLAGS fill_flag,
		const DOF_REAL_VEC *uh,
		const QUAD *quad, bool rel_err, REAL *(*rw_err_el)(EL *),
		REAL *max_h1_err2);
REAL L2_err_loc_dow(LOC_FCT_D_AT_QP u_at_qp, void *ud, FLAGS fill_flag,
		    const DOF_REAL_VEC_D *uh,
		    const QUAD *quad,
		    bool rel_err, bool mean_value_adjust,
		    REAL *(*rw_err_el)(EL *el), REAL *max_l2_err2);
REAL H1_err_loc_dow(GRD_LOC_FCT_D_AT_QP grd_u_at_qp, void *ud, FLAGS fill_flag,
		    const DOF_REAL_VEC_D *uh, const QUAD *quad,
		    bool rel_err,
		    REAL *(*rw_err_el)(EL *), REAL *max_h1_err2);
REAL deform_err_loc(GRD_LOC_FCT_D_AT_QP grd_u_loc, void *ud, FLAGS fill_flag,
		    const DOF_REAL_VEC_D *uh, const QUAD *quad,
		    bool rel_err, 
		    REAL *(*rw_err_el)(EL *), REAL *max_el_err2);

/* DOF_REAL_D_VEC versions */

static inline
REAL max_err_d_at_qp(FCT_D_AT_X u, const DOF_REAL_D_VEC *uh, const QUAD *quad)
{
  return max_err_dow_at_qp(u, (const DOF_REAL_VEC_D *)uh, quad);
}
static inline
REAL max_err_d_at_vert(FCT_D_AT_X u, const DOF_REAL_D_VEC *uh)
{
  return max_err_dow_at_vert(u, (const DOF_REAL_VEC_D *)uh);
}
static inline
REAL L2_err_d(FCT_D_AT_X u, const DOF_REAL_D_VEC *uh,
	      const QUAD *quad,
	      bool rel_err, bool mean_value_adjust,
	      REAL *(*rw_err_el)(EL *el), REAL *max_l2_err2)
{
  return L2_err_dow(u, (const DOF_REAL_VEC_D *)uh,
		    quad, rel_err, mean_value_adjust, rw_err_el, max_l2_err2);
}
static inline
REAL H1_err_d(GRD_FCT_D_AT_X grd_u, const DOF_REAL_D_VEC *uh,
	      const QUAD *quad, bool rel_err,
	      REAL *(*rw_err_el)(EL *el), REAL *max_h1_err2)
{
  return H1_err_dow(grd_u, (const DOF_REAL_VEC_D *)uh,
		    quad, rel_err, rw_err_el, max_h1_err2);
}
static inline
REAL max_err_d_at_qp_loc(LOC_FCT_D_AT_QP u_at_qp, void *ud, FLAGS fill_flag,
			 const DOF_REAL_D_VEC *uh, const QUAD *quad)
{
  return max_err_dow_at_qp_loc(u_at_qp, ud, fill_flag,
			       (const DOF_REAL_VEC_D *)uh, quad);
}
static inline
REAL max_err_d_at_vert_loc(LOC_FCT_D_AT_QP u_at_qp, void *ud, FLAGS fill_flag,
			   const DOF_REAL_D_VEC *uh)
{
  return max_err_dow_at_vert_loc(u_at_qp, ud, fill_flag,
				 (const DOF_REAL_VEC_D *)uh);
}
static inline
REAL L2_err_loc_d(LOC_FCT_D_AT_QP u_at_qp, void *ud, FLAGS fill_flag,
		  const DOF_REAL_D_VEC *uh,
		  const QUAD *quad,
		  bool rel_err, bool mean_value_adjust,
		  REAL *(*rw_err_el)(EL *el), REAL *max_l2_err2)
{
  return L2_err_loc_dow(u_at_qp, ud, fill_flag,
			(const DOF_REAL_VEC_D *)uh,
			quad,
			rel_err, mean_value_adjust, rw_err_el, max_l2_err2);
}
static inline
REAL H1_err_loc_d(GRD_LOC_FCT_D_AT_QP grd_u_at_qp, void *ud, FLAGS fill_flag,
		  const DOF_REAL_D_VEC *uh, const QUAD *quad,
		  bool rel_err,
		  REAL *(*rw_err_el)(EL *), REAL *max_h1_err2)
{
  return H1_err_loc_dow(grd_u_at_qp, ud, fill_flag,
			(const DOF_REAL_VEC_D *)uh,
			quad, rel_err, rw_err_el, max_h1_err2);
}

/***   file eval.c   **********************************************************/

/* evaluation routines are defined as inline functions in the file evaluate.h */

REAL H1_norm_uh(const QUAD *quad, const DOF_REAL_VEC *u_h);
REAL L2_norm_uh(const QUAD *quad, const DOF_REAL_VEC *u_h);
REAL L8_uh_at_qp(REAL *minp, REAL *maxp,
		 const QUAD *quad, const DOF_REAL_VEC *u_h);
REAL H1_norm_uh_dow(const QUAD *quad, const DOF_REAL_VEC_D *u_h);
REAL L2_norm_uh_dow(const QUAD *quad, const DOF_REAL_VEC_D *u_h);
REAL L8_uh_at_qp_dow(REAL *minp, REAL *maxp,
		     const QUAD *quad, const DOF_REAL_VEC_D *u_h);
static inline REAL L2_norm_uh_d(const QUAD *quad, const DOF_REAL_D_VEC *u_h)
{
  return L2_norm_uh_dow(quad, (const DOF_REAL_VEC_D *)u_h);
}
static inline REAL L8_uh_at_qp_d(REAL *minp, REAL *maxp,
				 const QUAD *quad, const DOF_REAL_D_VEC *u_h)
{
  return L8_uh_at_qp_dow(minp, maxp, quad, (const DOF_REAL_VEC_D *)u_h);
}
static inline REAL H1_norm_uh_d(const QUAD *quad, const DOF_REAL_D_VEC *u_h)
{
  return H1_norm_uh_dow(quad, (const DOF_REAL_VEC_D *)u_h);
}

void interpol(FCT_AT_X fct, DOF_REAL_VEC *u_h);
void interpol_dow(FCT_D_AT_X fct, DOF_REAL_VEC_D *uh);
void interpol_loc(DOF_REAL_VEC *vec,
		  LOC_FCT_AT_QP fct_at_qp, void *app_data,
		  FLAGS fill_flags);
void interpol_loc_dow(DOF_REAL_VEC_D *vec,
		      LOC_FCT_D_AT_QP fct_at_qp, void *app_data,
		      FLAGS fill_flags);

static inline
void interpol_d(FCT_D_AT_X f, DOF_REAL_D_VEC *u_h)
{
  interpol_dow(f, (DOF_REAL_VEC_D *)u_h);
}

static inline
void interpol_loc_d(DOF_REAL_D_VEC *vec,
		    LOC_FCT_D_AT_QP fct_at_qp, void *app_data,
		    FLAGS fill_flags)
{
  interpol_loc_dow((DOF_REAL_VEC_D *)vec, fct_at_qp, app_data, fill_flags);
}

/***   file graphXO.c   *******************************************************/
GRAPH_WINDOW graph_open_window(const char *title, const char *geometry,
			       REAL *world, MESH *mesh);
void graph_close_window(GRAPH_WINDOW win);
void graph_clear_window(GRAPH_WINDOW win, const GRAPH_RGBCOLOR c);

void graph_mesh(GRAPH_WINDOW win, MESH *mesh, const GRAPH_RGBCOLOR c,
		FLAGS flag);
void graph_drv(GRAPH_WINDOW win, const DOF_REAL_VEC *uh,
	       REAL min, REAL max, int refine);
void graph_drv_d(GRAPH_WINDOW win, const DOF_REAL_D_VEC *uh,
		 REAL min, REAL max, int refine);
void graph_el_est(GRAPH_WINDOW win, MESH *mesh, REAL (*get_el_est)(EL *el),
		  REAL min, REAL max);

void graph_point(GRAPH_WINDOW, const REAL [2], const GRAPH_RGBCOLOR, float);
void graph_points(GRAPH_WINDOW win, int np, REAL (*p)[2],
		  const GRAPH_RGBCOLOR c, float ps);
void graph_line(GRAPH_WINDOW, const REAL [2], const REAL [2],
		const GRAPH_RGBCOLOR, float);

void graph_fvalues_2d(GRAPH_WINDOW win, MESH *mesh,
		      REAL(*fct)(const EL_INFO *el_info, const REAL *lambda),
		      FLAGS flags, REAL min, REAL max, int refine);
void graph_level_2d(GRAPH_WINDOW win, const DOF_REAL_VEC *v, REAL level,
		    const GRAPH_RGBCOLOR c, int refine);
void graph_levels_2d(GRAPH_WINDOW win, const DOF_REAL_VEC *v,
		     int n, REAL const *levels, const GRAPH_RGBCOLOR *color,
		     int refine);
void graph_level_d_2d(GRAPH_WINDOW, const DOF_REAL_D_VEC *,
		      REAL, const GRAPH_RGBCOLOR, int);
void graph_levels_d_2d(GRAPH_WINDOW, const DOF_REAL_D_VEC *,
		       int, const REAL *, const GRAPH_RGBCOLOR *, int);

/* multigrid level display routines:  */
void graph_mesh_mg_2d(GRAPH_WINDOW win, MESH *mesh, const GRAPH_RGBCOLOR c,
		      FLAGS flags, int mg_level);
void graph_values_mg_2d(GRAPH_WINDOW win, const DOF_REAL_VEC *v,
			REAL min, REAL max, int refine,
			int mg_level, const FE_SPACE *fe_space,
			const int *sort_dof_invers);

/***   file l2scp.c   *********************************************************/
void L2scp_fct_bas(FCT_AT_X f, const QUAD *quad, DOF_REAL_VEC *fh);
void L2scp_fct_bas_dow(FCT_D_AT_X, const QUAD *quad, DOF_REAL_VEC_D *fhd);
void L2scp_fct_bas_loc(DOF_REAL_VEC *fh,
		       LOC_FCT_AT_QP f_at_qp, void *fct_data, FLAGS fill_flag,
		       const QUAD *quad);
void L2scp_fct_bas_loc_dow(DOF_REAL_VEC_D *fh,
			   LOC_FCT_D_AT_QP f_at_qp, void *ud, FLAGS fill_flag,
			   const QUAD *quad);

void H1scp_fct_bas(GRD_FCT_AT_X f, const QUAD *quad, DOF_REAL_VEC *fh);
void H1scp_fct_bas_dow(GRD_FCT_D_AT_X f, const QUAD *quad, DOF_REAL_VEC_D *fh);
void H1scp_fct_bas_loc(DOF_REAL_VEC *fh,
		       GRD_LOC_FCT_AT_QP f, void *fd, FLAGS fill_flag,
		       const QUAD *quad);
void H1scp_fct_bas_loc_dow(DOF_REAL_VEC_D *fh,
			   GRD_LOC_FCT_D_AT_QP f, void *fd, FLAGS fill_flag,
			   const QUAD *quad);

bool bndry_L2scp_fct_bas(DOF_REAL_VEC *fh,
			 REAL (*f)(const REAL_D x, const REAL_D normal),
			 const BNDRY_FLAGS bndry_seg,
			 const WALL_QUAD *quad);
bool bndry_L2scp_fct_bas_loc(DOF_REAL_VEC *fh,
			     LOC_FCT_AT_QP f_at_qp, void *ud, FLAGS fill_flag,
			     const BNDRY_FLAGS bndry_seg,
			     const WALL_QUAD *quad);
bool bndry_L2scp_fct_bas_dow(DOF_REAL_VEC_D *fh,
			     const REAL *(*f)(const REAL_D x,
					      const REAL_D normal,
					      REAL_D result),
			     const BNDRY_FLAGS bndry_seg,
			     const WALL_QUAD *quad);
bool bndry_L2scp_fct_bas_loc_dow(DOF_REAL_VEC_D *fh,
				 LOC_FCT_D_AT_QP f_at_qp, void *ud,
				 FLAGS fill_flag,
				 const BNDRY_FLAGS bndry_seg,
				 const WALL_QUAD *quad);

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

bool bndry_H1scp_fct_bas(DOF_REAL_VEC *fh,
			 const REAL *(*f)(REAL_D result,
					  const REAL_D x, const REAL_D normal),
			 const BNDRY_FLAGS scp_segment,
			 const WALL_QUAD *quad);
bool bndry_H1scp_fct_bas_loc(DOF_REAL_VEC *fh,
			     GRD_LOC_FCT_AT_QP f, void *fd, FLAGS fill_flag,
			     const BNDRY_FLAGS scp_segment,
			     const WALL_QUAD *quad);
bool bndry_H1scp_fct_bas_dow(DOF_REAL_VEC_D *fh,
			     const REAL_D *(*f)(REAL_DD result,
						const REAL_D x,
						const REAL_D normal),
			     const BNDRY_FLAGS bndry_seg,
			     const WALL_QUAD *quad);
bool bndry_H1scp_fct_bas_loc_dow(DOF_REAL_VEC_D *fh,
				 GRD_LOC_FCT_D_AT_QP f_loc,
				 void *fd, FLAGS fill_flag,
				 const BNDRY_FLAGS bndry_seg,
				 const WALL_QUAD *quad);

extern bool
dirichlet_bound(DOF_REAL_VEC *fh, DOF_REAL_VEC *uh,
		DOF_SCHAR_VEC *bound,
		const BNDRY_FLAGS dirichlet_segment,
		REAL (*g)(const REAL_D));
extern bool
dirichlet_bound_dow(DOF_REAL_VEC_D *fh, DOF_REAL_VEC_D *uh,
		    DOF_SCHAR_VEC *bound,
		    const BNDRY_FLAGS dirichlet_segment,
		    const REAL *(*g)(const REAL_D, REAL_D));
extern bool
dirichlet_bound_loc(DOF_REAL_VEC *fh,
		    DOF_REAL_VEC *uh,
		    DOF_SCHAR_VEC *bound,
		    const BNDRY_FLAGS dirichlet_segment,
		    LOC_FCT_AT_QP g_at_qp, void *ud, FLAGS fill_flags);
extern bool
dirichlet_bound_loc_dow(DOF_REAL_VEC_D *fh, DOF_REAL_VEC_D *uh,
			DOF_SCHAR_VEC *bound,
			const BNDRY_FLAGS dirichlet_segment,
			LOC_FCT_D_AT_QP g_at_qp, void *ud, FLAGS fill_flags);

REAL mean_value(MESH *mesh,
		REAL (*f)(const REAL_D),
		const DOF_REAL_VEC *fh, const QUAD *quad);
const REAL *mean_value_dow(MESH *mesh,
			   FCT_D_AT_X f,
			   const DOF_REAL_VEC_D *fh,
			   const QUAD *quad,
			   REAL_D mean);
REAL mean_value_loc(MESH *mesh,
		    LOC_FCT_AT_QP f_at_qp, void *ud, FLAGS fill_flags,
		    const DOF_REAL_VEC *fh, const QUAD *quad);
const REAL *mean_value_loc_dow(REAL_D mean,
			       MESH *mesh,
			       LOC_FCT_D_AT_QP f_at_qp,
			       void *ud, FLAGS fill_flag,
			       const DOF_REAL_VEC_D *fh,
			       const QUAD *quad);

void robin_bound_matrix_info(EL_MATRIX_INFO *robin_info,
			     const FE_SPACE *row_fe_space,
			     const FE_SPACE *col_fe_space,
			     const BNDRY_FLAGS robin_segment,
			     REAL alpha_r,
			     const WALL_QUAD *wall_quad,
			     REAL exponent);
void robin_bound(DOF_MATRIX *matrix,
		 const BNDRY_FLAGS robin_seg,
		 REAL alpha_r,
		 const WALL_QUAD *wall_quad,
		 REAL exponent);

bool boundary_conditions(DOF_MATRIX *matrix,
			 DOF_REAL_VEC *fh,
			 DOF_REAL_VEC *uh,
			 DOF_SCHAR_VEC *bound,
			 const BNDRY_FLAGS dirichlet_segment,
			 REAL (*g)(const REAL_D x),
			 REAL (*gn)(const REAL_D x, const REAL_D normal),
			 REAL alpha_r,
			 const WALL_QUAD *wall_quad);
bool boundary_conditions_loc(DOF_MATRIX *matrix,
			     DOF_REAL_VEC *fh,
			     DOF_REAL_VEC *uh,
			     DOF_SCHAR_VEC *bound,
			     const BNDRY_FLAGS dirichlet_segment,
			     LOC_FCT_AT_QP g_at_qp,
			     LOC_FCT_AT_QP gn_at_qp,
			     void *app_data, FLAGS fill_flags,
			     REAL alpha_r,
			     const WALL_QUAD *wall_quad);
bool boundary_conditions_dow(DOF_MATRIX *matrix,
			     DOF_REAL_VEC_D *fh,
			     DOF_REAL_VEC_D *uh,
			     DOF_SCHAR_VEC *bound,
			     const BNDRY_FLAGS dirichlet_segment,
			     const REAL *(*g)(const REAL_D x, REAL_D res),
			     const REAL *(*gn)(const REAL_D x,
					       const REAL_D normal,
					       REAL_D res),
			     REAL alpha_r,
			     const WALL_QUAD *wall_quad);
bool boundary_conditions_loc_dow(DOF_MATRIX *matrix,
				 DOF_REAL_VEC_D *fh,
				 DOF_REAL_VEC_D *uh,
				 DOF_SCHAR_VEC *bound,
				 const BNDRY_FLAGS dirichlet_segment,
				 LOC_FCT_D_AT_QP g_at_qp,
				 LOC_FCT_D_AT_QP gn_at_qp,
				 void *app_data, FLAGS fill_flags,
				 REAL alpha_r,
				 const WALL_QUAD *wall_quad);

/* ... _d versions */
static inline
void L2scp_fct_bas_d(FCT_D_AT_X fct, const QUAD *quad, DOF_REAL_D_VEC *fhd)
{
  L2scp_fct_bas_dow(fct, quad, (DOF_REAL_VEC_D *)fhd);
}
static inline
void L2scp_fct_bas_loc_d(DOF_REAL_D_VEC *fh,
			 LOC_FCT_D_AT_QP f_at_qp, void *fd, FLAGS fill_flag,
			 const QUAD *quad)
{
  L2scp_fct_bas_loc_dow((DOF_REAL_VEC_D *)fh, f_at_qp, fd, fill_flag, quad);
}
static inline
void H1scp_fct_bas_d(GRD_FCT_D_AT_X f, const QUAD *quad, DOF_REAL_D_VEC *fh)
{
  H1scp_fct_bas_dow(f, quad, (DOF_REAL_VEC_D *)fh);
}
static inline
void H1scp_fct_bas_loc_d(DOF_REAL_D_VEC *fh,
			 GRD_LOC_FCT_D_AT_QP f, void *fd, FLAGS fill_flag,
			 const QUAD *quad)
{
  H1scp_fct_bas_loc_dow((DOF_REAL_VEC_D *)fh, f, fd, fill_flag, quad);
}
static inline
bool bndry_L2scp_fct_bas_d(DOF_REAL_D_VEC *fh,
			   const REAL *(*f)(const REAL_D x,
					    const REAL_D normal,
					    REAL_D result),
			   const BNDRY_FLAGS bndry_seg,
			   const WALL_QUAD *quad)
{
  return bndry_L2scp_fct_bas_dow((DOF_REAL_VEC_D *)fh, f, bndry_seg, quad);
}
static inline
bool bndry_L2scp_fct_bas_loc_d(DOF_REAL_D_VEC *fh,
			       LOC_FCT_D_AT_QP f_at_qp, void *fd,
			       FLAGS fill_flag,
			       const BNDRY_FLAGS bndry_seg,
			       const WALL_QUAD *quad)
{
  return bndry_L2scp_fct_bas_loc_dow((DOF_REAL_VEC_D *)fh,
				     f_at_qp, fd, fill_flag, bndry_seg, quad);
}
static inline
bool dirichlet_bound_d(DOF_REAL_D_VEC *fh, DOF_REAL_D_VEC *uh,
		       DOF_SCHAR_VEC *bound,
		       const BNDRY_FLAGS dirichlet_segment,
		       const REAL *(*g)(const REAL_D, REAL_D))
{
  return dirichlet_bound_dow((DOF_REAL_VEC_D *)fh, (DOF_REAL_VEC_D *)uh,
			     bound, dirichlet_segment, g);
}
static inline
bool dirichlet_bound_loc_d(DOF_REAL_D_VEC *fh, DOF_REAL_D_VEC *uh,
			   DOF_SCHAR_VEC *bound,
			   const BNDRY_FLAGS dirichlet_segment,
			   LOC_FCT_D_AT_QP g_at_qp, void *gd, FLAGS fill_flags)
{
  return dirichlet_bound_loc_dow((DOF_REAL_VEC_D *)fh, (DOF_REAL_VEC_D *)uh,
				 bound, dirichlet_segment,
				 g_at_qp, gd, fill_flags);
}
static inline
const REAL *mean_value_d(MESH *mesh,
			 FCT_D_AT_X f,
			 const DOF_REAL_D_VEC *fh,
			 const QUAD *quad,
			 REAL_D mean)
{
  return mean_value_dow(mesh, f, (const DOF_REAL_VEC_D *)fh, quad, mean);
}
static inline
const REAL *mean_value_loc_d(REAL_D mean,
			     MESH *mesh,
			     LOC_FCT_D_AT_QP f_at_qp, void *fd, FLAGS fill_flag,
			     const DOF_REAL_D_VEC *fh,
			     const QUAD *quad)
{
  return mean_value_loc_dow(
    mean, mesh, f_at_qp, fd, fill_flag, (const DOF_REAL_VEC_D *)fh, quad);
}
static inline
bool boundary_conditions_d(DOF_MATRIX *matrix,
			   DOF_REAL_D_VEC *fh,
			   DOF_REAL_D_VEC *uh,
			   DOF_SCHAR_VEC *bound,
			   const BNDRY_FLAGS dirichlet_segment,
			   const REAL *(*g)(const REAL_D x, REAL_D res),
			   const REAL *(*gn)(const REAL_D x,
					     const REAL_D normal,
					     REAL_D res),
			   REAL alpha_r,
			   const WALL_QUAD *wall_quad)
{
  return boundary_conditions_dow(matrix,
				 (DOF_REAL_VEC_D *)fh, (DOF_REAL_VEC_D *)uh,
				 bound, dirichlet_segment, g, gn, alpha_r,
				 wall_quad);
}
static inline
bool boundary_conditions_loc_d(DOF_MATRIX *matrix,
			       DOF_REAL_D_VEC *fh,
			       DOF_REAL_D_VEC *uh,
			       DOF_SCHAR_VEC *bound,
			       const BNDRY_FLAGS dirichlet_segment,
			       LOC_FCT_D_AT_QP g_at_qp,
			       LOC_FCT_D_AT_QP gn_at_qp,
			       void *gdata, FLAGS fill_flags,
			       REAL alpha_r,
			       const WALL_QUAD *wall_quad)
{
  return boundary_conditions_loc_dow(matrix,
				     (DOF_REAL_VEC_D *)fh, (DOF_REAL_VEC_D *)uh,
				     bound, dirichlet_segment,
				     g_at_qp, gn_at_qp, gdata, fill_flags,
				     alpha_r, wall_quad);
}

/*  file oem_solve.c  *******************************************************/
extern OEM_MV_FCT get_oem_solver(OEM_SOLVER solver);
extern OEM_DATA *init_oem_solve(const DOF_MATRIX *A,
				const DOF_SCHAR_VEC *bound,
				REAL tol, const PRECON *precon,
				int restart, int max_iter, int info);
extern const PRECON *init_oem_precon(const DOF_MATRIX *A,
				     const DOF_SCHAR_VEC *mask,
				     int info, OEM_PRECON precon,
				     ... /* ssor_omega, ssor_n_iter etc. */);
extern const PRECON *vinit_oem_precon(const DOF_MATRIX *A,
				      const DOF_SCHAR_VEC *mask,
				      int info, OEM_PRECON,
				      va_list ap);
const PRECON *init_precon_from_type(const DOF_MATRIX *A,
				    const DOF_SCHAR_VEC *mask,
				    int info,
				    const PRECON_TYPE *prec_type);

extern void release_oem_solve(const OEM_DATA *oem);
extern int oem_mat_vec(void *ud, int dim, const REAL *x, REAL *y);
extern OEM_MV_FCT
init_oem_mat_vec(void **datap,
		 MatrixTranspose transpose, const DOF_MATRIX *A,
		 const DOF_SCHAR_VEC *bound);
extern void exit_oem_mat_vec(void *);

extern int call_oem_solve_s(const OEM_DATA *oem, OEM_SOLVER solver,
			    const DOF_REAL_VEC *f, DOF_REAL_VEC *u);
extern int
oem_solve_s(const DOF_MATRIX *A, const DOF_SCHAR_VEC *bound,
	    const DOF_REAL_VEC *f, DOF_REAL_VEC *u, OEM_SOLVER solver,
	    REAL tol, const PRECON *precon,
	    int restart, int max_iter, int info);

extern int call_oem_solve_dow(const OEM_DATA *oem, OEM_SOLVER solver,
			      const DOF_REAL_VEC_D *f, DOF_REAL_VEC_D *u);
extern int
oem_solve_dow(const DOF_MATRIX *A, const DOF_SCHAR_VEC *bound,
	      const DOF_REAL_VEC_D *f, DOF_REAL_VEC_D *u, OEM_SOLVER solver,
	      REAL tol, const PRECON *precon,
	      int restart, int max_iter, int info);

static inline
int call_oem_solve_d(const OEM_DATA *oem, OEM_SOLVER solver,
		     const DOF_REAL_D_VEC *f, DOF_REAL_D_VEC *u)
{
  return call_oem_solve_dow(oem, solver,
			    (const DOF_REAL_VEC_D *)f, (DOF_REAL_VEC_D *)u);
}
static inline
int oem_solve_d(const DOF_MATRIX *A, const DOF_SCHAR_VEC *bound,
		const DOF_REAL_D_VEC *f, DOF_REAL_D_VEC *u, OEM_SOLVER solver,
		REAL tol, const PRECON *precon,
		int restart, int max_iter, int info)
{
  return oem_solve_dow(A, bound,
		       (const DOF_REAL_VEC_D *)f, (DOF_REAL_VEC_D *)u, solver,
		       tol, precon, restart, max_iter, info);
}

/*  file oem_sp_solve.c  ******************************************************/
typedef struct sp_constraint
{
  const DOF_MATRIX    *B, *Bt;
  const DOF_SCHAR_VEC *bound;
  OEM_MV_FCT          project;
  void                *project_data; /* possibly "OEM_DATA *" */
  OEM_MV_FCT          precon;
  void                *precon_data; /* possibly "OEM_DATA *" */
  REAL                proj_factor, prec_factor;
} SP_CONSTRAINT;

SP_CONSTRAINT *init_sp_constraint(const DOF_MATRIX *B,
				  const DOF_MATRIX *Bt,
				  const DOF_SCHAR_VEC *bound,
				  REAL tol, int info,
				  const DOF_MATRIX *Yproj,
				  OEM_SOLVER Yproj_solver,
				  int Yproj_max_iter, const PRECON *Yproj_prec,
				  const DOF_MATRIX *Yprec,
				  OEM_SOLVER Yprec_solver,
				  int Yprec_max_iter, const PRECON *Yprec_prec,
				  REAL Yproj_frac, REAL Yprec_frac);
void release_sp_constraint(SP_CONSTRAINT *constr);
int oem_sp_schur_solve(OEM_SOLVER sp_solver,
		       REAL sp_tol, int sp_max_iter, int sp_info,
		       OEM_MV_FCT principal_inverse,
		       OEM_DATA *principal_data,
		       const DOF_REAL_VEC_D *f,
		       DOF_REAL_VEC_D *u,
		       SP_CONSTRAINT *constraint,
		       const DOF_REAL_VEC *g,
		       DOF_REAL_VEC *p,
		       ...);
int oem_sp_solve_dow_scl(OEM_SOLVER sp_solver,
			 REAL sp_tol, REAL tol_incr,
			 int sp_max_iter, int sp_info,
			 const DOF_MATRIX *A, const DOF_SCHAR_VEC *bound,
			 OEM_SOLVER A_solver,
			 int A_max_iter, const PRECON *A_precon,
			 DOF_MATRIX *B,
			 DOF_MATRIX *Bt,
			 DOF_MATRIX *Yproj,
			 OEM_SOLVER Yproj_solver,
			 int Yproj_max_iter, const PRECON *Yproj_precon,
			 DOF_MATRIX *Yprec,
			 OEM_SOLVER Yprec_solver,
			 int Yprec_max_iter, const PRECON *Yprec_precon,
			 REAL Yproj_frac, REAL Ypre_frac,
			 const DOF_REAL_VEC_D *f,
			 const DOF_REAL_VEC *g,
			 DOF_REAL_VEC_D *x,
			 DOF_REAL_VEC *y);
REAL sp_dirichlet_bound_dow_scl(MatrixTranspose transpose,
				const DOF_MATRIX *Bt,
				const DOF_SCHAR_VEC *bound,
				const DOF_REAL_VEC_D *u_h,
				DOF_REAL_VEC *g_h);
REAL sp_flux_adjust_dow_scl(MatrixTranspose transpose,
			    const DOF_MATRIX *Bt,
			    const DOF_SCHAR_VEC *bound,
			    const DOF_REAL_VEC_D *u_h,
			    DOF_REAL_VEC *g_h,
			    REAL initial_flux,
			    bool do_adjust);

static inline
int oem_sp_solve_ds(
  OEM_SOLVER sp_solver,
  REAL sp_tol, REAL tol_incr,
  int sp_max_iter, int sp_info,
  const DOF_MATRIX *A, const DOF_SCHAR_VEC *bound,
  OEM_SOLVER A_solver, int A_max_iter, const PRECON *A_prec,
  DOF_MATRIX *B,
  DOF_MATRIX *Bt,
  DOF_MATRIX *Yproj,
  OEM_SOLVER Yproj_solver, int Yproj_max_iter, const PRECON *Yproj_prec,
  DOF_MATRIX *Yprec,
  OEM_SOLVER Yprec_solver, int Yprec_max_iter, const PRECON *Yprec_prec,
  REAL Yproj_frac, REAL Yprec_frac,
  const DOF_REAL_D_VEC *f,
  const DOF_REAL_VEC *g,
  DOF_REAL_D_VEC *x,
  DOF_REAL_VEC *y)
{
  return oem_sp_solve_dow_scl(sp_solver, sp_tol, tol_incr, sp_max_iter, sp_info,
			      A, bound, A_solver, A_max_iter, A_prec,
			      B, Bt,
			      Yproj, Yproj_solver, Yproj_max_iter, Yproj_prec,
			      Yprec, Yprec_solver, Yprec_max_iter, Yprec_prec,
			      Yproj_frac, Yprec_frac,
			      (const DOF_REAL_VEC_D *)f, g,
			      (DOF_REAL_VEC_D *)x, y);
}
static inline
REAL sp_dirichlet_bound_ds(MatrixTranspose transpose,
			   const DOF_MATRIX *Bt,
			   const DOF_SCHAR_VEC *bound,
			   const DOF_REAL_D_VEC *u_h,
			   DOF_REAL_VEC *g_h)
{
  return sp_dirichlet_bound_dow_scl(transpose, Bt, bound,
				    (const DOF_REAL_VEC_D *)u_h, g_h);
}


/*  file parametric.c  ********************************************************/
void use_lagrange_parametric(MESH *mesh, int degree,
			     NODE_PROJ *n_proj,
			     FLAGS flags);
DOF_REAL_D_VEC *get_lagrange_coords(MESH *mesh);
DOF_PTR_VEC *get_lagrange_edge_projections(MESH *mesh);

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 tomesh);

/*--  file sor.c  *************************************************************/
int sor_d(DOF_MATRIX *a, const DOF_REAL_D_VEC *f, const DOF_SCHAR_VEC *b,
	  DOF_REAL_D_VEC *u, REAL omega, REAL tol, int max_iter, int info);
int sor_s(DOF_MATRIX *a, const DOF_REAL_VEC *f, const DOF_SCHAR_VEC *b,
	  DOF_REAL_VEC *u, REAL omega, REAL tol, int max_iter, int info);

/***  file ssor.c  ************************************************************/
int ssor_d(DOF_MATRIX *a, const DOF_REAL_D_VEC *f, const DOF_SCHAR_VEC *b,
	   DOF_REAL_D_VEC *u, REAL omega, REAL tol, int max_iter, int info);
int ssor_s(DOF_MATRIX *a, const DOF_REAL_VEC *f, const DOF_SCHAR_VEC *b,
	   DOF_REAL_VEC *u, REAL omega, REAL tol, int max_iter, int info);

/***   file traverse_r.c  *****************************************************/
extern void mesh_traverse(MESH *mesh, int level, FLAGS fill_flag,
			  void (*el_fct)(const EL_INFO *, void *data),
			  void *data);
extern void fill_macro_info(MESH *mesh, const MACRO_EL *mel, EL_INFO *elinfo);
extern void fill_elinfo(int ichild, FLAGS mask,
			const EL_INFO *parent_info, EL_INFO *elinfo);

/***   file traverse_nr.c *****************************************************/
extern TRAVERSE_STACK *get_traverse_stack(void);
extern void free_traverse_stack(TRAVERSE_STACK *stack);
extern const EL_INFO *traverse_first(TRAVERSE_STACK *stack,
				     MESH *mesh, int level, FLAGS fill_flag);
extern const EL_INFO *traverse_next(TRAVERSE_STACK *stack, const EL_INFO *);
extern const EL_INFO *traverse_neighbour(TRAVERSE_STACK *stack, const EL_INFO *,
					 int neighbour);
extern const EL_INFO *traverse_parent(const TRAVERSE_STACK *stack,
				      const EL_INFO *child);
extern const EL_INFO *subtree_traverse_first(TRAVERSE_STACK *stack,
					     const EL_INFO *local_root,
					     int level, FLAGS fill_flag);
void clear_traverse_mark(TRAVERSE_STACK *stack);

#define TRAVERSE_FIRST(mesh, level, fill_flag)				\
  {									\
    TRAVERSE_STACK *stack = get_traverse_stack();			\
    const EL_INFO  *el_info;						\
    if ((el_info = traverse_first(stack, (mesh), (level), (fill_flag)))) do
#define TRAVERSE_NEXT(/**/)						\
  while ((el_info = traverse_next(stack, el_info)));			\
  free_traverse_stack(stack);						\
  }

/* Compatibility */
#define TRAVERSE_START(mesh, level, fill_flag)	\
  TRAVERSE_FIRST(mesh, level, fill_flag)
#define TRAVERSE_STOP(/**/) TRAVERSE_NEXT()

#define TRAVERSE_NEIGHBOUR(el_info, neighbour)	\
  traverse_neighbour(stack, el_info, neighbour)

/*  file trav_xy.c  ***********************************************************/
extern int find_el_at_pt(MESH *mesh, const REAL_D xy,
			 EL_INFO **el_info_p, FLAGS flag, REAL_B bary,
			 const MACRO_EL *start_mel,
			 const REAL_D xy0, REAL *sp);

/***   file read_mesh.c  ******************************************************/
extern MESH *
read_mesh(const char *fn, REAL *timeptr,
	  NODE_PROJ *(*n_proj)(MESH *, MACRO_EL *, int),
	  MESH *master);
DOF_REAL_VEC   *read_dof_real_vec(const char *, MESH *, FE_SPACE *);
DOF_REAL_D_VEC *read_dof_real_d_vec(const char *, MESH *, FE_SPACE *);
DOF_REAL_VEC_D *read_dof_real_vec_d(const char *, MESH *, FE_SPACE *);
DOF_INT_VEC    *read_dof_int_vec(const char *, MESH *, FE_SPACE *);
DOF_SCHAR_VEC  *read_dof_schar_vec(const char *, MESH *, FE_SPACE *);
DOF_UCHAR_VEC  *read_dof_uchar_vec(const char *, MESH *, FE_SPACE *);

extern MESH *
fread_mesh(FILE *fp, REAL *timeptr,
	   NODE_PROJ *(*n_proj)(MESH *, MACRO_EL *, int),
	   MESH *master);
DOF_REAL_VEC   *fread_dof_real_vec(FILE *fp, MESH *, FE_SPACE *);
DOF_REAL_D_VEC *fread_dof_real_d_vec(FILE *fp, MESH *, FE_SPACE *);
DOF_REAL_VEC_D *fread_dof_real_vec_d(FILE *fp, MESH *, FE_SPACE *);
DOF_INT_VEC    *fread_dof_int_vec(FILE *fp, MESH *, FE_SPACE *);
DOF_SCHAR_VEC  *fread_dof_schar_vec(FILE *fp, MESH *, FE_SPACE *);
DOF_UCHAR_VEC  *fread_dof_uchar_vec(FILE *fp, MESH *, FE_SPACE *);

extern MESH *
read_mesh_xdr(const char *file, REAL *timeptr,
	      NODE_PROJ *(*)(MESH *, MACRO_EL *, int),
	      MESH *master);
DOF_REAL_VEC   *read_dof_real_vec_xdr(const char *, MESH *, FE_SPACE *);
DOF_REAL_D_VEC *read_dof_real_d_vec_xdr(const char *, MESH *, FE_SPACE *);
DOF_REAL_VEC_D *read_dof_real_vec_d_xdr(const char *, MESH *, FE_SPACE *);
DOF_INT_VEC    *read_dof_int_vec_xdr(const char *, MESH *, FE_SPACE *);
DOF_SCHAR_VEC  *read_dof_schar_vec_xdr(const char *, MESH *, FE_SPACE *);
DOF_UCHAR_VEC  *read_dof_uchar_vec_xdr(const char *, MESH *, FE_SPACE *);

extern MESH *
fread_mesh_xdr(FILE *fp, REAL *timeptr,
	       NODE_PROJ *(*)(MESH *, MACRO_EL *, int),
	       MESH *master);
DOF_REAL_VEC   *fread_dof_real_vec_xdr(FILE *fp, MESH *, FE_SPACE *);
DOF_REAL_D_VEC *fread_dof_real_d_vec_xdr(FILE *fp, MESH *, FE_SPACE *);
DOF_REAL_VEC_D *fread_dof_real_vec_d_xdr(FILE *fp, MESH *, FE_SPACE *);
DOF_INT_VEC    *fread_dof_int_vec_xdr(FILE *fp, MESH *, FE_SPACE *);
DOF_SCHAR_VEC  *fread_dof_schar_vec_xdr(FILE *fp, MESH *, FE_SPACE *);
DOF_UCHAR_VEC  *fread_dof_uchar_vec_xdr(FILE *fp, MESH *, FE_SPACE *);

#if 0
/* IFF format (multiple objects in one file) */

/* FORM
   length
   AFEM
   MESH
   length
   data
   AVEC
   length
   data
   AVEC
   length
   data

   where AVEC is one of DRV DRDV DUCV DSCV DINV
 */

#define IFF_TAG_ALBERTA    "AFEM"
#define IFF_TAG_MESH       "MESH"
#define IFF_TAG_REAL_VEC   "DRV "
#define IFF_TAG_REAL_D_VEC "DRDV"
#define IFF_TAG_INT_VEC    "DINV"
#define IFF_TAG_UCHAR_VEC  "DUCV"
#define IFF_TAG_SCHAR_VEC  "DSCV"

static inline bool is_iff_tag(const char tag1[4], const char tag2[4])
{
  return memcmp(tag1, tag2, 4) == 0;
}

/* fread_iff() just read in the tag and the size, afterwards the
 * normal read_xdr() routines can be used to suck in the following
 * data, depending on what TAG contains.
 */
bool fread_iff(FILE *fp, char tag[4], unsigned int *size);
FILE *read_iff(const char *filename, char tag[4], unsigned int *size);
#endif

/***   file write_mesh.c  *****************************************************/
bool write_mesh(MESH *, const char *, REAL);
bool write_dof_real_vec(const DOF_REAL_VEC *, const char *);
bool write_dof_real_vec_d(const DOF_REAL_VEC_D *, const char *);
bool write_dof_real_d_vec(const DOF_REAL_D_VEC *, const char *);
bool write_dof_int_vec(const DOF_INT_VEC *, const char *);
bool write_dof_schar_vec(const DOF_SCHAR_VEC *, const char *);
bool write_dof_uchar_vec(const DOF_UCHAR_VEC *, const char *);

bool fwrite_mesh(MESH *, FILE *fp, REAL);
bool fwrite_dof_real_vec(const DOF_REAL_VEC *, FILE *fp);
bool fwrite_dof_real_d_vec(const DOF_REAL_D_VEC *, FILE *fp);
bool fwrite_dof_int_vec(const DOF_INT_VEC *, FILE *fp);
bool fwrite_dof_schar_vec(const DOF_SCHAR_VEC *, FILE *fp);
bool fwrite_dof_uchar_vec(const DOF_UCHAR_VEC *, FILE *fp);

bool write_mesh_xdr(MESH *, const char *, REAL);
bool write_dof_real_vec_xdr(const DOF_REAL_VEC *, const char *);
bool write_dof_real_vec_d_xdr(const DOF_REAL_VEC_D *, const char *);
bool write_dof_real_d_vec_xdr(const DOF_REAL_D_VEC *, const char *);
bool write_dof_int_vec_xdr(const DOF_INT_VEC *, const char *);
bool write_dof_schar_vec_xdr(const DOF_SCHAR_VEC *, const char *);
bool write_dof_uchar_vec_xdr(const DOF_UCHAR_VEC *, const char *);

bool fwrite_mesh_xdr(MESH *, FILE *fp, REAL time);
bool fwrite_dof_real_vec_xdr(const DOF_REAL_VEC *, FILE *fp);
bool fwrite_dof_real_vec_d_xdr(const DOF_REAL_VEC_D *, FILE *fp);
bool fwrite_dof_real_d_vec_xdr(const DOF_REAL_D_VEC *, FILE *fp);
bool fwrite_dof_int_vec_xdr(const DOF_INT_VEC *, FILE *fp);
bool fwrite_dof_schar_vec_xdr(const DOF_SCHAR_VEC *, FILE *fp);
bool fwrite_dof_uchar_vec_xdr(const DOF_UCHAR_VEC *, FILE *fp);

bool write_dof_matrix_pbm(const DOF_MATRIX *matrix, const char *filename);
bool fwrite_dof_matrix_pbm(const DOF_MATRIX *matrix, FILE *file);

/* Writing of IFF-files, writing requires file-positioning support. */
#if 0
bool fwrite_iff(FILE *fp);
bool fterm_iff(FILE *fp);
FILE *fopen_iff(const char *filename, bool append);
bool fclose_iff(FILE *fp);
bool fwrite_mesh_iff(MESH *mesh, REAL time, FILE *fp);
bool fwrite_dof_real_vec_iff(const DOF_REAL_VEC *v, FILE *fp);
bool fwrite_dof_real_d_vec_iff(const DOF_REAL_D_VEC *v, FILE *fp);
bool fwrite_dof_int_vec_iff(const DOF_INT_VEC *v, FILE *fp);
bool fwrite_dof_schar_vec_iff(const DOF_SCHAR_VEC *v, FILE *fp);
bool fwrite_dof_uchar_vec_iff(const DOF_UCHAR_VEC *v, FILE *fp);
#endif

/***   file write_mesh_gmv.c  *************************************************/

bool write_mesh_gmv(MESH *mesh, const char *file_name, bool write_ascii,
		    bool use_refined_grid,
		    const int n_drv,
		    DOF_REAL_VEC **drv_ptr,
		    const int n_drdv,
		    DOF_REAL_VEC_D **drv_dow_ptr,
		    DOF_REAL_VEC_D *velocity,
		    REAL sim_time);

bool write_dof_vec_gmv(MESH *mesh,
		       const char *mesh_file,
		       const char *file_name, bool write_ascii,
		       bool use_refined_grid,
		       const int n_drv,
		       DOF_REAL_VEC **drv_ptr,
		       const int n_drdv,
		       DOF_REAL_VEC_D **drv_dow_ptr,
		       DOF_REAL_VEC_D *velocity,
		       REAL sim_time);

/*--  file write_mesh_ps.c  ***************************************************/
void write_mesh_ps(MESH *mesh, const char *filename, const char *title,
                   const REAL x[2], const REAL y[2], bool keepaspect,
                   bool draw_bound);

/*******************************************************************************
 *  interface for Lagrange elements for the gltools
 *  file  gltools.c
 ******************************************************************************/
typedef void*    GLTOOLS_WINDOW;

GLTOOLS_WINDOW open_gltools_window(const char *, const char *, const REAL *,
				   MESH *, int);
void close_gltools_window(GLTOOLS_WINDOW);

extern int gltools_get_next_dialog(void);
extern void gltools_set_next_dialog(int dialog);
extern void gltools_est(GLTOOLS_WINDOW, MESH *,
			REAL (*rw_est)(EL *el), REAL min, REAL max, REAL time);
extern void gltools_disp_mesh(GLTOOLS_WINDOW, MESH *,
			      bool mark, const DOF_REAL_VEC_D *, REAL time);
extern void gltools_mesh(GLTOOLS_WINDOW win, MESH *, bool mark, REAL time);
extern void gltools_disp_drv(GLTOOLS_WINDOW, const DOF_REAL_VEC *,
			     REAL min, REAL max, const DOF_REAL_VEC_D *,
			     REAL time);
extern void gltools_drv(GLTOOLS_WINDOW, const DOF_REAL_VEC *,
			REAL min, REAL max, REAL time);
extern void gltools_disp_drv_d(GLTOOLS_WINDOW, const DOF_REAL_VEC_D *,
			       REAL min, REAL max, const DOF_REAL_VEC_D *,
			       REAL time);
extern void gltools_drv_d(GLTOOLS_WINDOW, const DOF_REAL_VEC_D *,
			  REAL min, REAL max, REAL time);
extern void gltools_disp_vec(GLTOOLS_WINDOW, const DOF_REAL_VEC_D *,
			     REAL min, REAL max, const DOF_REAL_VEC_D *,
			     REAL time);
extern void gltools_vec(GLTOOLS_WINDOW, const DOF_REAL_VEC_D *,
			REAL min, REAL max, REAL time);

/*******************************************************************************
 *  interface for Lagrange elements for the dxtools
 *  file  dxtools.c
 ******************************************************************************/

typedef struct dxtools_window DXTOOLS_WINDOW;

extern DXTOOLS_WINDOW *open_dxtools_window(const char *title,
					   const char *geometry);

extern void close_dxtools_window(DXTOOLS_WINDOW *win);
extern void dxtools_mesh(DXTOOLS_WINDOW *win, MESH *mesh);
extern void dxtools_drv(DXTOOLS_WINDOW *win, const DOF_REAL_VEC *u);
extern void dxtools_drdv(DXTOOLS_WINDOW *win, const DOF_REAL_D_VEC *u);

#ifdef __cplusplus
} /* extern "C" */
#endif

/*******************************************************************************
 *
 * A couple of header files containing inline functions for various purposes
 *
 ******************************************************************************/

#include "alberta_inlines.h" /* DIM_OF_WORLD blas */

#include "dof_chains.h"      /* support for chains of objects */

#include "el_vec.h"          /* element vectors and matrices */

#include "evaluate.h"        /* evaluation of finite element functions */

#endif  /* !_ALBERTA_H_ */