// ............
verbosity=0;
border a0(t=2*pi,0) {x=0.5*cos(t); y=0.5*sin(t); label=1;}
border a1(t=20,-10) {x=t; y=10; label=2;}
border a2(t=10,-10) {x=-10; y=t; label=2;}
border a3(t=-10,20) {x=t; y=-10; label=2;}
border a4(t=-10,10) {x=20; y=t; label=3;}

// Ndelta controls the wall resolution
//int Ndelta=32;
int Ndelta=64;
//int Ndelta=16;
mesh Th= buildmesh(a0(Ndelta)+a1(30)+a2(20)+a3(30)+a4(20));

plot(Th,wait=0,ps="mesh.eps");

load "Element_P3";

fespace Vh(Th,P2);
fespace Vh1(Th,P1);

//fespace Vh(Th,P3);
//fespace Vh1(Th,P2);

real uinf=1.0; //free stream velocity
real nu = 0.005; // kinematic viscosity
real rho=1.0; //density
real dt = 0.025; //time step
int Ndt=4000;

real area=1.0;
real fDrag,fLift,CD,CL;
real volumeDomain,error,time;

volumeDomain=int2d(Th)(1.0);

{
  ofstream myfile("forces.dat");
}

Vh w, u = 0, v =0, uold, vold, vorticity;
Vh1 w1, p = 0, q=0;

time=0.0;

for(int n=0;n<Ndt;n++)
 {	
  uold = u;
  vold = v;
  
  solve pb4u(u,w,init=n,solver=GMRES,eps=1.e-6)
        =int2d(Th)(u*w/dt +nu*(dx(u)*dx(w)+dy(u)*dy(w)))
        -int2d(Th)((convect([uold,vold],-dt,uold)/dt-(1.0/rho)*dx(p))*w)
        + on(1,u = 0) + on(2,u = uinf);

  solve pb4v(v,w,init=n,solver=GMRES,eps=1.e-6)
        = int2d(Th)(v*w/dt +nu*(dx(v)*dx(w)+dy(v)*dy(w)))
        -int2d(Th)((convect([uold,vold],-dt,vold)/dt-(1.0/rho)*dy(p))*w)
        + on(1,v = 0) + on(2,v = 0);

  solve pb4p(q,w1,solver=GMRES,init=n,eps=1.e-6) 
        = int2d(Th)(dx(q)*dx(w1)+dy(q)*dy(w1))
        + int2d(Th)((dx(u)+ dy(v))*w1*rho/dt)+ on(3,q=0);

  p = p + q;
  u = u - dt*(1.0/rho)*dx(q);
  v = v - dt*(1.0/rho)*dy(q);

  vorticity = dx(v)-dy(u);

// calculate forces
  fDrag=int1d(Th,1)(-2.0*rho*nu*dx(u)*N.x-rho*nu*(dy(u)+dx(v))*N.y+p*N.x);
  fLift=int1d(Th,1)(-rho*nu*(dy(u)+dx(v))*N.x-2.0*rho*nu*dy(v)*N.y+p*N.y);
  CD=2.0*fDrag/(rho*uinf*uinf*area);
  CL=2.0*fLift/(rho*uinf*uinf*area);
  time=time+dt;

  {
    ofstream myfile("forces.dat", append );
    myfile << time << " " << CD << " " << CL << endl;
  }


  plot(p,wait=0,ps="pressure.eps",nbiso=20, fill=1);

  error=int2d(Th)(dx(u)+dy(v))*dt/volumeDomain;

  cout << " time step " << n  << " mass conservation error " << error << endl; 
  cout << " CD " << CD  << " CL " << CL << endl;
  cout << endl;
 }

//export results in vtk format for paraview
load"iovtk";
savevtk("results.vtk", Th, [u,v,0], p, vorticity, dataname=" velocity pressure vorticity");