// AmrCore functions
CactusAmrMesh::CactusAmrMesh() {}
CactusAmrMesh::CactusAmrMesh(const RealBox *rb, int max_level_in,
                             const Vector<int> &n_cell_in, int coord,
                             Vector<IntVect> ref_ratios, const int *is_per)
    : AmrMesh(rb, max_level_in, n_cell_in, coord, ref_ratios, is_per) {}
CactusAmrMesh::CactusAmrMesh(const RealBox &rb, int max_level_in,
                             const Vector<int> &n_cell_in, int coord,
                             Vector<IntVect> const &ref_ratios,
                             Array<int, AMREX_SPACEDIM> const &is_per)
    : AmrMesh(rb, max_level_in, n_cell_in, coord, ref_ratios, is_per) {}
CactusAmrMesh::~CactusAmrMesh() {}

void CactusAmrMesh::ErrorEst(int lev, TagBoxArray &tags, Real time, int ngrow) {
  // // refine everywhere
  // tags.setVal(boxArray(lev), TagBox::SET);
  // refine centre
  const BoxArray &ba = boxArray(lev);
  const Box &bx = ba.minimalBox();
  Box nbx;
  for (int d = 0; d < dim; ++d) {
    int sz = bx.bigEnd(d) - bx.smallEnd(d) + 1;
    nbx.setSmall(d, bx.smallEnd(d) + sz / 4 + 1);
    nbx.setBig(d, bx.bigEnd(d) - sz / 4 - 1);
  }
  cout << "nbx: " << nbx << "\n";
  tags.setVal(intersect(ba, nbx), TagBox::SET);
}
////////////////////////////////////////////////////////////////////////////////

  const RealBox domain({AMREX_D_DECL(xmin, ymin, zmin)},
                       {AMREX_D_DECL(xmax, ymax, zmax)});

  const RealBox domain({xmin, ymin, zmin}, {xmax, ymax, zmax});

  const int maxlevels = 1;

  const int maxnumlevels = 2;

  ghext->amrmesh = make_unique<AmrMesh>(domain, maxlevels, ncells, coord,
                                        reffacts, periodic);

  ghext->amrmesh = make_unique<CactusAmrMesh>(domain, maxnumlevels - 1, ncells,
                                              coord, reffacts, periodic);
#warning "TODO: increase blocking factor"
  // const int blocking_factor = 8;
  const int blocking_factor = 1;
  ghext->amrmesh->SetBlockingFactor(blocking_factor);

  ghext->amrmesh->MakeNewGrids(0.0);

  CCTK_VINFO("Geometry:");
  cout << ghext->amrmesh->Geom(0) << "\n";
  CCTK_VINFO("BoxArray:");
  cout << ghext->amrmesh->boxArray(0) << "\n";
  CCTK_VINFO("DistributionMap:");
  cout << ghext->amrmesh->DistributionMap(0) << "\n";

  for (int level = 0; level < maxnumlevels; ++level) {
    CCTK_VINFO("Geometry level %d:", level);
    cout << ghext->amrmesh->Geom(level) << "\n";
  }

  const int nlevels = 1;
  ghext->leveldata.resize(nlevels);
  for (int level = 0; level < nlevels; ++level) {

  const int numlevels = maxnumlevels;
  ghext->leveldata.resize(numlevels);
  for (int level = 0; level < numlevels; ++level) {

    CCTK_REAL time = 0.0; // dummy time

    // Create grid structure
    if (level == 0) {
      // Create coarse grid
      ghext->amrmesh->MakeNewGrids(time);
    } else {
      // Create refined grid
      int new_finest = -999;
      Vector<BoxArray> new_grids;
      ghext->amrmesh->MakeNewGrids(0, time, new_finest, new_grids);
      assert(new_finest == level);
    }
    // CCTK_VINFO("BoxArray level %d:", level);
    // cout << ghext->amrmesh->boxArray(level) << "\n";
    // CCTK_VINFO("DistributionMap level %d:", level);
    // cout << ghext->amrmesh->DistributionMap(level) << "\n";

}

} // namespace AMReX

