// Setup Problem Domain
real b = 1.0;  // 1, 0.5, 0.25, 0.125, 0.0625
real h = 2.0;
real t = 0.1;

real Xc = -3.*b^2/(h+6.*b);

// Mesh
// "Inside" is to the "left" of boundary
border C1(s=0., b+t) { x=b+t/2-Xc-s; y=h/2.+t/2.; }
border C2(s=0., h+t) { x=-t/2.-Xc; y=h/2.+t/2.-s; }
border C3(s=0., b+t) { x=-t/2.-Xc+s; y=-h/2.-t/2.; }
border C4(s=0., t) { x=b+t/2.-Xc; y=-h/2.-t/2.+s; }
border C5(s=0., b) { x=b+t/2.-Xc-s; y=-h/2.+t/2.; };
border C6(s=0., h-t) { x=t/2.-Xc; y=-h/2.+t/2.+s; }
border C7(s=0., b) { x=t/2.-Xc+s; y=h/2.-t/2.; }
border C8(s=0., t) { x=b+t/2.-Xc; y=h/2.-t/2.+s; }

int nNt = 12;
int nNb = 30;
int nNh = 60;

mesh Th=buildmesh( C1(nNb+nNt) + C2(nNh+nNt) + C3(nNb+nNt) + C4(nNt)
		   + C5(nNb) + C6(nNh) + C7(nNb) + C8(nNt) );

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

// Case 1
solve poisson(phi, v) =
  int2d(Th)( dx(phi)*dx(v) + dy(phi)*dy(v) )
  + int2d(Th)( flux*v )
  + on(C1,C2,C3,C4,C5,C6,C7,C8, phi=0);

solve poisson2(u, v) =
  int2d(Th)( dx(u)*dx(v) + dy(u)*dy(v))
  - int1d(Th, C1,C2,C3,C4,C5,C6,C7,C8)( dth1*(y*N.x - x*N.y)*v );

real umean = int2d(Th)(u)/int2d(Th)(1), phimean = int2d(Th)(phi)/int2d(Th)(1);
Vh upl = u - umean;
real urms = sqrt(int2d(Th)(upl^2)), phirms = sqrt(int2d(Th)((phi-phimean)^2));
Vh phipl = phi*(urms/phirms);

Vh Sx = dy(phi), Sy = -dx(phi);

cout << "Hey " << urms << endl << phirms << endl;

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

plot(Sx, value=true, fill=true, dim=3, wait=true);
plot(Sy, value=true, fill=true, dim=3, wait=true);