//  run with MPI:  ff-mpirun -np 2 script.edp
// NBPROC 2

/* to run : in a Shell.  cut and past:
cd /Users/hecht/work/freefem++/examples/mpi 
mpirun -np 2 FreeFem++-mpi schwarz-3.edp -glut ffglut
*/
if ( mpisize != 2 ) {
cout << " sorry number of processeur !=2 " << endl;
exit(1);}
verbosity=4;
real pi=4*atan(1);
int inside = 2;
int outside = 1;
real zmin =0, zmax=1;
border a(t=1,2){x=t;y=0;label=outside;};
border b(t=0,1){x=2;y=t;label=outside;};
border c(t=2,0){x=t ;y=1;label=outside;};
border d(t=1,0){x = 1-t; y = t;label=inside;};
border e(t=0, pi/2){ x= cos(t); y = sin(t);label=inside;};
border e1(t=pi/2, 2*pi){ x= cos(t); y = sin(t);label=outside;};
int n=1;
mesh th2,TH2;
mesh3 th,TH;
int[int] rup=[0,0], rdown=[0,0], rmid=[outside,outside,inside,inside];;
if (mpirank == 0)
 {
 th2 = buildmesh( a(5*n) + b(5*n) + c(10*n) + d(5*n));
 th =  buildlayers(th2,n,
	     zbound=[zmin,zmax],  labelmid=rmid,
	     reffaceup = rup,     reffacelow = rdown);

 cout << " end th  " << endl;
}
else
 {
 TH2 = buildmesh ( e(5*n) + e1(25*n) );
 TH =  buildlayers(TH2,n,
	     zbound=[zmin,zmax],  labelmid=rmid,
	     reffaceup = rup,     reffacelow = rdown);

 cout << " end TH  " << endl;
 }
broadcast(processor(0),th);
broadcast(processor(1),TH);


fespace vh(th,P1);
fespace VH(TH,P1);
fespace ph(th,P1);
fespace PH(TH,P1);
vh u=0,v; VH U,V;
ph chi=1;
PH CHI=1;
int i=0;

problem PB(U,V,init=i,solver=CG) =
  int3d(TH)( dx(U)*dx(V)+dy(U)*dy(V) +dz(U)*dz(V))
  + int3d(TH)( -V) + on(inside,U = u)  +    on(outside,U= 0 ) ;
problem pb(u,v,init=i,solver=CG) =
    int3d(th)( dx(u)*dx(v)+dy(u)*dy(v)+dz(u)*dz(v) )
  + int3d(th)( -v) + on(inside ,u = U) +    on(outside,u = 0 ) ;


for ( i=0 ;i< 10; i++)
{
  cout << mpirank << " looP " << i << endl;
  if (mpirank == 0)
   PB;
  else
    pb;
   broadcast(processor(0),U[]);
   broadcast(processor(1),u[]);
 
  if (mpirank==0)     plot(U,u,wait=1,fill=1);

  if(mpirank==0)
    {
      real ERR= int3d(TH)( square(U-u)*CHI*chi) ;
      cout << " ERR= " << ERR<< endl;
    }
  else
    {
      real err= int3d(th)( square(U-u)*CHI*chi) ;
      cout << " err= " << err<< endl;
    }

};
 if (mpirank==0)
    plot(U,u,ps="Uu.eps",fill=1);
