// remark: the sign of p is correct 
real s0=clock();
real pi=4*atan(1);
border a(t=0,8){x=t;y=0;label=1;};
border b(t=0,2){x=8;y=t;label=4;};
border c(t=0,8){x=8-t;y=2;label=1;};
border d(t=0,2){x=0;y=2-t;label=2;};

real inSquare=0.15;
real inSquareX=3;
real inSquareY=1;

border e(t=0,2){x=inSquareX+inSquare*cos(pi*t);y=inSquareY+inSquare*sin(pi*t);label=3;};

mesh Th = buildmesh ( a(20)+ b(10) + c(20) +d(10) + e(-10));
mesh[int] Thg(50);

fespace Vh2(Th,P2);
fespace Vh(Th,P1);
Vh2 u2,v2,up1,up2,up3,up4,up5;
Vh2 u1,v1; 
Vh2 w1,w2,wp1,wp2,t1,t2;
//Vh  u1x=0,u1y,u2x,u2y, vv;

real reylnods=200;
real Cs=0.1;
real delta=0.0;

//cout << " Enter the reynolds number :"; cin >> reylnods;
assert(reylnods>1 && reylnods < 100000); 
up1=0;
up2=0; 
func g=(y)*(2-y)*3; 
Vh p=0,q,pw,qw,uno;
real alpha=0;
real  nu=1;
int i=0,iter=0;
real dt=0;

problem Euler ([u1,u2,p],[v1,v2,q],solver=Crout,init=0) =
    int2d(Th)(
             alpha*( u1*v1 + u2*v2) 
            + nu * ( dx(u1)*dx(v1) + dy(u1)*dy(v1)
            +        dx(u2)*dx(v2) + dy(u2)*dy(v2) )
            + p*q*(0.000001) 
            - p*dx(v1) - p*dy(v2)
            - dx(u1)*q - dy(u2)*q
           )
  + int2d(Th) ( -alpha*convect([up1,up2],-dt,up1)*v1 -alpha*convect([up1,up2],-dt,up2)*v2 )
  + on(2,u2=0,u1=g) + on(1,u1=0,u2=0)+on(3,u2=0,u1=0);

problem NS ([u1,u2,p],[v1,v2,q],solver=Crout,init=0) =
    int2d(Th)(
             alpha*( u1*v1 + u2*v2) 
            + nu * ( dx(u1)*dx(v1) + dy(u1)*dy(v1)
            +        dx(u2)*dx(v2) + dy(u2)*dy(v2) )
	    + Cs*Cs*delta*delta*(abs(dx(up1))+abs(dy(up2))+0.5*abs(dx(up2)+dy(up1)))*( dx(u1)*dx(v1) + 0.5*(dx(u2)+dy(u1))*(dx(v2)+dy(v1))
            +        dy(u2)*dy(v2) )
            + p*q*(0.000001) 
            - p*dx(v1) - p*dy(v2)
            - dx(u1)*q - dy(u2)*q
           )
  + int2d(Th) ( -alpha*convect([up1,up2],-dt,up1)*v1 -alpha*convect([up1,up2],-dt,up2)*v2 )
  + on(2,u2=0,u1=g) + on(1,u1=0,u2=0)+on(3,u2=0,u1=0);

