tree.h

Go to the documentation of this file.

/** @file tree.h
 */
typedef double real;
 
#include "mempool.h"
 
#define TWO_ONE 1 // enforce 2:1 refinement ratio
#define GHOSTS  2
 
struct _Point {
  /* the current cell index and level */
  int i;
#if dimension >= 2
  int j;
#endif
#if dimension >= 3
  int k;
#endif
  int level;
#if LAYERS
  int l;
  @define _BLOCK_INDEX , point.l
#else
  @define _BLOCK_INDEX
#endif
};
 
//#include "memindex/range.h"
#include "memindex/virtual.h"
 
#if LAYERS
# include "grid/layers.h"
#endif
 
/* By default only one layer of ghost cells is used on the boundary to
   optimise the cost of boundary conditions. */
 
#ifndef BGHOSTS
@ define BGHOSTS 1
#endif
 
#define _I     (point.i - GHOSTS)
#if dimension >= 2
# define _J    (point.j - GHOSTS)
#endif
#if dimension >= 3
# define _K    (point.k - GHOSTS)
#endif
#define _DELTA (1./(1 << point.level))
 
typedef struct {
  unsigned short flags;
  // number of refined neighbors in a 3^dimension neighborhood
  unsigned short neighbors;
  int pid; // process id
} Cell;
 
enum {
  active    = 1 << 0,
  leaf      = 1 << 1,
  border    = 1 << 2,
  vertex    = 1 << 3,
  user      = 4,
 
  face_x = 1 << 0
#if dimension >= 2
  , face_y = 1 << 1
#endif
#if dimension >= 3
  , face_z = 1 << 2
#endif
};
 
@define is_active(cell)  ((cell).flags & active)
@define is_leaf(cell)    ((cell).flags & leaf)
@define is_coarse()      ((cell).neighbors > 0)
@define is_border(cell)  ((cell).flags & border)
@define is_local(cell)   ((cell).pid == pid())
@define is_vertex(cell)  ((cell).flags & vertex)
 
// Caches
 
typedef struct {
  int i;
#if dimension >= 2
  int j;
#endif
#if dimension >= 3
  int k;
#endif  
} IndexLevel;
 
typedef struct {
  IndexLevel * p;
  int n, nm;
} CacheLevel;
 
typedef struct {
  int i;
#if dimension >= 2
  int j;
#endif
#if dimension >= 3
  int k;
#endif  
  int level, flags;
} Index;
 
typedef struct {
  Index * p;
  int n, nm;
} Cache;
 
// Layer
 
typedef struct {
  Memindex m; // the structure indexing the data
  Mempool * pool; // the memory pool actually holding the data
  long nc;     // the number of allocated elements
  int len;    // the (1D) size of the array
} Layer;
 
static size_t _size (size_t depth)
{
  return (1 << depth) + 2*GHOSTS;
}
 
static size_t poolsize (size_t depth, size_t size)
{
  // the maximum amount of data at a given level
#if dimension == 1
  return _size(depth)*size;  
#elif dimension == 2
  return sq(_size(depth))*size;
#else
  return cube(_size(depth))*size;
#endif
}
 
static Layer * new_layer (int depth)
{
  Layer * l = qmalloc (1, Layer);
  l->len = _size (depth);
  if (depth == 0)
    l->pool = NULL; // the root layer does not use a pool
  else {
    size_t size = sizeof(Cell) + datasize;
    // the block size is 2^dimension*size because we allocate
    // 2^dimension children at a time
    l->pool = mempool_new (poolsize (depth, size), (1 << dimension)*size);
  }
  l->m = mem_new (l->len);
  l->nc = 0;
  return l;
}
 
static void destroy_layer (Layer * l)
{
  if (l->pool)
    mempool_destroy (l->pool);
  mem_destroy (l->m, l->len);
  free (l);
}
 
// Tree
 
typedef struct {
  Grid g;
  Layer ** L; /* the grids at each level */
 
  Cache        leaves;   /* leaf indices */
  Cache        faces;    /* face indices */
  Cache        vertices; /* vertex indices */
  Cache        refined;  /* refined cells */
  CacheLevel * active;   /* active cells indices for each level */
  CacheLevel * prolongation; /* halo prolongation indices for each level */
  CacheLevel * boundary;  /* boundary indices for each level */
  /* indices of boundary cells with non-boundary parents */
  CacheLevel * restriction;
  
  bool dirty;       /* whether caches should be updated */
} Tree;
 
#define tree ((Tree *)grid)
 
static Point last_point;
 
#define BSIZE 128
 
static void cache_level_append (CacheLevel * c, Point p)
{
  if (c->n >= c->nm) {
    c->nm += BSIZE;
    qrealloc (c->p, c->nm, IndexLevel);
  }
  c->p[c->n].i = p.i;
#if dimension >= 2
  c->p[c->n].j = p.j;
#endif
#if dimension >= 3
  c->p[c->n].k = p.k;
#endif
  c->n++;
}
 
static void cache_level_shrink (CacheLevel * c)
{
  if (c->nm > (c->n/BSIZE + 1)*BSIZE) {
    c->nm = (c->n/BSIZE + 1)*BSIZE;
    assert (c->nm > c->n);
    c->p = (IndexLevel *) realloc (c->p, sizeof (Index)*c->nm);
  }
}
 
static void cache_append (Cache * c, Point p, unsigned short flags)
{
  if (c->n >= c->nm) {
    c->nm += BSIZE;
    qrealloc (c->p, c->nm, Index);
  }
  c->p[c->n].i = p.i;
#if dimension >= 2
  c->p[c->n].j = p.j;
#endif
#if dimension >= 3
  c->p[c->n].k = p.k;
#endif  
  c->p[c->n].level = p.level;
  c->p[c->n].flags = flags;
  c->n++;
}
 
void cache_shrink (Cache * c)
{
  cache_level_shrink ((CacheLevel *)c);
}
 
#undef BSIZE
 
/* low-level memory management */
#if dimension == 1
@define allocated(k,l,n) (mem_allocated (tree->L[point.level]->m, point.i+k))
@define NEIGHBOR(k,l,n) (mem_data (tree->L[point.level]->m, point.i+k))
@def PARENT(k,l,n) (mem_data (tree->L[point.level-1]->m,
            (point.i+GHOSTS)/2+k))
@
@def allocated_child(k,l,n) (level < depth() &&
           mem_allocated(tree->L[point.level+1]->m,
             2*point.i-GHOSTS+k))
@
@define CHILD(k,l,n)  (mem_data (tree->L[point.level+1]->m, 2*point.i-GHOSTS+k))
#elif dimension == 2
@def allocated(k,l,n) (mem_allocated (tree->L[point.level]->m,
              point.i+k, point.j+l))
@
@def NEIGHBOR(k,l,n)  (mem_data (tree->L[point.level]->m,
           point.i+k, point.j+l))
@
@def PARENT(k,l,n) (mem_data (tree->L[point.level-1]->m,
            (point.i+GHOSTS)/2+k, (point.j+GHOSTS)/2+l))
@
@def allocated_child(k,l,n)  (level < depth() &&
            mem_allocated (tree->L[point.level+1]->m,
               2*point.i-GHOSTS+k,
               2*point.j-GHOSTS+l))
@
@def CHILD(k,l,n)  (mem_data (tree->L[point.level+1]->m,
            2*point.i-GHOSTS+k, 2*point.j-GHOSTS+l))
@
#else // dimension == 3
@def allocated(a,l,n) (mem_allocated (tree->L[point.level]->m,
              point.i+a, point.j+l, point.k+n))
