/*

 mesh.hh - This file is part of MUSIC -
 a code to generate multi-scale initial conditions
 for cosmological simulations

 Copyright (C) 2010  Oliver Hahn

*/

#ifndef __MESH_HH
#define __MESH_HH

#include <iomanip>
#include <iostream>
#include <stdexcept>
#include <vector>

#include <math.h>

#include "config_file.hh"
#include "log.hh"

#include "region_generator.hh"

class refinement_mask {
protected:
  std::vector<short> mask_;
  size_t nx_, ny_, nz_;

public:
  refinement_mask(void) : nx_(0), ny_(0), nz_(0) {}

  refinement_mask(size_t nx, size_t ny, size_t nz, short value = 0.) : nx_(nx), ny_(ny), nz_(nz) {
    mask_.assign(nx_ * ny_ * nz_, value);
  }

  refinement_mask(const refinement_mask &r) {
    nx_ = r.nx_;
    ny_ = r.ny_;
    nz_ = r.nz_;
    mask_ = r.mask_;
  }

  refinement_mask &operator=(const refinement_mask &r) {
    nx_ = r.nx_;
    ny_ = r.ny_;
    nz_ = r.nz_;
    mask_ = r.mask_;

    return *this;
  }

  void init(size_t nx, size_t ny, size_t nz, short value = 0.) {
    nx_ = nx;
    ny_ = ny;
    nz_ = nz;
    mask_.assign(nx_ * ny_ * nz_, value);
  }

  const short &operator()(size_t i, size_t j, size_t k) const { return mask_[(i * ny_ + j) * nz_ + k]; }

  short &operator()(size_t i, size_t j, size_t k) { return mask_[(i * ny_ + j) * nz_ + k]; }

  size_t count_flagged(void) {
    size_t count = 0;
    for (size_t i = 0; i < mask_.size(); ++i)
      if (mask_[i])
        ++count;
    return count;
  }

  size_t count_notflagged(void) {
    size_t count = 0;
    for (size_t i = 0; i < mask_.size(); ++i)
      if (!mask_[i])
        ++count;
    return count;
  }
};

//! base class for all things that have rectangular mesh structure
template <typename T> class Meshvar {
public:
  typedef T real_t;

  size_t m_nx, //!< x-extent of the rectangular mesh
      m_ny,    //!< y-extent of the rectangular mesh
      m_nz;    //!< z-extent of the rectangular mesh

  int m_offx, //!< x-offset of the grid (just as a helper, not used inside the class)
      m_offy, //!< y-offset of the grid (just as a helper, not used inside the class)
      m_offz; //!< z-offset of the grid (just as a helper, not used inside the class)

  real_t *m_pdata; //!< pointer to the dynamic data array

  //! constructor for cubic mesh
  explicit Meshvar(size_t n, int offx, int offy, int offz)
      : m_nx(n), m_ny(n), m_nz(n), m_offx(offx), m_offy(offy), m_offz(offz) {
    m_pdata = new real_t[m_nx * m_ny * m_nz];
  }

  //! constructor for rectangular mesh
  Meshvar(size_t nx, size_t ny, size_t nz, int offx, int offy, int offz)
      : m_nx(nx), m_ny(ny), m_nz(nz), m_offx(offx), m_offy(offy), m_offz(offz) {
    m_pdata = new real_t[m_nx * m_ny * m_nz];
  }

  //! variant copy constructor with optional copying of the actual data
  Meshvar(const Meshvar<real_t> &m, bool copy_over = true) {
    m_nx = m.m_nx;
    m_ny = m.m_ny;
    m_nz = m.m_nz;

    m_offx = m.m_offx;
    m_offy = m.m_offy;
    m_offz = m.m_offz;

    m_pdata = new real_t[m_nx * m_ny * m_nz];

    if (copy_over)
      for (size_t i = 0; i < m_nx * m_ny * m_nz; ++i)
        m_pdata[i] = m.m_pdata[i];
  }

  //! standard copy constructor
  explicit Meshvar(const Meshvar<real_t> &m) {
    m_nx = m.m_nx;
    m_ny = m.m_ny;
    m_nz = m.m_nz;

    m_offx = m.m_offx;
    m_offy = m.m_offy;
    m_offz = m.m_offz;

    m_pdata = new real_t[m_nx * m_ny * m_nz];

    for (size_t i = 0; i < m_nx * m_ny * m_nz; ++i)
      m_pdata[i] = m.m_pdata[i];
  }

  //! destructor
  ~Meshvar() {
    if (m_pdata != NULL)
      delete[] m_pdata;
  }

  //! deallocate the data, but keep the structure
  inline void deallocate(void) {
    if (m_pdata != NULL)
      delete[] m_pdata;
    m_pdata = NULL;
  }

  //! get extent of the mesh along a specified dimension (const)
  inline size_t size(unsigned dim) const {
    if (dim == 0)
      return m_nx;
    if (dim == 1)
      return m_ny;
    return m_nz;
  }

  //! get offset of the mesh along a specified dimension  (const)
  inline int offset(unsigned dim) const {
    if (dim == 0)
      return m_offx;
    if (dim == 1)
      return m_offy;
    return m_offz;
  }

  //! get extent of the mesh along a specified dimension
  inline int &offset(unsigned dim) {
    if (dim == 0)
      return m_offx;
    if (dim == 1)
      return m_offy;
    return m_offz;
  }

  //! set all the data to zero values
  void zero(void) {
    for (size_t i = 0; i < m_nx * m_ny * m_nz; ++i)
      m_pdata[i] = 0.0;
  }

  //! direct array random acces to the data block
  inline real_t *operator[](const size_t i) { return &m_pdata[i]; }

  //! direct array random acces to the data block (const)
  inline const real_t *operator[](const size_t i) const { return &m_pdata[i]; }

  //! 3D random access to the data block via index 3-tuples
  inline real_t &operator()(const int ix, const int iy, const int iz) {
#ifdef DEBUG
    if (ix < 0 || ix >= (int)m_nx || iy < 0 || iy >= (int)m_ny || iz < 0 || iz >= (int)m_nz)
      LOGERR("Array index (%d,%d,%d) out of bounds", ix, iy, iz);
#endif

    return m_pdata[((size_t)ix * m_ny + (size_t)iy) * m_nz + (size_t)iz];
  }

  //! 3D random access to the data block via index 3-tuples (const)
  inline const real_t &operator()(const int ix, const int iy, const int iz) const {
#ifdef DEBUG
    if (ix < 0 || ix >= (int)m_nx || iy < 0 || iy >= (int)m_ny || iz < 0 || iz >= (int)m_nz)
      LOGERR("Array index (%d,%d,%d) out of bounds", ix, iy, iz);
#endif

    return m_pdata[((size_t)ix * m_ny + (size_t)iy) * m_nz + (size_t)iz];
  }

  //! direct multiplication of the whole data block with a number
  Meshvar<real_t> &operator*=(real_t x) {
    for (size_t i = 0; i < m_nx * m_ny * m_nz; ++i)
      m_pdata[i] *= x;
    return *this;
  }

  //! direct addition of a number to the whole data block
  Meshvar<real_t> &operator+=(real_t x) {
    for (size_t i = 0; i < m_nx * m_ny * m_nz; ++i)
      m_pdata[i] += x;
    return *this;
  }

  //! direct element-wise division of the whole data block by a number
  Meshvar<real_t> &operator/=(real_t x) {
    for (size_t i = 0; i < m_nx * m_ny * m_nz; ++i)
      m_pdata[i] /= x;
    return *this;
  }

  //! direct subtraction of a number from the whole data block
  Meshvar<real_t> &operator-=(real_t x) {
    for (size_t i = 0; i < m_nx * m_ny * m_nz; ++i)
      m_pdata[i] -= x;
    return *this;
  }

  //! direct element-wise multiplication with another compatible mesh
  Meshvar<real_t> &operator*=(const Meshvar<real_t> &v) {
    if (v.m_nx * v.m_ny * v.m_nz != m_nx * m_ny * m_nz) {
      LOGERR("Meshvar::operator*= : attempt to operate on incompatible data");
      throw std::runtime_error("Meshvar::operator*= : attempt to operate on incompatible data");
    }
    for (size_t i = 0; i < m_nx * m_ny * m_nz; ++i)
      m_pdata[i] *= v.m_pdata[i];

    return *this;
  }

  //! direct element-wise division with another compatible mesh
  Meshvar<real_t> &operator/=(const Meshvar<real_t> &v) {
    if (v.m_nx * v.m_ny * v.m_nz != m_nx * m_ny * m_nz) {
      LOGERR("Meshvar::operator/= : attempt to operate on incompatible data");
      throw std::runtime_error("Meshvar::operator/= : attempt to operate on incompatible data");
    }

    for (size_t i = 0; i < m_nx * m_ny * m_nz; ++i)
      m_pdata[i] /= v.m_pdata[i];

    return *this;
  }

  //! direct element-wise addition of another compatible mesh
  Meshvar<real_t> &operator+=(const Meshvar<real_t> &v) {
    if (v.m_nx * v.m_ny * v.m_nz != m_nx * m_ny * m_nz) {
      LOGERR("Meshvar::operator+= : attempt to operate on incompatible data");
      throw std::runtime_error("Meshvar::operator+= : attempt to operate on incompatible data");
    }
    for (size_t i = 0; i < m_nx * m_ny * m_nz; ++i)
      m_pdata[i] += v.m_pdata[i];

    return *this;
  }

  //! direct element-wise subtraction of another compatible mesh
  Meshvar<real_t> &operator-=(const Meshvar<real_t> &v) {
    if (v.m_nx * v.m_ny * v.m_nz != m_nx * m_ny * m_nz) {
      LOGERR("Meshvar::operator-= : attempt to operate on incompatible data");
      throw std::runtime_error("Meshvar::operator-= : attempt to operate on incompatible data");
    }
    for (size_t i = 0; i < m_nx * m_ny * m_nz; ++i)
      m_pdata[i] -= v.m_pdata[i];

    return *this;
  }

  //! assignment operator for rectangular meshes
  Meshvar<real_t> &operator=(const Meshvar<real_t> &m) {
    m_nx = m.m_nx;
    m_ny = m.m_ny;
    m_nz = m.m_nz;

    m_offx = m.m_offx;
    m_offy = m.m_offy;
    m_offz = m.m_offz;

    if (m_pdata != NULL)
      delete m_pdata;

    m_pdata = new real_t[m_nx * m_ny * m_nz];

    for (size_t i = 0; i < m_nx * m_ny * m_nz; ++i)
      m_pdata[i] = m.m_pdata[i];

    return *this;
  }

  real_t *get_ptr(void) { return m_pdata; }
};

