// Setup Problem Domain
real a = 0.5;  // 1, 0.5, 0.25, 0.125, 0.0625
real b = 1.0;

real t = 0.025;
real ah = a-t;  // 1, 0.5, 0.25, 0.125, 0.0625
real bh = b-t;

int N1 = 20;
int N2 = round(N1*b/a);

// Mesh
// "Inside" is to the "left" of boundary
border C1(t=-a/2., a/2.) { x=t; y=-b/2.; }
border C2(t=-b/2., b/2.) { x=a/2.; y=t; }
border C3(t=-a/2., a/2.) { x=-t; y=b/2; }
border C4(t=-b/2., b/2.) { x=-a/2.; y=-t; }

border Ch1(t=-ah/2., ah/2.) { x=t; y=-bh/2.; }
border Ch2(t=-bh/2., bh/2.) { x=ah/2.; y=t; }
border Ch3(t=-ah/2., ah/2.) { x=-t; y=bh/2; }
border Ch4(t=-bh/2., bh/2.) { x=-ah/2.; y=-t; }

mesh Th=buildmesh( C1(N1) + C2(N2) + C3(N1) + C4(N2) +
		   Ch1(-N1) + Ch2(-N2) + Ch3(-N1) + Ch4(-N2) ); 

// fespace
fespace Vh(Th, P1);
Vh phi,v,u;
real dth1 = 1, G = 1;
func flux = -2*G*dth1;

// Stress Function Formulation. Wrong.
solve poisson(phi, v) =
  int2d(Th)( dx(phi)*dx(v) + dy(phi)*dy(v) )
  + int2d(Th)( flux*v )
  + on(C1,C2,C3,C4,Ch1,Ch2,Ch3,Ch4, phi=0);

// Warping Formulation
solve poisson2(u, v) =
  int2d(Th)( dx(u)*dx(v) + dy(u)*dy(v))
  - intalledges(Th)( dth1*(y*N.x - x*N.y)*v );      
  // - int1d(Th, C1,C2,C3,C4,Ch1,Ch2,Ch3,Ch4 )( dth1*(y*N.x - x*N.y)*v );  

real umean = int2d(Th)(u)/int2d(Th)(1);
Vh upl;
upl = u - umean;

cout << endl << umean << endl << endl;

// Visualizations
// plot(phi, value=true, fill=true, dim=3, wait=true);
 plot(upl, value=true, fill=true, dim=3);