@
@def NEIGHBOR(a,l,n)  (mem_data (tree->L[point.level]->m,
           point.i+a, point.j+l, point.k+n))
@     
@def PARENT(a,l,n) (mem_data (tree->L[point.level-1]->m,
            (point.i+GHOSTS)/2+a,
            (point.j+GHOSTS)/2+l,
            (point.k+GHOSTS)/2+n))
@
@def allocated_child(a,l,n)  (level < depth() &&
            mem_allocated (tree->L[point.level+1]->m,
               2*point.i-GHOSTS+a,
               2*point.j-GHOSTS+l,
               2*point.k-GHOSTS+n))
@
@def CHILD(a,l,n)  (mem_data (tree->L[point.level+1]->m,
            2*point.i-GHOSTS+a,
            2*point.j-GHOSTS+l,
            2*point.k-GHOSTS+n))
@
#endif // dimension == 3
@define CELL(m) (*((Cell *)(m)))
 
/***** Multigrid macros *****/
@define depth()        (grid->depth)
@define aparent(k,l,n) CELL(PARENT(k,l,n))
@define child(k,l,n)   CELL(CHILD(k,l,n))
 
/***** Tree macros ****/
@define cell     CELL(NEIGHBOR(0,0,0))
@define neighbor(k,l,n)  CELL(NEIGHBOR(k,l,n))
@def neighborp(l,m,n) (Point) {
    point.i + l,
#if dimension >= 2
    point.j + m,
#endif
#if dimension >= 3
    point.k + n,
#endif
    point.level
    _BLOCK_INDEX
}
@
 
/***** Data macros *****/
@define data(k,l,n)     ((double *) (NEIGHBOR(k,l,n) + sizeof(Cell)))
@define fine(a,k,p,n)   ((double *) (CHILD(k,p,n) + sizeof(Cell)))[_index(a,n)]
@define coarse(a,k,p,n) ((double *) (PARENT(k,p,n) + sizeof(Cell)))[_index(a,n)]
 
macro POINT_VARIABLES (Point point = point) {
  VARIABLES();
  int level = point.level; NOT_UNUSED(level);
#if dimension == 1
  struct { int x; } child = { 2*((point.i+GHOSTS)%2)-1 };
#elif dimension == 2
  struct { int x, y; } child = {
    2*((point.i+GHOSTS)%2)-1, 2*((point.j+GHOSTS)%2)-1
  };
#else
  struct { int x, y, z; } child = {
    2*((point.i+GHOSTS)%2)-1, 2*((point.j+GHOSTS)%2)-1, 2*((point.k+GHOSTS)%2)-1
  };
#endif
  NOT_UNUSED(child);
  Point parent = point; NOT_UNUSED(parent);
  parent.level--;
  parent.i = (point.i + GHOSTS)/2;
#if dimension >= 2
  parent.j = (point.j + GHOSTS)/2;
#endif
#if dimension >= 3
  parent.k = (point.k + GHOSTS)/2;
#endif
#if TRASH
  Cell * cellp = point.level <= depth() && allocated(0,0,0) ?
    (Cell *) NEIGHBOR(0,0,0) : NULL;
  NOT_UNUSED(cellp);
#endif
}
 
#include "foreach_cell.h"
 
#if dimension == 1
macro1 foreach_child (Point point = point, break = (_k = 2)) {
  {
    int _i = 2*point.i - GHOSTS;
    point.level++;
    for (int _k = 0; _k < 2; _k++) {
      point.i = _i + _k;
      POINT_VARIABLES();
      {...}
    }
    point.i = (_i + GHOSTS)/2;
    point.level--;
  }
}
#elif dimension == 2
macro1 foreach_child (Point point = point, break = (_k = _l = 2)) {
  {
    int _i = 2*point.i - GHOSTS, _j = 2*point.j - GHOSTS;
    point.level++;
    for (int _k = 0; _k < 2; _k++) {
      point.i = _i + _k;
      for (int _l = 0; _l < 2; _l++) {
  point.j = _j + _l;
  POINT_VARIABLES();
  {...}
      }
    }
    point.i = (_i + GHOSTS)/2; point.j = (_j + GHOSTS)/2;
    point.level--;
  }
}
#else // dimension == 3
macro1 foreach_child (Point point = point, break = (_l = _m = _n = 2)) {
  {
    int _i = 2*point.i - GHOSTS, _j = 2*point.j - GHOSTS, _k = 2*point.k - GHOSTS;
    point.level++;
    for (int _l = 0; _l < 2; _l++) {
      point.i = _i + _l;
      for (int _m = 0; _m < 2; _m++) {
  point.j = _j + _m;
  for (int _n = 0; _n < 2; _n++) {
    point.k = _k + _n;
    POINT_VARIABLES();
    {...}
  }
      }
    }
    point.i = (_i + GHOSTS)/2;point.j = (_j + GHOSTS)/2;point.k = (_k + GHOSTS)/2;
    point.level--;
  }
}
#endif // dimension == 3
  
#define update_cache() { if (tree->dirty) update_cache_f(); }
 
#define is_refined(cell)      (!is_leaf (cell) && cell.neighbors && cell.pid >= 0)
#define is_prolongation(cell) (!is_leaf(cell) && !cell.neighbors && cell.pid >= 0)
#define is_boundary(cell)     (cell.pid < 0)
  
@def is_refined_check() (is_refined(cell) &&
       point.i > 0 && point.i < (1 << level) + 2*GHOSTS - 1
#if dimension > 1
       && point.j > 0 && point.j < (1 << level) + 2*GHOSTS - 1
#endif
#if dimension > 2
       && point.k > 0 && point.k < (1 << level) + 2*GHOSTS - 1
#endif
       )
@
 
macro2 for (int _i = 0; _i < 1; _i++) /* cache Cache cache, Reduce reductions = None */
{
  OMP_PARALLEL (reductions) {
    int ig = 0, jg = 0, kg = 0; NOT_UNUSED(ig); NOT_UNUSED(jg); NOT_UNUSED(kg);
    Point point = {0}; NOT_UNUSED (point);
    point.i = GHOSTS;
#if dimension > 1
    point.j = GHOSTS;
#endif
#if dimension > 2
    point.k = GHOSTS;
#endif
    int _k; unsigned short _flags; NOT_UNUSED(_flags);
    OMP(omp for schedule(static))
      for (_k = 0; _k < cache.n; _k++) {
  point.i = cache.p[_k].i;
#if dimension >= 2
  point.j = cache.p[_k].j;
#endif
#if dimension >= 3
  point.k = cache.p[_k].k;
#endif
  point.level = cache.p[_k].level;
  _flags = cache.p[_k].flags;
  {...}
      }
  }
}
 
macro2 foreach_cache_level (Cache cache, int _l, Reduce reductions = None)
{
  OMP_PARALLEL (reductions) {
    int ig = 0, jg = 0, kg = 0; NOT_UNUSED(ig); NOT_UNUSED(jg); NOT_UNUSED(kg);
    Point point = {0}; NOT_UNUSED (point);
    point.i = GHOSTS;
#if dimension > 1
    point.j = GHOSTS;
#endif
#if dimension > 2
    point.k = GHOSTS;
#endif
    point.level = _l;
    int _k;
    OMP(omp for schedule(static))
      for (_k = 0; _k < cache.n; _k++) {
  point.i = cache.p[_k].i;
#if dimension >= 2
  point.j = cache.p[_k].j;
#endif
#if dimension >= 3
  point.k = cache.p[_k].k;
#endif
  {...}
      }
  }
}
 
