// example 5
// Poisson equation: : a toy problem to check the error estimate
// with P1 or P1b element 
// COE-tutorial 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);
fespace Vh2(Th2,P1b);
// higher polynomial to avoid errors form numerical integration
fespace Vh0(Th0,P2);	

// external force
func ff = (-5.0)/8.0 * pi * pi * sin(pi * x) * sin(pi * y / 2.0); 

// finite element solutions and test functions
Vh1 u1,v1;
Vh2 u2,v2;
// heigher order interpolation to avoid error caused by numerical integration
Vh0 xx, ww;

// Dirichlet boundary conditions
Vh1 w1;
Vh2 w2;
w1 = sin(pi * x) * sin(pi * y / 2.0) / (-2.0);
w2 = sin(pi * x) * sin(pi * y / 2.0) / (-2.0);

//  Poisson problem on finite element space Th1
solve poisson1(u1,v1) = 
int2d(Th1)( dx(u1)*dx(v1)+dy(u1)*dy(v1) ) 
-  int2d(Th1)( ff*v1 ) //
+  on(1,2,3,4,u1=w1);

//  Poisson problem on finite element space Th2
solve poisson2(u2,v2) = 
int2d(Th2)( dx(u2)*dx(v2)+dy(u2)*dy(v2) ) 
-  int2d(Th2)( ff*v2 ) //
+  on(1,2,3,4,u2=w2);

// exact solution
// heigher order interpolation to avoid error caused by numerical integration
ww = sin(pi * x) * sin(pi * y / 2.0) / (-2.0);

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

xx = u1 - ww;
// H1-norm
err1 = int2d(Th0)( dx(xx)*dx(xx) + dy(xx)*dy(xx) + xx*xx );
err1 = sqrt(err1);

xx = u2 - ww;
// H1-norm
err2 = int2d(Th0)( dx(xx)*dx(xx) + dy(xx)*dy(xx) + xx*xx );
err2 = sqrt(err2);

err0 = int2d(Th0)( dx(ww)*dx(ww) + dy(ww)*dy(ww) + ww*ww );
err0 = sqrt(err0);


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