problem Sensitivity ([w1,w2,pw],[t1,t2,qw],solver=LU,init=0) =
    int2d(Th)(
             alpha*( w1*t1 + w2*t2) 
            + nu * ( dx(w1)*dx(t1) + dy(w1)*dy(t1)
            +        dx(w2)*dx(t2) + dy(w2)*dy(t2) )
	    + Cs*Cs*delta*delta*(abs(dx(u1))+abs(dy(u2))+0.5*abs(dx(u2)+dy(u1)))*( dx(w1)*dx(t1) + 0.5*(dx(w2)+dy(w1))*(dx(t2)+dy(t1))
            +        dy(w2)*dy(t2) )
	    + Cs*Cs*delta*delta*(dx(w1)*dx(u1)/(abs(dx(u1))+dt*dt)+dy(w2)*dy(u2)/(abs(dy(u2))+dt*dt)+
		0.5*(dx(w2)+dy(w1))*(dx(u2)+dy(u1))/(abs(dx(u2)+dy(u1))+dt*dt))*
                (dx(u1)*dx(t1) + 0.5*(dx(u2)+dy(u1))*(dx(t2)+dy(t1))+dy(u2)*dy(t2) )
            + pw*qw*(0.000001) 
            - pw*dx(t1) - pw*dy(t2)
            - dx(w1)*qw- dy(w2)*qw
           )
  - int2d(Th) (alpha*(wp1*t1+wp2*t2))
  + int2d(Th) ( t1*(u1*dx(w1)+u2*dy(w2))+t2*(u1*dx(w2)+u2*dy(w2))) 
  + int2d(Th) ( t1*(w1*dx(u1)+w2*dy(u2))+t2*(w1*dx(u2)+w2*dy(u2))) 
  + int2d(Th) (Cs*Cs*(abs(dx(u1))+abs(dy(u2))+0.5*abs(dx(u2)+dy(u1)))*( dx(u1)*dx(t1) + 0.5*(dx(u2)+dy(u1))*(dx(t2)+dy(t1))
            +        dy(u2)*dy(t2) ))
  + on(2,w2=0,w1=0) + on(1,w1=0,w2=0)+on(3,w2=0,w1=0);

//plot(coef=0.2,cmm=" [u1,u2] et p  ",p,[u1,u2],ps="StokesP2P1.eps",value=1,wait=1);
{
  real[int] xx(21),yy(21),pp(21);
  for (int i=0;i<21;i++)
   {
     yy[i]=i/20.;
     xx[i]=u1(0.5,i/20.);
     pp[i]=p(i/20.,0.999);
    }
      cout << " " << yy << endl;
     plot([xx,yy],wait=1,cmm="u1 x=0.5 cup");
     plot([yy,pp],wait=1,cmm="pressure y=0.999 cup");
}

Euler;
up1=u1;	
up2=u2;

dt = 0.015;
alpha=1/dt;
nu=1./reylnods;   

// int nbiter = 5;
int nbiter=2;
// int nbiter=1;
// int nrpasses = 1;
int nrpasses = 5;
real coefdt = 0.25^(1./nbiter);
real coefcut = 0.25^(1./nbiter);
real cut=0.01;
real tol=0.07,coeftol = 0.25^(1./nbiter);
alpha=1/dt;
nu=1./reylnods;   
int nrdeltas=10;
real deltaDelta=0.05;

Thg[0]=Th;
int nrMeshes=1;


int stepCounter=0;
for (iter=1;iter<=nrpasses*nbiter;iter++)
{
  cout << " dt = " << dt << " ------------------------ " << endl;
  Th=adaptmesh(Th,[dx(u1),dy(u1),dx(u1),dy(u2)],
             abserror=0,cutoff=cut, err=tol, inquire=0,ratio=1.5,hmin=1./1000);
//  Thg[nrMeshes]=Th;
// old version, in which I was using adaptation
  nrMeshes=nrMeshes+1; 	
  alpha=1/dt;
for (i=0;i<=50;i++)
// for (i=0;i<=20;i++)
// for (i=0;i<=5;i++)
// for (i=0;i<=2;i++)
   {
  cout << " step is" << stepCounter << endl;
  stepCounter++;
     up1=u1;
     up2=u2;
     NS;
     up3=sqrt(u1*u1+u2*u2);
     up4=sqrt(abs(dx(u2)-dy(u1)));	
//     if ( !(i % 5 )) 
//	{ 
//     	plot(coef=1, bb=[[2.5,0.5],[3.5,1.5]], cmm="Norm of Velocity ",up3,ps="plotNS_"+(iter+10)+"_"+i+".eps",value=1);
//     	plot(coef=0.2,cmm="Square root of the Vorticity ",up4,ps="plotNSVor_"+(iter+10)+"_"+i+".eps",value=1);  
//     	plot(coef=1, bb=[[2.5,0.5],[3.5,1.5]], cmm="Dif ",up5,ps="plotNSDif_"+(iter+10)+"_"+i+".eps",value=1);  
//	}
     cout << "CPU " << clock()-s0 << "s " << endl;     
   } 
 // tol=tol/2; 	
// tol=tol/1.02;
  if (iter>= nrpasses*nbiter) break;
//  dt = dt*coefdt;
//  tol = tol *coeftol;
//  cut = cut *coefcut;
}

