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

int N=10;	// Nb of edges on the boundary
func f=1;	// volumic flux
real k1=1,k2=10;// conductivity in Omega1 and Omega2

macro vh() P1	// finite element space used

mesh Sh;
int GammaD;
int Omega1,Omega2;
{
	real r=0.2;		// radius of inner inclusion
	real cx=0.5,cy=0;	// center of the inclusion
	GammaD=1;
	border a(t=0,2.*pi){x=cos(t);y=sin(t);label=GammaD;}
	border b(t=0,2.*pi){x=r*cos(t)+cx;y=r*sin(t)+cy;label=2;}
	Sh=buildmesh(a(N*2*pi)+b(N*2*pi*r));
	Omega1=Sh(0,1).region;
	Omega2=Sh(cx,cy).region;
}

fespace Vh(Sh,vh);	// finite element space
Vh u,v;			// state and test functions

macro grad(u) [dx(u),dy(u)] //	gradient operator
problem Laplacian(u,v)=int2d(Sh,Omega1)(k1*grad(u)'*grad(v))+int2d(Sh,Omega2)(k2*grad(u)'*grad(v))-int2d(Sh)(f*v)+on(GammaD,u=0);

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