load "Morley"

// 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; }

border C41(t=0., b/4) { x=0; y=b-t; label=41; }
border C42(t=b/4, 3.*b/4.) { x=0; y=b-t; label=42; }
border C43(t=3.*b/4., b) { x=0; y=b-t; label=43; }

border C21(t=0., b/4) { x=a; y=t; label=21; }
border C22(t=b/4, 3.*b/4.) { x=a; y=t; label=22; }
border C23(t=3.*b/4., b) { x=a; y=t; label=23; }


int N1 = 80;
int N2 = 40;
// mesh Th=buildmesh( C1(N1) + C2(N2) + C3(N1) + C4(N2));
mesh Th=buildmesh( C1(N1) +
		   C21(N2/4) + C22(N2/2) + C23(N2/4) + 
		   C3(N1) +
		   C41(N2/4) + C42(N2/2) + C43(N2/4) );

// fespace
fespace Vh(Th, P2Morley);
Vh [u1, u1x, u1y], [u2, u2x, u2y], [v, vx, vy];

solve biharm([u1, u1x, u1y], [v, vx, vy]) =
  int2d(Th)( dxx(v)*dxx(u1)+2*dxy(v)*dxy(u1)+dyy(v)*dyy(u1) )
  + on(21,22,23, u1=4, u1x=0, u1y=0)  
  + on(41,42,43, u1=0, u1x=0, u1y=0);

solve biharm2([u2, u2x, u2y], [v, vx, vy]) =
  int2d(Th)( dxx(v)*dxx(u2)+2*dxy(v)*dxy(u2)+dyy(v)*dyy(u2) )
  + on(21,22,23, u2=4, u2x=0, u2y=0)
  + on(41,43, u2=-2, u2x=0, u2y=0)
  + on(42, u2=2, u2x=0, u2y=0);

plot(u1, value=true, fill=true, wait=true);
plot(u2, value=true, fill=true, wait=true);
plot(u1, u2, fill=true, dim=3, wait=true);