// for the make check:
// NBPROC 4
// PARAM -raspart -ffddm_schwarz_method oras -waven 20 -nsuper 2 -nw -ns

macro dimension 2// 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 = P1;

real k = getARGV("-waven",20.);
func f = 100*exp(-((x-.5)^2+(y-.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[int] chlab=[1, Robin, 2, Robin, 3, Robin, 4, Robin];
//int[int] chlab=[1, Dirichlet, 2, Dirichlet, 3, Dirichlet, 4, Dirichlet];

real npplo = 15.;

mesh Th = square(npplo*1./lambda,npplo*1./lambda);    // global mesh

Th=change(Th,refe=chlab);

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

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

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

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

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

int nsuper = getARGV("-nsuper",2);

macro Hsuperwithchunk()1 // EOM
ffddmbuildSuper(H,Th,nsuper,complex,def,init,Pk)

ffddmsetupsuperinit(H,H)

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

ffddmsetupOperator(H,H,Varf)

Hchunkverbosity = 1;

ffddmsetupPrecondchunk(H,VarfOpt)

ffddmsetupPrecond(Hchunk,VarfOpt)

complex[int] rhs(1);

ffddmbuildrhs(H,Varfrhs,rhs)

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

HVhi<complex> def(u), def(err);

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

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

Hwritesummary

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