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

/*******************************************************************************
 * ALBERTA:  an Adaptive multi Level finite element toolbox using
 *           Bisectioning refinement and Error control by Residual
 *           Techniques for scientific Applications
 *
 * file:     oem_solve.h
 *
 * description: Some inline support routines and definitions for chains of
 *              objects, linked via a doubly linked list. The list-nodes are
 *              assumed to have the name `chain', `row_chain' and `col_chain'.
 *
 *******************************************************************************
 *
 *  authors:   Claus-Justus Heine
 *             Abteilung fuer Angewandte Mathematik
 *             Albert-Ludwigs-Universitaet Freiburg
 *             Hermann-Herder-Str. 10
 *             D-79104 Freiburg im Breisgau, Germany
 *
 *  http://www.alberta-fem.de
 *
 *  (c) by C.-J. Heine (2009)
 *
 ******************************************************************************/
#ifndef _ALBERTA_DOF_CHAINS_H_
#define _ALBERTA_DOF_CHAINS_H_

#include "alberta.h"

#ifndef CHAINED_BASIS_FUNCTIONS
# define CHAINED_BASIS_FUNCTIONS 1
#endif

static inline
int __chain_length(const DBL_LIST_NODE *head) 
{
  const DBL_LIST_NODE *ptr;
  int len;
  
  if (!CHAINED_BASIS_FUNCTIONS) {
    return 1;
  }

  for (len = 1, ptr = head->next; ptr != head; ptr = ptr->next, ++len);

  return len;
}

/* <<< CHAIN_...() */

#define CHAIN_INIT(elem)  DBL_LIST_INIT(&(elem)->chain)
#define CHAIN_INITIALIZER(name) DBL_LIST_INITIALIZER((name).chain)
#define CHAIN_LENGTH(head) __chain_length(&(head)->chain)
#define CHAIN_SINGLE(var)			\
  (!CHAINED_BASIS_FUNCTIONS || dbl_list_empty(&(var)->chain))
#define CHAIN_NEXT(var, type) dbl_list_entry((var)->chain.next, type, chain)
#define CHAIN_PREV(var, type) dbl_list_entry((var)->chain.prev, type, chain)
#define CHAIN_ADD_HEAD(head, elem)			\
  dbl_list_add_head(&(head)->chain, &(elem)->chain)
#define CHAIN_ADD_TAIL(head, elem)			\
  dbl_list_add_tail(&(head)->chain, &(elem)->chain)
#define CHAIN_DEL(elem)			\
  dbl_list_del(&(elem)->chain)
#define CHAIN_ENTRY(node, type)			\
  dbl_list_entry((node)x, type, chain)

/* A loop over all elements of the chain _AFTER_ the first one */
#define CHAIN_FOREACH(ptr, head, type)				\
  if (!CHAIN_SINGLE((head)))					\
    dbl_list_for_each_entry((ptr), &(head)->chain, type, chain)
/* A loop over all elements safe the first, elements may be deleted */
#define CHAIN_FOREACH_SAFE(ptr, next, head, type)			\
  if (!CHAIN_SINGLE((head)))						\
    dbl_list_for_each_entry_safe((ptr), (next), &(head)->chain, type, chain)

/* A loop over all elements of the chain _AFTER_ the first one, in
 * reverse direction.
 */
#define CHAIN_FOREACH_REV(ptr, head, type)				\
  if (!CHAIN_SINGLE((head)))						\
    dbl_list_for_each_entry_rev((ptr), &(head)->chain, type, chain)
/* A loop over all elements safe the first, elements may be deleted,
 * loop goes in reverse direction.
 */
#define CHAIN_FOREACH_REV_SAFE(ptr, next, head, type)			\
  if (!CHAIN_SINGLE((head)))						\
    dbl_list_for_each_entry_rev_safe((ptr), (next), &(head)->chain, type, chain)

/* A do-while-loop over the entire cyclic list. NEVER delete elements
 * from the list.
 */
#if CHAINED_BASIS_FUNCTIONS
# define CHAIN_DO(list, type)    dbl_list_do_cyclic(list, type, chain)
# define CHAIN_WHILE(list, type) dbl_list_while_cyclic(list, type, chain)
# define CHAIN_DO_REV(list, type)		\
  dbl_list_do_cyclic_rev((list), type, chain)
# define CHAIN_WHILE_REV(list, type)		\
  dbl_list_while_cyclic_rev((list), type, chain)
#else
# define CHAIN_DO(list, type)    do
# define CHAIN_WHILE(list, type) while (false)
# define CHAIN_DO_REV(list, type)    do
# define CHAIN_WHILE_REV(list, type) while (false)
#endif

/* >>> */

/* <<< ROW_CHAIN_...() */

#define ROW_CHAIN_INIT(elem)  DBL_LIST_INIT(&(elem)->row_chain)
#define ROW_CHAIN_INITIALIZER(name) DBL_LIST_INITIALIZER((name).col_chain)
#define ROW_CHAIN_LENGTH(head) __chain_length(&(head)->row_chain)
#define ROW_CHAIN_SINGLE(var)						\
  (!CHAINED_BASIS_FUNCTIONS || dbl_list_empty(&(var)->row_chain))
#define ROW_CHAIN_NEXT(var, type)			\
  dbl_list_entry((var)->row_chain.next, type, row_chain)
#define ROW_CHAIN_PREV(var, type)			\
  dbl_list_entry((var)->row_chain.prev, type, row_chain)
#define ROW_CHAIN_ADD_HEAD(head, elem)				\
  dbl_list_add_head(&(head)->row_chain, &(elem)->row_chain)
