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

// Mesh
// "Inside" is to the "left" of boundary
border C(t=0., 2.*pi) { x=a*cos(t); y = b*sin(t); }

int Nn = 100;
mesh Th=buildmesh( C(Nn) );

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

// Stress function
solve poisson(phi, v) =
  int2d(Th)( dx(phi)*dx(v) + dy(phi)*dy(v) )
  + int2d(Th)( flux*v )
  + on(C, phi=0);

// Axial Displacement
solve poisson2(u, v) =
  int2d(Th)( dx(u)*dx(v) + dy(u)*dy(v) )
  + int1d(Th, C)( 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, wait=true);

// Warping function