//  Mortar  (4 sub domain)
//  with matrix ---
//  ------------
verbosity=2;
real cpu=clock();

macro laplacien(u,v) (dx(u)*dx(v)+dy(u)*dy(u)) //

// --- begin  meshes  building --------------
int nbsd=4;
int labext= nbsd+1;
real[int] theta(nbsd+1),cost(nbsd),sint(nbsd);

for (int i=0;i<nbsd;i++)
 {
  real t=i*2*pi/nbsd;
  theta[i]= t;
  theta[i+1]= (i+1)*2*pi/nbsd;
  cost[i]=cos(t);
  sint[i]=sin(t);
 }


border g1(t=0,1){x=cost[0]*t;y=sint[0]*t;label=1;};
border g2(t=0,1){x=cost[1]*t;y=sint[1]*t;label=1;};
border g3(t=0,1){x=cost[2]*t;y=sint[2]*t;label=1;};
border g4(t=0,1){x=cost[3]*t;y=sint[3]*t;label=1;};

border e12(t=theta[0],theta[1]){x=cos(t);y=sin(t);label=labext;};
border e23(t=theta[1],theta[2]){x=cos(t);y=sin(t);label=labext;};
border e34(t=theta[2],theta[3]){x=cos(t);y=sin(t);label=labext;};
border e41(t=theta[3],theta[4]){x=cos(t);y=sin(t);label=labext;};

 int Ng = 10;
 int Ne = 10;

plot(g1(Ng)+g2(Ng)+g3(Ng)+g4(Ng) + e12(Ne) + e23(Ne)+ e34(Ne) + e41(Ne) ,wait=1);

mesh Tha = buildmesh( g1(Ng)+g2(Ng)+g3(Ng)+g4(Ng) + e12(Ne) + e23(Ne)+ e34(Ne) + e41(Ne) );

int [int] regi(4);
for (int i=0;i<4;i++)
  {
    real tt = (theta[i]+theta[i+1])*0.5;
    regi[i]=Tha(0.5*cos(tt),0.5*sin(tt)).region;
  }
int reg1= regi[0];
int reg2= regi[1];
int reg3= regi[2];
int reg4= regi[3];
cout << " reg " << regi << endl;

real meshsize=0.1;
Tha=adaptmesh(Tha,0.1,IsMetric=1);
Tha=adaptmesh(Tha,meshsize,IsMetric=1);
mesh Thaa=splitmesh(Tha,4);


fespace Pha(Tha,P0);
Pha reg=region;
real[int] VVV=[1,2,3,4,5,6,7,8,9];
plot(reg,fill=1,grey=1,viso=VVV);
//  definition des espaces de joints
Pha UNN; // defini l'element des triangles touchant les joints 1
varf vNNN(u,v) = int1d(Tha,1)( v); 
UNN[]=vNNN(0,Pha);
UNN= abs(UNN)>0.01; 
plot(UNN,fill=1,wait=1,cmm="UNN");
mesh Thm=trunc(Tha,UNN!=0,label=0); 
mesh Thmm=trunc(Thaa,UNN!=0,label=0,split=1); 
plot(Thmm,wait=1);
verbosity=2;

plot(Thm,wait=1);
mesh Th0=trunc(Tha,region==reg1,split=1);// les sous domaines
mesh Th1=trunc(Tha,region==reg2,split=1);// 
mesh Th2=trunc(Tha,region==reg3,split=1);// 
mesh Th3=trunc(Tha,region==reg4,split=1);// 

Th0=adaptmesh(Th0,meshsize*0.5,IsMetric=1);
Th1=adaptmesh(Th1,meshsize*0.6,IsMetric=1);
Th2=adaptmesh(Th2,meshsize*0.7,IsMetric=1);
Th3=adaptmesh(Th3,meshsize*0.8,IsMetric=1);

fespace Vha(Tha,P1);
//Vhaa ssd0=region==reg1,ssd1=region==reg2,ssd2=region==reg3,ssd3=region==reg4;