static void update_cache_f (void);
 
macro2 foreach_boundary_level (int _l, Reduce reductions = None)
{
  if (_l <= depth()) {
    update_cache();
    CacheLevel _boundary = tree->boundary[_l];
    foreach_cache_level (_boundary, _l, reductions)
      {...}
  }
}
 
#define bid(cell) (- cell.pid - 1)
 
macro2 foreach_halo_prolongation (int _l)
{
  if (_l <= depth()) {
    update_cache();
    CacheLevel _cache = tree->prolongation[_l];
    foreach_cache_level (_cache, _l)
      {...}
  }
}
 
macro2 foreach_halo_restriction (int _l)
{
  if (_l <= depth()) {
    update_cache();
    CacheLevel _cache = tree->restriction[_l];
    foreach_cache_level (_cache, _l)
      {...}
  }
}
 
#define foreach_halo(type, l) foreach_halo_##type(l)
 
#include "neighbors.h"
 
static inline bool has_local_children (Point point)
{
  for (int _c = 0; _c < 4; _c++) /* foreach_child */
    if (is_local(cell))
      return true;
  return false;
}
 
static inline void cache_append_face (Point point, unsigned short flags)
{
  Tree * q = tree;
  cache_append (&q->faces, point, flags);
#if dimension == 2
  if (!is_vertex(cell)) {
    cache_append (&q->vertices, point, 0);
    cell.flags |= vertex;
  }
  for (int _d = 0; _d < dimension; _d++)
    if ((flags & face_y) && !is_vertex(neighbor(1))) {
      cache_append (&q->vertices, neighborp(1), 0);
      neighbor(1).flags |= vertex;
    }
#elif dimension == 3
  for (int _d = 0; _d < dimension; _d++)
    if (flags & face_x)
      for (int i = 0; i <= 1; i++)
  for (int j = 0; j <= 1; j++)
    if (!is_vertex(neighbor(0,i,j))) {
      cache_append (&q->vertices, neighborp(0,i,j), 0);
      neighbor(0,i,j).flags |= vertex;
    }
#endif
}
 
#define FBOUNDARY 1 // fixme: this should work with zero
 