//! MeshvarBnd derived class adding boundary ghost cell functionality
template <typename T> class MeshvarBnd : public Meshvar<T> {
  using Meshvar<T>::m_nx;
  using Meshvar<T>::m_ny;
  using Meshvar<T>::m_nz;
  using Meshvar<T>::m_pdata;

public:
  typedef T real_t;

  //! number of boundary (ghost) cells
  int m_nbnd;

  //! most general constructor
  MeshvarBnd(int nbnd, size_t nx, size_t ny, size_t nz, size_t xoff, size_t yoff, size_t zoff)
      : Meshvar<real_t>(nx + 2 * nbnd, ny + 2 * nbnd, nz + 2 * nbnd, xoff, yoff, zoff), m_nbnd(nbnd) {}

  //! zero-offset constructor
  MeshvarBnd(size_t nbnd, size_t nx, size_t ny, size_t nz)
      : Meshvar<real_t>(nx + 2 * nbnd, ny + 2 * nbnd, nz + 2 * nbnd, 0, 0, 0), m_nbnd(nbnd) {}

  //! constructor for cubic meshes
  MeshvarBnd(size_t nbnd, size_t n, size_t xoff, size_t yoff, size_t zoff)
      : Meshvar<real_t>(n + 2 * nbnd, xoff, yoff, zoff), m_nbnd(nbnd) {}

  //! constructor for cubic meshes with zero offset
  MeshvarBnd(size_t nbnd, size_t n) : Meshvar<real_t>(n + 2 * nbnd, 0, 0, 0), m_nbnd(nbnd) {}

  //! modified copy constructor, allows to avoid copying actual data
  MeshvarBnd(const MeshvarBnd<real_t> &v, bool copyover) : Meshvar<real_t>(v, copyover), m_nbnd(v.m_nbnd) {}

  //! copy constructor
  explicit MeshvarBnd(const MeshvarBnd<real_t> &v) : Meshvar<real_t>(v, true), m_nbnd(v.m_nbnd) {}

  //! get extent of the mesh along a specified dimension
  inline size_t size(unsigned dim = 0) const {
    if (dim == 0)
      return m_nx - 2 * m_nbnd;
    if (dim == 1)
      return m_ny - 2 * m_nbnd;
    return m_nz - 2 * m_nbnd;
  }

  //! 3D random access to the data block via index 3-tuples
  inline real_t &operator()(const int ix, const int iy, const int iz) {
    size_t iix(ix + m_nbnd), iiy(iy + m_nbnd), iiz(iz + m_nbnd);
    return m_pdata[(iix * m_ny + iiy) * m_nz + iiz];
  }

  //! 3D random access to the data block via index 3-tuples (const)
  inline const real_t &operator()(const int ix, const int iy, const int iz) const {
    size_t iix(ix + m_nbnd), iiy(iy + m_nbnd), iiz(iz + m_nbnd);
    return m_pdata[(iix * m_ny + iiy) * m_nz + iiz];
  }

  //! assignment operator for rectangular meshes with ghost zones
  MeshvarBnd<real_t> &operator=(const MeshvarBnd<real_t> &m) {
    if (this->m_nx != m.m_nx || this->m_ny != m.m_ny || this->m_nz != m.m_nz) {
      this->m_nx = m.m_nx;
      this->m_ny = m.m_ny;
      this->m_nz = m.m_nz;

      if (m_pdata != NULL)
        delete[] m_pdata;

      m_pdata = new real_t[m_nx * m_ny * m_nz];
    }

    for (size_t i = 0; i < m_nx * m_ny * m_nz; ++i)
      this->m_pdata[i] = m.m_pdata[i];

    return *this;
  }

  //! sets the value of all ghost zones to zero
  void zero_bnd(void) {

    int nx, ny, nz;
    nx = this->size(0);
    ny = this->size(1);
    nz = this->size(2);

    for (int j = -m_nbnd; j < ny + m_nbnd; ++j)
      for (int k = -m_nbnd; k < nz + m_nbnd; ++k) {
        for (int i = -m_nbnd; i < 0; ++i) {
          (*this)(i, j, k) = 0.0;
          (*this)(nx - 1 - i, j, k) = 0.0;
        }
      }

    for (int i = -m_nbnd; i < nx + m_nbnd; ++i)
      for (int k = -m_nbnd; k < nz + m_nbnd; ++k) {
        for (int j = -m_nbnd; j < 0; ++j) {
          (*this)(i, j, k) = 0.0;
          (*this)(i, ny - j - 1, k) = 0.0;
        }
      }

    for (int i = -m_nbnd; i < nx + m_nbnd; ++i)
      for (int j = -m_nbnd; j < ny + m_nbnd; ++j) {
        for (int k = -m_nbnd; k < 0; ++k) {
          (*this)(i, j, k) = 0.0;
          (*this)(i, j, nz - k - 1) = 0.0;
        }
      }
  }

  //! outputs the data, for debugging only, not practical for large datasets
  void print(void) const {
    int nbnd = m_nbnd;

    std::cout << "size is [" << this->size(0) << ", " << this->size(1) << ", " << this->size(2) << "]\n";
    std::cout << "ghost region has length of " << nbnd << std::endl;

    std::cout.precision(3);
    for (int i = -nbnd; i < (int)this->size(0) + nbnd; ++i) {
      std::cout << "ix = " << i << ": \n";

      for (int j = -nbnd; j < (int)this->size(1) + nbnd; ++j) {
        for (int k = -nbnd; k < (int)this->size(2) + nbnd; ++k) {
          if (i < 0 || i >= this->size(0) || j < 0 || j >= this->size(1) || k < 0 || k >= this->size(2))
            std::cout << "[" << std::setw(6) << (*this)(i, j, k) << "] ";
          else
            std::cout << std::setw(8) << (*this)(i, j, k) << " ";
        }
        std::cout << std::endl;
      }

      std::cout << std::endl;
    }
  }
};