class CactusAmrMesh : public AmrMesh {
public:
  CactusAmrMesh();
  CactusAmrMesh(const RealBox *rb, int max_level_in,
                const Vector<int> &n_cell_in, int coord = -1,
                Vector<IntVect> ref_ratios = Vector<IntVect>(),
                const int *is_per = nullptr);
  CactusAmrMesh(const RealBox &rb, int max_level_in,
                const Vector<int> &n_cell_in, int coord,
                Vector<IntVect> const &ref_ratios,
                Array<int, AMREX_SPACEDIM> const &is_per);
  CactusAmrMesh(const AmrMesh &rhs) = delete;
  CactusAmrMesh &operator=(const AmrMesh &rhs) = delete;
  virtual ~CactusAmrMesh();
  // virtual void MakeNewLevelFromScratch(int lev, Real time, const BoxArray
  // &ba,
  //                                      const DistributionMapping &dm)
  //                                      override;
  virtual void ErrorEst(int lev, TagBoxArray &tags, Real time,
                        int ngrow) override;
  // virtual void ManualTagsPlacement(int lev, TagBoxArray &tags,
  //                                  const Vector<IntVect> &bf_lev) override;
  // virtual BoxArray GetAreaNotToTag(int lev) override;
};

  unique_ptr<AmrMesh> amrmesh;

  unique_ptr<CactusAmrMesh> amrmesh;