static void update_cache_f (void)
{
  Tree * q = tree;
 
  for (int _i = 0; _i < 1; _i++) /* cache q->vertices */
    if (level <= depth() && allocated(0))
      cell.flags &= ~vertex;
  
  /* empty caches */
  q->leaves.n = q->faces.n = q->vertices.n = 0;
  for (int l = 0; l <= depth(); l++)
    q->active[l].n = q->prolongation[l].n =
      q->boundary[l].n = q->restriction[l].n = 0;
#if FBOUNDARY
  const unsigned short fboundary = 1 << user;
  for (int _i = 0; _i < _N; _i++) /* foreach_cell */ {
#else    
  foreach_cell_all() {
#endif
    if (is_local(cell) && is_active(cell)) {
      // active cells
      //      assert (is_active(cell));
      cache_level_append (&q->active[level], point);
    }
#if !FBOUNDARY
    if (is_boundary(cell)) {
      // boundary conditions
      bool has_neighbors = false;
      for (int _n = 0; _n < BGHOSTS; _n++)
  if (allocated(0) && !is_boundary(cell)) {
    has_neighbors = true; break;
  }
      if (has_neighbors)
  cache_level_append (&q->boundary[level], point);
      // restriction for masked cells
      if (level > 0 && is_local(aparent(0)))
  cache_level_append (&q->restriction[level], point);
    }
#else
    // boundaries
    if (!is_boundary(cell)) {
      // look in a 5x5 neighborhood for boundary cells
      for (int _n = 0; _n < BGHOSTS; _n++)
  if (allocated(0) && is_boundary(cell) && !(cell.flags & fboundary)) {
    cache_level_append (&q->boundary[level], point);
    cell.flags |= fboundary;
  }
    }
    // restriction for masked cells
    else if (level > 0 && is_local(aparent(0)))
      cache_level_append (&q->restriction[level], point);
#endif
    if (is_leaf (cell)) {
      if (is_local(cell)) {
  cache_append (&q->leaves, point, 0);
  // faces
  unsigned short flags = 0;
  for (int _d = 0; _d < dimension; _d++)
    if (is_boundary(neighbor(-1)) || is_prolongation(neighbor(-1)) ||
        is_leaf(neighbor(-1)))
      flags |= face_x;
  if (flags)
    cache_append (&q->faces, point, flags);
  for (int _d = 0; _d < dimension; _d++)
    if (is_boundary(neighbor(1)) || is_prolongation(neighbor(1)) ||
        (!is_local(neighbor(1)) && is_leaf(neighbor(1))))
      cache_append (&q->faces, neighborp(1), face_x);
  // vertices
  for (int i = 0; i <= 1; i++)
        #if dimension >= 2
    for (int j = 0; j <= 1; j++)
        #endif
          #if dimension >= 3
      for (int k = 0; k <= 1; k++)
    #endif
        if (!is_vertex(neighbor(i,j,k))) {
    cache_append (&q->vertices, neighborp(i,j,k), 0);
    neighbor(i,j,k).flags |= vertex;
        }
  // halo prolongation
        if (cell.neighbors > 0)
    cache_level_append (&q->prolongation[level], point);
      }
      else if (!is_boundary(cell) || is_local(aparent(0))) { // non-local
  // faces
  unsigned short flags = 0;
  for (int _d = 0; _d < dimension; _d++)
    if (allocated(-1) &&
        is_local(neighbor(-1)) && is_prolongation(neighbor(-1)))
      flags |= face_x;
  if (flags)
    cache_append_face (point, flags);
  for (int _d = 0; _d < dimension; _d++)
    if (allocated(1) && is_local(neighbor(1)) &&
        is_prolongation(neighbor(1)))
      cache_append_face (neighborp(1), face_x);
      }
#if FBOUNDARY // fixme: this should always be included
      continue; 
#endif
    }
  }
 
  /* optimize caches */
  cache_shrink (&q->leaves);
  cache_shrink (&q->faces);
  cache_shrink (&q->vertices);
  for (int l = 0; l <= depth(); l++) {
    cache_level_shrink (&q->active[l]);
    cache_level_shrink (&q->prolongation[l]);
    cache_level_shrink (&q->boundary[l]);
    cache_level_shrink (&q->restriction[l]);
}
  
  q->dirty = false;
 
#if FBOUNDARY
  for (int l = depth(); l >= 0; l--)
    foreach_boundary_level (l)
      cell.flags &= ~fboundary;
#endif
  
  // mesh size
  grid->n = q->leaves.n;
  // for MPI the reduction operation over all processes is done by balance()
@if !_MPI
  grid->tn = grid->n;
  grid->maxdepth = grid->depth;
@endif
}
 
macro2 foreach (char flags = 0, Reduce reductions = None) {
  update_cache();
  for (int _i = 0; _i < 1; _i++) /* cache tree->leaves, reductions */
    {...}
}
 
macro2 foreach_face_generic (char flags = 0, Reduce reductions = None,
        const char * order = "xyz")
{
  update_cache();
  for (int _i = 0; _i < 1; _i++) /* cache tree->faces, reductions */
    {...}
}
 
macro1 is_face_x (unsigned short _f = _flags) {
  if (_f & face_x) {
    int ig = -1; NOT_UNUSED(ig);
    {...}
  }
}
 
#if dimension >= 2
macro1 is_face_y (unsigned short _f = _flags) {
  if (_f & face_y) {
    int jg = -1; NOT_UNUSED(jg);
    {...}
  }
}
#endif
#if dimension >= 3
macro1 is_face_z (unsigned short _f = _flags) {
  if (_f & face_z) {
    int kg = -1; NOT_UNUSED(kg);
    {...}
  }
}
#endif
 
#if dimension == 3
# define foreach_edge(...)      \
  for (int _i = 0; _i < _N; _i++) /* foreach_vertex */      \
    for (int _d = 0; _d < dimension; _d++)        \
    if (is_vertex(neighbor(1)))
#else // dimension < 3
# define foreach_edge(...) for (int _i = 0; _i < _N; _i++) /* foreach_face */
#endif
 
macro2 for (int _l = 0; _l < int l, char flags = 0, Reduce reductions = None; _l++) {
  if (l <= depth()) {
    update_cache();
    CacheLevel _active = tree->active[l];
    foreach_cache_level (_active, l, reductions)
      {...}
  }
}
 
macro2 for (int _l = 0; _l < int l, char flags = 0, Reduce reductions = None; _l++) {
  for (int _l = 0; _l < l, flags, reductions; _l++)
    if (!is_leaf(cell))
      {...}
}
 
macro2 for (int _i = 0; _i < int l, char flags = 0, Reduce reductions = None; _i++) {
  for (int _l1 = l; _l1 >= 0; _l1--)
    for (int _l = 0; _l < _l1, flags, reductions; _l++)
      if (_l1 == l || is_leaf (cell))
  {...}
}
 
@if TRASH
@ undef trash
@ define trash(list) reset(list, undefined)
@endif
 
void reset (void * alist, double val)
{
  scalar * list = (scalar *) alist;
  Tree * q = tree;
  /* low-level memory management */
  for (int l = 0; l <= depth(); l++) {
    Layer * L = q->L[l];
    foreach_mem (L->m, L->len, 1) {
      point.level = l;
      for (int _s = 0; _s < 1; _s++) /* scalar in list */ {
  if (!is_constant(s))
    for (int b = 0; b < s.block; b++)
      data(0,0,0)[s.i + b] = val;
      }
    }
  }
}
 
static CacheLevel * cache_level_resize (CacheLevel * name, int a)
{
  for (int i = 0; i <= depth() - a; i++)
    free (name[i].p);
  free (name);
  return qcalloc (depth() + 1, CacheLevel);
}
  
static void update_depth (int inc)
{
  Tree * q = tree;
  grid->depth += inc;
  q->L = &(q->L[-1]);
  qrealloc (q->L, grid->depth + 2, Layer *);
  q->L = &(q->L[1]);
  if (inc > 0)
    q->L[grid->depth] = new_layer (grid->depth);
  q->active = cache_level_resize (q->active, inc);
  q->prolongation = cache_level_resize (q->prolongation, inc);
  q->boundary = cache_level_resize (q->boundary, inc);
  q->restriction = cache_level_resize (q->restriction, inc);
}
 
#if dimension == 1
typedef void (* PeriodicFunction) (Memindex, int, int, void *);
 
static void periodic_function (Memindex m, int i, int len, void * b,
             PeriodicFunction f)
{
  f(m, i, len, b);
  if (Period.x) {
    int nl = len - 2*GHOSTS;
    for (int l = - 1; l <= 1; l += 2)
      for (int n = i + l*nl; n >= 0 && n < len; n += l*nl)
  f(m, n, len, b);
  }
}
             
static void assign_periodic (Memindex m, int i, int len, void * b)
{
  periodic_function (m, i, len, b, mem_assign);
}
 
static void free_periodic (Memindex m, int i, int len)
{
  periodic_function (m, i, len, NULL, (PeriodicFunction) mem_free);
}
#elif dimension == 2
typedef void (* PeriodicFunction) (Memindex, int, int, int, void *);
 
static void periodic_function (Memindex m, int i, int j, int len, void * b,
             PeriodicFunction f)
{
  f(m, i, j, len, b);
  if (Period.x) {
    int nl = len - 2*GHOSTS;
    for (int l = - 1; l <= 1; l += 2)
      for (int n = i + l*nl; n >= 0 && n < len; n += l*nl)
  f(m, n, j, len, b);
    if (Period.y)
      for (int l = - 1; l <= 1; l += 2)
  for (int n = j + l*nl; n >= 0 && n < len; n += l*nl) {
    f(m, i, n, len, b);
    for (int o = - 1; o <= 1; o += 2)
      for (int p = i + o*nl; p >= 0 && p < len; p += o*nl)
        f(m, p, n, len, b);
  }
  }
  else if (Period.y) {
    int nl = len - 2*GHOSTS;
    for (int l = - 1; l <= 1; l += 2)
      for (int n = j + l*nl; n >= 0 && n < len; n += l*nl)
  f(m, i, n, len, b);
  }
}
             
static void assign_periodic (Memindex m, int i, int j, int len, void * b)
{
  periodic_function (m, i, j, len, b, mem_assign);
}
 
static void free_periodic (Memindex m, int i, int j, int len)
{
  periodic_function (m, i, j, len, NULL, mem_free);
}
#else // dimension == 3
typedef void (* PeriodicFunction) (Memindex, int, int, int, int, void *);
 
static void periodic_function (Memindex m, int i, int j, int k, int len,
             void * b, PeriodicFunction f)
{
  f(m, i, j, k, len, b);
  if (Period.x) {
    int nl = len - 2*GHOSTS;
    for (int l = - 1; l <= 1; l += 2)
      for (int n = i + l*nl; n >= 0 && n < len; n += l*nl)
  f(m, n, j, k, len, b);
    if (Period.y) {
      for (int l = - 1; l <= 1; l += 2)
  for (int n = j + l*nl; n >= 0 && n < len; n += l*nl) {
    f(m, i, n, k, len, b);
    for (int o = - 1; o <= 1; o += 2)
      for (int p = i + o*nl; p >= 0 && p < len; p += o*nl)
        f(m, p, n, k, len, b);
  }
      if (Period.z)
  for (int l = - 1; l <= 1; l += 2)
    for (int n = k + l*nl; n >= 0 && n < len; n += l*nl) {
      f(m, i, j, n, len, b);
      for (int q = - 1; q <= 1; q += 2)
        for (int r = j + q*nl; r >= 0 && r < len; r += q*nl)
    f(m, i, r, n, len, b);
      for (int o = - 1; o <= 1; o += 2)
        for (int p = i + o*nl; p >= 0 && p < len; p += o*nl) {
    f(m, p, j, n, len, b);
    for (int q = - 1; q <= 1; q += 2)
      for (int r = j + q*nl; r >= 0 && r < len; r += q*nl)
        f(m, p, r, n, len, b);
        }
    }
    }
    else if (Period.z)
      for (int l = - 1; l <= 1; l += 2)
  for (int n = k + l*nl; n >= 0 && n < len; n += l*nl) {
    f(m, i, j, n, len, b);
    for (int o = - 1; o <= 1; o += 2)
      for (int p = i + o*nl; p >= 0 && p < len; p += o*nl)
        f(m, p, j, n, len, b);
  }
  }
  else if (Period.y) {
    int nl = len - 2*GHOSTS;
    for (int l = - 1; l <= 1; l += 2)
      for (int n = j + l*nl; n >= 0 && n < len; n += l*nl)
  f(m, i, n, k, len, b);
    if (Period.z)
      for (int l = - 1; l <= 1; l += 2)
  for (int n = k + l*nl; n >= 0 && n < len; n += l*nl) {
    f(m, i, j, n, len, b);
    for (int o = - 1; o <= 1; o += 2)
      for (int p = j + o*nl; p >= 0 && p < len; p += o*nl)
        f(m, i, p, n, len, b);
  }
  }
  else if (Period.z) {
    int nl = len - 2*GHOSTS;
    for (int l = - 1; l <= 1; l += 2)
      for (int n = k + l*nl; n >= 0 && n < len; n += l*nl)
  f(m, i, j, n, len, b);
  }
}
 
static void assign_periodic (Memindex m, int i, int j, int k, int len, void * b)
{
  periodic_function (m, i, j, k, len, b, mem_assign);
}
 
static void free_periodic (Memindex m, int i, int j, int k, int len)
{
  periodic_function (m, i, j, k, len, NULL, (PeriodicFunction) mem_free);
}
#endif // dimension == 3
 
static void alloc_children (Point point)
{
  if (point.level == grid->depth)
    update_depth (+1);
  else if (allocated_child(0,0,0))
    return;
  
  /* low-level memory management */
  Layer * L = tree->L[point.level + 1];
  L->nc++;
  size_t len = sizeof(Cell) + datasize;
  char * b = (char *) mempool_alloc0 (L->pool);
  int i = 2*point.i - GHOSTS;
  for (int k = 0; k < 2; k++, i++) {
#if dimension == 1
    assign_periodic (L->m, i, L->len, b);
    b += len;
#elif dimension == 2
    int j = 2*point.j - GHOSTS;
    for (int l = 0; l < 2; l++, j++) {
      assign_periodic (L->m, i, j, L->len, b);
      b += len;
    }
#else // dimension == 3
    int j = 2*point.j - GHOSTS;
    for (int l = 0; l < 2; l++, j++) {
      int m = 2*point.k - GHOSTS;
      for (int n = 0; n < 2; n++, m++) {
  assign_periodic (L->m, i, j, m, L->len, b);
  b += len;
      }
    }
#endif
  }
  
  int pid = cell.pid;
  for (int _c = 0; _c < 4; _c++) /* foreach_child */ {
    cell.pid = pid;
@if TRASH
    for (int _s = 0; _s < 1; _s++) /* scalar in all */
      s[] = undefined;
@endif    
  }
}
 
#if dimension == 1 
static void free_children (Point point)
{
  /* low-level memory management */
  Layer * L = tree->L[point.level + 1];
  int i = 2*point.i - GHOSTS;
  assert (mem_data (L->m,i));
  mempool_free (L->pool, mem_data (L->m,i));
  for (int k = 0; k < 2; k++, i++)
    free_periodic (L->m, i, L->len);
  if (--L->nc == 0) {
    destroy_layer (L);
    assert (point.level + 1 == grid->depth);
    update_depth (-1);
  }
}
#elif dimension == 2
static void free_children (Point point)
{
  /* low-level memory management */
  Layer * L = tree->L[point.level + 1];
  int i = 2*point.i - GHOSTS, j = 2*point.j - GHOSTS;
  assert (mem_data (L->m,i,j));
  mempool_free (L->pool, mem_data (L->m,i,j));
  for (int k = 0; k < 2; k++)
    for (int l = 0; l < 2; l++)
      free_periodic (L->m, i + k, j + l, L->len);
  if (--L->nc == 0) {
    destroy_layer (L);
    assert (point.level + 1 == grid->depth);
    update_depth (-1);
  }
}
#else // dimension == 3
static void free_children (Point point)
{
  /* low-level memory management */
  Layer * L = tree->L[point.level + 1];
  int i = 2*point.i - GHOSTS;
  assert (mem_data (L->m,i,2*point.j - GHOSTS,2*point.k - GHOSTS));
  mempool_free (L->pool, mem_data (L->m,
           i,2*point.j - GHOSTS,2*point.k - GHOSTS));
  for (int k = 0; k < 2; k++, i++) {
    int j = 2*point.j - GHOSTS;
    for (int l = 0; l < 2; l++, j++) {
      int m = 2*point.k - GHOSTS;
      for (int n = 0; n < 2; n++, m++)
  free_periodic (L->m, i, j, m, L->len);
    }
  }
  if (--L->nc == 0) {
    destroy_layer (L);
    assert (point.level + 1 == grid->depth);
    update_depth (-1);
  }
}
#endif // dimension == 3
 
void increment_neighbors (Point point)
{
  tree->dirty = true;
  if (cell.neighbors++ == 0)
    alloc_children (point);
  for (int _n = 0; _n < GHOSTS/2; _n++)
    if (cell.neighbors++ == 0)
      alloc_children (point);
  cell.neighbors--;
}
 
void decrement_neighbors (Point point)
{
  tree->dirty = true;
  for (int _n = 0; _n < GHOSTS/2; _n++)
    if (allocated(0)) {
      cell.neighbors--;
      if (cell.neighbors == 0)
  free_children (point);
    }
  if (cell.neighbors) {
    int pid = cell.pid;
    for (int _c = 0; _c < 4; _c++) /* foreach_child */ {
      cell.flags = 0;
      cell.pid = pid;
    }
  }
}
 
void realloc_scalar (int size)
{
  /* low-level memory management */
  Tree * q = tree;
  size_t oldlen = sizeof(Cell) + datasize;
  size_t newlen = oldlen + size;
  datasize += size;
  /* the root level is allocated differently */
  Layer * L = q->L[0];
  foreach_mem (L->m, L->len, 1) {
#if dimension == 1
    char * p = (char *) realloc (mem_data (L->m, point.i), newlen*sizeof(char));
    assign_periodic (L->m, point.i, L->len, p);
#elif dimension == 2
    char * p = (char *) realloc (mem_data (L->m, point.i, point.j),
         newlen*sizeof(char));
    assign_periodic (L->m, point.i, point.j, L->len, p);
#else
    char * p = (char *) realloc (mem_data (L->m, point.i, point.j, point.k),
         newlen*sizeof(char));
    assign_periodic (L->m, point.i, point.j, point.k, L->len, p);
#endif
  }
  /* all other levels */
  for (int l = 1; l <= depth(); l++) {
    Layer * L = q->L[l];
    Mempool * oldpool = L->pool;
    L->pool = mempool_new (poolsize (l, newlen), (1 << dimension)*newlen);
    foreach_mem (L->m, L->len, 2) {
      char * new = (char *) mempool_alloc (L->pool);
#if dimension == 1
      for (int k = 0; k < 2; k++) {
  memcpy (new, mem_data (L->m, point.i + k), oldlen);
  assign_periodic (L->m, point.i + k, L->len, new);
  new += newlen;
      }
#elif dimension == 2
      for (int k = 0; k < 2; k++)
  for (int o = 0; o < 2; o++) {
    memcpy (new, mem_data (L->m, point.i + k,point.j + o), oldlen);
    assign_periodic (L->m, point.i + k, point.j + o, L->len, new);
    new += newlen;
  }
#else // dimension == 3
      for (int l = 0; l < 2; l++)
  for (int m = 0; m < 2; m++)
    for (int n = 0; n < 2; n++) {
      memcpy (new, mem_data (L->m, point.i + l, point.j + m, point.k + n),
        oldlen);
      assign_periodic (L->m, point.i + l, point.j + m, point.k + n,
         L->len, new);
      new += newlen;
    }
#endif // dimension == 3
    }
    mempool_destroy (oldpool);
  }
}
 
/* Boundaries */
 
@define VN v.x
@define VT v.y
@define VR v.z
 
#define is_neighbor(...) (allocated(__VA_ARGS__) &&   \
        !is_boundary(neighbor(__VA_ARGS__)))
 