#define ROW_CHAIN_ADD_TAIL(head, elem)				\
  dbl_list_add_tail(&(head)->row_chain, &(elem)->row_chain)
#define ROW_CHAIN_DEL(elem)						\
  dbl_list_del(&(elem)->row_chain)
#define ROW_CHAIN_FOREACH(ptr, head, type)				\
  if (!ROW_CHAIN_SINGLE((head)))					\
    dbl_list_for_each_entry((ptr), &(head)->row_chain, type, row_chain)
#define ROW_CHAIN_FOREACH_SAFE(ptr, next, head, type)			\
  if (!ROW_CHAIN_SINGLE((head)))					\
    dbl_list_for_each_entry_safe(					\
      (ptr), (next), &(head)->row_chain, type, row_chain)
#define ROW_CHAIN_FOREACH_REV(ptr, head, type)				\
  if (!ROW_CHAIN_SINGLE((head)))					\
    dbl_list_for_each_entry_rev((ptr), &(head)->row_chain, type, row_chain)
#define ROW_CHAIN_FOREACH_REV_SAFE(ptr, next, head, type)		\
  if (!ROW_CHAIN_SINGLE((head)))					\
    dbl_list_for_each_entry_rev_safe(					\
      (ptr), (next), &(head)->row_chain, type, row_chain)
#if CHAINED_BASIS_FUNCTIONS
# define ROW_CHAIN_DO(list, type)		\
  dbl_list_do_cyclic((list), type, row_chain)
# define ROW_CHAIN_WHILE(list, type)		\
  dbl_list_while_cyclic((list), type, row_chain)
# define ROW_CHAIN_DO_REV(list, type)		\
  dbl_list_do_cyclic_rev((list), type, row_chain)
# define ROW_CHAIN_WHILE_REV(list, type)		\
  dbl_list_while_cyclic_rev((list), type, row_chain)
#else
# define ROW_CHAIN_DO(list, type)    do
# define ROW_CHAIN_WHILE(list, type) while (false)
# define ROW_CHAIN_DO_REV(list, type)    do
# define ROW_CHAIN_WHILE_REV(list, type) while (false)
#endif

/* >>> */

/* <<< COL_CHAIN...() */

#define COL_CHAIN_INIT(elem)  DBL_LIST_INIT(&(elem)->col_chain)
#define COL_CHAIN_INITIALIZER(name) DBL_LIST_INITIALIZER((name).col_chain)
#define COL_CHAIN_LENGTH(head) __chain_length(&(head)->col_chain)
#define COL_CHAIN_SINGLE(var)			\
  (!CHAINED_BASIS_FUNCTIONS || dbl_list_empty(&(var)->col_chain))
#define COL_CHAIN_NEXT(var, type) \
  dbl_list_entry((var)->col_chain.next, type, col_chain)
#define COL_CHAIN_PREV(var, type) \
  dbl_list_entry((var)->col_chain.prev, type, col_chain)
#define COL_CHAIN_ADD_HEAD(head, elem)				\
  dbl_list_add_head(&(head)->col_chain, &(elem)->col_chain)
#define COL_CHAIN_ADD_TAIL(head, elem)				\
  dbl_list_add_tail(&(head)->col_chain, &(elem)->col_chain)
#define COL_CHAIN_DEL(elem)						\
  dbl_list_del(&(elem)->col_chain)
#define COL_CHAIN_FOREACH(ptr, head, type)				\
  if (!COL_CHAIN_SINGLE((head)))					\
    dbl_list_for_each_entry((ptr), &(head)->col_chain, type, col_chain)
#define COL_CHAIN_FOREACH_SAFE(ptr, next, head, type)			\
  if (!COL_CHAIN_SINGLE((head)))					\
    dbl_list_for_each_entry_safe(					\
      (ptr), (next), &(head)->col_chain, type, col_chain)
#define COL_CHAIN_FOREACH_REV(ptr, head, type)				\
  if (!COL_CHAIN_SINGLE((head)))					\
    dbl_list_for_each_entry_rev((ptr), &(head)->col_chain, type, col_chain)
#define COL_CHAIN_FOREACH_REV_SAFE(ptr, next, head, type)		\
  if (!COL_CHAIN_SINGLE((head)))					\
    dbl_list_for_each_entry_rev_safe(					\
      (ptr), (next), &(head)->col_chain, type, col_chain)
#if CHAINED_BASIS_FUNCTIONS
# define COL_CHAIN_DO(list, type)		\
  dbl_list_do_cyclic((list), type, col_chain)
# define COL_CHAIN_WHILE(list, type)		\
  dbl_list_while_cyclic((list), type, col_chain)
# define COL_CHAIN_DO_REV(list, type)		\
  dbl_list_do_cyclic_rev((list), type, col_chain)
# define COL_CHAIN_WHILE_REV(list, type)		\
  dbl_list_while_cyclic_rev((list), type, col_chain)
#else
# define COL_CHAIN_DO(list, type)    do
# define COL_CHAIN_WHILE(list, type) while (false)
# define COL_CHAIN_DO_REV(list, type)    do
# define COL_CHAIN_WHILE_REV(list, type) while (false)
#endif

/* >>> */

/* <<< DEP_CHAIN...() */

/* "depth" chain, this is currently only use for quadrature "tensors" */

