// Resolution of
// - \Delta u = f 	in \Omega
//	u=0		on \Gamma_D
//	du/dn=g 	on \Gamma_N

int N=3;	// Nb of edges on the boundary
func f=1;	// volumic flux
func g=-1;	// boundary flux

macro vh() P1	// finite element space used

mesh Sh;
int GammaN,GammaD;
{
	GammaN=1;GammaD=2;
	int[int] labels=[GammaN,GammaD,GammaD,GammaN];
	Sh=square(N,N,label=labels);
}

fespace Vh(Sh,vh);	// finite element space

macro grad(u) [dx(u),dy(u)] //	gradient operator
varf a(u,v)=int2d(Sh)(grad(u)'*grad(v))+on(GammaD,u=0);
varf L(u,v)=int1d(Sh,GammaN)(g*v);

matrix A=a(Vh,Vh);	// Building matrix
real[int] b=L(0,Vh);	// Building second member

cout<<"The Matrix="<<A<<endl;
cout<<"The second member="<<b<<endl;

Vh u=0;			// FE and allow memory space
u[]=A^-1*b;		// solving system

plot(u,cmm="potential",wait=1);