@if _MPI // fixme
@ define disable_fpe_for_mpi() disable_fpe (FE_DIVBYZERO|FE_INVALID)
@ define enable_fpe_for_mpi()  enable_fpe (FE_DIVBYZERO|FE_INVALID)
@else
@ define disable_fpe_for_mpi()
@ define enable_fpe_for_mpi()
@endif
 
static inline void no_restriction (Point point, scalar s);
 
static bool normal_neighbor (Point point, scalar * scalars, vector * vectors)
{
  for (int k = 1; k <= BGHOSTS; k++)
    for (int _d = 0; _d < dimension; _d++)
      for (int i = -k; i <= k; i += 2*k)
  if (is_neighbor(i)) {
    int id = bid(cell);
    for (int _s = 0; _s < 1; _s++) /* scalar in scalars */
      /* foreach_blockf */
        s[] = s.boundary[id](neighborp(i), point, s, NULL);
    for (int _v = 0; _v < 1; _v++) /* vector in vectors */ {
      scalar vn = VN;
      /* foreach_blockf */
        v.x[] = vn.boundary[id](neighborp(i), point, v.x, NULL);
#if dimension >= 2
      scalar vt = VT;
      /* foreach_blockf */
        v.y[] = vt.boundary[id](neighborp(i), point, v.y, NULL);
#endif
#if dimension >= 3
      scalar vr = VR;
      /* foreach_blockf */
        v.z[] = vr.boundary[id](neighborp(i), point, v.z, NULL);
#endif
    }
    return true;
  }
  return false;
}
 
