if (HaveUMFPACK) //  to heavy computation if not 
{
load "BernardiRaugel"

real s0=clock();

// Parameters
real reylnods = 400.;
func g = (x)*(1-x)*4;
real nu = 1;
real dt = 0.1;
real penalty = 1.e-6;
int nbiter = 3;
real coefdt = 0.25^(1./nbiter);
real coefcut = 0.25^(1./nbiter);
real cut = 0.01;
real tol = 0.3;
real coeftol = 0.25^(1./nbiter);

// Mesh
mesh Th = square(2, 2);

// Fespace
fespace Vh2(Th, P2BR);
fespace Vh(Th, P0);
fespace Wh(Th, [P2BR, P0]);
Wh [u1, u2, p], [v1, v2, q], [up1, up2, pp];

// Macro
macro div(u1, u2) (dx(u1) + dy(u2)) //

// Problem
real alpha = 0;
int i = 0;

varf NS ([u1, u2, p], [v1, v2, q], init=i)
	= int2d(Th)(
		  alpha*(u1*v1 + u2*v2)
		+ nu*(
			  dx(u1)*dx(v1)
			+ dy(u1)*dy(v1)
			+ dx(u2)*dx(v2)
			+ dy(u2)*dy(v2)
		)
		- penalty * p*q
		- p*div(v1, v2)
		- div(u1, u2)*q
	)
	- int2d(Th)(
		- alpha*convect([up1, up2], -dt, up1)*v1
		- alpha*convect([up1, up2], -dt, up2)*v2
	)
	+ on(3, u1=g, u2=0)
	+ on(1, 2, 4, u1=0, u2=0)
	;

// Initialization
[up1, up2, pp] = [0., 0., 0.];

// Solve (Stokes)
matrix A = NS(Wh, Wh, solver="SPARSESOLVER");
real[int] b = NS(0, Wh);
u1[] = A^-1 * b;

// Plot
plot([u1,u2], p, coef=0.2, cmm="[u1, u2] and p", value=true, wait=true);

// Convergence loop
int iter = 0;
nu = 1./reylnods;
for (iter = 1; iter <= nbiter; iter++){
	cout << "dt = " << dt << endl;
	alpha = 1./dt;
	
	A = NS(Wh, Wh, solver="SPARSESOLVER");
	
	// Time loop
	for (i = 0; i <= 10; i++){
		// Update
		up1[] = u1[];
		real[int] b = NS(0, Wh);
		
		// Solve
		u1[] = A^-1 * b;
		
		// Plot
		if (!(i % 10))
			plot([u1,u2], p, coef=0.2, cmm="[u1,u2] and p");
		
		// Display
		cout << "CPU " << clock()-s0 << "s " << endl;
	} 
	
	if (iter >= nbiter) break;
	
	// Mesh adaptation
	Th = adaptmesh(Th, [u1, u2], iso=0, abserror=0, cutoff=cut, err=tol, inquire=0, ratio=1.5, hmin=1./1000);
	plot(Th);
	
	// Update
	dt = dt*coefdt;
	tol = tol *coeftol;
	cut = cut *coefcut;
	
	// Re-interpolation
	[u1, u2, p] = [u1, u2, p];
	[up1, up2, pp] = [u1, u2, p];
}

// Display
cout << "CPU " << clock()-s0 << "s " << endl;
}