#define DEP_CHAIN_INIT(elem)  DBL_LIST_INIT(&(elem)->dep_chain)
#define DEP_CHAIN_INITIALIZER(name) DBL_LIST_INITIALIZER((name).dep_chain)
#define DEP_CHAIN_LENGTH(head) __chain_length(&(head)->dep_chain)
#define DEP_CHAIN_SINGLE(var)			\
  (!CHAINED_BASIS_FUNCTIONS || dbl_list_empty(&(var)->dep_chain))
#define DEP_CHAIN_NEXT(var, type) \
  dbl_list_entry((var)->dep_chain.next, type, dep_chain)
#define DEP_CHAIN_PREV(var, type) \
  dbl_list_entry((var)->dep_chain.prev, type, dep_chain)
#define DEP_CHAIN_ADD_HEAD(head, elem)				\
  dbl_list_add_head(&(head)->dep_chain, &(elem)->dep_chain)
#define DEP_CHAIN_ADD_TAIL(head, elem)				\
  dbl_list_add_tail(&(head)->dep_chain, &(elem)->dep_chain)
#define DEP_CHAIN_DEL(elem)						\
  dbl_list_del(&(elem)->dep_chain)
#define DEP_CHAIN_FOREACH(ptr, head, type)				\
  if (!DEP_CHAIN_SINGLE((head)))					\
    dbl_list_for_each_entry((ptr), &(head)->dep_chain, type, dep_chain)
#define DEP_CHAIN_FOREACH_SAFE(ptr, next, head, type)			\
  if (!DEP_CHAIN_SINGLE((head)))					\
    dbl_list_for_each_entry_safe(					\
      (ptr), next, &(head)->dep_chain, type, dep_chain)
#define DEP_CHAIN_FOREACH_REV(ptr, head, type)				\
  if (!DEP_CHAIN_SINGLE((head)))					\
    dbl_list_for_each_entry_rev((ptr), &(head)->dep_chain, type, dep_chain)
#define DEP_CHAIN_FOREACH_REV_SAFE(ptr, next, head, type)		\
  if (!DEP_CHAIN_SINGLE((head)))					\
    dbl_list_for_each_entry_rev_safe(					\
      (ptr), (next), &(head)->dep_chain, type, dep_chain)
#if CHAINED_BASIS_FUNCTIONS
# define DEP_CHAIN_DO(list, type)		\
  dbl_list_do_cyclic((list), type, dep_chain)
# define DEP_CHAIN_WHILE(list, type)		\
  dbl_list_while_cyclic((list), type, dep_chain)
# define DEP_CHAIN_DO_REV(list, type)		\
  dbl_list_do_cyclic_rev((list), type, dep_chain)
# define DEP_CHAIN_WHILE_REV(list, type)		\
  dbl_list_while_cyclic_rev((list), type, dep_chain)
#else
# define DEP_CHAIN_DO(list, type)    do
# define DEP_CHAIN_WHILE(list, type) while (false)
# define DEP_CHAIN_DO_REV(list, type)    do
# define DEP_CHAIN_WHILE_REV(list, type) while (false)
#endif

/* >>> */

/* <<< FOREACH_DOF...() */

#define FOREACH_DOF(fe_space, todo, next)		\
  {							\
    const FE_SPACE *_AI_fe_space = (fe_space);		\
    CHAIN_DO(_AI_fe_space, const FE_SPACE) {		\
      FOR_ALL_DOFS(_AI_fe_space->admin, { todo; });	\
      next;						\
    } CHAIN_WHILE(_AI_fe_space, const FE_SPACE);	\
  }
    
  
#define FOREACH_DOF_DOW(fe_space, todo, todo_cart, next)		\
  {									\
    const FE_SPACE *_AI_fe_space = (fe_space);				\
    CHAIN_DO(_AI_fe_space, const FE_SPACE) {				\
      if (_AI_fe_space->rdim != 1 && _AI_fe_space->bas_fcts->rdim == 1) { \
	FOR_ALL_DOFS(_AI_fe_space->admin, { todo_cart; });		\
      } else {								\
	FOR_ALL_DOFS(_AI_fe_space->admin, { todo; });			\
      }									\
      next;								\
    } CHAIN_WHILE(_AI_fe_space, const FE_SPACE);			\
  }
    

#define FOREACH_FREE_DOF(fe_space, todo, next)		\
  {							\
    const FE_SPACE *_AI_fe_space = (fe_space);		\
    CHAIN_DO(_AI_fe_space, const FE_SPACE) {		\
      FOR_ALL_FREE_DOFS(_AI_fe_space->admin, todo);	\
      next;						\
    } CHAIN_WHILE(_AI_fe_space, const FE_SPACE);	\
  }

#define FOREACH_FREE_DOF_DOW(fe_space, todo, todo_cart, next)		\
  {									\
    const FE_SPACE *_AI_fe_space = (fe_space);				\
    CHAIN_DO(_AI_fe_space, const FE_SPACE) {				\
      if (_AI_fe_space->rdim != 1 && _AI_fe_space->bas_fcts->rdim == 1) { \
	FOR_ALL_FREE_DOFS(_AI_fe_space->admin, todo_cart);		\
      } else {								\
	FOR_ALL_FREE_DOFS(_AI_fe_space->admin, todo);			\
      }									\
      next;								\
    } CHAIN_WHILE(_AI_fe_space, const FE_SPACE);			\
  } 

/* >>> */

/* <<< generate a DOF-vec skeleton for a possibly chained fe-space */