//! class that subsumes a nested grid collection
template <typename T> class GridHierarchy {
public:
  //! number of ghost cells on boundary
  size_t m_nbnd;

  //! highest level without adaptive refinement
  unsigned m_levelmin;

  //! vector of pointers to the underlying rectangular mesh data for each level
  std::vector<MeshvarBnd<T> *> m_pgrids;

  std::vector<int> m_xoffabs, //!< vector of x-offsets of a level mesh relative to the coarser level
      m_yoffabs,              //!< vector of x-offsets of a level mesh relative to the coarser level
      m_zoffabs;              //!< vector of x-offsets of a level mesh relative to the coarser level

  std::vector<refinement_mask *> m_ref_masks;
  bool bhave_refmask;

protected:
  //! check whether a given grid has identical hierarchy, dimensions to this
  bool is_consistent(const GridHierarchy<T> &gh) {
    if (gh.levelmax() != levelmax())
      return false;

    if (gh.levelmin() != levelmin())
      return false;

    for (unsigned i = levelmin(); i <= levelmax(); ++i)
      for (int j = 0; j < 3; ++j) {
        if (size(i, j) != gh.size(i, j))
          return false;
        if (offset(i, j) != gh.offset(i, j))
          return false;
      }

    return true;
  }

public:
  //! return a pointer to the MeshvarBnd object representing data for one level
  MeshvarBnd<T> *get_grid(unsigned ilevel) {

    if (ilevel >= m_pgrids.size()) {
      LOGERR("Attempt to access level %d but maxlevel = %d", ilevel, m_pgrids.size() - 1);
      throw std::runtime_error("Fatal: attempt to access non-existent grid");
    }
    return m_pgrids[ilevel];
  }

  //! return a pointer to the MeshvarBnd object representing data for one level (const)
  const MeshvarBnd<T> *get_grid(unsigned ilevel) const {
    if (ilevel >= m_pgrids.size()) {
      LOGERR("Attempt to access level %d but maxlevel = %d", ilevel, m_pgrids.size() - 1);
      throw std::runtime_error("Fatal: attempt to access non-existent grid");
    }

    return m_pgrids[ilevel];
  }

  //! constructor for a collection of rectangular grids representing a multi-level hierarchy
  /*! creates an empty hierarchy, levelmin is initially zero, no grids are stored
   * @param nbnd number of ghost zones added at the boundary
   */
  explicit GridHierarchy(size_t nbnd) : m_nbnd(nbnd), m_levelmin(0), bhave_refmask(false) { m_pgrids.clear(); }

  //! copy constructor
  explicit GridHierarchy(const GridHierarchy<T> &gh) {
    for (unsigned i = 0; i <= gh.levelmax(); ++i)
      m_pgrids.push_back(new MeshvarBnd<T>(*gh.get_grid(i)));

    m_nbnd = gh.m_nbnd;
    m_levelmin = gh.m_levelmin;

    m_xoffabs = gh.m_xoffabs;
    m_yoffabs = gh.m_yoffabs;
    m_zoffabs = gh.m_zoffabs;

    // ref_mask   = gh.ref_mask;
    bhave_refmask = gh.bhave_refmask;

    if (bhave_refmask) {
      for (size_t i = 0; i < gh.m_ref_masks.size(); ++i)
        m_ref_masks.push_back(new refinement_mask(*(gh.m_ref_masks[i])));
    }
  }

  //! destructor
  ~GridHierarchy() { this->deallocate(); }

  //! free all memory occupied by the grid hierarchy
  void deallocate() {
    for (unsigned i = 0; i < m_pgrids.size(); ++i)
      delete m_pgrids[i];
    m_pgrids.clear();
    std::vector<MeshvarBnd<T> *>().swap(m_pgrids);

    m_xoffabs.clear();
    m_yoffabs.clear();
    m_zoffabs.clear();
    m_levelmin = 0;

    for (size_t i = 0; i < m_ref_masks.size(); ++i)
      delete m_ref_masks[i];
    m_ref_masks.clear();
  }

  // meaning of the mask:
  //  -1  =  outside of mask
  //  0.5 =  in mask and refined (i.e. cell exists also on finer level)
  //  1   =  in mask and not refined (i.e. cell exists only on this level)

  void add_refinement_mask(const double *shift) {
    bhave_refmask = false;

    //! generate a mask
    if (m_levelmin != levelmax()) {
      for (int ilevel = (int)levelmax(); ilevel >= (int)levelmin(); --ilevel) {
        double xq[3], dx = 1.0 / (1ul << ilevel);

        m_ref_masks[ilevel]->init(size(ilevel, 0), size(ilevel, 1), size(ilevel, 2), 0);

        for (size_t i = 0; i < size(ilevel, 0); i += 2) {
          xq[0] = (offset_abs(ilevel, 0) + i) * dx + 0.5 * dx + shift[0];
          for (size_t j = 0; j < size(ilevel, 1); j += 2) {
            xq[1] = (offset_abs(ilevel, 1) + j) * dx + 0.5 * dx + shift[1];
            for (size_t k = 0; k < size(ilevel, 2); k += 2) {
              xq[2] = (offset_abs(ilevel, 2) + k) * dx + 0.5 * dx + shift[2];

              short mask_val = -1; // outside mask
              if (the_region_generator->query_point(xq, ilevel) || ilevel == (int)levelmin())
                mask_val = 1; // inside mask

              (*m_ref_masks[ilevel])(i + 0, j + 0, k + 0) = mask_val;
              (*m_ref_masks[ilevel])(i + 0, j + 0, k + 1) = mask_val;
              (*m_ref_masks[ilevel])(i + 0, j + 1, k + 0) = mask_val;
              (*m_ref_masks[ilevel])(i + 0, j + 1, k + 1) = mask_val;
              (*m_ref_masks[ilevel])(i + 1, j + 0, k + 0) = mask_val;
              (*m_ref_masks[ilevel])(i + 1, j + 0, k + 1) = mask_val;
              (*m_ref_masks[ilevel])(i + 1, j + 1, k + 0) = mask_val;
              (*m_ref_masks[ilevel])(i + 1, j + 1, k + 1) = mask_val;
            }
          }
        }
      }

      bhave_refmask = true;

      for (int ilevel = (int)levelmin(); ilevel < (int)levelmax(); ++ilevel) {
        for (size_t i = 0; i < size(ilevel, 0); i++)
          for (size_t j = 0; j < size(ilevel, 1); j++)
            for (size_t k = 0; k < size(ilevel, 2); k++) {
              bool fine_is_flagged = false;

              int ifine[] = {
                  2 * (int)i - 2 * (int)offset(ilevel + 1, 0), 2 * (int)j - 2 * (int)offset(ilevel + 1, 1),
                  2 * (int)k - 2 * (int)offset(ilevel + 1, 2),
              };

              if (ifine[0] >= 0 && ifine[0] < (int)size(ilevel + 1, 0) && ifine[1] >= 0 &&
                  ifine[1] < (int)size(ilevel + 1, 1) && ifine[2] >= 0 && ifine[2] < (int)size(ilevel + 1, 2)) {
                fine_is_flagged |= (*m_ref_masks[ilevel + 1])(ifine[0] + 0, ifine[1] + 0, ifine[2] + 0) > 0;
                fine_is_flagged |= (*m_ref_masks[ilevel + 1])(ifine[0] + 0, ifine[1] + 0, ifine[2] + 1) > 0;
                fine_is_flagged |= (*m_ref_masks[ilevel + 1])(ifine[0] + 0, ifine[1] + 1, ifine[2] + 0) > 0;
                fine_is_flagged |= (*m_ref_masks[ilevel + 1])(ifine[0] + 0, ifine[1] + 1, ifine[2] + 1) > 0;
                fine_is_flagged |= (*m_ref_masks[ilevel + 1])(ifine[0] + 1, ifine[1] + 0, ifine[2] + 0) > 0;
                fine_is_flagged |= (*m_ref_masks[ilevel + 1])(ifine[0] + 1, ifine[1] + 0, ifine[2] + 1) > 0;
                fine_is_flagged |= (*m_ref_masks[ilevel + 1])(ifine[0] + 1, ifine[1] + 1, ifine[2] + 0) > 0;
                fine_is_flagged |= (*m_ref_masks[ilevel + 1])(ifine[0] + 1, ifine[1] + 1, ifine[2] + 1) > 0;

                if (fine_is_flagged) {
                  (*m_ref_masks[ilevel])(i, j, k) = 2; // cell is refined

                  (*m_ref_masks[ilevel + 1])(ifine[0] + 0, ifine[1] + 0, ifine[2] + 0) = 1;
                  (*m_ref_masks[ilevel + 1])(ifine[0] + 0, ifine[1] + 0, ifine[2] + 1) = 1;
                  (*m_ref_masks[ilevel + 1])(ifine[0] + 0, ifine[1] + 1, ifine[2] + 0) = 1;
                  (*m_ref_masks[ilevel + 1])(ifine[0] + 0, ifine[1] + 1, ifine[2] + 1) = 1;
                  (*m_ref_masks[ilevel + 1])(ifine[0] + 1, ifine[1] + 0, ifine[2] + 0) = 1;
                  (*m_ref_masks[ilevel + 1])(ifine[0] + 1, ifine[1] + 0, ifine[2] + 1) = 1;
                  (*m_ref_masks[ilevel + 1])(ifine[0] + 1, ifine[1] + 1, ifine[2] + 0) = 1;
                  (*m_ref_masks[ilevel + 1])(ifine[0] + 1, ifine[1] + 1, ifine[2] + 1) = 1;
                }
              }
            }
      }
    }
  }

  //! get offset of a grid at specified refinement level
  /*! the offset describes the shift of a refinement grid with respect to its coarser parent grid
   *  @param ilevel the level for which the offset is to be determined
   *  @param idim the dimension along which the offset is to be determined
   *  @return integer value denoting the offset in units of coarse grid cells
   *  @sa offset_abs
   */
  int offset(int ilevel, int idim) const { return m_pgrids[ilevel]->offset(idim); }

  //! get size of a grid at specified refinement level
  /*! the size describes the number of cells along one dimension of a grid
   *  @param ilevel the level for which the size is to be determined
   *  @param idim the dimension along which the size is to be determined
   *  @return integer value denoting the size of refinement grid at level ilevel along dimension idim
   */
  size_t size(int ilevel, int idim) const { return m_pgrids[ilevel]->size(idim); }

  //! get the absolute offset of a grid at specified refinement level
  /*! the absolute offset describes the shift of a refinement grid with respect to the simulation volume
   *  @param ilevel the level for which the offset is to be determined
   *  @param idim the dimension along which the offset is to be determined
   *  @return integer value denoting the absolute offset in units of fine grid cells
   *  @sa offset
   */
  int offset_abs(int ilevel, int idim) const {
    if (idim == 0)
      return m_xoffabs[ilevel];
    if (idim == 1)
      return m_yoffabs[ilevel];
    return m_zoffabs[ilevel];
  }

  //! get the coordinate posisition of a grid cell
  /*! returns the position of a grid cell at specified level relative to the simulation volume
   *  @param ilevel the refinement level of the grid cell
   *  @param i the x-index of the cell in the level grid
   *  @param j the y-index of the cell in the level grid
   *  @param k the z-index of the cell in the level grid
   *  @param ppos pointer to a double[3] array to which the coordinates are written
   *  @return none
   */
  void cell_pos(unsigned ilevel, int i, int j, int k, double *ppos) const {
    double h = 1.0 / (1 << ilevel); //, htop = h*2.0;
    ppos[0] = h * ((double)offset_abs(ilevel, 0) + (double)i + 0.5);
    ppos[1] = h * ((double)offset_abs(ilevel, 1) + (double)j + 0.5);
    ppos[2] = h * ((double)offset_abs(ilevel, 2) + (double)k + 0.5);

    if (ppos[0] >= 1.0 || ppos[1] >= 1.0 || ppos[2] >= 1.0)
      std::cerr << " - Cell seems outside domain! : (" << ppos[0] << ", " << ppos[1] << ", " << ppos[2] << "\n";
  }

  //! get the bounding box of a grid in code units
  /*! returns the bounding box of a grid at specified level in code units
   *  @param ilevel the refinement level of the grid
   *  @param left pointer to a double[3] array to which the left corner is written
   *  @param right pointer to a double[3] array to which the right corner is written
   *  @return none
   */
  void grid_bbox(unsigned ilevel, double *left, double *right) const {
    double h = 1.0 / (1 << ilevel); //, htop = h*2.0;
    left[0] = h * ((double)offset_abs(ilevel, 0));
    left[1] = h * ((double)offset_abs(ilevel, 1));
    left[2] = h * ((double)offset_abs(ilevel, 2));

    right[0] = left[0] + h * ((double)size(ilevel, 0));
    right[1] = left[1] + h * ((double)size(ilevel, 1));
    right[2] = left[2] + h * ((double)size(ilevel, 2));
  }

  //! checks whether a given grid cell is refined
  /*! a grid cell counts as refined if it is divided into 8 cells at the next higher level
   *  @param ilevel the refinement level of the grid cell
   *  @param i the x-index of the cell in the level grid
   *  @param j the y-index of the cell in the level grid
   *  @param k the z-index of the cell in the level grid
   *  @return true if cell is refined, false otherwise
   */
  bool is_refined(unsigned ilevel, int i, int j, int k) const {
    // meaning of the mask:
    //  -1  =  outside of mask
    //  2 =  in mask and refined (i.e. cell exists also on finer level)
    //  1   =  in mask and not refined (i.e. cell exists only on this level)

    if (bhave_refmask) {
      return (*m_ref_masks[ilevel])(i, j, k) == 2;
    }

    if (!bhave_refmask && ilevel == levelmax())
      return false;

    if (i < offset(ilevel + 1, 0) || i >= offset(ilevel + 1, 0) + (int)size(ilevel + 1, 0) / 2 ||
        j < offset(ilevel + 1, 1) || j >= offset(ilevel + 1, 1) + (int)size(ilevel + 1, 1) / 2 ||
        k < offset(ilevel + 1, 2) || k >= offset(ilevel + 1, 2) + (int)size(ilevel + 1, 2) / 2)
      return false;

    return true;
  }

  bool is_in_mask(unsigned ilevel, int i, int j, int k) const {
    // meaning of the mask:
    //  -1  =  outside of mask
    //  2 =  in mask and refined (i.e. cell exists also on finer level)
    //  1   =  in mask and not refined (i.e. cell exists only on this level)

    if (bhave_refmask) {
      return ((*m_ref_masks[ilevel])(i, j, k) >= 0);
    }

    return true;
  }

  //! sets the values of all grids on all levels to zero
  void zero(void) {
    for (unsigned i = 0; i < m_pgrids.size(); ++i)
      m_pgrids[i]->zero();
  }

  //! count the number of cells that are not further refined (=leafs)
  /*! for allocation purposes it is useful to query the number of cells to be expected
   *  @param lmin the minimum refinement level to consider
   *  @param lmax the maximum refinement level to consider
   *  @return the integer number of cells between lmin and lmax that are not further refined
   */
  size_t count_leaf_cells(unsigned lmin, unsigned lmax) const {
    size_t npcount = 0;

    for (int ilevel = lmax; ilevel >= (int)lmin; --ilevel)
      for (unsigned i = 0; i < get_grid(ilevel)->size(0); ++i)
        for (unsigned j = 0; j < get_grid(ilevel)->size(1); ++j)
          for (unsigned k = 0; k < get_grid(ilevel)->size(2); ++k)
            if (is_in_mask(ilevel, i, j, k) && !is_refined(ilevel, i, j, k))
              ++npcount;

    return npcount;
  }

  //! count the number of cells that are not further refined (=leafs)
  /*! for allocation purposes it is useful to query the number of cells to be expected
   *  @return the integer number of cells in the hierarchy that are not further refined
   */
  size_t count_leaf_cells(void) const { return count_leaf_cells(levelmin(), levelmax()); }

  //! creates a hierarchy of coextensive grids, refined by factors of 2
  /*! creates a hierarchy of lmax grids, each extending over the whole simulation volume with
   *  grid length 2^n for level 0<=n<=lmax
   *  @param lmax the maximum refinement level to be added (sets the resolution to 2^lmax for each dim)
   *  @return none
   */
  void create_base_hierarchy(unsigned lmax) {
    size_t n = 1;

    this->deallocate();

    m_pgrids.clear();

    m_xoffabs.clear();
    m_yoffabs.clear();
    m_zoffabs.clear();

    for (unsigned i = 0; i <= lmax; ++i) {
      // std::cout << "....adding level " << i << " (" << n << ", " << n << ", " << n << ")" << std::endl;
      m_pgrids.push_back(new MeshvarBnd<T>(m_nbnd, n, n, n, 0, 0, 0));
      m_pgrids[i]->zero();
      n *= 2;

      m_xoffabs.push_back(0);
      m_yoffabs.push_back(0);
      m_zoffabs.push_back(0);
    }

    m_levelmin = lmax;

    for (unsigned i = 0; i <= lmax; ++i)
      m_ref_masks.push_back(new refinement_mask(size(i, 0), size(i, 1), size(i, 2), (short)(i != lmax)));
  }

  //! multiply entire grid hierarchy by a constant
  GridHierarchy<T> &operator*=(T x) {
    for (unsigned i = 0; i < m_pgrids.size(); ++i)
      (*m_pgrids[i]) *= x;
    return *this;
  }

  //! divide entire grid hierarchy by a constant
  GridHierarchy<T> &operator/=(T x) {
    for (unsigned i = 0; i < m_pgrids.size(); ++i)
      (*m_pgrids[i]) /= x;
    return *this;
  }

  //! add a constant to the entire grid hierarchy
  GridHierarchy<T> &operator+=(T x) {
    for (unsigned i = 0; i < m_pgrids.size(); ++i)
      (*m_pgrids[i]) += x;
    return *this;
  }

  //! subtract a constant from the entire grid hierarchy
  GridHierarchy<T> &operator-=(T x) {
    for (unsigned i = 0; i < m_pgrids.size(); ++i)
      (*m_pgrids[i]) -= x;
    return *this;
  }

  //! multiply (element-wise) two grid hierarchies
  GridHierarchy<T> &operator*=(const GridHierarchy<T> &gh) {
    if (!is_consistent(gh)) {
      LOGERR("GridHierarchy::operator*= : attempt to operate on incompatible data");
      throw std::runtime_error("GridHierarchy::operator*= : attempt to operate on incompatible data");
    }
    for (unsigned i = 0; i < m_pgrids.size(); ++i)
      (*m_pgrids[i]) *= *gh.get_grid(i);
    return *this;
  }

  //! divide (element-wise) two grid hierarchies
  GridHierarchy<T> &operator/=(const GridHierarchy<T> &gh) {
    if (!is_consistent(gh)) {
      LOGERR("GridHierarchy::operator/= : attempt to operate on incompatible data");
      throw std::runtime_error("GridHierarchy::operator/= : attempt to operate on incompatible data");
    }
    for (unsigned i = 0; i < m_pgrids.size(); ++i)
      (*m_pgrids[i]) /= *gh.get_grid(i);
    return *this;
  }

  //! add (element-wise) two grid hierarchies
  GridHierarchy<T> &operator+=(const GridHierarchy<T> &gh) {
    if (!is_consistent(gh))
      throw std::runtime_error("GridHierarchy::operator+= : attempt to operate on incompatible data");

    for (unsigned i = 0; i < m_pgrids.size(); ++i)
      (*m_pgrids[i]) += *gh.get_grid(i);
    return *this;
  }

  //! subtract (element-wise) two grid hierarchies
  GridHierarchy<T> &operator-=(const GridHierarchy<T> &gh) {
    if (!is_consistent(gh)) {
      LOGERR("GridHierarchy::operator-= : attempt to operate on incompatible data");
      throw std::runtime_error("GridHierarchy::operator-= : attempt to operate on incompatible data");
    }
    for (unsigned i = 0; i < m_pgrids.size(); ++i)
      (*m_pgrids[i]) -= *gh.get_grid(i);
    return *this;
  }

  //! assign (element-wise) two grid hierarchies
  GridHierarchy<T> &operator=(const GridHierarchy<T> &gh) {
    bhave_refmask = gh.bhave_refmask;

    if (bhave_refmask) {
      for (unsigned i = 0; i <= gh.levelmax(); ++i)
        m_ref_masks.push_back(new refinement_mask(*(gh.m_ref_masks[i])));
    }

    if (!is_consistent(gh)) {
      for (unsigned i = 0; i < m_pgrids.size(); ++i)
        delete m_pgrids[i];
      m_pgrids.clear();

      for (unsigned i = 0; i <= gh.levelmax(); ++i)
        m_pgrids.push_back(new MeshvarBnd<T>(*gh.get_grid(i)));
      m_levelmin = gh.levelmin();
      m_nbnd = gh.m_nbnd;

      m_xoffabs = gh.m_xoffabs;
      m_yoffabs = gh.m_yoffabs;
      m_zoffabs = gh.m_zoffabs;

      return *this;
    } // throw std::runtime_error("GridHierarchy::operator= : attempt to operate on incompatible data");

    for (unsigned i = 0; i < m_pgrids.size(); ++i)
      (*m_pgrids[i]) = *gh.get_grid(i);
    return *this;
  }

  /*
  //! assignment operator
  GridHierarchy& operator=( const GridHierarchy<T>& gh )
  {
          for( unsigned i=0; i<m_pgrids.size(); ++i )
                  delete m_pgrids[i];
          m_pgrids.clear();

          for( unsigned i=0; i<=gh.levelmax(); ++i )
                  m_pgrids.push_back( new MeshvarBnd<T>( *gh.get_grid(i) ) );
          m_levelmin = gh.levelmin();
          m_nbnd = gh.m_nbnd;

          m_xoffabs = gh.m_xoffabs;
          m_yoffabs = gh.m_yoffabs;
          m_zoffabs = gh.m_zoffabs;


          return *this;
  }
  */

  /*! add a new refinement patch to the so-far finest level
   * @param xoff x-offset in units of the coarser grid (finest grid before adding new patch)
   * @param yoff y-offset in units of the coarser grid (finest grid before adding new patch)
   * @param zoff z-offset in units of the coarser grid (finest grid before adding new patch)
   * @param nx x-extent in fine grid cells
   * @param ny y-extent in fine grid cells
   * @param nz z-extent in fine grid cells
   */
  void add_patch(unsigned xoff, unsigned yoff, unsigned zoff, unsigned nx, unsigned ny, unsigned nz) {
    m_pgrids.push_back(new MeshvarBnd<T>(m_nbnd, nx, ny, nz, xoff, yoff, zoff));
    m_pgrids.back()->zero();

    //.. add absolute offsets (in units of current level grid cells)
    m_xoffabs.push_back(2 * (m_xoffabs.back() + xoff));
    m_yoffabs.push_back(2 * (m_yoffabs.back() + yoff));
    m_zoffabs.push_back(2 * (m_zoffabs.back() + zoff));

    m_ref_masks.push_back(new refinement_mask(nx, ny, nz, 0));
  }

  /*! cut a refinement patch to the specified size
   * @param ilevel grid level on which to perform the size adjustment
   * @param xoff new x-offset in units of the coarser grid (finest grid before adding new patch)
   * @param yoff new y-offset in units of the coarser grid (finest grid before adding new patch)
   * @param zoff new z-offset in units of the coarser grid (finest grid before adding new patch)
   * @param nx new x-extent in fine grid cells
   * @param ny new y-extent in fine grid cells
   * @param nz new z-extent in fine grid cells
   */
  void cut_patch(unsigned ilevel, unsigned xoff, unsigned yoff, unsigned zoff, unsigned nx, unsigned ny, unsigned nz) {
    unsigned dx, dy, dz, dxtop, dytop, dztop;

    dx = xoff - m_xoffabs[ilevel];
    dy = yoff - m_yoffabs[ilevel];
    dz = zoff - m_zoffabs[ilevel];

    assert(dx % 2 == 0 && dy % 2 == 0 && dz % 2 == 0);

    dxtop = m_pgrids[ilevel]->offset(0) + dx / 2;
    dytop = m_pgrids[ilevel]->offset(1) + dy / 2;
    dztop = m_pgrids[ilevel]->offset(2) + dz / 2;

    MeshvarBnd<T> *mnew = new MeshvarBnd<T>(m_nbnd, nx, ny, nz, dxtop, dytop, dztop);

    //... copy data
    for (unsigned i = 0; i < nx; ++i)
      for (unsigned j = 0; j < ny; ++j)
        for (unsigned k = 0; k < nz; ++k)
          (*mnew)(i, j, k) = (*m_pgrids[ilevel])(i + dx, j + dy, k + dz);

    //... replace in hierarchy
    delete m_pgrids[ilevel];
    m_pgrids[ilevel] = mnew;

    //... update offsets
    m_xoffabs[ilevel] += dx;
    m_yoffabs[ilevel] += dy;
    m_zoffabs[ilevel] += dz;

    if (ilevel < levelmax()) {
      m_pgrids[ilevel + 1]->offset(0) -= dx;
      m_pgrids[ilevel + 1]->offset(1) -= dy;
      m_pgrids[ilevel + 1]->offset(2) -= dz;
    }

    find_new_levelmin();
  }

  void cut_patch_enforce_top_density(unsigned ilevel, unsigned xoff, unsigned yoff, unsigned zoff, unsigned nx,
                                     unsigned ny, unsigned nz) {
    unsigned dx, dy, dz, dxtop, dytop, dztop;

    dx = xoff - m_xoffabs[ilevel];
    dy = yoff - m_yoffabs[ilevel];
    dz = zoff - m_zoffabs[ilevel];

    assert(dx % 2 == 0 && dy % 2 == 0 && dz % 2 == 0);

    dxtop = m_pgrids[ilevel]->offset(0) + dx / 2;
    dytop = m_pgrids[ilevel]->offset(1) + dy / 2;
    dztop = m_pgrids[ilevel]->offset(2) + dz / 2;

    MeshvarBnd<T> *mnew = new MeshvarBnd<T>(m_nbnd, nx, ny, nz, dxtop, dytop, dztop);

    double coarsesum = 0.0, finesum = 0.0;
    size_t coarsecount = 0, finecount = 0;

    //... copy data
    for (unsigned i = 0; i < nx; ++i)
      for (unsigned j = 0; j < ny; ++j)
        for (unsigned k = 0; k < nz; ++k) {
          (*mnew)(i, j, k) = (*m_pgrids[ilevel])(i + dx, j + dy, k + dz);
          finesum += (*mnew)(i, j, k);
          finecount++;
        }

    //... replace in hierarchy
    delete m_pgrids[ilevel];
    m_pgrids[ilevel] = mnew;

    //... update offsets
    m_xoffabs[ilevel] += dx;
    m_yoffabs[ilevel] += dy;
    m_zoffabs[ilevel] += dz;

    if (ilevel < levelmax()) {
      m_pgrids[ilevel + 1]->offset(0) -= dx;
      m_pgrids[ilevel + 1]->offset(1) -= dy;
      m_pgrids[ilevel + 1]->offset(2) -= dz;
    }

    //... enforce top mean density over same patch
    if (ilevel > levelmin()) {
      int ox = m_pgrids[ilevel]->offset(0);
      int oy = m_pgrids[ilevel]->offset(1);
      int oz = m_pgrids[ilevel]->offset(2);

      for (unsigned i = 0; i < nx / 2; ++i)
        for (unsigned j = 0; j < ny / 2; ++j)
          for (unsigned k = 0; k < nz / 2; ++k) {
            coarsesum += (*m_pgrids[ilevel - 1])(i + ox, j + oy, k + oz);
            coarsecount++;
          }

      coarsesum /= (double)coarsecount;
      finesum /= (double)finecount;

      for (unsigned i = 0; i < nx; ++i)
        for (unsigned j = 0; j < ny; ++j)
          for (unsigned k = 0; k < nz; ++k)
            (*mnew)(i, j, k) += (coarsesum - finesum);

      LOGINFO("level %d : corrected patch overlap mean density by %f", ilevel, coarsesum - finesum);
    }

    find_new_levelmin();
  }

  /*! determine level for which grid extends over entire domain */
  void find_new_levelmin(void) {
    for (unsigned i = 0; i <= levelmax(); ++i) {
      unsigned n = 1 << i;
      if (m_pgrids[i]->size(0) == n && m_pgrids[i]->size(1) == n && m_pgrids[i]->size(2) == n) {
        m_levelmin = i;
      }
    }
  }

  //! return maximum level in refinement hierarchy
  unsigned levelmax(void) const { return m_pgrids.size() - 1; }

  //! return minimum level in refinement hierarchy (the level which extends over the entire domain)
  unsigned levelmin(void) const { return m_levelmin; }
};