reduction<T> reduce_array(const T *restrict var, const int *restrict ash,
                          const int *restrict lsh) {

reduction<T> reduce_array(const Geometry &geom, const T *restrict var,
                          const int *restrict ash, const int *restrict lsh) {
  const CCTK_REAL *restrict dx = geom.CellSize();
  CCTK_REAL dV = 1.0;
  for (int d = 0; d < dim; ++d)
    dV *= dx[d];

        red = reduction<T>(red, reduction<T>(var[idx]));

        red = reduction<T>(red, reduction<T>(dV, var[idx]));

      const auto &restrict geom = ghext->amrmesh->Geom(leveldata.level);

        red += reduce_array(var, ash, lsh);

        red += reduce_array(geom, var, ash, lsh);

  T count, maxabs, sumabs, sum2abs;

  T vol, maxabs, sumabs, sum2abs;

  reduction(const T &x);

  reduction(const T &V, const T &x);

  T avg() const { return sum / vol; }
  T sdv() const { return sqrt(fmax(T(0.0), vol * sum2 - sum * sum)); }
  T norm0() const { return vol; }
  T norm1() const { return sumabs / vol; }
  T norm2() const { return sqrt(sum2abs / vol); }
  T norm_inf() const { return maxabs; }

    : min(1.0 / 0.0), max(-1.0 / 0.0), sum(0.0), sum2(0.0), count(0.0),

    : min(1.0 / 0.0), max(-1.0 / 0.0), sum(0.0), sum2(0.0), vol(0.0),

reduction<T>::reduction(const T &x)
    : min(x), max(x), sum(x), sum2(x * x), count(1.0), maxabs(fabs(x)),
      sumabs(fabs(x)), sum2abs(fabs(x) * fabs(x)) {}

reduction<T>::reduction(const T &V, const T &x)
    : min(x), max(x), sum(V * x), sum2((V * x) * (V * x)), vol(V),
      maxabs(fabs(x)), sumabs(fabs(V * x)), sum2abs(fabs(V * x) * fabs(V * x)) {
}

      sum2(x.sum2 + y.sum2), count(x.count + y.count),

      sum2(x.sum2 + y.sum2), vol(x.vol + y.vol),

  int ncells[dim] = {ncells_x, ncells_y, ncells_z};
  CCTK_REAL x0[dim] = {xmin, ymin, zmin};
  CCTK_REAL x1[dim] = {xmax, ymax, zmax};
  CCTK_REAL dx[dim];
  for (int d = 0; d < dim; ++d)
    dx[d] = (x1[d] - x0[d]) / ncells[d];
  CCTK_REAL mindx = 1.0 / 0.0;
  for (int d = 0; d < dim; ++d)
    mindx = fmin(mindx, dx[d]);

  const Geometry &geom = ghext->amrmesh->Geom(0);
  const CCTK_REAL *restrict x0 = geom.ProbLo();
  const CCTK_REAL *restrict dx = geom.CellSize();

    cctkGH->cctk_origin_space[d] = x0[d] + 0.5 * dx[d];

    cctkGH->cctk_origin_space[d] = x0[d] + 0.0 * dx[d];

  CCTK_REAL mindx = 1.0 / 0.0;
  const int numlevels = ghext->amrmesh->finestLevel() + 1;
  for (int level = 0; level <numlevels; ++level) {
    const Geometry &geom = ghext->amrmesh->Geom(level);
    const CCTK_REAL *restrict dx = geom.CellSize();
    for (int d = 0; d < dim; ++d)
      mindx = fmin(mindx, dx[d]);
  }

  // TODO: Get this from mfab
#warning "TODO"
  int ncells[dim] = {ncells_x, ncells_y, ncells_z};

  const Box &domain = ghext->amrmesh->Geom(leveldata.level).Domain();

    cctkGH->cctk_gsh[d] = ncells[d];

    cctkGH->cctk_gsh[d] = domain[orient(d, 1)] - domain[orient(d, 0)] + 1;

    // The refinement factor over the top level (coarsest) grid
#warning "TODO"

  // The refinement factor over the top level (coarsest) grid

    cctkGH->cctk_levfac[d] = 1; // TODO

    cctkGH->cctk_levfac[d] = 1 << leveldata.level;

    // Offset between this level's and the coarsest level's origin
#warning "TODO"
  for (int d = 0; d < dim; ++d)
    cctkGH->cctk_levoff[d] = 0; // TODO
#warning "TODO"
  for (int d = 0; d < dim; ++d)
    cctkGH->cctk_levoffdenom[d] = 0; // TODO

  // Offset between this level's and the coarsest level's origin
  for (int d = 0; d < dim; ++d) {
    cctkGH->cctk_levoff[d] = 1;
    cctkGH->cctk_levoffdenom[d] = 2;
  }

    int ncells[dim];
    for (int d = 0; d < dim; ++d)
      ncells[d] = cctkGH->cctk_gsh[d];

    enter_level_mode(cctkGH, ghext->leveldata.at(0));
    const int *restrict gsh = cctkGH->cctk_gsh;

      x0[d] = cctkGH->cctk_origin_space[d] - 0.5 * dx[d];

      x0[d] = cctkGH->cctk_origin_space[d] - 0.0 * dx[d];

      x1[d] = x0[d] + ncells[d] * dx[d];

      x1[d] = x0[d] + gsh[d] * dx[d];

    CCTK_VINFO("  gsh=[%d,%d,%d]", ncells[0], ncells[1], ncells[2]);

    CCTK_VINFO("  gsh=[%d,%d,%d]", gsh[0], gsh[1], gsh[2]);

    leave_level_mode(cctkGH, ghext->leveldata.at(0));

  const CCTK_REAL dt = cctk_delta_time;
  const CCTK_REAL *restrict const x0 = cctk_origin_space;
  const CCTK_REAL *restrict const dx = cctk_delta_space;

  const CCTK_REAL dt = CCTK_DELTA_TIME;
  CCTK_REAL x0[dim], dx[dim];
  for (int d = 0; d < dim; ++d) {
    x0[d] = CCTK_ORIGIN_SPACE(d);
    dx[d] = CCTK_DELTA_SPACE(d);
  }

  const CCTK_REAL *restrict const dx = cctk_delta_space;
  const CCTK_REAL dt = cctk_delta_time;

  const CCTK_REAL dt = CCTK_DELTA_TIME;
  CCTK_REAL dx[dim];
  for (int d = 0; d < dim; ++d) {
    dx[d] = CCTK_DELTA_SPACE(d);
  }

  const CCTK_REAL *restrict const x0 = cctk_origin_space;
  const CCTK_REAL *restrict const dx = cctk_delta_space;

  CCTK_REAL x0[dim], dx[dim];
  for (int d = 0; d < dim; ++d) {
    x0[d] = CCTK_ORIGIN_SPACE(d);
    dx[d] = CCTK_DELTA_SPACE(d);
  }

AMReX::ncells_x = 128
AMReX::ncells_y = 128
AMReX::ncells_z = 128