static inline
DOF_REAL_VEC *init_dof_real_vec_skel(DOF_REAL_VEC vecs[],
				     const char *name,
				     const FE_SPACE *fe_space)
{
  DOF_REAL_VEC *head = vecs;
  const FE_SPACE *fe_chain;
  
  memset(head, 0, sizeof(*head));
  head->fe_space = fe_space;
  head->name     = name;
  head->size     = fe_space->admin->size_used;
  head->reserved = 1;
  CHAIN_INIT(head);
  CHAIN_FOREACH(fe_chain, fe_space, const FE_SPACE) {
    memset(++vecs, 0, sizeof(*vecs));
    vecs->fe_space = fe_chain;
    vecs->name     = name;
    vecs->size     = fe_chain->admin->size_used;
    vecs->reserved = 1;
    CHAIN_ADD_TAIL(head, vecs);
  }

  return head;
}

static inline
DOF_REAL_D_VEC *init_dof_real_d_vec_skel(DOF_REAL_D_VEC vecs[],
					 const char *name,
					 const FE_SPACE *fe_space)
{
  DOF_REAL_D_VEC *head = vecs;
  const FE_SPACE *fe_chain;

  memset(head, 0, sizeof(*head));
  head->fe_space = fe_space;
  head->name     = name;
  head->size     = fe_space->admin->size_used;
  head->reserved = DIM_OF_WORLD;
  CHAIN_INIT(head);
  CHAIN_FOREACH(fe_chain, fe_space, const FE_SPACE) {
    memset(++vecs, 0, sizeof(*vecs));
    vecs->fe_space = fe_chain;
    vecs->name     = name;
    vecs->size     = fe_chain->admin->size_used;
    vecs->reserved = DIM_OF_WORLD;
    CHAIN_ADD_TAIL(head, vecs);
  }

  return head;
}

static inline
DOF_REAL_VEC_D *init_dof_real_vec_d_skel(DOF_REAL_VEC_D vecs[],
					 const char *name,
					 const FE_SPACE *fe_space)
{
  DOF_REAL_VEC_D *head = vecs;
  const FE_SPACE *fe_chain;
  
  memset(head, 0, sizeof(*head));
  head->fe_space = fe_space;
  head->name     = name;
  head->size     = fe_space->admin->size_used;
  head->stride   = 
    (fe_space->rdim == fe_space->bas_fcts->rdim) ? 1 : DIM_OF_WORLD;
  CHAIN_INIT(head);
  CHAIN_FOREACH(fe_chain, fe_space, const FE_SPACE) {
    memset(++vecs, 0, sizeof(*vecs));
    vecs->fe_space = fe_chain;
    vecs->name     = name;
    vecs->size     = fe_chain->admin->size_used;
    vecs->stride   =
      (fe_chain->rdim == fe_chain->bas_fcts->rdim) ? 1 : DIM_OF_WORLD;
    CHAIN_ADD_TAIL(head, vecs);
  }

  return head;
}

static inline
DOF_SCHAR_VEC *init_dof_schar_vec_skel(DOF_SCHAR_VEC vecs[],
				       const char *name,
				       const FE_SPACE *fe_space)
{
  DOF_SCHAR_VEC *head = vecs;
  const FE_SPACE *fe_chain;
  
  memset(head, 0, sizeof(*head));
  head->fe_space = fe_space;
  head->name     = name;
  head->size     = fe_space->admin->size_used;
  head->reserved = 1;
  CHAIN_INIT(head);
  CHAIN_FOREACH(fe_chain, fe_space, const FE_SPACE) {
    memset(++vecs, 0, sizeof(*vecs));
    vecs->fe_space = fe_chain;
    vecs->name     = name;
    vecs->size     = fe_chain->admin->size_used;
    vecs->reserved = 1;
    CHAIN_ADD_TAIL(head, vecs);
  }

  return head;
}


static inline
DOF_REAL_VEC *get_dof_real_vec_skel(const char *name,
				    const FE_SPACE *fe_space,
				    SCRATCH_MEM scr)
{
  DOF_REAL_VEC *vecs;
  
  vecs = SCRATCH_MEM_ALLOC(scr, CHAIN_LENGTH(fe_space), DOF_REAL_VEC);

  return init_dof_real_vec_skel(vecs, name, fe_space);
}

static inline
DOF_REAL_D_VEC *get_dof_real_d_vec_skel(const char *name,
					const FE_SPACE *fe_space,
					SCRATCH_MEM scr)
{
  DOF_REAL_D_VEC *vecs;
  
  vecs = SCRATCH_MEM_ALLOC(scr, CHAIN_LENGTH(fe_space), DOF_REAL_D_VEC);

  return init_dof_real_d_vec_skel(vecs, name, fe_space);
}

static inline
DOF_REAL_VEC_D *get_dof_real_vec_d_skel(const char *name,
					const FE_SPACE *fe_space,
					SCRATCH_MEM scr)
{
  DOF_REAL_VEC_D *vecs;
  
  vecs = SCRATCH_MEM_ALLOC(scr, CHAIN_LENGTH(fe_space), DOF_REAL_VEC_D);

  return init_dof_real_vec_d_skel(vecs, name, fe_space);
}

static inline
DOF_SCHAR_VEC *get_dof_schar_vec_skel(const char *name,
				      const FE_SPACE *fe_space,
				      SCRATCH_MEM scr)
{
  DOF_SCHAR_VEC *vecs;
  
  vecs = SCRATCH_MEM_ALLOC(scr, CHAIN_LENGTH(fe_space), DOF_SCHAR_VEC);

  return init_dof_schar_vec_skel(vecs, name, fe_space);
}

/* >>> */

/* <<< distribute a contiguous vector to a DOF_REAL[_D]_VEC[_D] */