//! class that computes the refinement structure given parameters
class refinement_hierarchy {
  std::vector<double> x0_, //!< x-coordinates of grid origins (in [0..1[)
      y0_,                 //!< y-coordinates of grid origins (in [0..1[)
      z0_,                 //!< z-coordinates of grid origins (in [0..1[)
      xl_,                 //!< x-extent of grids (in [0..1[)
      yl_,                 //!< y-extent of grids (in [0..1[)
      zl_;                 //!< z-extent of grids (in [0..1[)

  std::vector<unsigned> ox_, //!< relative x-coordinates of grid origins (in coarser grid cells)
      oy_,                   //!< relative y-coordinates of grid origins (in coarser grid cells)
      oz_,                   //!< relative z-coordinates of grid origins (in coarser grid cells)
      oax_,                  //!< absolute x-coordinates of grid origins (in fine grid cells)
      oay_,                  //!< absolute y-coordinates of grid origins (in fine grid cells)
      oaz_,                  //!< absolute z-coordinates of grid origins (in fine grid cells)
      nx_,                   //!< x-extent of grids (in fine grid cells)
      ny_,                   //!< y-extent of grids (in fine grid cells)
      nz_;                   //!< z-extent of grids (in fine grid cells)

  unsigned levelmin_, //!< minimum grid level for Poisson solver
      levelmax_,      //!< maximum grid level for all operations
      levelmin_tf_,   //!< minimum grid level for density calculation
      padding_,       //!< padding in number of coarse cells between refinement levels
      blocking_factor_, gridding_unit_;

