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

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

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

// fespace
fespace Vh(Th, P1);
Vh u1,u2,v;
func flux2 = .5*(y-b/2.)^2;

// Case 1
func T1 = 1;
solve laplacian1(u1, v) =
  int2d(Th)( dx(u1)*dx(v)+dy(u1)*dy(v) )
  + int2d(Th)( v*0 ) - int1d(Th, C1,C3)( 0*v )
  - int1d(Th, C2)( flux2*v )
  + on(C4, u1=T1);

// Case 2
func T2 = (y>b/4.)*(y<3.*b/4.)*2;
// func T2 = 2.*y/b;
// func T2 = 3.*pi/2.*sin(3.*pi*y/2./b);
solve laplacian2(u2, v) =
  int2d(Th)( dx(u2)*dx(v)+dy(u2)*dy(v) )
  + int2d(Th)( v*0 ) - int1d(Th, C1,C3)( 0*v )
  - int1d(Th, C2)( flux2*v )
  + on(C4, u2=T2);

// Visualizations
plot(u1, u2, value=true, fill=true, dim=3, WindowIndex=0);
plot(u1, value=true, fill=true, WindowIndex=1);
plot(u2, value=true, fill=true, WindowIndex=2);