// Setup Problem Domain
real a = 4.0;
real b = 4.0;

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

int Nn = 50;
mesh Th=buildmesh( C1(Nn) + C2(Nn) + C3(Nn) ); 

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

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

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

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

// Warping function