  config_file &cf_; //!< reference to config_file

  bool align_top_,   //!< bool whether to align all grids with coarsest grid cells
      equal_extent_; //!< bool whether the simulation code requires squared refinement regions (e.g. RAMSES)

  double x0ref_[3], //!< coordinates of refinement region origin (in [0..1[)
      lxref_[3];    //!< extent of refinement region (int [0..1[)

  size_t lnref_[3];
  bool bhave_nref;

  int xshift_[3]; //!< shift of refinement region in coarse cells (in order to center it in the domain)
  double rshift_[3];


  //! calculates the greatest common divisor
  int gcd(int a, int b) const {
    return b == 0 ? a : gcd(b, a % b);
  }

  // calculates the cell shift in units of levelmin grid cells if there is an additional constraint to be
  // congruent with another grid that partitions the same space in multiples of "base_unit"
  int get_shift_unit( int base_unit, int levelmin ) const {
    /*int Lp = 0;
    while( base_unit * (1<<Lp) < 1<<(levelmin+1) ){
      ++Lp;
    }
    int U = base_unit * (1<<Lp);

    return std::max<int>( 1, (1<<(levelmin+1)) / (2*gcd(U,1<<(levelmin+1) )) );*/

    int level_m = 0;
    while( base_unit * (1<<level_m) < (1<<levelmin) )
      ++level_m;

    return std::max<int>( 1, (1<<levelmin)/gcd(base_unit * (1<<level_m),(1<<levelmin)) );

  }

public:
  //! copy constructor
  refinement_hierarchy(const refinement_hierarchy &rh) : cf_(rh.cf_) { *this = rh; }