/* Given a possibly chained DOF_VEC skeleton, distribute the given
 * contiguous vector to the DOF_VEC skeleton. This is meant to be used
 * inside a matrix-vector multiplication routine. SKEL is a fake-chain
 * of DOF-vectors, e.g. obtained by get_dof_real_vec_skel().
 */

static inline
size_t distribute_to_dof_real_vec_skel(DOF_REAL_VEC *skel, const REAL *data)
{
  size_t length = 0;

  CHAIN_DO(skel, DOF_REAL_VEC) {
    skel->vec   = (REAL *)data;
    skel->size  = skel->fe_space->admin->size_used;
    data       += skel->size;
    length     += skel->size;
  } CHAIN_WHILE(skel, DOF_REAL_VEC);

  return length;
}

static inline
size_t distribute_to_dof_real_d_vec_skel(DOF_REAL_D_VEC *skel, const REAL *ptr)
{
  size_t length = 0;

  REAL_D *data = (REAL_D *)ptr;
  CHAIN_DO(skel, DOF_REAL_VEC) {
    skel->vec   = data;
    skel->size  = skel->fe_space->admin->size_used;
    data       += skel->size;
    length     += skel->size;
  } CHAIN_WHILE(skel, DOF_REAL_D_VEC);

  return length * DIM_OF_WORLD;
}

static inline
size_t distribute_to_dof_real_vec_d_skel(DOF_REAL_VEC_D *skel, const REAL *data)
{
  size_t length = 0, l;

  CHAIN_DO(skel, DOF_REAL_VEC_D) {
    skel->vec   = (REAL *)data;
    skel->size  = skel->fe_space->admin->size_used;
    l           = skel->size * (skel->stride != 1 ? DIM_OF_WORLD : 1);
    data       += l;
    length     += l;
  } CHAIN_WHILE(skel, DOF_REAL_VEC_D);

  return length;
}

static inline
size_t distribute_to_dof_schar_vec_skel(DOF_SCHAR_VEC *skel, const S_CHAR *data)
{
  size_t length = 0;

  CHAIN_DO(skel, DOF_SCHAR_VEC) {
    skel->vec   = (S_CHAR *)data;
    skel->size  = skel->fe_space->admin->size_used;
    data       += skel->size;
    length     += skel->size;
  } CHAIN_WHILE(skel, DOF_SCHAR_VEC);

  return length;
}

/* >>> */

/* <<< copy a contiguous vector to a DOF_REAL[_D]_VEC[_D] */

/* This is meant to be called after some solver has done its job to
 * copy the contiguous result vector back to the DOF-vector chain.
 */

static inline
size_t copy_to_dof_real_vec(DOF_REAL_VEC *vecs, const REAL *data)
{
  size_t length = 0;
  
  CHAIN_DO(vecs, DOF_REAL_VEC) {
    DOF size_used = vecs->fe_space->admin->size_used;
    memcpy(vecs->vec, data, size_used * sizeof(REAL));
    data   += size_used;
    length += size_used;
  } CHAIN_WHILE(vecs, DOF_REAL_VEC);

  return length;
}

static inline
size_t copy_to_dof_real_d_vec(DOF_REAL_D_VEC *vecs, const REAL *_data)
{
  size_t length = 0;
  const REAL_D *data = (const REAL_D *)_data;

  CHAIN_DO(vecs, DOF_REAL_D_VEC) {
    DOF size_used = vecs->fe_space->admin->size_used;
    memcpy(vecs->vec, data, size_used * sizeof(REAL_D));
    data   += size_used;
    length += size_used;
  } CHAIN_WHILE(vecs, DOF_REAL_D_VEC);

  return length;
}

static inline
size_t copy_to_dof_real_vec_d(DOF_REAL_VEC_D *vecs, const REAL *data)
{
  size_t length = 0;

  CHAIN_DO(vecs, DOF_REAL_VEC_D) {
    DOF size_used = vecs->fe_space->admin->size_used;
    if (vecs->stride != 1) {
      size_used *= DIM_OF_WORLD;
    }
    memcpy(vecs->vec, data, size_used * sizeof(REAL));
    data   += size_used;
    length += size_used;
  } CHAIN_WHILE(vecs, DOF_REAL_VEC_D);

  return length;
}

static inline
size_t copy_to_dof_schar_vec(DOF_SCHAR_VEC *vecs, const S_CHAR *data)
{
  size_t length = 0;

  CHAIN_DO(vecs, DOF_SCHAR_VEC) {
    DOF size_used = vecs->fe_space->admin->size_used;
    memcpy(vecs->vec, data, size_used * sizeof(S_CHAR));
    data   += size_used;
    length += size_used;
  } CHAIN_WHILE(vecs, DOF_SCHAR_VEC);

  return length;
}

/* >>> */

/* <<< copy a DOF_REAL[_D]_VEC[_D] to a contiguous vector */

/* This is meant to be called before some solver expecting a
 * contiguous vector tries to do its job. This takes the data of a
 * DOF-vector chain and copies it to the given contiguous memory
 * portion.
 */

static inline
size_t copy_from_dof_real_vec(REAL *data, const DOF_REAL_VEC *vecs)
{
  size_t length = 0;

  CHAIN_DO(vecs, DOF_REAL_VEC) {
    DOF size_used = vecs->fe_space->admin->size_used;
    memcpy(data, vecs->vec, size_used * sizeof(REAL));
    FOR_ALL_USED_FREE_DOFS(vecs->fe_space->admin, data[dof] = 0.0);
    data   += size_used;
    length += size_used;
  } CHAIN_WHILE(vecs, DOF_REAL_VEC);

  return length;
}

