struct Adapt_Field {
  scalar f;       // The field with leaf-filled values
  double max;     // refinement criterion
  double min;     // coarseing criterion 
  int maxlevel;   // maximum level of refinement
  int minlevel;   // minimum level of refinement (default 1)
  scalar * list;  // list of fields to update (default all)
};

astats adapt_field (struct Adapt_Field p);

extern vector u;

struct Adapt_Flow {
  double ue;
  int maxlevel;
  double p;
  double cfac;
};

attribute {
  void (* interpolant) (Point, scalar);
}

void wavelet_interpolant (scalar s, scalar w)
{
  restriction ({s});
  for (int l = grid->maxdepth - 1; l >= 0; l--) {
    foreach_coarse_level (l) {
      foreach_child()
        w[] = s[];
      s.interpolant (point, s);
      foreach_child() {
        double sp = s[];
        s[] = w[];
        /* difference between fine value and its prolongation */
        w[] -= sp;
      }
    }
    boundary_level ({w}, l + 1);
  }
  /* root cell */
  foreach_level(0) 
    w[] = s[];
  boundary_level ({w}, 0);
}

astats adapt_flow (struct Adapt_Flow pa) {
  double ue = pa.ue;
  int maxlevel = pa.maxlevel;
  double p = pa.p;
  double uec = ue/1.5;
  if (pa.cfac)
    uec = ue/pa.cfac; 
  scalar w[];
  vector wv[];
  foreach_dimension() 
    wavelet_interpolant (u.x, wv.x);
    
  for (int l = 1; l < grid->maxdepth; l++)
    foreach_coarse_level (l) { 
      double maxw = 0.;
      foreach_child() {
	double a = 0;
	foreach_dimension()
	  a += fabs(wv.x[]);
	maxw = max(maxw, a);
      }
      w[] = maxw;
    }
#if _MPI
  w.restriction = no_restriction;
  w.prolongation = no_restriction;
  boundary ({w});
#endif
  for (int l = grid->maxdepth; l > 1; l--)  
    foreach_level(l)  {// 3x3(x3) Gaussian kernel
#if (dimension == 2)
      w[] = pow(Delta, p)*(1.  * coarse(w,0,0,0)   +
			   .5  *(coarse(w,1,0,0)   + coarse(w,0,1,0)   +
				 coarse(w,-1,0,0)  + coarse(w,0,-1,0)) +
			   0.25*(coarse(w,1,1,0)   + coarse(w,1,-1,0)  +
				 coarse(w,-1,-1,0) + coarse(w,-1,1,0)))/4.;
#else //(dimension == 3) 
      w[] = pow(Delta, p)*(1.   * coarse(w,0,0,0)   +
			   .5   *(coarse(w,1,0,0)   + coarse(w,-1,0,0)  +
				  coarse(w,0,1,0)   + coarse(w,0,-1,0)  +
				  coarse(w,0,0,1)   + coarse(w,0,0,-1)) +
			   0.25 *(coarse(w,1,1,0)   + coarse(w,1,-1,0)  +
				  coarse(w,-1,-1,0) + coarse(w,-1,1,0)  +
				  coarse(w,0,1,1)   + coarse(w,0,-1,1)  +
				  coarse(w,0,-1,-1) + coarse(w,0,1,-1)  +
				  coarse(w,1,0,1)   + coarse(w,1,0,-1)  +
				  coarse(w,-1,0,-1) + coarse(w,-1,0,1)) +
			   0.125*(coarse(w,1,1,1)   + coarse(w,-1,1,1)  +
				  coarse(w,1,-1,1)  + coarse(w,1,1,-1)  +
				  coarse(w,-1,-1,1) + coarse(w,-1,1,-1) +
				  coarse(w,1,-1,-1) + coarse(w,-1,-1,-1)))/8.;
#endif
    }
  return adapt_field (w, ue, uec, maxlevel);
}



struct Adapt_List {
  scalar * list;
  double * ues;
  int maxlevel;
  double p;
  double fac;
};

