 //ff-mpirun -np 4 Helmholtz-3d-simple.edp -wg -raspart -ffddm_schwarz_method oras  -ffddm_medit

// If you have openmpi you may need to add the option --oversubscribe to allow more processes than the number of cores available on your computer

// for the make check:
// NBPROC 4
// PARAM -raspart -ffddm_schwarz_method oras -waven 6

macro dimension 3// EOM            // 2D or 3D

include "ffddm.idp"


macro def(i)i// EOM                         // scalar field definition
macro init(i)i// EOM                        // scalar field initialization
func Pk = P2;  

real k = getARGV("-waven",10.);
func f = exp(-((x-.5)^2+(y-.5)^2+(z-.5)^2)*k*10);

real lambda = 2*pi/k;

real epsilonprob = 0;
real be = getARGV("-betaEprec",1.);
real epsilonEprec = 0;//k^be;

int Dirichlet = 1;
int Robin = 2;

int mysplit = 3;

real npplo = 10./mysplit;

//int[int] LL = [Dirichlet,Dirichlet, Dirichlet,Dirichlet, Dirichlet,Dirichlet];
int[int] LL = [Robin,Robin, Robin,Robin, Robin,Robin];

mesh3 Thc = cube(npplo*1./lambda,npplo*1./lambda,npplo*1./lambda, [x, y, z], label = LL);

mesh3 Th = trunc(Thc,1,split=mysplit);

macro Grad(u) [dx(u),dy(u),dz(u)] // EOM

macro Varf(varfName, meshName, PhName)
    varf varfName(u,v) =
  int3d(meshName)(-(k^2-1i*epsilonprob)*u*v+Grad(u)'*Grad(v))
				+ int2d(meshName,Robin)(1i*k*u*v)
				+ on(Dirichlet,u=0); // EOM

// for the preconditioner
macro VarfEprec(varfName, meshName, PhName)
    varf varfName(u,v) =
  int3d(meshName)(-(k^2-1i*epsilonEprec)*u*v+Grad(u)'*Grad(v))
				+ int2d(meshName,Robin)(1i*k*u*v)
				+ on(Dirichlet,u=0); // EOM

// for the preconditioner
macro VarfOpt(varfName, meshName, PhName)
    varf varfName(u,v) =
  int3d(meshName)(-(k^2-1i*epsilonEprec)*u*v+Grad(u)'*Grad(v))
				+ int2d(meshName,Robin,10)(1i*k*u*v)
				+ on(Dirichlet,u=0); // EOM

macro Varfrhs(varfName, meshName, PhName)
    varf varfName(u,v) =
                -int3d(meshName)(f*v)
                + on(Dirichlet,u=0); // EOM

macro Hdefmplot(u)real(u)//

vsym = 2; // symmetric, not hermitian
vtgv = 1.e+30;
vtgvelim = 1.e+30;

ffddmbuild(H,Th,complex,def,init,Pk,mpiCommWorld)

ffddmsetup(H,H,Varf,VarfOpt)

complex[int] rhs(1);

ffddmbuildrhs(H,Varfrhs,rhs)

complex[int] x0(rhs.n);
x0 = 0;

HVhi<complex> def(u);

if (mpirank == 0) cout << "ORAS :" << endl;

u[] = HfGMRES(x0, rhs, 1.e-6, 200, "right");

Hwritesummary

if (mpirank == 0) cout <<  endl << "ORAS TWO-LEVEL :" << endl;

ffddmcoarsemeshsetup(H,Thc,VarfEprec,VarfEprec)

u[] = HfGMRES(x0, rhs, 1.e-6, 200, "right");

Hwritesummary

ffddmplot(H,real(u), "Global solution");
