load "MUMPS_seq"
//load "MUMPS"

mesh Th = square(10, 10);
int[int] Wall = [1, 2, 3, 4];

fespace Uh(Th, P1);
Uh u, uh;

real f = 1.;

macro grad(u) [dx(u), dy(u)] //

solve Laplacian(u, uh, solver=sparsesolver)
	= int2d(Th)(
		grad(u)' * grad(uh)
	)
	- int2d(Th)(
		f * uh
	)
	+ on(Wall, u=0)
	;

plot(u);