vector n = (15,15,15);
vertex A = (0,0,0);
vertex B = (1,1,1);

mesh M = structured(n,A,B);

double eps = 1e-3;
double mu  = 1e-1; // viscosity
double dt = 0.5;

solve(u1,u2,u3,p) in M
    method(type=eliminate),krylov(type=bicg),
    bicg(epsilon=1e-3, maxiter=2000)
{
   test(U1,U2,U3,P)
     int(grad(u1)*grad(U1) - p*dx(U1)
        +grad(u2)*grad(U2) - p*dy(U2)
        +grad(u3)*grad(U3) - p*dz(U3)
        +eps*grad(p)*grad(P) + dx(u1)*P + dy(u2)*P + dz(u3)*P)
     = 0;

   u1 = 1 on M ymax;
   u1 = 0 on M xmin, M xmax,         M ymax, M zmin, M zmax;
   u2 = 0 on M xmin, M xmax, M ymin, M ymax, M zmin, M zmax;
   u3 = 0 on M xmin, M xmax, M ymin, M ymax, M zmin, M zmax;
};

function un1 = u1;
function un2 = u2;
function un3 = u3;

double i=0;
for (i=0; i<10; i++)
{
  solve(u1,u2,u3,p) in M
    method(type=eliminate),krylov(type=bicg),
    bicg(epsilon=1e-5, maxiter=2000) {
       test(U1,U2,U3,P)
	    int(1/dt*u1*U1 + mu*grad(u1)*grad(U1) - p*dx(U1)
	       +1/dt*u2*U2 + mu*grad(u2)*grad(U2) - p*dy(U2)
	       +1/dt*u3*U3 + mu*grad(u3)*grad(U3) - p*dz(U3)
	       +eps*grad(p)*grad(P) + eps*p*P + dx(u1)*P + dy(u2)*P + dz(u3)*P)
	   = int(1/dt*convect(un1, un2, un3, dt, un1)*U1
		+1/dt*convect(un1, un2, un3, dt, un2)*U2
		+1/dt*convect(un1, un2, un3, dt, un3)*U3);

	u1 = 1 on M ymax;
	u1 = 0 on M xmin, M xmax,         M ymax, M zmin, M zmax;
	u2 = 0 on M xmin, M xmax, M ymin, M ymax, M zmin, M zmax;
	u3 = 0 on M xmin, M xmax, M ymin, M ymax, M zmin, M zmax;
    };
    
  un1 = u1;
  un2 = u2;
  un3 = u3;
}

save(medit,"u",M);
save(medit,"u",[u1,u2,u3],M);

save(medit,"p",M);
save(medit,"p",p,M);