// example 6
// Stokes equations : a toy problem to check the error estimate
// P1b/P1 element
// COE-tutorail 2007, Atsushi Suzuki 2007/03/08

// finite element meshes
mesh Th1=square(20,20); 
mesh Th2=square(40,40);
mesh Th0=square(160,160);

// finite element spaces
fespace Vh1(Th1,P1b),Qh1(Th1,P1);
fespace Vh2(Th2,P1b),Qh2(Th2,P1);
fespace Vh0(Th0,P2),Qh0(Th0,P2);

// external force
func f1 = 5.0/8.0 * pi * pi * sin(pi * x) * sin(pi * y / 2.0) + 2.0 * x;
func f2 = 5.0/4.0 * pi * pi * cos(pi * x) * cos(pi * y / 2.0) + 2.0 * y;

// finite element solutions and test functions
Vh1 u11,u12,v11,v12;
Qh1 p1,q1;
Vh2 u21,u22,v21,v22;
Qh2 p2,q2;
// heigher order interpolation to avoid error caused by numerical integration 
Vh0 uu1, uu2, vv1, vv2;
Qh0 pp, qq;

// Dirichlet boundary conditions
Vh1 w11, w12;
Vh2 w21, w22;

w11 = sin(pi * x) * sin(pi * y / 2.0) / 2.0;
w12 = cos(pi * x) * cos(pi * y / 2.0);

w21 = sin(pi * x) * sin(pi * y / 2.0) / 2.0;
w22 = cos(pi * x) * cos(pi * y / 2.0);

real epsln = 1.0e-6; 
// definitions of macros for strain rate tensor
macro d11(u1)     dx(u1)// 
macro d22(u2)     dy(u2)//
macro d12(u1,u2)  (dy(u1) + dx(u2))/2.0//

// Stokes problem on finite element space Th1 
solve stokes1(u11,u12,p1, v11,v12,q1) = 
int2d(Th1)(
	2.0*(d11(u11)*d11(v11)+2.0*d12(u11,u12)*d12(v11,v12)+d22(u12)*d22(v12))
	- p1*dx(v11) - p1*dy(v12) 
	- dx(u11)*q1 - dy(u12)*q1
    	- p1*q1*epsln
) 
-  int2d(Th1)( f1 * v11 + f2 * v12 ) // 
+  on(1,2,3,4,u11=w11,u12=w12);

real area, meanp;
area = int2d(Th1)(1.0);
meanp = int2d(Th1)(p1) / area;
p1 = p1 - meanp;

plot([u11,u12],p1,wait=1,value=true,coef=0.1);

// Stokes problem on finite element space Th2
solve stokes2(u21,u22,p2, v21,v22,q2) = 
int2d(Th2)(
	2.0*(d11(u21)*d11(v21)+2.0*d12(u21,u22)*d12(v21,v22)+d22(u22)*d22(v22))
	- p2*dx(v21) - p2*dy(v22) 
	- dx(u21)*q2 - dy(u22)*q2
    	- p2*q2*epsln
) 
-  int2d(Th2)( f1 * v21 + f2 * v22) // 
+  on(1,2,3,4,u21=w21,u22=w22);

area = int2d(Th2)(1.0);
meanp = int2d(Th2)(p2) / area;
p2 = p2 - meanp;

plot([u21,u22],p2,wait=1,value=true,coef=0.1);

// exact solution
// heigher order interpolation to avoid error caused by numerical integration 
uu1 = sin(pi * x) * sin(pi * y / 2.0) / 2.0;
uu2 = cos(pi * x) * cos(pi * y / 2.0);
pp = x*x+y*y - 2.0/3.0;

real err0,err1,err2,hh1,hh2;
Vh1 h1 = hTriangle;
Vh2 h2 = hTriangle;
hh1 = h1[].max;
hh2 = h2[].max;

vv1 = u11 - uu1;
vv2 = u12 - uu2;
qq = p1 - pp;
// (H1-norm for velocity)^2 + (L2-norm for pressure)^2 
err1 = int2d(Th0)( dx(vv1)*dx(vv1) + dy(vv1)*dy(vv1) + vv1*vv1 +
       		   dx(vv2)*dx(vv2) + dy(vv2)*dy(vv2) + vv2*vv2 +
		   qq*qq);
err1 = sqrt(err1);

vv1 = u21 - uu1;
vv2 = u22 - uu2;
qq = p2 - pp;
// (H1-norm for velocity)^2 + (L2-norm for pressure)^2
err2 = int2d(Th0)( dx(vv1)*dx(vv1) + dy(vv1)*dy(vv1) + vv1*vv1 +
       		   dx(vv2)*dx(vv2) + dy(vv2)*dy(vv2) + vv2*vv2 +
		   qq*qq);
err2 = sqrt(err2);

cout << "coarse mesh: h=" << hh1 << " err-H1=" << err1 << endl;
cout << "fine mesh:   h=" << hh2 << " err-H1=" << err2 << endl;
cout << "O(h)=" << log(err1/err2)/log(hh1/hh2) << endl;