diff --git a/include/grid.h b/include/grid.h
new file mode 100644
index 0000000000000000000000000000000000000000..576499613c6134134a6eb13f58227bb4a901c547
--- /dev/null
+++ b/include/grid.h
@@ -0,0 +1,159 @@
+/** @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