static inline
size_t copy_from_dof_real_d_vec(REAL_D *data, const DOF_REAL_D_VEC *vecs)
{
  size_t length = 0;

  CHAIN_DO(vecs, DOF_REAL_D_VEC) {
    DOF size_used = vecs->fe_space->admin->size_used;
    memcpy(data, vecs->vec, size_used * sizeof(REAL_D));
    FOR_ALL_USED_FREE_DOFS(vecs->fe_space->admin, SET_DOW(0.0, data[dof]));
    data   += size_used;
    length += size_used;
  } CHAIN_WHILE(vecs, DOF_REAL_D_VEC);

  return length;
}

static inline
size_t copy_from_dof_real_vec_d(REAL *data, const DOF_REAL_VEC_D *vecs)
{
  size_t length = 0;

  CHAIN_DO(vecs, DOF_REAL_VEC_D) {
    DOF size_used = vecs->fe_space->admin->size_used;
    if (vecs->stride != 1) {
      size_used *= DIM_OF_WORLD;
      memcpy(data, vecs->vec, size_used * sizeof(REAL));
      FOR_ALL_USED_FREE_DOFS(vecs->fe_space->admin,
			     SET_DOW(0.0, ((REAL_D *)data)[dof]));
    } else {
      memcpy(data, vecs->vec, size_used * sizeof(REAL));
      FOR_ALL_USED_FREE_DOFS(vecs->fe_space->admin, data[dof] = 0.0);
    }
    data   += size_used;
    length += size_used;
  } CHAIN_WHILE(vecs, DOF_REAL_VEC_D);

  return length;
}

static inline
size_t copy_from_dof_schar_vec(S_CHAR *data, const DOF_SCHAR_VEC *vecs)
{
  size_t length = 0;

  CHAIN_DO(vecs, DOF_SCHAR_VEC) {
    DOF size_used = vecs->fe_space->admin->size_used;
    memcpy(data, vecs->vec, size_used * sizeof(S_CHAR));
    FOR_ALL_USED_FREE_DOFS(vecs->fe_space->admin, data[dof] = 0);
    data   += size_used;
    length += size_used;
  } CHAIN_WHILE(vecs, DOF_SCHAR_VEC);

  return length;
}

/* >>> */

/* <<< Length (in REAL's) of a vector */

/* Compute the size of the given DOF_REAL[_D]_VEC[_D] chain in terms
 * of REAL components.
 */
static inline size_t dof_real_vec_d_length(const FE_SPACE *fe_space)
{
  size_t len = 0;

  CHAIN_DO(fe_space, const FE_SPACE) {
    if (fe_space->bas_fcts->rdim == 1 && fe_space->rdim != 1) {
      len += fe_space->admin->size_used * DIM_OF_WORLD;
    } else {
      len += fe_space->admin->size_used;
    }
  } CHAIN_WHILE(fe_space, const FE_SPACE);

  return len;
}

static inline size_t dof_real_d_vec_length(const FE_SPACE *fe_space)
{
  return dof_real_vec_d_length(fe_space);
}

static inline size_t dof_real_vec_length(const FE_SPACE *fe_space)
{
  return dof_real_vec_d_length(fe_space);
}

/* >>> */

/* <<< refine_inter hooks etc. */

/* Install the standard hooks from the basis function implementations
 * into each component of a possibly chained discrete function.
 */
/* <<< refine_inter */

static inline void set_refine_inter(DOF_REAL_VEC *uh)
{
  CHAIN_DO(uh,DOF_REAL_VEC) {
    uh->refine_interpol = uh->fe_space->bas_fcts->real_refine_inter;
  } CHAIN_WHILE(uh, DOF_REAL_VEC);
}

static inline void set_refine_inter_d(DOF_REAL_D_VEC *uh)
{
  CHAIN_DO(uh,DOF_REAL_D_VEC) {
    uh->refine_interpol = uh->fe_space->bas_fcts->real_d_refine_inter;
  } CHAIN_WHILE(uh, DOF_REAL_D_VEC);
}

static inline void set_refine_inter_dow(DOF_REAL_VEC_D *uh)
{
  CHAIN_DO(uh,DOF_REAL_VEC_D) {
    uh->refine_interpol = uh->fe_space->bas_fcts->real_refine_inter_d;
  } CHAIN_WHILE(uh, DOF_REAL_VEC_D);
}

/* >>> */

/* <<< coarse_inter */

static inline void set_coarse_inter(DOF_REAL_VEC *uh)
{
  CHAIN_DO(uh,DOF_REAL_VEC) {
    uh->coarse_restrict = uh->fe_space->bas_fcts->real_coarse_inter;
  } CHAIN_WHILE(uh, DOF_REAL_VEC);
}

static inline void set_coarse_inter_d(DOF_REAL_D_VEC *uh)
{
  CHAIN_DO(uh,DOF_REAL_D_VEC) {
    uh->coarse_restrict = uh->fe_space->bas_fcts->real_d_coarse_inter;
  } CHAIN_WHILE(uh, DOF_REAL_D_VEC);
}

static inline void set_coarse_inter_dow(DOF_REAL_VEC_D *uh)
{
  CHAIN_DO(uh,DOF_REAL_VEC_D) {
    uh->coarse_restrict = uh->fe_space->bas_fcts->real_coarse_inter_d;
  } CHAIN_WHILE(uh, DOF_REAL_VEC_D);
}

/* >>> */

/* <<< coarse_restrict */