static bool diagonal_neighbor_2D (Point point,
          scalar * scalars, vector * vectors)
{
#if dimension >= 2
  for (int k = 1; k <= BGHOSTS; k++)
#if dimension == 3
    for (int _d = 0; _d < dimension; _d++)
#endif
      for (int i = -k; i <= k; i += 2*k)
  for (int j = -k; j <= k; j += 2*k)
    if (allocated(i,j) && is_neighbor(i,j) &&
        allocated(i,0) && is_boundary(neighbor(i,0)) &&
        allocated(0,j) && is_boundary(neighbor(0,j))) {
      int id1 = bid(neighbor(i,0)), id2 = bid(neighbor(0,j));
      for (int _s = 0; _s < 1; _s++) /* scalar in scalars */
        /* foreach_blockf */
    s[] = (s.boundary[id1](neighborp(i,j), neighborp(i,0), s, NULL) +
           s.boundary[id2](neighborp(i,j), neighborp(0,j), s,NULL) -
           s[i,j]);
      for (int _v = 0; _v < 1; _v++) /* vector in vectors */ {
        scalar vt = VT, vn = VN;
        /* foreach_blockf */
    v.x[] = (vt.boundary[id1](neighborp(i,j), neighborp(i,0), v.x, NULL) +
       vn.boundary[id2](neighborp(i,j), neighborp(0,j), v.x,NULL) -
       v.x[i,j]);
        /* foreach_blockf */
    v.y[] = (vn.boundary[id1](neighborp(i,j), neighborp(i,0), v.y, NULL) +
       vt.boundary[id2](neighborp(i,j), neighborp(0,j), v.y,NULL) -
       v.y[i,j]);
#if dimension == 3
        scalar vr = VR;
        /* foreach_blockf */
    v.z[] = (vr.boundary[id1](neighborp(i,j), neighborp(i,0), v.z, NULL) +
       vr.boundary[id2](neighborp(i,j), neighborp(0,j), v.z,NULL) -
       v.z[i,j]);
#endif
      }
      return true;
    }
#endif // dimension >= 2
  return false;
}
 
static bool diagonal_neighbor_3D (Point point,
          scalar * scalars, vector * vectors)
{
#if dimension == 3
  for (int n = 1; n <= BGHOSTS; n++)
    for (int i = -n; i <= n; i += 2*n)
      for (int j = -n; j <= n; j += 2*n)
  for (int k = -n; k <= n; k += 2*n)
    if (is_neighbor(i,j,k) &&
        is_boundary(neighbor(i,j,0)) &&
        is_boundary(neighbor(i,0,k)) &&
        is_boundary(neighbor(0,j,k))) {
      Point
        n0 = neighborp(i,j,k),
        n1 = neighborp(i,j,0),
        n2 = neighborp(i,0,k),
        n3 = neighborp(0,j,k);
      int
        id1 = bid(neighbor(i,j,0)),
        id2 = bid(neighbor(i,0,k)),
        id3 = bid(neighbor(0,j,k));
      for (int _s = 0; _s < 1; _s++) /* scalar in scalars */
        s[] = (s.boundary[id1](n0,n1,s,NULL) +
         s.boundary[id2](n0,n2,s,NULL) +
         s.boundary[id3](n0,n3,s,NULL) -
         2.*s[i,j,k]);
      for (int _v = 0; _v < 1; _v++) /* vector in vectors */ {
        scalar vt = VT, vn = VN, vr = VR;
        v.x[] = (vt.boundary[id1](n0,n1,v.x,NULL) +
           vt.boundary[id2](n0,n2,v.x,NULL) +
           vn.boundary[id3](n0,n3,v.x,NULL) -
           2.*v.x[i,j,k]);
        v.y[] = (vt.boundary[id1](n0,n1,v.y,NULL) +
           vn.boundary[id2](n0,n2,v.y,NULL) +
           vt.boundary[id3](n0,n3,v.y,NULL) -
           2.*v.y[i,j,k]);
        v.z[] = (vn.boundary[id1](n0,n1,v.z,NULL) +
           vr.boundary[id2](n0,n2,v.z,NULL) +
           vr.boundary[id3](n0,n3,v.z,NULL) -
           2.*v.z[i,j,k]);
      }
      return true;
    }
#endif
  return false;
}
 
#if dimension > 1
for (int _d = 0; _d < dimension; _d++)
static Point tangential_neighbor_x (Point point, bool * zn)
{
  for (int k = 1; k <= BGHOSTS; k++)
    for (int j = -k; j <= k; j += 2*k) {
      if (is_neighbor(0,j) || is_neighbor(-1,j)) {
  *zn = false;
  return neighborp(0,j);
      }
#if dimension == 3
      // fixme: what about diagonals?
      if (is_neighbor(0,0,j) || is_neighbor(-1,0,j)) {
  *zn = true;
  return neighborp(0,0,j);
      }
#endif // dimension == 3
    }
  return (Point){.level = -1};
}
#endif // dimension > 1
 
static inline bool is_boundary_point (Point point) {
  return is_boundary (cell);
}
 
