#include "grid/quadtree.h"
#include "navier-stokes/centered.h"

int main(){
  int jmax=9;
  int n=4;
  int nm=((n-1)*pow(4,jmax-1))+1;
  double xn[nm],yn[nm];
  double xk[nm];
  double yk[nm];
  xk[0]=0;
  xk[1]=1;
  xk[2]=0.5;
  xk[3]=0;
  yk[0]=0;
  yk[1]=0;
  yk[2]=-sqrt(3.)/2;
  yk[3]=0;
  fprintf(ferr,"\nyk[n-1]=%g\n\n",yk[n-1]);
  for (int j=1;j<jmax;j++){
    for (int k=0;k<(n-1);k++){
      xn[(k*4)]=xk[(k)];
      xn[(k*4)+1]=(2*xk[k]+xk[k+1])/3;
      xn[(k*4)+2]=((xk[k]+xk[k+1])/2)+((yk[k]-yk[k+1])*sqrt(3.)/6);      
      xn[(k*4)+3]=(xk[k]+2*xk[k+1])/3;
      yn[(k*4)]=yk[k];
      yn[(k*4)+1]=(2*yk[k]+yk[k+1])/3;
      yn[(k*4)+2]=((yk[k]+yk[k+1])/2)-((xk[k]-xk[k+1])*sqrt(3.)/6);
      yn[(k*4)+3]=(yk[k]+2*yk[k+1])/3;
    }
    yn[4*(n-1)]=yk[n-1];
    xn[4*(n-1)]=xk[n-1];
    FILE * fpkoch = fopen("koch.dat","w");
    for (int l=0;l<(4*(n-1)+1);l++){
      yk[l]=yn[l];
      xk[l]=xn[l];
    }
    if (j==jmax-1){
      for (int l=0;l<(4*(n-1)+1);l++){
	fprintf(fpkoch,"%g %g\n",xk[l],yk[l]);
      }
    }
    n=(n-1)*4+1;
  }
  
  
  FILE * fpit = fopen("cells.dat","w");
  X0=-0.5;
  Y0=-1.2;
  L0=2.0;
  init_grid(32);
  
  scalar c[];
  for (int m=0;m<7;m++){
    foreach()
      c[]=0;
    for(int i=0;i<n;i+=64/pow(2,m)){
      Point point = locate(xn[i],yn[i]);
      c[]+=1.0;
    }
    int gcc=0;
    int gcn=0;
    foreach(reduction(+:gcc) reduction(+:gcn)){
      gcn++;
      if(c[]>0.5)
	gcc++;
    }

    refine((c[]>0.05)&&(level<7+m));
      
    fprintf(fpit,"%d %d %d\n",m,gcc,gcn);
    fflush(fpit);
    int i=m;
    scalar lev[];
    foreach(){
      lev[]=level;
    }
    static FILE * fp1 = popen ("ppm2gif --delay 200 > f.gif", "w");
    output_ppm (lev, fp1,512,min=4,max=15);
  }
}