dt=dt/100;
alpha=1/dt;

// save the mesh 
// savemesh(Th,"toto.msh");


// start the data loop 
// define the 

real[int] loopDelta(4);
real[int] scalingP(100);
real areaDomain;
real[int] drag(100);
loopDelta[0]=0.0;
loopDelta[1]=sqrt(0.04);
loopDelta[2]=sqrt(0.041);
loopDelta[3]=sqrt(0.042);
// the drag in the x direction (since due to symmetry, it will be 0 in the other direction)
real[int] err(5);

// compute the pressure and stresses at the boundary of the body 
	real[int] sa11(1001),sa22(1001),sa12(1001),pa(1001);
	real[int] sb11(1001),sb22(1001),sb12(1001),pb(1001);
	real[int] sc11(1001),sc22(1001),sc12(1001),pc(1001);
        real[int] indices(1001);

// need to do this because freefem does not support a matrix mode
uno=1;
areaDomain=int2d(Th)(uno);
int nrDeltaSteps=1;
assert(nrDeltaSteps<=49);

for(int iDelta=0; iDelta<=nrDeltaSteps; iDelta++){
	drag[iDelta]=0;
//	delta=loopDelta[iDelta];
	delta=sqrt(0.04*iDelta/nrDeltaSteps);
	NS;
	scalingP[iDelta]=int2d(Th)(p);
        cout << "scaling is " << scalingP[iDelta] << "\n" ;
	for (int jj=0;jj<1001;jj++)
   	{
        indices[jj]=1/(1001.00)*jj;
	sa11[jj]=dx(u1)(3+cos(jj*2*pi/1001)*0.151,1+0.151*sin(jj*2*pi/1001));
	sa22[jj]=dy(u2)(3+cos(jj*2*pi/1001)*0.151,1+0.151*sin(jj*2*pi/1001));
	sa12[jj]=dx(u2)(3+cos(jj*2*pi/1001)*0.151,1+0.151*sin(jj*2*pi/1001));
	sa12[jj]=0.5*(sa12[jj]+dy(u1)(3+cos(jj*2*pi/1001)*0.151,1+0.151*sin(jj*2*pi/1001)));
	pa[jj]=p(3+cos(jj*2*pi/1001)*0.151,1+0.151*sin(jj*2*pi/1001))-scalingP[iDelta]/areaDomain;
	drag[iDelta]=drag[iDelta]+cos(jj*2*pi/1001)*(pa[jj]-2*nu*sa11[jj])-2*nu*sin(jj*2*pi/1001)*sa12[jj];
	}
	drag[iDelta]=drag[iDelta]*2*pi*0.151/1001;
}
Sensitivity;
scalingP[4]=int2d(Th)(pw);
real sdrag=0;
	for (int jj=0;jj<1001;jj++)
   	{
        indices[jj]=1/(1001.00)*jj;
	sa11[jj]=dx(w1)(3+cos(jj*2*pi/1001)*0.151,1+0.151*sin(jj*2*pi/1001));
	sa22[jj]=dy(w2)(3+cos(jj*2*pi/1001)*0.151,1+0.151*sin(jj*2*pi/1001));
	sa12[jj]=dx(w2)(3+cos(jj*2*pi/1001)*0.151,1+0.151*sin(jj*2*pi/1001));
	sa12[jj]=0.5*(sa12[jj]+dy(w1)(3+cos(jj*2*pi/1001)*0.151,1+0.151*sin(jj*2*pi/1001)));
	pa[jj]=pw(3+cos(jj*2*pi/1001)*0.151,1+0.151*sin(jj*2*pi/1001))-scalingP[4]/areaDomain;
	sdrag=sdrag+cos(jj*2*pi/1001)*(pa[jj]-2*nu*sa11[jj])-2*nu*sin(jj*2*pi/1001)*sa12[jj];
	}
	sdrag=sdrag*2*pi*0.151/1001;

ofstream f("dragMihai5.txt");
for (int iout=0; iout<nrDeltaSteps; iout++){
	cout << iout << " " << drag[iout] << endl;
	f << iout << " " << drag[iout] << endl;
}

