// -*- c++ -*-


function uin = -4;
double itsteps = 100;

vector n2 = (60,15,30);
vector n1 = 2*n2;

vertex A = (0,0,0);
vertex B = (4,1,2);
mesh M = structured(n1,A,B);
mesh M2= structured(n2,A,B);
scene S = pov("final.pov");
domain Omega = domain(S,outside(<1,0,0>));

mesh obstacle = surface(Omega,M);

save (vtk,"o",obstacle);

double eps = 1E-3;
function P = 1/eps*one(<1,0,0>)+1;

double area = int[M](1);

function mu = 0.1;
femfunction u0(M) = 0;
femfunction v0(M) = 0;
femfunction w0(M) = 0;
femfunction p0(M2) = 0;

function ki = one(<1,0,0>); 

double dt = 0.01;
double pp;
double qq;

double i = 1;

function u = 0;
function v = 0;
function w = 0;

function p1 = 0;

do{
  function deltaP = 2*p0-p1;

   cout << "========== Navier-Stokes step " << i << "\n";

   // Reinitialise pour le gradient conjugue
   u = 0;
   v = 0;
   w = 0;

   solve(u) in M
     krylov(type=cg,precond=diagonal),
     memory(matrix=none)
   {
     test(U)
	int(P/dt*u*U)

       +int(mu*grad(u)*grad(U))
       =
        int(1/dt*convect([u0,v0,w0],dt,u0)*U)

       -int(dx(deltaP)*U)
      ;

      u = uin on 0;// M xmin;
      u = 0   on 2,3,4,5; // other boundaries but xmax
    };

   solve(v) in M
     krylov(type=cg,precond=diagonal),
     memory(matrix=none)
   {
     test(V)
	int(P/dt*v*V)

       +int(mu*grad(v)*grad(V))
       =
        int(1/dt*convect([u0,v0,w0],dt,v0)*V)

       -int(dy(deltaP)*V)
      ;

      v = 0 on 0,2,3,4,5;
    };

   solve(w) in M
     krylov(type=cg,precond=diagonal),
     memory(matrix=none)
   {
     test(W)
	int(P/dt*w*W)

       +int(mu*grad(w)*grad(W))
       =
        int(1/dt*convect([u0,v0,w0],dt,w0)*W)

       -int(dz(deltaP)*W)
      ;

      w = 0 on 0,2,3,4,5;
    };

    solve(p) in M2
      krylov(precond=ichol)
    {
      test(q)
	 int(grad(p)*grad(q))
       = int(grad(p0)*grad(q))
       - int(1/dt*dx(u)*q)
       - int(1/dt*dy(v)*q)
       - int(1/dt*dz(w)*q);
       p = 0 on 1;
    }

   u0 = u;
   v0 = v;
   w0 = w;
   p1 = p0;
   p0 = p;

   if (i<10) {
     save(vtk,"ns.0".i,{[u,v,w],p},M);
   } else {
     save(vtk,"ns.".i,{[u,v,w],p},M);
   }
    i = i+1;

   double fx = int[M](ki/eps*u);
   double fy = int[M](ki/eps*v);
   double fz = int[M](ki/eps*w);
   
   cout << "<#i=" << i << "\n";
   cout << "<#Fx=" << fx << "\n";
   cout << "<#Fy=" << fy << "\n";
   cout << "<#Fz=" << fz << "\n";

} while(i <= itsteps);


mesh T = tetrahedrize(Omega,M);

save(medit,"velocity",T);
save(medit,"velocity",[u,v,w],T);
save(medit,"pressure",T);
save(medit,"pressure",p0,T);