static inline void set_coarse_restrict(DOF_REAL_VEC *uh)
{
  CHAIN_DO(uh,DOF_REAL_VEC) {
    uh->coarse_restrict = uh->fe_space->bas_fcts->real_coarse_restr;
  } CHAIN_WHILE(uh, DOF_REAL_VEC);
}

static inline void set_coarse_restrict_d(DOF_REAL_D_VEC *uh)
{
  CHAIN_DO(uh,DOF_REAL_D_VEC) {
    uh->coarse_restrict = uh->fe_space->bas_fcts->real_d_coarse_restr;
  } CHAIN_WHILE(uh, DOF_REAL_D_VEC);
}

static inline void set_coarse_restrict_dow(DOF_REAL_VEC_D *uh)
{
  CHAIN_DO(uh,DOF_REAL_VEC_D) {
    uh->coarse_restrict = uh->fe_space->bas_fcts->real_coarse_restr_d;
  } CHAIN_WHILE(uh, DOF_REAL_VEC_D);
}

/* >>> */

/* >>> */

/* <<< sub-chains */

/* A sub-chain is a shallow-copy of the real objects (shallow means,
 * e.g., that DOF_REAL_VEC::vec is not cloned, but just the pointer is
 * installed into the copy).
 *
 * A sub-chain has to be updated explicitly when the backing "real"
 * object changes (e.g. because of mesh-refinement etc., matrix
 * update).
 */

/* <<< clone basis functions */

static inline
BAS_FCTS *bas_fcts_sub_chain(SCRATCH_MEM scr,
			     const BAS_FCTS *bfcts, FLAGS which)
{
  BAS_FCTS *first = NULL, *chain;

  CHAIN_DO(bfcts, const BAS_FCTS) {
    if ((which & 1)) {
      chain = SCRATCH_MEM_ALLOC(scr, 1, BAS_FCTS);
      *chain = *bfcts;
      if (first == NULL) {
	first = chain;
	CHAIN_INIT(first);
      } else {
	CHAIN_ADD_TAIL(first, chain);
      }
    }
    which >>= 1;
  } CHAIN_WHILE(bfcts, const BAS_FCTS);

  return first;
}

/* >>> */

/* <<< clone fe-spaces */

static inline
FE_SPACE *fe_space_sub_chain(SCRATCH_MEM scr,
			     const FE_SPACE *fe_space, FLAGS which)
{
  FE_SPACE *first = NULL, *chain;
  const BAS_FCTS *bfcts;

  bfcts = bas_fcts_sub_chain(scr, fe_space->bas_fcts, which);
  
  CHAIN_DO(fe_space, const FE_SPACE) {
    if ((which & 1)) {
      chain = SCRATCH_MEM_ALLOC(scr, 1, FE_SPACE);
      *chain = *fe_space;
      chain->bas_fcts = bfcts;
      if (first == NULL) {
	first = chain;
	CHAIN_INIT(first);
      } else {
	CHAIN_ADD_TAIL(first, chain);
      }
      bfcts = CHAIN_NEXT(bfcts, const BAS_FCTS);
    }
    which >>= 1;
  } CHAIN_WHILE(fe_space, const FE_SPACE);

  return first;
}

/* >>> */

/* <<< sub-chains for DOF_..._VECs */