astats adapt_list (struct Adapt_List pa) {
  scalar * list = pa.list;
  double * ues = pa.ues;
  int maxlevel = pa.maxlevel;
  double p = pa.p;
  if (!pa.fac)
    pa.fac = 1.5;
  scalar w[], * wl = list_clone (list);
 scalar s, ws;
  for (s,ws in list,wl)
    wavelet (s, ws);
  for (int l = 1; l < depth() - 1; l++)
    foreach_coarse_level (l) { 
      double maxw = 0.;
      foreach_child() {
	double a = 0;
	scalar s, ws;
	NOT_UNUSED(s);
	int is = 0;
	for (s, ws in list, wl) {
	  a += fabs(ws[])/ues[is]; //scale w with ue
	  is++;
	}
	maxw = max(maxw, a);
      }
      w[] = maxw;
    }
#if _MPI
  w.restriction = no_restriction;
  w.prolongation = no_restriction;
  boundary ({w});
#endif
  for (int l = depth(); l > 1; l--)  
    foreach_level(l)  {// 3x3(x3) Gaussian kernel
#if (dimension == 2)
      w[] = pow(Delta, p)*(1.  * coarse(w,0,0,0)   +
			   .5  *(coarse(w,1,0,0)   + coarse(w,0,1,0)   +
				 coarse(w,-1,0,0)  + coarse(w,0,-1,0)) +
			   0.25*(coarse(w,1,1,0)   + coarse(w,1,-1,0)  +
				 coarse(w,-1,-1,0) + coarse(w,-1,1,0)))/4.;
#else //(dimension == 3) 
      w[] = pow(Delta, p)*(1.   * coarse(w,0,0,0)   +
			   .5   *(coarse(w,1,0,0)   + coarse(w,-1,0,0)  +
				  coarse(w,0,1,0)   + coarse(w,0,-1,0)  +
				  coarse(w,0,0,1)   + coarse(w,0,0,-1)) +
			   0.25 *(coarse(w,1,1,0)   + coarse(w,1,-1,0)  +
				  coarse(w,-1,-1,0) + coarse(w,-1,1,0)  +
				  coarse(w,0,1,1)   + coarse(w,0,-1,1)  +
				  coarse(w,0,-1,-1) + coarse(w,0,1,-1)  +
				  coarse(w,1,0,1)   + coarse(w,1,0,-1)  +
				  coarse(w,-1,0,-1) + coarse(w,-1,0,1)) +
			   0.125*(coarse(w,1,1,1)   + coarse(w,-1,1,1)  +
				  coarse(w,1,-1,1)  + coarse(w,1,1,-1)  +
				  coarse(w,-1,-1,1) + coarse(w,-1,1,-1) +
				  coarse(w,1,-1,-1) + coarse(w,-1,-1,-1)))/8.;
#endif
    }
  return adapt_field (w, 1., 1./pa.fac, maxlevel); //w is already scaled
}



trace
astats adapt_field (struct Adapt_Field p) {
  scalar f = p.f;
  if (p.list == NULL)
    p.list = all;
  scalar * listc = NULL;
  for (scalar s in p.list)
    if (!is_constant(s) && s.restriction != no_restriction)
      listc = list_add (listc, s);
  astats st = {0, 0};
  // refinement
  if (p.minlevel < 1)
    p.minlevel = 1;
  tree->refined.n = 0;
  static const int refined = 1 << user, too_fine = 1 << (user + 1);
  foreach_cell() {
    if (is_active(cell)) {
      static const int too_coarse = 1 << (user + 2);
      if (is_leaf (cell)) {
  	if (cell.flags & too_coarse) {
  	  cell.flags &= ~too_coarse;
  	  refine_cell (point, listc, refined, &tree->refined);
  	  st.nf++;
  	}
  	continue;
      }
      else { // !is_leaf (cell)
  	if (cell.flags & refined) {
  	  // cell has already been refined, skip its children
  	  cell.flags &= ~too_coarse;
  	  continue;
  	}
  	// check whether the cell or any of its children is local
  	bool local = is_local(cell);
  	if (!local)
  	  foreach_child()
  	    if (is_local(cell)) {
  	      local = true;
	      break;
	    }
  	if (local) {
  	  static const int just_fine = 1 << (user + 3);
	  double max = p.max, min = p.min;
	  foreach_child() {
	    double e = fabs(f[]);
	    if (e > max && level < p.maxlevel) {
	      cell.flags &= ~too_fine;
	      cell.flags |= too_coarse;
	    }
	    else if ((e <= min || level > p.maxlevel) &&
		     !(cell.flags & (too_coarse|just_fine))) {
	      if (level >= p.minlevel)
		cell.flags |= too_fine;
	    }
	    else if (!(cell.flags & too_coarse)) {
	      cell.flags &= ~too_fine;
	      cell.flags |= just_fine;
	    }
	  }
	  foreach_child() {
  	    cell.flags &= ~just_fine;
  	    if (!is_leaf(cell)) {
  	      cell.flags &= ~too_coarse;
  	      if (level >= p.maxlevel)
  		cell.flags |= too_fine;
  	    }
  	    else if (!is_active(cell))
  	      cell.flags &= ~too_coarse;
  	  }
  	}
      }
    }
    else // inactive cell
      continue;
  }
  mpi_boundary_refine (listc);
  
  // coarsening
  // the loop below is only necessary to ensure symmetry of 2:1 constraint
  for (int l = depth(); l >= 0; l--) {
    foreach_cell()
      if (!is_boundary(cell)) {
  	if (level == l) {
  	  if (!is_leaf(cell)) {
  	    if (cell.flags & refined)
  	      // cell was refined previously, unset the flag
  	      cell.flags &= ~(refined|too_fine);
  	    else if (cell.flags & too_fine) {
  	      if (is_local(cell) && coarsen_cell (point, listc))
  		st.nc++;
  	      cell.flags &= ~too_fine; // do not coarsen parent
  	    }
  	  }
  	  if (cell.flags & too_fine)
  	    cell.flags &= ~too_fine;
  	  else if (level > 0 && (aparent(0).flags & too_fine))
  	    aparent(0).flags &= ~too_fine;
  	  continue;
  	}
  	else if (is_leaf(cell))
  	  continue;
      }
    mpi_boundary_coarsen (l, too_fine);
    }
  free (listc);
  mpi_all_reduce (st.nf, MPI_INT, MPI_SUM);
  mpi_all_reduce (st.nc, MPI_INT, MPI_SUM);
  if (st.nc || st.nf)
    mpi_boundary_update (p.list);
  return st;
}