cout << " " << sdrag << endl;
f << " " << sdrag << endl;

real[int] meshDepSensDrag(50);
real[int] meshDepVelDrag(50);
// open file for delta dependent output
ofstream gDelta("sensDeltaMihai5.txt");
// new stuff, I am just using split to generate the meshes
nrMeshes=4;
assert(nrMeshes<50);
/* for(int iMeshes=0; iMeshes<(nrMeshes-1); iMeshes++){
Thg[iMeshes+1]=trunc(Thg[iMeshes],1,split=2);
string FileName="trunc"+iMeshes+".ps";
plot(Thg[iMeshes],ps=FileName);
}*/
Th=Thg[0];
for (int iMeshes=0; iMeshes<nrMeshes; iMeshes++){	
  ofstream velValue("velValue"+iMeshes+".txt");	
  NS;
  scalingP[0]=int2d(Th)(p);
  // initialize drag variable
  meshDepVelDrag[iMeshes]=0;
  for (int jj=0;jj<1001;jj++)
    {
      sa11[jj]=dx(u1)(3+cos(jj*2*pi/1001)*0.151,1+0.151*sin(jj*2*pi/1001));
      sa22[jj]=dy(u2)(3+cos(jj*2*pi/1001)*0.151,1+0.151*sin(jj*2*pi/1001));
      sa12[jj]=dx(u2)(3+cos(jj*2*pi/1001)*0.151,1+0.151*sin(jj*2*pi/1001));
      sa12[jj]=0.5*(sa12[jj]+dy(u1)(3+cos(jj*2*pi/1001)*0.151,1+0.151*sin(jj*2*pi/1001)));
      pa[jj]=p(3+cos(jj*2*pi/1001)*0.151,1+0.151*sin(jj*2*pi/1001))-scalingP[0]/areaDomain;
      velValue << jj << " " << sa11[jj] << " " << sa22[jj] << " " << sa12[jj] << " " << pa[jj] << endl;
      meshDepVelDrag[iMeshes]=meshDepVelDrag[iMeshes]+cos(jj*2*pi/1001)*(pa[jj]-2*nu*sa11[jj])-2*nu*sin(jj*2*pi/1001)*sa12[jj];
    }	
  meshDepVelDrag[iMeshes]=meshDepVelDrag[iMeshes]*2*pi*0.151/1001;  
  // now compute sensitivity information for the drag  
  Sensitivity;
  scalingP[0]=int2d(Th)(pw);
  meshDepSensDrag[iMeshes]=0;
  ofstream sensValue("sensValue"+iMeshes+".txt");	
  for (int jj=0;jj<1001;jj++)
    {
      sa11[jj]=dx(w1)(3+cos(jj*2*pi/1001)*0.151,1+0.151*sin(jj*2*pi/1001));
      sa22[jj]=dy(w2)(3+cos(jj*2*pi/1001)*0.151,1+0.151*sin(jj*2*pi/1001));
      sa12[jj]=dx(w2)(3+cos(jj*2*pi/1001)*0.151,1+0.151*sin(jj*2*pi/1001));
      sa12[jj]=0.5*(sa12[jj]+dy(w1)(3+cos(jj*2*pi/1001)*0.151,1+0.151*sin(jj*2*pi/1001)));
      pa[jj]=pw(3+cos(jj*2*pi/1001)*0.151,1+0.151*sin(jj*2*pi/1001))-scalingP[0]/areaDomain;
      sensValue << jj << " " << sa11[jj] << " " << sa22[jj] << " " << sa12[jj] << " " << pa[jj] << endl;
      meshDepSensDrag[iMeshes]=meshDepSensDrag[iMeshes]+cos(jj*2*pi/1001)*(pa[jj]-2*nu*sa11[jj])-2*nu*sin(jj*2*pi/1001)*sa12[jj];
    }
  meshDepSensDrag[iMeshes]=meshDepSensDrag[iMeshes]*2*pi*0.151/1001;
  gDelta << iMeshes << " " << meshDepVelDrag[iMeshes] << " " << meshDepSensDrag[iMeshes] << endl;
  //split mesh
  Th=trunc(Th,1,split=2);
}