load "Morley"

// Parameter
real f = 1;

// Mesh
mesh Th = square(40, 40); //mesh definition of Omega

// Fespaces
fespace Wh(Th, P2);
fespace Vh(Th, P2Morley);	// HCT finite element space
Vh [u, ux, uy], [v, vx, vy];

// Macro
macro bilaplacien(u, v) (dxx(u)*dxx(v) + dyy(u)*dyy(v) + 2.*dxy(u)*dxy(v))	// end of macro

// Problem
solve bilap ([u, ux, uy], [v, vx, vy])
	= int2d(Th)(bilaplacien(u, v))
	- int2d(Th)(f*v)
	+ on(1, 2, 3, 4, u=0, ux=0, uy=0)
	;

// Plot
plot(u, cmm="u", wait=1, fill=1);
plot(ux, wait=1, cmm="u_x");
plot(uy, wait=1, cmm="u_y");

// Max & Error
Wh uu = u;
real umax = uu[].max;
int err = (abs(umax-0.0012782) > 1e-4);
cout << " uu max = " << umax << " ~ 0.0012782, err = " << err << endl;

// Plot
int n = 100, nn = n+10;
real[int] xx(nn), yy(nn);
for (int i = 0; i <= n; i++) {
	xx[i] = real(i)/n;
	yy[i] = u(0.5, real(i)/n);	// value of uh at point (0.5, i/10.)
}
plot([xx(0:n), yy(0:n)], wait=1);

// End
assert(err == 0);

