#include "grid/octree.h"
#include "utils.h"
#include "additional_attributes.h"

#define func (1.+x+y+z+x*y+x*z+y*z+x*y*z)

scalar c[],err[];
c[top]   = dirichlet(func); // On coarse grid we limit the effect of the boundaries 
c[bottom]= dirichlet(func);
c[left]  = dirichlet(func);
c[right] = dirichlet(func);
c[front] = dirichlet(func);
c[back]  = dirichlet(func);
FILE * fp1;
int cells = 0;
stats s;
int maxlevel = 8;
int main(){
  fp1=fopen("data.dat","w");
  L0=10.;
  X0=Y0=-L0/2;
  init_grid(1<<maxlevel);
  foreach()
    c[]=func;

  for (int n=1;n<5;n++){
    cells=0;
    foreach(reduction(+:cells))
      cells++;
    fprintf(ferr,"%d\t %d\n",n,cells);
    fprintf(fp1,"%d\t%g\t",n,L0/(pow(2.,(double)(maxlevel-n+1.))));
    for (int j=1;j<6;j++){
      if (j==1) c.prolongation=refine_injection;
      if (j==2) c.prolongation=refine_bilinear;
      if (j==3) c.prolongation=refine_linear;
      if (j==4) c.prolongation=refine_quadratic; //Third order?
      if (j==5) c.prolongation=refine_3x3x3_linear;
      wavelet(c,err);
      foreach(){// We are not interested in cells at the boundary. 
	if (abs(x)<3. && abs(y)<3. &&abs(z)<3.) 
	  err[]=fabs(err[]);
	else
	  err[]=0.;
      }
      s=statsf(err);
      fprintf(fp1,"%g\t",s.sum+0.00000001);
    }
    fprintf(fp1,"\n");
    unrefine(level>=(maxlevel-n)); 
  }
}