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

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

int N1 = 100;
int N2 = 50;
mesh Th=buildmesh( C1(N1) + C2(N2) + C3(N1) + C4(N2)); 

// 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, phi=0);

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

real umean = int2d(Th)(u)/int2d(Th)(1);
Vh upl;
upl = u - umean*2;
real uplmean = (int2d(Th)(u^4)/int2d(Th)(1))^(1./4.);

cout << "Hey " << umean << endl << uplmean << endl;

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

// plot(phi, value=true, fill=true, dim=2, wait=true);

// Warping function