  //! constructor from a config_file holding information about the desired refinement geometry
  explicit refinement_hierarchy(config_file &cf) : cf_(cf) {
    //... query the parameter data we need
    levelmin_ = cf_.getValue<unsigned>("setup", "levelmin");
    levelmax_ = cf_.getValue<unsigned>("setup", "levelmax");
    levelmin_tf_ = cf_.getValueSafe<unsigned>("setup", "levelmin_TF", levelmin_);
    align_top_ = cf_.getValueSafe<bool>("setup", "align_top", false);
    equal_extent_ = cf_.getValueSafe<bool>("setup", "force_equal_extent", false);
    blocking_factor_ = cf.getValueSafe<unsigned>("setup", "blocking_factor", 0);
    gridding_unit_ = cf.getValueSafe<unsigned>("setup", "gridding_unit", 2);

    if (gridding_unit_ != 2) {
      blocking_factor_ = gridding_unit_;
    }

    bool bnoshift = cf_.getValueSafe<bool>("setup", "no_shift", false);
    bool force_shift = cf_.getValueSafe<bool>("setup", "force_shift", false);

    //... call the region generator
    if (levelmin_ != levelmax_) {

      double x1ref[3];
      the_region_generator->get_AABB(x0ref_, x1ref, levelmax_);
      for (int i = 0; i < 3; ++i)
        lxref_[i] = x1ref[i] - x0ref_[i];
      bhave_nref = false;

      std::string region_type = cf.getValueSafe<std::string>("setup", "region", "box");

      LOGINFO("refinement region is \'%s\', w/ bounding box\n        left = [%f,%f,%f]\n       right = [%f,%f,%f]",
              region_type.c_str(), x0ref_[0], x0ref_[1], x0ref_[2], x1ref[0], x1ref[1], x1ref[2]);

      bhave_nref = the_region_generator->is_grid_dim_forced(lnref_);
    }

    // if not doing any refinement levels, set extent to full box
    if (levelmin_ == levelmax_) {
      x0ref_[0] = 0.0;
      x0ref_[1] = 0.0;
      x0ref_[2] = 0.0;

      lxref_[0] = 1.0;
      lxref_[1] = 1.0;
      lxref_[2] = 1.0;
    }

    unsigned ncoarse = 1 << levelmin_;

    //... determine shift

    double xc[3];
    xc[0] = fmod(x0ref_[0] + 0.5 * lxref_[0], 1.0);
    xc[1] = fmod(x0ref_[1] + 0.5 * lxref_[1], 1.0);
    xc[2] = fmod(x0ref_[2] + 0.5 * lxref_[2], 1.0);

    if ((levelmin_ != levelmax_) && (!bnoshift || force_shift)) {
      int random_base_grid_unit = cf.getValueSafe<int>("random","base_unit",1);
      int shift_unit = get_shift_unit( random_base_grid_unit, levelmin_ );
      if( shift_unit != 1 ){
        LOGINFO("volume can only be shifted by multiples of %d coarse cells.",shift_unit);
      }

      xshift_[0] = (int)((0.5-xc[0]) * (double)ncoarse / shift_unit + 0.5) * shift_unit;//ARJ(int)((0.5 - xc[0]) * ncoarse);
      xshift_[1] = (int)((0.5-xc[1]) * (double)ncoarse / shift_unit + 0.5) * shift_unit;//ARJ(int)((0.5 - xc[1]) * ncoarse);
      xshift_[2] = (int)((0.5-xc[2]) * (double)ncoarse / shift_unit + 0.5) * shift_unit;//ARJ(int)((0.5 - xc[2]) * ncoarse);
    } else {
      xshift_[0] = 0;
      xshift_[1] = 0;
      xshift_[2] = 0;
    }

    char strtmp[32];
    sprintf(strtmp, "%d", xshift_[0]);
    cf_.insertValue("setup", "shift_x", strtmp);
    sprintf(strtmp, "%d", xshift_[1]);
    cf_.insertValue("setup", "shift_y", strtmp);
    sprintf(strtmp, "%d", xshift_[2]);
    cf_.insertValue("setup", "shift_z", strtmp);

    rshift_[0] = -(double)xshift_[0] / ncoarse;
    rshift_[1] = -(double)xshift_[1] / ncoarse;
    rshift_[2] = -(double)xshift_[2] / ncoarse;

    x0ref_[0] += (double)xshift_[0] / ncoarse;
    x0ref_[1] += (double)xshift_[1] / ncoarse;
    x0ref_[2] += (double)xshift_[2] / ncoarse;

    //... initialize arrays
    x0_.assign(levelmax_ + 1, 0.0);
    xl_.assign(levelmax_ + 1, 1.0);
    y0_.assign(levelmax_ + 1, 0.0);
    yl_.assign(levelmax_ + 1, 1.0);
    z0_.assign(levelmax_ + 1, 0.0);
    zl_.assign(levelmax_ + 1, 1.0);
    ox_.assign(levelmax_ + 1, 0);
    nx_.assign(levelmax_ + 1, 0);
    oy_.assign(levelmax_ + 1, 0);
    ny_.assign(levelmax_ + 1, 0);
    oz_.assign(levelmax_ + 1, 0);
    nz_.assign(levelmax_ + 1, 0);

    oax_.assign(levelmax_ + 1, 0);
    oay_.assign(levelmax_ + 1, 0);
    oaz_.assign(levelmax_ + 1, 0);

    nx_[levelmin_] = ncoarse;
    ny_[levelmin_] = ncoarse;
    nz_[levelmin_] = ncoarse;

    // set up base hierarchy sizes
    for (unsigned ilevel = 0; ilevel <= levelmin_; ++ilevel) {
      unsigned n = 1 << ilevel;

      xl_[ilevel] = yl_[ilevel] = zl_[ilevel] = 1.0;
      nx_[ilevel] = ny_[ilevel] = nz_[ilevel] = n;
    }

    // if no refinement, we can exit here
    if (levelmax_ == levelmin_)
      return;

    //... determine the position of the refinement region on the finest grid
    int il, jl, kl, ir, jr, kr;
    int nresmax = 1 << levelmax_;

    il = (int)(x0ref_[0] * nresmax);
    jl = (int)(x0ref_[1] * nresmax);
    kl = (int)(x0ref_[2] * nresmax);
    ir = (int)((x0ref_[0] + lxref_[0]) * nresmax); //+ 1.0);
    jr = (int)((x0ref_[1] + lxref_[1]) * nresmax); //+ 1.0);
    kr = (int)((x0ref_[2] + lxref_[2]) * nresmax); //+ 1.0);

    //... align with coarser grids ...
    if (align_top_) {
      //... require alignment with top grid
      unsigned nref = 1 << (levelmax_ - levelmin_ + 1);

      if (bhave_nref) {
        if (lnref_[0] % (1ul << (levelmax_ - levelmin_)) != 0 || lnref_[1] % (1ul << (levelmax_ - levelmin_)) != 0 ||
            lnref_[2] % (1ul << (levelmax_ - levelmin_)) != 0) {
          LOGERR("specified ref_dims and align_top=yes but cannot be aligned with coarse grid!");
          throw std::runtime_error("specified ref_dims and align_top=yes but cannot be aligned with coarse grid!");
        }
      }

      il = (int)((double)il / nref) * nref;
      jl = (int)((double)jl / nref) * nref;
      kl = (int)((double)kl / nref) * nref;

      int irr = (int)((double)ir / nref) * nref;
      int jrr = (int)((double)jr / nref) * nref;
      int krr = (int)((double)kr / nref) * nref;

      if (irr < ir)
        ir = (int)((double)ir / nref + 1.0) * nref;
      else
        ir = irr;

      if (jrr < jr)
        jr = (int)((double)jr / nref + 1.0) * nref;
      else
        jr = jrr;

      if (krr < kr)
        kr = (int)((double)kr / nref + 1.0) * nref;
      else
        kr = krr;

    } else {
      //... require alignment with coarser grid
      LOGINFO("Internal refinement bounding box error: [%d,%d]x[%d,%d]x[%d,%d]", il, ir, jl, jr, kl, kr);

      il -= il % gridding_unit_;
      jl -= jl % gridding_unit_;
      kl -= kl % gridding_unit_;

      ir = ((ir%gridding_unit_)!=0)? (ir/gridding_unit_ + 1)*gridding_unit_ : ir;
      jr = ((jr%gridding_unit_)!=0)? (jr/gridding_unit_ + 1)*gridding_unit_ : jr;
      kr = ((kr%gridding_unit_)!=0)? (kr/gridding_unit_ + 1)*gridding_unit_ : kr;

    }

    // if doing unigrid, set region to whole box
    if (levelmin_ == levelmax_) {
      il = jl = kl = 0;
      ir = jr = kr = nresmax - 1;
    }
    if (bhave_nref) {
      ir = il + lnref_[0];
      jr = jl + lnref_[1];
      kr = kl + lnref_[2];
    }

    //... make sure bounding box lies in domain
    il = (il + nresmax) % nresmax;
    ir = (ir + nresmax) % nresmax;
    jl = (jl + nresmax) % nresmax;
    jr = (jr + nresmax) % nresmax;
    kl = (kl + nresmax) % nresmax;
    kr = (kr + nresmax) % nresmax;

    if (il >= ir || jl >= jr || kl >= kr) {
      LOGERR("Internal refinement bounding box error: [%d,%d]x[%d,%d]x[%d,%d]", il, ir, jl, jr, kl, kr);
      throw std::runtime_error("refinement_hierarchy: Internal refinement bounding box error 1");
    }
    //... determine offsets
    if (levelmin_ != levelmax_) {

      oax_[levelmax_] = (il + nresmax) % nresmax;
      oay_[levelmax_] = (jl + nresmax) % nresmax;
      oaz_[levelmax_] = (kl + nresmax) % nresmax;
      nx_[levelmax_] = ir - il;
      ny_[levelmax_] = jr - jl;
      nz_[levelmax_] = kr - kl;

      if (equal_extent_) {

        if (bhave_nref && (lnref_[0] != lnref_[1] || lnref_[0] != lnref_[2])) {
          LOGERR("Specified equal_extent=yes conflicting with ref_dims which are not equal.");
          throw std::runtime_error("Specified equal_extent=yes conflicting with ref_dims which are not equal.");
        }
        size_t ilevel = levelmax_;
        size_t nmax = std::max(nx_[ilevel], std::max(ny_[ilevel], nz_[ilevel]));
        int dx = (int)((double)(nmax - nx_[ilevel]) * 0.5);
        int dy = (int)((double)(nmax - ny_[ilevel]) * 0.5);
        int dz = (int)((double)(nmax - nz_[ilevel]) * 0.5);

        oax_[ilevel] -= dx;
        oay_[ilevel] -= dy;
        oaz_[ilevel] -= dz;
        nx_[ilevel] = nmax;
        ny_[ilevel] = nmax;
        nz_[ilevel] = nmax;

        il = oax_[ilevel];
        jl = oay_[ilevel];
        kl = oaz_[ilevel];
        ir = il + nmax;
        jr = jl + nmax;
        kr = kl + nmax;
      }
    }

    padding_ = cf_.getValueSafe<unsigned>("setup", "padding", 8);

    //... determine position of coarser grids
    for (unsigned ilevel = levelmax_ - 1; ilevel > levelmin_; --ilevel) {
      il = (int)((double)il * 0.5 - padding_);
      jl = (int)((double)jl * 0.5 - padding_);
      kl = (int)((double)kl * 0.5 - padding_);

      ir = (int)((double)ir * 0.5 + padding_);
      jr = (int)((double)jr * 0.5 + padding_);
      kr = (int)((double)kr * 0.5 + padding_);

      //... align with coarser grids ...
      if (align_top_) {
        //... require alignment with top grid
        unsigned nref = 1 << (ilevel - levelmin_);

        il = (int)((double)il / nref) * nref;
        jl = (int)((double)jl / nref) * nref;
        kl = (int)((double)kl / nref) * nref;

        ir = (int)((double)ir / nref + 1.0) * nref;
        jr = (int)((double)jr / nref + 1.0) * nref;
        kr = (int)((double)kr / nref + 1.0) * nref;

      } else {
        //... require alignment with coarser grid
        /*il -= il % 2;
        jl -= jl % 2;
        kl -= kl % 2;
        ir += ir % 2;
        jr += jr % 2;
        kr += kr % 2;*/

        il -= il % gridding_unit_;
        jl -= jl % gridding_unit_;
        kl -= kl % gridding_unit_;

        ir = ((ir%gridding_unit_)!=0)? (ir/gridding_unit_ + 1)*gridding_unit_ : ir;
        jr = ((jr%gridding_unit_)!=0)? (jr/gridding_unit_ + 1)*gridding_unit_ : jr;
        kr = ((kr%gridding_unit_)!=0)? (kr/gridding_unit_ + 1)*gridding_unit_ : kr;
      }

      if (il >= ir || jl >= jr || kl >= kr || il < 0 || jl < 0 || kl < 0) {
        LOGERR("Internal refinement bounding box error: [%d,%d]x[%d,%d]x[%d,%d], level=%d", il, ir, jl, jr, kl, kr,
               ilevel);
        throw std::runtime_error("refinement_hierarchy: Internal refinement bounding box error 2");
      }
      oax_[ilevel] = il;
      oay_[ilevel] = jl;
      oaz_[ilevel] = kl;
      nx_[ilevel] = ir - il;
      ny_[ilevel] = jr - jl;
      nz_[ilevel] = kr - kl;

      if (blocking_factor_ > 0) {
        LOGINFO("applying blocking factor to enforce grid sizes");
        nx_[ilevel] += nx_[ilevel] % blocking_factor_;
        ny_[ilevel] += ny_[ilevel] % blocking_factor_;
        nz_[ilevel] += nz_[ilevel] % blocking_factor_;
      }

      if (equal_extent_) {
        size_t nmax = std::max(nx_[ilevel], std::max(ny_[ilevel], nz_[ilevel]));
        int dx = (int)((double)(nmax - nx_[ilevel]) * 0.5);
        int dy = (int)((double)(nmax - ny_[ilevel]) * 0.5);
        int dz = (int)((double)(nmax - nz_[ilevel]) * 0.5);

        oax_[ilevel] -= dx;
        oay_[ilevel] -= dy;
        oaz_[ilevel] -= dz;
        nx_[ilevel] = nmax;
        ny_[ilevel] = nmax;
        nz_[ilevel] = nmax;

        il = oax_[ilevel];
        jl = oay_[ilevel];
        kl = oaz_[ilevel];
        ir = il + nmax;
        jr = jl + nmax;
        kr = kl + nmax;
      }
    }

    //... determine relative offsets between grids
    for (unsigned ilevel = levelmax_; ilevel > levelmin_; --ilevel) {
      ox_[ilevel] = (oax_[ilevel] / 2 - oax_[ilevel - 1]);
      oy_[ilevel] = (oay_[ilevel] / 2 - oay_[ilevel - 1]);
      oz_[ilevel] = (oaz_[ilevel] / 2 - oaz_[ilevel - 1]);
    }

    //... do a forward sweep to ensure that absolute offsets are also correct now
    for (unsigned ilevel = levelmin_ + 1; ilevel <= levelmax_; ++ilevel) {
      oax_[ilevel] = 2 * oax_[ilevel - 1] + 2 * ox_[ilevel];
      oay_[ilevel] = 2 * oay_[ilevel - 1] + 2 * oy_[ilevel];
      oaz_[ilevel] = 2 * oaz_[ilevel - 1] + 2 * oz_[ilevel];
    }

    for (unsigned ilevel = levelmin_ + 1; ilevel <= levelmax_; ++ilevel) {
      double h = 1.0 / (1ul << ilevel);

      x0_[ilevel] = h * (double)oax_[ilevel];
      y0_[ilevel] = h * (double)oay_[ilevel];
      z0_[ilevel] = h * (double)oaz_[ilevel];

      xl_[ilevel] = h * (double)nx_[ilevel];
      yl_[ilevel] = h * (double)ny_[ilevel];
      zl_[ilevel] = h * (double)nz_[ilevel];
    }

    // do a consistency check that largest subgrid in zoom is not larger than half the box size
    for (unsigned ilevel = levelmin_ + 1; ilevel <= levelmax_; ++ilevel) {
      if (nx_[ilevel] > (1ul << (ilevel - 1)) || ny_[ilevel] > (1ul << (ilevel - 1)) ||
          nz_[ilevel] > (1ul << (ilevel - 1))) {
        LOGERR("On level %d, subgrid is larger than half the box. This is not allowed!", ilevel);
        throw std::runtime_error("Fatal: Subgrid larger than half boxin zoom.");
      }
    }

    // update the region generator with what has been actually created
    double left[3] = {x0_[levelmax_] + rshift_[0], y0_[levelmax_] + rshift_[1], z0_[levelmax_] + rshift_[2]};
    double right[3] = {left[0] + xl_[levelmax_], left[1] + yl_[levelmax_], left[2] + zl_[levelmax_]};
    the_region_generator->update_AABB(left, right);
  }