fespace Lh(Thm,P1);
fespace RTh(Tha,[P0edge,P0edge]);
 RTh [Nmx,Nmy]=[N.x,N.y]; // ne marche pas car la normal 
//                             n'est definie que une un bord

varf  vNN([ux,uy],[nx,ny]) = int1d(Tha,1)(( nx*N.x + ny*N.y)/lenEdge);
Nmx[]= vNN(0,RTh);


plot([Nmx,Nmy],wait=1,cmm="Nmx,Nmy");

// les joint P0 sur le squelette 
// -----   \int [q] l + \int[p] m 
Lh  lh,rhsl;
lh=0; 
rhsl=0;
macro AA(u,v)  (dx(u)*dx(v)+dy(u)*dy(v)) // 

// remark: operator # is  the concatenation operator in macro 
macro defspace(i)

   cout << " Domaine " << i<< " --------" << endl;
   fespace Vh#i(Th#i,P1);
   fespace Eh#i(Th#i,P0edge);


   Vh#i u#i;
   Vh#i rhs#i;

varf veps#i(u,v)= int1d(Th#i,1,qforder=5)( (Nmx*N.x + Nmy*N.y)*v/lenEdge);
   Eh#i eps#i = 0;
   eps#i[]= veps#i(0,Eh#i);
   plot(eps#i,wait=1);
   eps#i = -real(eps#i <-0.01) + real(eps#i >0.01);

   varf vLapM#i([u#i],[v#i]) =
       int2d(Th#i)(  AA(u#i,v#i) )
    +  int2d(Th#i) (f*v#i)
    
    +  on(labext,u#i=g)
   ;  

   varf  vrhsM#i(u#i,v#i) =   on(labext,u#i=g);

   varf cc#i([l],[u]) = int1d(Thmm,1)(l*u*eps#i);
   cout << "Build mat C and A "<< endl;
   matrix C#i = cc#i(Lh,Vh#i);
   matrix A#i = vLapM#i(Vh#i,Vh#i,solver=GMRES);
   cout << " ... done "<< endl;
   rhs#i[]=vLapM#i(0,Vh#i);
//  fin macro defspace(i) 


varf vfmortar(l,m) = int1d(Thm,1)(l*m);
matrix Mm = vfmortar(Lh,Lh);




func f=1+x+y;
real g=1; 


defspace(0)
defspace(1)
defspace(2)
defspace(3)

  plot(eps0,eps1,eps2,eps3,cmm="eps 0,1,2,3",wait=1,value=1);


lh[]=0;
varf  vDD(u,v) = int2d(Thm)(u*v*1e-10);

verbosity=4;
int iter=0;
//

matrix DD=vDD(Lh,Lh);

matrix M=[ 
  [ A0 ,0  ,0  ,0  ,C0 ],
  [ 0  ,A1 ,0  ,0  ,C1 ],
  [ 0  ,0  ,A2 ,0 , C2 ],
  [ 0  ,0  ,0  ,A3, C3 ],
  [ C0',C1',C2',C3',DD ] 
 ];

real[int] xx(M.n);
real[int] bb(M.n);

real[int] bbb=[rhs0[], rhs1[],rhs2[],rhs3[],rhsl[] ];
bb=[rhs0[], rhs1[],rhs2[],rhs3[],rhsl[] ];
set(M,solver=UMFPACK);
xx = M^-1 * bb;
[u0[],u1[],u2[],u3[],lh[]] = xx; // buf



Vha vah,uah;
solve vLapMM([uah],[vah]) =
       int2d(Tha)(  AA(uah,vah) )
    -  int2d(Tha) (f*vah)
    +  on(labext,uah=g)
   ;
verbosity =3;
plot(uah,cmm="uah",wait=1); 
plot(u0,u1,u2,u3,cmm="u1,u2,u3,u4",wait=1);