//ff-mpirun -np 8 Maxwell-3d-simple.edp -wg -raspart -ffddm_schwarz_method oras -ffddm_medit

// for the make check:
// NBPROC 4
// PARAM -raspart -ffddm_schwarz_method oras -global 10

load "Element_Mixte3d"

macro dimension 3// EOM

include "ffddm.idp"


macro def(i)[i, i#y, i#z]// EOM             // vector field definition
macro init(i)[i, i, i]// EOM                // vector field initialization
macro defPart(u)u// EOM                     // partition of unity definition
macro initPart(u)u// EOM                    // partition of unity initialization
macro Curl(ux, uy, uz)[dy(uz)-dz(uy), dz(ux)-dx(uz), dx(uy)-dy(ux)]// EOM
macro CrossN(ux, uy, uz)[uy*N.z-uz*N.y, uz*N.x-ux*N.z, ux*N.y-uy*N.x]// EOM
macro Curlabs(ux, uy, uz)[abs(dy(uz)-dz(uy)), abs(dz(ux)-dx(uz)), abs(dx(uy)-dy(ux))]// EOM
func Pk = Edge03d;
func PkPart = Edge03ds0;

int Dirichlet = 1, Robin = 2;

int ng = getARGV("-global", 30);

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

mesh3 Th = cube(ng, ng, ng, [x,y,z], label = labs);
func k = 6 * pi;

func f = exp(-((x-0.5)^2+(y-0.5)^2+(z-0.5)^2)*60);

int mysplit = 3;

mesh3 Thc = cube(ng/mysplit, ng/mysplit, ng/mysplit, [x,y,z], label = labs);

real epsilonprob = 0;//0 k k^2;                  // epsilon of the problem
real be = getARGV("-betaEprec",1.); // beta for epsilonEprec  
real epsilonEprec = 0;//k^be; 

//macro Mdefplot(u)sqrt(real(u)^2+real(u#y)^2+real(u#z)^2)//

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

ffddmbuildDmesh(M,Th,mpiCommWorld)
ffddmbuildDfespaceEdge(M,M,complex,def,init,Pk,defPart,initPart,PkPart)

macro Varf(varfName, meshName, PhName)
    varf varfName([Ex,Ey,Ez],[vx,vy,vz]) =
  int3d(meshName)(Curl(vx,vy,vz)'*Curl(Ex,Ey,Ez))
                + int3d(meshName)(-(k^2-1i*epsilonprob)*[vx,vy,vz]'*[Ex,Ey,Ez])
                + int2d(meshName,Robin)(1i*k*CrossN(vx,vy,vz)'*CrossN(Ex,Ey,Ez))
                + on(Dirichlet,Ex=0,Ey=0,Ez=0);    // EOM

// for the preconditioner
macro VarfEprec(varfName, meshName, PhName)
    varf varfName([Ex,Ey,Ez],[vx,vy,vz]) =
  int3d(meshName)(Curl(vx,vy,vz)'*Curl(Ex,Ey,Ez))
                + int3d(meshName)(-(k^2-1i*epsilonEprec)*[vx,vy,vz]'*[Ex,Ey,Ez])
                + int2d(meshName,Robin)(1i*k*CrossN(vx,vy,vz)'*CrossN(Ex,Ey,Ez))
                + on(Dirichlet,Ex=0,Ey=0,Ez=0);    // EOM
                
// for the preconditioner
macro VarfOpt(varfName, meshName, PhName)

    varf varfName([Ex,Ey,Ez],[vx,vy,vz]) =
  int3d(meshName)(Curl(vx,vy,vz)'*Curl(Ex,Ey,Ez))
                + int3d(meshName)(-(k^2-1i*epsilonEprec)*[vx,vy,vz]'*[Ex,Ey,Ez])
                + int2d(meshName,Robin)(1i*k*CrossN(vx,vy,vz)'*CrossN(Ex,Ey,Ez))
                + int2d(meshName,10)(1i*k*CrossN(vx,vy,vz)'*CrossN(Ex,Ey,Ez))
                + on(Dirichlet,Ex=0,Ey=0,Ez=0);    // EOM

macro Varfrhs(varfName, meshName, PhName)
    varf varfName([Ex,Ey,Ez],[vx,vy,vz]) =
        -int3d(meshName)([vx,vy,vz]'*[0,0,f])
       + on(Dirichlet,Ex=0,Ey=0,Ez=0);  // EOM

ffddmsetup(M,M,Varf,VarfOpt)

complex[int] rhs(1);

ffddmbuildrhs(M,Varfrhs,rhs)

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

MVhi<complex> def(u);

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

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

Mwritesummary

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

ffddmcoarsemeshsetup(M,Thc,VarfEprec,VarfEprec)

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

Mwritesummary

ffddmplot(M,sqrt(real(u)^2+real(uy)^2+real(uz)^2), "Global solution");