  //! asignment operator
  refinement_hierarchy &operator=(const refinement_hierarchy &o) {
    levelmin_ = o.levelmin_;
    levelmax_ = o.levelmax_;
    padding_ = o.padding_;
    cf_ = o.cf_;
    align_top_ = o.align_top_;
    for (int i = 0; i < 3; ++i) {
      x0ref_[i] = o.x0ref_[i];
      lxref_[i] = o.lxref_[i];
      xshift_[i] = o.xshift_[i];
      rshift_[i] = o.rshift_[i];
    }

    x0_ = o.x0_;
    y0_ = o.y0_;
    z0_ = o.z0_;
    xl_ = o.xl_;
    yl_ = o.yl_;
    zl_ = o.zl_;
    ox_ = o.ox_;
    oy_ = o.oy_;
    oz_ = o.oz_;
    oax_ = o.oax_;
    oay_ = o.oay_;
    oaz_ = o.oaz_;
    nx_ = o.nx_;
    ny_ = o.ny_;
    nz_ = o.nz_;

    return *this;
  }

  /*! cut a grid level to the specified size
   * @param ilevel grid level on which to perform the size adjustment
   * @param nx new x-extent in fine grid cells
   * @param ny new y-extent in fine grid cells
   * @param nz new z-extent in fine grid cells
   * @param oax new x-offset in units fine grid units
   * @param oay new y-offset in units fine grid units
   * @param oaz new z-offset in units fine grid units

   */
  void adjust_level(unsigned ilevel, int nx, int ny, int nz, int oax, int oay, int oaz) {
    double h = 1.0 / (1 << ilevel);

    int dx, dy, dz;

    dx = oax_[ilevel] - oax;
    dy = oay_[ilevel] - oay;
    dz = oaz_[ilevel] - oaz;

    ox_[ilevel] -= dx / 2;
    oy_[ilevel] -= dy / 2;
    oz_[ilevel] -= dz / 2;

    oax_[ilevel] = oax;
    oay_[ilevel] = oay;
    oaz_[ilevel] = oaz;

    nx_[ilevel] = nx;
    ny_[ilevel] = ny;
    nz_[ilevel] = nz;

    x0_[ilevel] = h * oax;
    y0_[ilevel] = h * oay;
    z0_[ilevel] = h * oaz;

    xl_[ilevel] = h * nx;
    yl_[ilevel] = h * ny;
    zl_[ilevel] = h * nz;

    if (ilevel < levelmax_) {
      ox_[ilevel + 1] += dx;
      oy_[ilevel + 1] += dy;
      oz_[ilevel + 1] += dz;
    }

    find_new_levelmin();
  }