static void box_boundary_level (const Boundary * b, scalar * list, int l)
{
  disable_fpe_for_mpi();
  scalar * scalars = NULL;
  vector * vectors = NULL, * faces = NULL;
  for (int _s = 0; _s < 1; _s++) /* scalar in list */
    if (!is_constant(s) && s.refine != no_restriction) {
      if (s.v.x.i == s.i) {
  if (s.face)
    faces = vectors_add (faces, s.v);
  else
    vectors = vectors_add (vectors, s.v);
      }
      else if (s.v.x.i < 0 && s.boundary[0])
  scalars = list_add (scalars, s);
    }
  
  foreach_boundary_level (l) {
    if (!normal_neighbor (point, scalars, vectors) &&
  !diagonal_neighbor_2D (point, scalars, vectors) &&
  !diagonal_neighbor_3D (point, scalars, vectors)) {
      // no neighbors
      for (int _s = 0; _s < 1; _s++) /* scalar in scalars */
  /* foreach_blockf */
    s[] = undefined;
      for (int _v = 0; _v < 1; _v++) /* vector in vectors */
  for (int _d = 0; _d < dimension; _d++)
    /* foreach_blockf */
      v.x[] = undefined;
    }
    if (faces) {
      int id = bid(cell);
      for (int _d = 0; _d < dimension; _d++)
  for (int i = -1; i <= 1; i += 2) {
    // normal neighbor for faces
    if (is_neighbor(i)) {
      Point neighbor = neighborp(i);
      for (int _v = 0; _v < 1; _v++) /* vector in faces */ {
        scalar vn = VN;
        if (vn.boundary[id])
    /* foreach_blockf */
      v.x[(i + 1)/2] = vn.boundary[id](neighbor, point, v.x, NULL);
      }
    }
#if dimension > 1
    else if (i == -1) {
      // tangential neighbor
      bool zn;
      Point neighbor = tangential_neighbor_x (point, &zn);
      if (neighbor.level >= 0) {
        int id = is_boundary_point (neighbor) ?
    bid(neighbor(-1)) : bid(cell);
        for (int _v = 0; _v < 1; _v++) /* vector in faces */ {
#if dimension == 2
    scalar vt = VT;
#else // dimension == 3
    scalar vt = zn ? VT : VR;
#endif
    /* foreach_blockf */
      v.x[] = vt.boundary[id](neighbor, point, v.x, NULL);
        }
      }
      else
        // no neighbor
        for (int _v = 0; _v < 1; _v++) /* vector in faces */
    /* foreach_blockf */
      v.x[] = 0.;
    }
#endif // dimension > 1
  }
    }
  }
 
  free (scalars);
  free (vectors);
  free (faces);
  enable_fpe_for_mpi();
}
 
#undef is_neighbor
 
@undef VN
@undef VT
@define VN _attribute[s.i].v.x
@define VT _attribute[s.i].v.y
 
static double masked_average (Point point, scalar s)
{
  double sum = 0., n = 0.;
  for (int _c = 0; _c < 4; _c++) /* foreach_child */
    if (!is_boundary(cell) && s[] != nodata)
      sum += s[], n++;
  return n ? sum/n : nodata;
}
 
for (int _d = 0; _d < dimension; _d++)
static double masked_average_x (Point point, scalar s)
{
  double sum = 0., n = 0.;
  for (int _c = 0; _c < 4; _c++) /* foreach_child */
    if (child.x < 0 && (!is_boundary(cell) || !is_boundary(neighbor(1))) &&
  s[1] != nodata)
      sum += s[1], n++;
  return n ? sum/n : nodata;
}
 
static void masked_boundary_restriction (const Boundary * b,
           scalar * list, int l)
{
  scalar * scalars = NULL;
  vector * faces = NULL;
  for (int _s = 0; _s < 1; _s++) /* scalar in list */
    if (!is_constant(s) && s.refine != no_restriction) {
      if (s.v.x.i == s.i && s.face)
  faces = vectors_add (faces, s.v);
      else
  scalars = list_add (scalars, s);
    }
  
  foreach_halo (restriction, l) /* foreach_block */ { // fixme: should use /* foreach_blockf */
    for (int _s = 0; _s < 1; _s++) /* scalar in scalars */
      s[] = masked_average (parent, s);
    for (int _v = 0; _v < 1; _v++) /* vector in faces */
      for (int _d = 0; _d < dimension; _d++) {
  double average = masked_average_x (parent, v.x);
  if (is_boundary(neighbor(-1)))
    v.x[] = average;
  if (is_boundary(neighbor(1)))
    v.x[1] = average;
      }
  }
 
  free (scalars);
  free (faces);  
}
 
macro mask (double func) {
  foreach_cell_post(!is_leaf(cell)) {
    if (is_leaf(cell)) {
      int bid = (func);
      if (bid >= 0)
  cell.pid = - bid - 1;
    }
    else { /* not a leaf */
      int pid = -1;
      for (int _c = 0; _c < 4; _c++) /* foreach_child */
  if (cell.pid >= 0 || pid < 0)
    pid = cell.pid;
      cell.pid = pid;
      if (pid < 0) {
  /* fixme: call coarsen_cell()? */
  cell.flags |= leaf;
  decrement_neighbors (point);
      }
    }
  }
  tree->dirty = true;
}
 
static void free_cache (CacheLevel * c)
{
  for (int l = 0; l <= depth(); l++)
    free (c[l].p);
  free (c);
}
 
void free_grid (void)
{
  if (!grid)
    return;
  free_boundaries();
  Tree * q = tree;
  free (q->leaves.p);
  free (q->faces.p);
  free (q->vertices.p);
  free (q->refined.p);
  /* low-level memory management */
  /* the root level is allocated differently */
  Layer * L = q->L[0];
  foreach_mem (L->m, L->len, 1) {
#if dimension == 1
    free (mem_data (L->m, point.i));
#elif dimension == 2
    free (mem_data (L->m, point.i, point.j));
#else // dimension == 3
    free (mem_data (L->m, point.i, point.j, point.k));
#endif // dimension == 3
  }
  for (int l = 0; l <= depth(); l++)
    destroy_layer (q->L[l]);
  q->L = &(q->L[-1]);
  free (q->L);
  free_cache (q->active);
  free_cache (q->prolongation);
  free_cache (q->boundary);
  free_cache (q->restriction);
  free (q);
  grid = NULL;
}
 
static void refine_level (int depth);
 
