// solve the Poisson equation in 69 lines (comments and blank lines included)
//  -Delta u=f 		on Omega
//  u=0			on partial Omega
// usinng Matrix Syntax + Domain Decomposition
//
// Note that such an approach is usefull if you want to solve
// multiphysic problems (for instance fluid/structure interaction)

real R=1,r=0.5,dn=10;
int Dirichlet=1,Interface=2;
func f=1;
mesh Th,Th1,Th2,Sh;
int reg1,reg2;
{
	border a(t=0,2.*pi){x=R*cos(t);y=R*sin(t);label=Dirichlet;};
	border b(t=0,2.*pi){x=r*cos(t);y=r*sin(t);label=Interface;};
	mesh Th=buildmesh(a(2.*pi*R*dn)+b(2.*pi*r*dn));
	// retrive region labels
		reg2=Th(0,0).region;
		reg1=Th(0,(r+R)/2.).region;
	// Extract submeshed Th1 and Th2
		Th1=trunc(Th,region==reg1,label=Interface);
		Th2=trunc(Th,region==reg2,label=Interface);
	// emptyMesh for the interface
		Sh=emptymesh(Th2);
}
plot(Th1,Th2,wait=1,cmm="Mesh of both domain");
plot(Th1,Sh,wait=1,cmm="The empty mesh of Th2");

// Definition of the variationnal formulation (to build matrices)
	macro grad(u) [dx(u),dy(u)] //
	varf a1(u,v)=int2d(Th1)(grad(u)'*grad(v))+on(Dirichlet,u=0);
	varf a2(u,v)=int2d(Th2)(grad(u)'*grad(v));
	varf c1(p,u)=int1d(Sh,Interface)(p*u);
	varf c2(p,u)=int1d(Sh,Interface)(-p*u);

// Definition of the variationnal formulation (to build second member)
	varf L1(u,v)=int2d(Th1)(f*v);
	varf L2(u,v)=int2d(Th2)(f*v);

// Definition of the fespace
fespace Vh1(Th1,P1),Vh2(Th2,P1),Rh(Sh,P1);

// Build Matrix of the system
matrix A;
{
	matrix A1=a1(Vh1,Vh1),A2=a2(Vh2,Vh2),C1=c1(Rh,Vh1),C2=c2(Rh,Vh2);
	A=[[A1 ,0  ,C1],
	   [0  ,A2 ,C2],
	   [C1',C2',0]];
	set(A,solver=UMFPACK);
}

// Build second member
real[int] b(A.n);
{
	real[int] b1=L1(0,Vh1), b2=L2(0,Vh2);
	real[int] b3(Rh.ndof);b3=0;
	b=[b1,b2,b3];
}

// sovle the system
Vh1 u1; Vh2 u2;
{
	real[int] sol=A^-1*b;
	u1=0;u1[]=sol(0:Vh1.ndof-1);
	u2=0;u2[]=sol(Vh1.ndof:Vh2.ndof+Vh1.ndof-1);
}
plot(u1,u2,cmm="The solution");