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

load "PETSc"                  // SLEPc plugin
include "macro_ddm.idp"       // additional DDM functions

macro def2(u)[u, u#B]// EOM
macro init2(u)[u, u]// EOM
macro def1(u)u// EOM
macro init1(u)u// EOM
macro grad(u)[dx(u), dy(u)]// EOM // two-dimensional gradient
real Sqrt = sqrt(2.);
macro epsilon(u)[dx(u), dy(u#B), (dy(u) + dx(u#B)) / Sqrt]// EOM
macro div(u)(dx(u) + dy(u#B))// EOM
func PkV = [P2, P2];          // finite element space for the velocities
func PkP = P1;                // finite element space for the pressure

int s = getARGV("-split", 1); // refinement factor

mesh Th = square(getARGV("-global", 40), getARGV("-global", 40), [x, y]);
fespace VhV(Th, PkV);         // local finite element space for the velocities
fespace VhP(Th, PkP);         // local finite element space for the pressure
DmeshCreate(Th);
Mat A, M;
{
    macro def(u)def2(u)//
    macro init(u)init2(u)//
    MatCreate(Th, A, PkV);
}
{
    macro def(u)def1(u)//
    macro init(u)init1(u)//
    MatCreate(Th, M, PkP);
}

varf vPbA([u, uB], [v, vB]) = int2d(Th)(grad(u)' * grad(v) + grad(uB)' * grad(vB)) + on(1, 2, 4, u = 0, uB = 0) + on(3, u = 1.0, uB = 0);
varf vPbB([q], [u, uB]) = int2d(Th)(-div(u) * q);
varf vPbM(p, q) = int2d(Th)(-p * q);
A = vPbA(VhV, VhV, tgv = -2);
M = vPbM(VhP, VhP);
Mat B(A, M);
Mat Z(A);
{
    B = vPbB(VhP, VhV);
    real[int] diag(VhV.ndof);
    matrix zeros = diag;
    Z = zeros;
    varf onGamma([u, uB], [v, vB]) = on(1, 3, 5, u = 1, uB = 1) + on(2, u = 1, uB = 1);
    VhV def2(g);
    g[] = onGamma(0, VhV, tgv = -1);
    real[int] gPETSc;
    ChangeNumbering(A, g[], gPETSc);
    MatZeroRows(B, gPETSc);
}
Mat LHS = [[A , B],
           [B', 0]];
Mat RHS = [[Z, 0],
           [0, M]];
int nV, nP;
{
    real[int] in(VhV.ndof), out;
    ChangeNumbering(A, in, out);
    nV = out.n;
    in.resize(VhP.ndof);
    ChangeNumbering(M, in, out);
    nP = out.n;
}
real[int] list(nV + nP);
list(0:nV - 1) = 1.0;
list(nV:nV + nP - 1) = 2.0;
real[int] nev(1);
EPSSolve(LHS, RHS, sparams = "-eps_gen_non_hermitian -st_pc_type fieldsplit -st_fieldsplit_1_pc_type jacobi -st_fieldsplit_0_pc_type hypre -st_ksp_type gmres -st_ksp_rtol 1e-1 " + " -eps_error_relative ::ascii_info_detail -st_ksp_converged_reason -eps_target_magnitude -st_type sinvert -eps_target 1.0e-5 -eps_nev 1", fields = list, values = nev);
assert(nev[0] > 0);