trace
void init_grid (int n)
{
  // check 64 bits structure alignment
  assert (sizeof(Cell) % 8 == 0);
  
  free_grid();
  int depth = 0;
  while (n > 1) {
    if (n % 2) {
      fprintf (stderr, "tree: N must be a power-of-two\n");
      exit (1);
    }
    n /= 2;
    depth++;
  }
  Tree * q = qcalloc (1, Tree);
  grid = (Grid *) q;
  grid->depth = 0;
 
  /* low-level memory management */
  q->L = qmalloc (2, Layer *);
  /* make sure we don't try to access level -1 */
  q->L[0] = NULL; q->L = &(q->L[1]);
  /* initialise the root cell */
  Layer * L = new_layer (0);
  q->L[0] = L;
#if dimension == 1
  for (int i = Period.x*GHOSTS; i < L->len - Period.x*GHOSTS; i++)
    assign_periodic (L->m, i, L->len, 
         (char *) calloc (1, sizeof(Cell) + datasize));
  CELL(mem_data (L->m,GHOSTS)).flags |= leaf;
  if (pid() == 0)
    CELL(mem_data (L->m,GHOSTS)).flags |= active;
  for (int k = -GHOSTS*(1 - Period.x); k <= GHOSTS*(1 - Period.x); k++)
    CELL(mem_data (L->m,GHOSTS+k)).pid = (k < 0 ? -1 - left :
            k > 0 ? -1 - right :
            0);
  CELL(mem_data (L->m,GHOSTS)).pid = 0;
#elif dimension == 2
  for (int i = Period.x*GHOSTS; i < L->len - Period.x*GHOSTS; i++)
    for (int j = Period.y*GHOSTS; j < L->len - Period.y*GHOSTS; j++)
      assign_periodic (L->m, i, j, L->len,
           (char *) calloc (1, sizeof(Cell) + datasize));
  CELL(mem_data (L->m,GHOSTS,GHOSTS)).flags |= leaf;
  if (pid() == 0)
    CELL(mem_data (L->m,GHOSTS,GHOSTS)).flags |= active;
  for (int k = - GHOSTS*(1 - Period.x); k <= GHOSTS*(1 - Period.x); k++)
    for (int l = -GHOSTS*(1 - Period.y); l <= GHOSTS*(1 - Period.y); l++)
      CELL(mem_data (L->m,GHOSTS+k,GHOSTS+l)).pid =
  (k < 0 ? -1 - left :
   k > 0 ? -1 - right :
   l > 0 ? -1 - top :
   l < 0 ? -1 - bottom :
   0);
  CELL(mem_data (L->m,GHOSTS,GHOSTS)).pid = 0;
#else // dimension == 3
  for (int i = Period.x*GHOSTS; i < L->len - Period.x*GHOSTS; i++)
    for (int j = Period.y*GHOSTS; j < L->len - Period.y*GHOSTS; j++)
      for (int k = Period.z*GHOSTS; k < L->len - Period.z*GHOSTS; k++)
  assign_periodic (L->m, i, j, k, L->len,
       (char *) calloc (1, sizeof(Cell) + datasize));
  CELL(mem_data (L->m,GHOSTS,GHOSTS,GHOSTS)).flags |= leaf;
  if (pid() == 0)
    CELL(mem_data (L->m,GHOSTS,GHOSTS,GHOSTS)).flags |= active;
  for (int k = - GHOSTS*(1 - Period.x); k <= GHOSTS*(1 - Period.x); k++)
    for (int l = -GHOSTS*(1 - Period.y); l <= GHOSTS*(1 - Period.y); l++)
      for (int n = -GHOSTS*(1 - Period.z); n <= GHOSTS*(1 - Period.z); n++)
  CELL(mem_data (L->m,GHOSTS+k,GHOSTS+l,GHOSTS+n)).pid =
    (k > 0 ? -1 - right :
     k < 0 ? -1 - left :
     l > 0 ? -1 - top :
     l < 0 ? -1 - bottom :
     n > 0 ? -1 - front :
     n < 0 ? -1 - back :
     0);
  CELL(mem_data (L->m,GHOSTS,GHOSTS,GHOSTS)).pid = 0;
#endif // dimension == 3
  q->active = qcalloc (1, CacheLevel);
  q->prolongation = qcalloc (1, CacheLevel);
  q->boundary = qcalloc (1, CacheLevel);
  q->restriction = qcalloc (1, CacheLevel);
  q->dirty = true;
  N = 1 << depth;
@if _MPI
  void mpi_boundary_new();
  mpi_boundary_new();
@endif
  // boundaries
  Boundary * b = qcalloc (1, Boundary);
  b->level = box_boundary_level;
  b->restriction = masked_boundary_restriction;
  add_boundary (b);
  refine_level (depth);
  reset (all, 0.);
  update_cache();
}
 
#if dimension == 2
void check_two_one (void)
{
  for (int _i = 0; _i < _N; _i++) /* foreach_leaf */
    if (level > 0)
      for (int k = -1; k <= 1; k++)
  for (int l = -1; l <= 1; l++) {
    /* fixme: all this mess is just to ignore ghost cells */
    int i = (point.i + GHOSTS)/2 + k;
    int j = (point.j + GHOSTS)/2 + l;
    double Delta = 1./(1 << point.level);
    double x = ((i - GHOSTS + 0.5)*Delta*2. - 0.5);
    double y = ((j - GHOSTS + 0.5)*Delta*2. - 0.5);
    if (x > -0.5 && x < 0.5 && y > -0.5 && y < 0.5 && 
        !(aparent(k,l).flags & active)) {
      FILE * fp = fopen("check_two_one_loc", "w");
      fprintf (fp,
         "# %d %d\n"
         "%g %g\n%g %g\n",
         k, l,
         ((point.i - GHOSTS + 0.5)* - 0.5),
         ((point.j - GHOSTS + 0.5)*Delta - 0.5),
         x, y);
      fclose (fp);
#if 0
      fp = fopen("check_two_one", "w");
      output_cells (fp);
      fclose (fp);
#endif
      assert (false);
    }
  }
}
#endif
 
Point locate (double xp = 0., double yp = 0., double zp = 0.)
{
  for (int l = depth(); l >= 0; l--) {
    Point point = {0};
    point.level = l;
    int n = 1 << point.level;
    point.i = (xp - X0)/L0*n + GHOSTS;
#if dimension >= 2
    point.j = (yp - Y0)/L0*n + GHOSTS;
#endif
#if dimension >= 3
    point.k = (zp - Z0)/L0*n + GHOSTS;
#endif
    if (point.i >= 0 && point.i < n + 2*GHOSTS
#if dimension >= 2
  && point.j >= 0 && point.j < n + 2*GHOSTS
#endif
#if dimension >= 3
  && point.k >= 0 && point.k < n + 2*GHOSTS
#endif
  ) {
      if (allocated(0) && is_local(cell) && is_leaf(cell))
  return point;
    }
    else
      break;
  }
  Point point = {0};
  point.level = -1;
  return point;
}
 
// return true if the tree is "full" i.e. all the leaves are at the
// same level
bool tree_is_full()
{
  update_cache();
  return (grid->tn == 1L << grid->maxdepth*dimension);
}
 
#include "variables.h"
 
macro2 for (int _b = 0; _b < 1; _b++) /* boundary int _b, Reduce reductions = None */ {
  for (int _l = depth(); _l >= 0; _l--)
    foreach_boundary_level (_l, reductions) {
      POINT_VARIABLES();
      if (bid(cell) == _b)
  for (int _d = 0; _d < dimension; _d++) {
    for (int _i = -1; _i <= 1; _i += 2) {
      if (_d == 0) ig = _i; else if (_d == 1) jg = _i; else kg = _i;
      if (allocated(-ig,-jg,-kg) &&
    is_leaf (neighbor(-ig,-jg,-kg)) &&
    !is_boundary(neighbor(-ig,-jg,-kg)) &&
    is_local(neighbor(-ig,-jg,-kg))) {
        point.i -= ig; x -= ig*Delta/2.; 
#if dimension >= 2
        point.j -= jg; y -= jg*Delta/2.; 
#endif
#if dimension >= 3
        point.k -= kg; z -= kg*Delta/2.;
#endif
        {...}
        point.i += ig; x += ig*Delta/2.;   
#if dimension >= 2
        point.j += jg; y += jg*Delta/2.; 
#endif
#if dimension >= 3
        point.k += kg; z += kg*Delta/2.;
#endif
            }
    }
    ig = jg = kg = 0;
  }
    }
}
 
macro2 for (int _i = 0; _i < _N; _i++) /* foreach_vertex */
{
  update_cache();
  for (int _i = 0; _i < 1; _i++) /* cache tree->vertices, reductions */ {
    int ig = -1; NOT_UNUSED (ig);
#if dimension >= 2
    int jg = -1; NOT_UNUSED (jg);
#endif
#if dimension >= 3
    int kg = -1; NOT_UNUSED (kg);
#endif    
    {...}
  }
}
 
#include "tree-common.h"
 
// overload the default periodic() function
void tree_periodic (int dir)
{
  int depth = grid ? depth() : -1;
  if (grid)
    free_grid();
  periodic (dir);
  if (depth >= 0)
    init_grid (1 << depth);
}
#define periodic(dir) tree_periodic(dir)
 
@if _MPI
#include "tree-mpi.h"
#include "balance.h"
@else // !_MPI
void mpi_boundary_refine  (scalar * list){}
void mpi_boundary_coarsen (int a, int b){}
void mpi_boundary_update  (scalar * list) {
  for (int _s = 0; _s < 1; _s++) /* scalar in list */
    set_dirty_stencil (s);
  boundary (list);
}
@endif // !_MPI