// NBPROC 4
assert(version>=2.23);
if ( mpisize != 4 )
  {
    cout << " sorry number of processeur !=4 " << endl;
    exit(1);
  }

//  Mortar  (4 sub domain)
//  with matrix -et Precon Conjugade Gradient --
//    Neuman -> Dirichlet .
//  -------------------------------

func f=1+x+y;
real g=1;
int withprecon=1;
macro Grad(u) [ dx(u), dy(u) ] //
int nbsd=4;

macro Psd(U) U[0],U[1],U[2],U[3] //

int labext= nbsd+1;
real meshsize=0.03;
real meshsizem=meshsize*1.5;
bool noconforme=0;

mesh Tha,Thm,Thmm;
int [int] regsd(4);


if(mpirank==0)
  {

    mesh Thacopy; //hack to dcl Tha ouside of mortar-msh.idp
    {
      include "mortar-msh.idp"
	Thacopy=Tha;
      regsd=regi;
    }
    Tha=Thacopy;
    cout << "mortar : " << endl;
    Thm=Tha;
    Thm=adaptmesh(Thm,meshsizem,IsMetric=1,thetamax=60);
    Thm=emptymesh(Thm);       Thmm=Thm;
    Thmm=trunc(Thm,split=4,1); // for fine integration
    Thmm=emptymesh(Thmm);
    plot(Thm,wait=0,ps="mortar-Thm.eps");
  }

// send meshes
verbosity=1;
broadcast(processor(0),Tha);
broadcast(processor(0),Thm);
broadcast(processor(0),Thmm);
broadcast(processor(0),regsd);

// build sub domaine mesh :
mesh Thi=trunc(Tha,region==regsd[mpirank],split=1);// les sous domaines

if(noconforme)
 Thi=adaptmesh(Thi,meshsize*(1+0.1*mpirank),IsMetric=1,nbvx=100000,thetamax=60);// les sous domaines


fespace Lh(Thm,P1);
fespace RTh(Thm,[P0edge,P0edge]);
 RTh [Nmx,Nmy]; // ne marche pas car la normal
//                             n'est definie que une un bord
varf  vNN([ux,uy],[nx,ny]) = int1d(Thm,1)(( nx*N.x + ny*N.y)/lenEdge);
Nmx[]= vNN(0,RTh);

// les joint P0 sur le squelette
// -----   \int [q] l + \int[p] m
Lh  lh=0,rhsl=0;

fespace Vhi(Thi,P1);
fespace Ehi(Thi,P0);
matrix Asd,Csd,PAsd,PIsd,PJsd;
Vhi usd,vsd,rhssd, pusd,bcsd;
Ehi epssd;

real tgv=1e30;



varf cci([l],[u]) = int1d(Thmm,1,qforder=3)(l*u*epssd);
varf vepsi(u,v)= int1d(Thi,1,qforder=10)( (Nmx*N.x + Nmy*N.y)*v/lenEdge);

varf vLapMi([ui],[vi],tgv=tgv) =
			   int2d(Thi)(  Grad(ui)'*Grad(vi)  )
			   // + int1d(Thi,1,qfe=qf1pElump)(alpha*ui*vi)
			   +  int2d(Thi) (f*vi)   +  on(labext,ui=g);

varf vPLapMi([ui],[vi],tgv=tgv) =
			   int2d(Thi)(  Grad(ui)'*Grad(vi)  )
			   // + int1d(Thi,1,qfe=qf1pElump)(alphap*ui*vi)
			   + on(labext,1,ui=0);
;

varf  vrhsMi(ui,vi) =   on(labext,ui=g);


usd=0;
vsd=0;
epssd[]= vepsi(0,Ehi);
epssd = -real(epssd <-0.00001) + real(epssd >0.00001);

Csd  = cci(Lh,Vhi);
Asd  = vLapMi(Vhi,Vhi,solver="SPARSESOLVER");
PAsd = vPLapMi(Vhi,Vhi,solver="SPARSESOLVER");
matrix IVL=interpolate(Vhi,Lh,inside=1);
//   v = IVL*l
varf vonext(u,v)=on(labext,u=1);
varf von1(u,v)=on(1,u=1);
real[int]  onext=vonext(0,Vhi);
real[int]  on1=von1(0,Vhi);
on1= on1 ? 1 : 0;
on1 = onext ? 0 : on1; //  remove df of ext
matrix I1(on1);//   matrix    tgv $i\in Gamma_1 \ Gamma_e $ , 0 otherwise
PIsd=  I1*IVL;//  remove of line not on $Gamma_1 \ Gamma_e $
// so PIsd*l  =  tgv * Interpole l on $Gamma_1 \ Gamma_e $
I1.diag=on1;
matrix AA=I1*Asd;//  remove line not on lab 1
PJsd= IVL'*AA;
rhssd[]=vLapMi(0,Vhi);

varf  vML(u,v) = int2d(Thm)(u*v*1e-10)+int1d(Thm,1)(u*v);
matrix ML=vML(Lh,Lh);


lh[]=0;
int itera=0;

varf  vbc(u,v) =   int1d(Thm,labext)(v);
real[int] lbc(Lh.ndof),lbc0(Lh.ndof);
lbc=vbc(0,Lh);
lbc = lbc ? 0 : 1 ;

int what; //  to choose which  funct call 1  SkPb 2:PSkPb, 3:end
// the preconditionner
func real[int] PSkPb(real[int] &l)
{
  if(withprecon)
    {
      if(mpirank==0)
	{
	  what=2;
	  broadcast(processor(0),what); // get from SkBk routine
	}
      if(what!=2) return l;

      int verb=verbosity;   verbosity=0;   itera++;
      real[int] ll= ML^-1*l;
      broadcast(processor(0),ll);
      ll= lbc .* ll;
      ll *= tgv;
      pusd[] = PAsd^-1*(vsd[]= PIsd* ll);
      ll  = PJsd*pusd[];
      if(mpirank==0)
	{
	  for (int i=1;i<4;++i)
	    {
	      processor(i) >> l;
	      ll += l;
	    }
	  l = ML^-1*ll;
	  l= lbc .* l;
	}
      else
	processor(0) << ll;
      verbosity=verb;
    }
  return l ;
};

func real[int] SkPb(real[int] &l)
{
  int verb=verbosity;   verbosity=0;   itera++;
  if(mpirank==0 && what!=3) what=1;
  broadcast(processor(0),what);
  if(what==2) return PSkPb(l);
  else if (what !=1) return l;
  broadcast(processor(0),l);

  vsd[]  = rhssd[];
  vsd[] += Csd* l;
  usd[]  = Asd^-1*vsd[];
  l  = Csd'*usd[];
  l= lbc .* l;

  if(mpirank==0)
   {
    real[int] ll(l.n);
    for (int i=1;i<4;++i)
      {
        processor(i) >> ll;
        l += ll;
      }
    }
  else
    processor(0) << l;
  verbosity=verb;
  return l ;
};





if(mpirank==0)
  {
    verbosity=100;
    lh[]=0;
     LinearCG(SkPb,lh[],eps=1.e-5,nbiter=100,precon=PSkPb);
    what=3;
    SkPb(lh[]);
  }
 else
   while(what!=3)
     SkPb(lh[]);

plot(usd,bb=[[-1,-1],[1,1]],ps="mortar-"+mpirank+".eps");
cout << "Fin CG  " << mpirank <<  endl;

/*
Brochet:examples/mpi hecht$ (grep -vh showpage  mortar-?.eps;echo showpage) > mortar.eps
Brochet:examples/mpi hecht$ gv mortar.eps
*/
