/** @file grid.h
 *  @brief Definitions for working with general grids.
 */
#ifndef _GRID_H_

#ifndef _TREE_H_
# include "tree.h"
#endif

#define RECT_COORDS  int r0, r1, z0, z1, ntheta
/*!< Coordinate names for the grids */

#define BND_AT_R        0
#define BND_AT_Z        2
#define BND_MIN         0
#define BND_MAX         1

#define BND_LEFT        (BND_AT_R | BND_MIN)
#define BND_RIGHT       (BND_AT_R | BND_MAX)
#define BND_BOTTOM      (BND_AT_Z | BND_MIN)
#define BND_TOP         (BND_AT_Z | BND_MAX)

#define BND_MASK(B_)    (1 << (B_))
#define BND_NONE             0
#define BND_MASK_LEFT        BND_MASK(BND_LEFT)
#define BND_MASK_RIGHT       BND_MASK(BND_RIGHT)
#define BND_MASK_BOTTOM      BND_MASK(BND_BOTTOM)
#define BND_MASK_TOP         BND_MASK(BND_TOP)
#define BND_MASK_ALL         (BND_MASK_LEFT | BND_MASK_RIGHT\
			      | BND_MASK_BOTTOM | BND_MASK_TOP)

#define GRID_INSIDE    (1 << (BND_TOP + 1))

/* Locations of the cell centers. */
#define r_at(r_, level_) (((double) (r_) + 0.5) * dr[level_])
#define z_at(z_, level_) (((double) (z_) + 0.5) * dz[level_])
#define theta_at(theta_) ((double) (theta_) * dtheta)

/* Locations of the electric field. */
#define er_r_at(i_, level_) ((i_ + 1) * dr[level_])
#define er_z_at(j_, level_) (((double) (j_) + 0.5) * dz[level_])
#define ez_r_at(i_, level_) (((double) (i_) + 0.5) * dr[level_])
#define ez_z_at(j_, level_) ((j_ + 1) * dz[level_])
#define etheta_theta_at(i_) (((double) (i_) + 0.5) * dtheta)


/** We can produce easily a true 2d code by replacing r when it appears
 *  due to the cylindrical coordinates by this macros.  This has a rather
 *  large performance penalty, so it is better to do that in compile time
 */

/* Locations of the electric field. */
#ifdef TRUE2D
# define cyl_q(X_) 1
# define cyl_er_r_at(i_, level_) 1
# define cyl_r_at(r_, level_) 1
#else
# define cyl_q(X_) (X_)
# define cyl_er_r_at(i_, level_) er_r_at(i_, level_) 
# define cyl_r_at(r_, level_) r_at(r_, level_)
#endif

/* These are shortcuts to make the iteration over grid cells easier.
 *  Note that it is important that the inner loop has the smallest
 *  stride to minimize the cache faults.
 */
#define iter_grid_z(grid_, iz_) \
  for(iz_ = grid_->z0; iz_ < grid_->z1; iz_++) 

#define iter_grid_r(grid_, ir_) \
  for(ir_ = grid_->r0; ir_ < grid_->r1; ir_++) 

#define iter_grid_theta(grid_, it_) \
  for(it_ = 0; it_ < grid_->ntheta; it_++) 

#define iter_grid(grid_, ir_, iz_) \
  iter_grid_z(grid_, iz_)	   \
       iter_grid_r(grid_, ir_)

#define iter_grid_3d(grid_, ir_, iz_, it_)		\
  iter_grid_theta(grid_, it_)				\
    iter_grid_z(grid_, iz_)				\
       iter_grid_r(grid_, ir_)


/* Functions to iterate over the grids _including_ 1 extra boundary */
#define iter_grid_z_1(grid_, iz_)			\
  for(iz_ = grid_->z0 - 1; iz_ < grid_->z1 + 1; iz_++) 

#define iter_grid_r_1(grid_, ir_) \
  for(ir_ = grid_->r0 - 1; ir_ < grid_->r1 + 1; ir_++) 

#define iter_grid_1(grid_, ir_, iz_) \
  iter_grid_z_1(grid_, iz_)	   \
       iter_grid_r(grid_, ir_)

#define iter_grid_3d_1(grid_, ir_, iz_, it_)		\
  iter_grid_theta(grid_, it_)				\
    iter_grid_z_1(grid_, iz_)				\
       iter_grid_r_1(grid_, ir_)


/* These include n extra cells in the boundary. We avoid that when 
 * ntheta == 1 because we are in the 2D case and there are no boundaries
 * in theta*/
#define iter_grid_theta_n(grid_, it_, n_)			\
  for(it_ = (grid_->ntheta == 1? 0: -n_); it_ < grid_->ntheta	\
	+ (grid_->ntheta == 1? 0: n_); it_++) 

#define iter_grid_z_n(grid_, iz_, n_)			\
  for(iz_ = grid_->z0 - n_; iz_ < grid_->z1 + n_; iz_++) 

#define iter_grid_r_n(grid_, ir_, n_)			\
  for(ir_ = grid_->r0 - n_; ir_ < grid_->r1 + n_; ir_++) 

#define iter_grid_n(grid_, ir_, iz_, n_)		\
  iter_grid_z_n(grid_, iz_, n_)				\
       iter_grid_r_n(grid_, ir_, n_)

#define iter_grid_3d_n(grid_, ir_, iz_, it_, n_)	\
  iter_grid_theta_n(grid_, it_, n_)			\
     iter_grid_z_n(grid_, iz_, n_)			\
       iter_grid_r_n(grid_, ir_, n_)


/* Iter through the parent indices that fill one grid. */
#define iter_grid_parent(grid_, ir_, iz_) \
  for(ir_ = (grid_->r0 >> 1) + (grid_->r0 % 2);				\
      ir_ < (grid_->r1 >> 1); ir_++)					\
    for(iz_ = (grid_->z0 >> 1) + (grid_->z0 % 2);			\
	iz_ < (grid_->z1 >> 1); iz_++)


/** @brief The cdr grids and poisson grids types are cdr_grid_t and
 *  pois_grid_t, defined in cdr.h and poisson.h.
 *
 *  This structure contains the common parts of both structures and
 *  thus can be used in general routines that accept any kind of grid.
 *  The C standard promises that we can du this.
 */
typedef struct grid_t grid_t;
struct grid_t {
  RECT_COORDS;
  LEAF_FIELDS(grid_t);
  /* ext_bound indicates which boundaries are external boundaries.
     e.g. ext_bound & BND_LEFT_TOP tells us whether the left-top boundary
     is an external one. */
  int ext_bound;
};

/** @brief To ease typing, you can write e.g.
 * printf ("Solving grid (" grid_printf_str ")\n", grid_printf_args(grid));
 */
#define grid_printf_str "{r0 = %d, z0 = %d, r1 = %d, z1 = %d, level = %d}"
#define grid_printf_args(G_) (G_)->r0, (G_)->z0, (G_)->r1, (G_)->z1,	\
    (G_)->level

#define _GRID_H_
#endif