#define DEFUN_DOF_VEC_SUB_CHAIN(TYPE, type)			\
  static inline							\
  DOF_##TYPE *dof_##type##_sub_chain(				\
    SCRATCH_MEM scr, const DOF_##TYPE *vec, FLAGS which)	\
  {								\
    DOF_##TYPE *first = NULL, *chain;				\
    const FE_SPACE *fecp;					\
								\
    fecp = fe_space_sub_chain(scr, vec->fe_space, which);	\
								\
    CHAIN_DO(vec, const DOF_##TYPE) {				\
      if ((which & 1)) {					\
	chain = SCRATCH_MEM_ALLOC(scr, 1, DOF_##TYPE);		\
	*chain = *vec;						\
	chain->unchained = vec;					\
	chain->fe_space = fecp;					\
	chain->mem_info = NULL;					\
	chain->next = NULL;					\
	if (first == NULL) {					\
	  first = chain;					\
	  CHAIN_INIT(first);					\
	} else {						\
	  CHAIN_ADD_TAIL(first, chain);				\
	}							\
	fecp = CHAIN_NEXT(fecp, const FE_SPACE);		\
      }								\
      which >>= 1;						\
    } CHAIN_WHILE(vec, const DOF_##TYPE);			\
    return first;						\
  }								\
  struct _AI_semicolon_dummy

/* Shallow update, i.e. make sure the sub_vec->vec points to
 * sub_vec->vec->unchained, sub_vec->fe_space is _NOT_ updated, nor
 * the underlying basis functions.
 */
#define DEFUN_UPDATE_DOF_VEC_SUB_CHAIN(TYPE, type)		\
  static inline							\
  void update_dof_##type##_sub_chain(const DOF_##TYPE *sub_vec)	\
  {								\
    CHAIN_DO(sub_vec, const DOF_REAL_VEC) {			\
      ((DOF_##TYPE *)sub_vec)->vec  = sub_vec->unchained->vec;	\
      ((DOF_##TYPE *)sub_vec)->size = sub_vec->unchained->size;	\
    } CHAIN_WHILE(sub_vec, const DOF_##TYPE);			\
  }								\
  struct _AI_semicolon_dummy

DEFUN_DOF_VEC_SUB_CHAIN(REAL_VEC, real_vec);
DEFUN_UPDATE_DOF_VEC_SUB_CHAIN(REAL_VEC, real_vec);

DEFUN_DOF_VEC_SUB_CHAIN(REAL_D_VEC, real_d_vec);
DEFUN_UPDATE_DOF_VEC_SUB_CHAIN(REAL_D_VEC, real_d_vec);

DEFUN_DOF_VEC_SUB_CHAIN(REAL_VEC_D, real_vec_d);
DEFUN_UPDATE_DOF_VEC_SUB_CHAIN(REAL_VEC_D, real_vec_d);

DEFUN_DOF_VEC_SUB_CHAIN(DOF_VEC, dof_vec);
DEFUN_UPDATE_DOF_VEC_SUB_CHAIN(DOF_VEC, dof_vec);

DEFUN_DOF_VEC_SUB_CHAIN(INT_VEC, int_vec);
DEFUN_UPDATE_DOF_VEC_SUB_CHAIN(INT_VEC, int_vec);

DEFUN_DOF_VEC_SUB_CHAIN(SCHAR_VEC, schar_vec);
DEFUN_UPDATE_DOF_VEC_SUB_CHAIN(SCHAR_VEC, schar_vec);

DEFUN_DOF_VEC_SUB_CHAIN(UCHAR_VEC, uchar_vec);
DEFUN_UPDATE_DOF_VEC_SUB_CHAIN(UCHAR_VEC, uchar_vec);

DEFUN_DOF_VEC_SUB_CHAIN(PTR_VEC, ptr_vec);
DEFUN_UPDATE_DOF_VEC_SUB_CHAIN(PTR_VEC, ptr_vec);

/* >>> */

/* <<< sub-chains for DOF_MATRIXes */

/* helper-function: flat dummy clone of a DOF-matrix */
static inline
void _AI_unchain_dof_matrix(DOF_MATRIX *to, const DOF_MATRIX *from,
			    const FE_SPACE *row_fe_space,
			    const FE_SPACE *col_fe_space)
{
  *to = *from;
  to->unchained = from;
  COL_CHAIN_INIT(to);
  ROW_CHAIN_INIT(to);
  to->row_fe_space = row_fe_space;
  to->col_fe_space = col_fe_space;
  to->mem_info = NULL;
  to->next = NULL;
}

/* a helper-function for extracting one row */
static inline
DOF_MATRIX *_AI_dm_row_sub_chain(SCRATCH_MEM scr,
				 const DOF_MATRIX *A, FLAGS col_which,
				 const FE_SPACE *row_fe_space,
				 const FE_SPACE *col_fe_space)
{
  DOF_MATRIX *first = NULL, *chain;

  ROW_CHAIN_DO(A, DOF_MATRIX) {
    if ((col_which & 1)) {
      chain = SCRATCH_MEM_ALLOC(scr, 1, DOF_MATRIX);
      _AI_unchain_dof_matrix(chain, A, row_fe_space, col_fe_space);
      if (first != NULL) {
	ROW_CHAIN_ADD_TAIL(first, chain);
      } else {
	first = chain;
      }
      col_fe_space = CHAIN_NEXT(col_fe_space, const FE_SPACE);
    }
    col_which >>= 1;
  } ROW_CHAIN_WHILE(A, DOF_MATRIX);

  return first;
}

static inline
DOF_MATRIX *dof_matrix_sub_chain(SCRATCH_MEM scr,
				 const DOF_MATRIX *A,
				 FLAGS row_which, FLAGS col_which)
{
  DOF_MATRIX *first = NULL, *chain;
  const FE_SPACE *row_fe_space;
  const FE_SPACE *col_fe_space;
  
  row_fe_space = fe_space_sub_chain(scr, A->row_fe_space, row_which);
  if (A->row_fe_space != A->col_fe_space || row_which != col_which) {
    col_fe_space = fe_space_sub_chain(scr, A->col_fe_space, col_which);
  } else {
    col_fe_space = row_fe_space;
  }

  COL_CHAIN_DO(A, DOF_MATRIX) {
    if (row_which & 1) {
      chain =
	_AI_dm_row_sub_chain(scr, A, col_which, row_fe_space, col_fe_space);
      if (first != NULL) {
	/* concatenate with first row */
	ROW_CHAIN_DO(first, DOF_MATRIX) {
	  COL_CHAIN_ADD_TAIL(first, chain);
	  chain = ROW_CHAIN_NEXT(chain, DOF_MATRIX);
	} ROW_CHAIN_WHILE(first, DOF_MATRIX);
      } else {
	first = chain;
      }
      row_fe_space = CHAIN_NEXT(row_fe_space, const FE_SPACE);
    }
    row_which >>= 1;
  } COL_CHAIN_WHILE(A, const DOF_MATRIX);

  return first;
}

/* Update a dummy sub-matrix chain, i.e. see that subM->matrix_row and
 * subM->size are up to date. NOTE: subM->{row,col}_fe_space is NOT
 * updated.
 */
static inline
void update_dof_matrix_sub_chain(DOF_MATRIX *subM)
{
  COL_CHAIN_DO(subM, DOF_MATRIX) {
    ROW_CHAIN_DO(subM, DOF_MATRIX) {
      subM->matrix_row = subM->unchained->matrix_row;
      subM->size       = subM->unchained->size;      
    } ROW_CHAIN_WHILE(subM, DOF_MATRIX);
  } COL_CHAIN_WHILE(subM, DOF_MATRIX);
}

/* >>> */

/* >>> */

#endif /* _ALBERTA_DOF_CHAINS_H_ */