  /*! determine level for which grid extends over entire domain */
  void find_new_levelmin(bool print = false) {
    unsigned old_levelmin(levelmin_);

    for (unsigned i = 0; i <= levelmax(); ++i) {
      unsigned n = 1 << i;

      if (oax_[i] == 0 && oay_[i] == 0 && oaz_[i] == 0 && nx_[i] == n && ny_[i] == n && nz_[i] == n) {
        levelmin_ = i;
      }
    }

    if ((old_levelmin != levelmin_) && print)
      LOGINFO("refinement_hierarchy: set new levelmin to %d", levelmin_);
  }

  //! get absolute grid offset for a specified level along a specified dimension (in fine grid units)
  unsigned offset_abs(unsigned ilevel, int dim) const {
    if (dim == 0)
      return oax_.at(ilevel);
    if (dim == 1)
      return oay_.at(ilevel);
    return oaz_.at(ilevel);
  }

  //! get relative grid offset for a specified level along a specified dimension (in coarser grid units)
  int offset(unsigned ilevel, int dim) const {
    if (dim == 0)
      return ox_.at(ilevel);
    if (dim == 1)
      return oy_.at(ilevel);
    return oz_.at(ilevel);
  }

  //! get grid size for a specified level along a specified dimension
  size_t size(unsigned ilevel, int dim) const {
    if (dim == 0)
      return nx_.at(ilevel);
    if (dim == 1)
      return ny_.at(ilevel);
    return nz_.at(ilevel);
  }

  //! get minimum grid level (the level for which the grid covers the entire domain)
  unsigned levelmin(void) const { return levelmin_; }

  //! get maximum grid level
  unsigned levelmax(void) const { return levelmax_; }

  //! get the total shift of the coordinate system in units of coarse cells
  int get_shift(int idim) const { return xshift_[idim]; }

  //! get the total shift of the coordinate system in box coordinates
  const double *get_coord_shift(void) const { return rshift_; }

  //! write refinement hierarchy to stdout
  void output(void) const {
    std::cout << "-------------------------------------------------------------\n";

    if (xshift_[0] != 0 || xshift_[1] != 0 || xshift_[2] != 0)
      std::cout << " - Domain will be shifted by (" << xshift_[0] << ", " << xshift_[1] << ", " << xshift_[2] << ")\n"
                << std::endl;

    std::cout << " - Grid structure:\n";

    for (unsigned ilevel = levelmin_; ilevel <= levelmax_; ++ilevel) {
      std::cout << "     Level " << std::setw(3) << ilevel << " :   offset = (" << std::setw(5) << ox_[ilevel] << ", "
                << std::setw(5) << oy_[ilevel] << ", " << std::setw(5) << oz_[ilevel] << ")\n"
                << "               offset_abs = (" << std::setw(5) << oax_[ilevel] << ", " << std::setw(5)
                << oay_[ilevel] << ", " << std::setw(5) << oaz_[ilevel] << ")\n"
                << "                   size   = (" << std::setw(5) << nx_[ilevel] << ", " << std::setw(5) << ny_[ilevel]
                << ", " << std::setw(5) << nz_[ilevel] << ")\n";
    }
    std::cout << "-------------------------------------------------------------\n";
  }

  void output_log(void) const {
    LOGUSER("   Domain shifted by      (%5d,%5d,%5d)", xshift_[0], xshift_[1], xshift_[2]);
    for (unsigned ilevel = levelmin_; ilevel <= levelmax_; ++ilevel) {
      LOGUSER("   Level %3d :   offset = (%5d,%5d,%5d)", ilevel, ox_[ilevel], oy_[ilevel], oz_[ilevel]);
      LOGUSER("                   size = (%5d,%5d,%5d)", nx_[ilevel], ny_[ilevel], nz_[ilevel]);
    }
  }
};

typedef GridHierarchy<real_t> grid_hierarchy;
typedef MeshvarBnd<real_t> meshvar_bnd;
typedef Meshvar<real_t> meshvar;

#endif