// NBPROC 4
// ff-mpirun -np 4 DDM-Schwarz-Lap-3d.edp -glut ffglut  -n 11 -k 1  -d 1 -ns -gmres 1
/*
  a first true parallele example fisrt freefem++ 
  Ok up to 200 proc for a Poisson equation.. 
  See the Doc for full explaiantion

  F Hecht Dec. 2010. 
  -------------------
usage :
ff-mpirun [mpi parameter] MPIGMRES2d.edp  [-glut ffglut]  [-n N] [-k K]  [-d D] [-ns] [-gmres [0|1]
 argument: 
   -glut ffglut : to see graphicaly the process
   -n N:  set the mesh cube split NxNxN
   -d D:  set debug flag D must be one for mpiplot 
   -k K:  to refined by K all  elemnt
   -ns: reomove script dump
   -gmres 0   : use iterative schwarz algo.  
          1   :  Algo GMRES on residu of schwarz algo.
          2   :  DDM GMRES 
          3   :  DDM GMRES with coarse grid preconditionner (Good one)  
*/
load "MPICG"  load "medit"  load "metis"
include "getARGV.idp"
include "DDM-Schwarz-macro.idp"
//include "AddLayer3d.idp"
include "DDM-funcs-v2.idp"
include "cube.idp" 


searchMethod=0; // more safe seach algo (warning can be very expensive in case lot of ouside point) 
// 0 by default
assert(version >=3.11);
real[int] ttt(10);int ittt=0;
macro settt {ttt[ittt++]=mpiWtime();}//


verbosity=getARGV("-vv",0);
int vdebug=getARGV("-d",1);
int ksplit=getARGV("-k",1);
int nloc = getARGV("-n",5);
string sff=getARGV("-p","");
int gmres=getARGV("-gmres",3); 
int nC = getARGV("-N" ,max(nloc/10,5)); 
int sizeoverlaps=1; // size of overlap
bool RAS=1; 

if(mpirank==0 && verbosity)
  cout << " vdebug: " << vdebug << " kspilt "<< ksplit << " nloc "<< nloc << " sff "<< sff <<"."<< endl;

string sPk="P2-Lame-3gd";     
func Pk=[P2,P2,P2];

func bool  plotMPIall(mesh3 &Th,real[int] & u,string  cm)
{if(vdebug) PLOTMPIALLU(mesh3,Pk, defPk3, Th, u, allu2, { cmm=cm,nbiso=20,fill=1,dim=3,value=1}); return 1;}


mpiComm comm(mpiCommWorld,0,0);// trick : make a no split mpiWorld 


int ipart= mpiRank(comm); // current partition number 

if(ipart==0)  cout << " Final N=" << ksplit*nloc << " nloc =" << nloc << " split =" << ksplit <<  endl;
int[int] l111=[1,1,1,1]; 
settt 


int[int,int] LL=[[1,1],[2,1],[1,1]];
real[int,int] BB=[[0,1],[0,5],[0,1]];
int[int] NN=[nloc,nloc*5,nloc]; 
int[int] NNC=[nC,nC*5,nC]; 


settt 
mesh3 Thg=Cube(NN,BB,LL);
mesh3 ThC=Cube(NNC,BB,LL);

fespace VhC(ThC,[P1,P1,P1]); // of the coarse problem.. 


BuildPartitioning(sizeoverlaps,mesh3,Thg,Thi,aThij,RAS,pii,jpart,comm,vdebug)

if(ksplit>1)
{
for(int jp=0;jp<jpart.n;++jp)
  aThij[jp]  = trunc(aThij[jp],1,split=ksplit);
Thi =   trunc(Thi,1,split=ksplit);
}

BuildTransferMat(ipart,mesh3,Pk,3,[0,1,2],Thi,Whi,Whij,Thij,aThij,Usend,Vrecv,jpart,vdebug)





/* the definition of the Problem .... */


// the definition of the Problem ....
real E = 21.5e4;
real sigma = 0.29;
real mu = E/(2*(1+sigma));
real lambda = E*sigma/((1+sigma)*(1-2*sigma));
real gravity = -0.05;

real sqrt2=sqrt(2.);
macro epsilon(u1,u2,u3)  [dx(u1),dy(u2),dz(u3),(dz(u2)+dy(u3))/sqrt2,(dz(u1)+dx(u3))/sqrt2,(dy(u1)+dx(u2))/sqrt2] // EOM
macro div(u1,u2,u3) ( dx(u1)+dy(u2)+dz(u3) ) // EOM
  
varf vPb([u1,u2,u3],[v1,v2,v3])=
  int3d(Thi)(  
	    lambda*div(u1,u2,u3)*div(v1,v2,v3)	
	    +2.*mu*( epsilon(u1,u2,u3)'*epsilon(v1,v2,v3) ) //')
	      )
  + int3d(Thi) (gravity*v3)
  + on(2,10,u1=0,u2=0,u3=0)
  ;
  
varf vPbC([u1,u2,u3],[v1,v2,v3])=
  int3d(ThC)(  
	    lambda*div(u1,u2,u3)*div(v1,v2,v3)	
	    +2.*mu*( epsilon(u1,u2,u3)'*epsilon(v1,v2,v3) ) //')
	      )
  + int3d(ThC) (gravity*v3)
  + on(2,u1=0,u2=0,u3=0)
  ;

varf vPbon10([u1,u2,u3],[v1,v2,v3])=on(10,u1=1,u2=1,u3=1)+on(2,u1=0,u2=0,u3=0);
varf vPBC(U,V)=on(2,U=0);



real[int] onG10 = vPbon10(0,Whi); // on 1 
real[int] Bi=vPb(0,Whi);


matrix Ai = vPb(Whi,Whi,solver=sparsesolver); 

DMMDeffuncAndGlobals(Lame3,comm,jpart,Whi,Vhc,3,Ai,vPbC,onG10,Pii,Usend,Vrecv,[0,1,2])

Lame3CheckUpdate();
  
Whi [u,u1,u2],[v,v1,v2];
 

u[]=vPBC(0,Whi,tgv=1); 
real eps=1e-10;
Lame3DDMSolver(Bi,u,v,gmres,eps,vdebug)


real errg =1,umaxg;
{ 
  real umax = u[].max,umaxg;
  real[int] aa=[umax], bb(1);
  mpiAllReduce(aa,bb,comm,mpiMAX);
  errg=bb[0];
  if(ipart==0)
    cout << " umax global  = " << bb[0] << " Wtime = " << (ttt[ittt-1]-ttt[ittt-2])  << " s " <<  " " << Lame3kiter <<  endl;
}

Lame3Saveff(sff,eps,ksplit,nloc,sizeoverlaps); 
cout << " On A Fin"<<endl; 
