// Parameters
func f = 1;
func g = 0;

int nMeshAdapt = 5;	//Warning: higher value is time consuming and require `nbvx`>default value in `adaptmesh`

// Mesh
border aaa(t=0, 1){x=t; y=0; label=1;}
border bbb(t=0, 0.5){x=1; y=t; label=2;}
border ccc(t=0, 0.5){x=1-t; y=0.5; label=3;}
border ddd(t=0.5, 1){x=0.5; y=t; label=4;}
border eee(t=0.5, 1){x=1-t; y=1; label=5;}
border fff(t=0, 1){x=0; y=1-t; label=6;}
mesh Th = buildmesh(aaa(6) + bbb(4) + ccc(4) +ddd(4) + eee(4) + fff(6));
mesh Th2 = movemesh(Th,[x+1.1, y]);

// Fespace
fespace Vh(Th, P1);
Vh u, v, zero;

fespace Vh2(Th2, P2);
Vh2 u2, v2;

fespace Vh21(Th2, P1);
Vh21 dxu2, dyu2;

// Initialization
u = 0;
u2 = 0;
zero = 0;

// Problem
int i = 0;
problem Probem1(u, v, solver=CG, init=i, eps=-1.0e-6)
	= int2d(Th)(
		  dx(u)*dx(v)
		+ dy(u)*dy(v)
	) 
	+ int2d(Th)(
		  v*f
	)
	+ on(1, 2, 3, 4, 5, 6, u=g)
	;

problem Probem2(u2, v2, solver=CG, init=i, eps=-1.0e-6)
	= int2d(Th2)(
		  dx(u2)*dx(v2)
		+ dy(u2)*dy(v2)
	)
	+ int2d(Th2)(
		  v2*f
	)
	+ on(1, 2, 3, 4, 5, 6, u2=g)
	;

// Mesh adaptation loop
real cpu=clock();
real error = 0.1;
real coef = 0.1^(1./5.);
for (i = 0; i < 5; i++) { 
	cout << "Iteration " << i << endl;
	real d = clock();
	
	// Solve
	Probem1;
	Probem2;
	
	// Plot
	plot(Th, u, Th2, u2, wait=true);
	
	// Mesh adpatation
	Th = adaptmesh(Th, u, inquire=1, err=error);
	
	dxu2 = dx(u2);
	dyu2 = dy(u2);
	Th2 = adaptmesh(Th2, dxu2, dyu2, inquire=1, err=error, ratio=1.5, hmin=1e-5);

	cout << " CPU = " << clock()-d << endl;
	error = error * coef;
}

cout << " CPU = " << clock()-cpu << endl;

