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

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

macro def(i)[i, i#B, i#C]// EOM     // vector field definition
macro init(i)[i, i, i]// EOM        // vector field initialization
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 Pk = [P2, P2, P1];             // finite element space

mesh Th;
{
    mesh ThGlobal = square(getARGV("-global", 40), getARGV("-global", 40), [x, y]);    // global mesh
    ThGlobal = trunc(ThGlobal, (x < 0.5) || (y < 0.5), label = 5);
    Th = movemesh(ThGlobal, [-x, y]);
    Th = ThGlobal + Th;
}
Mat A;
macro ThRefinementFactor()getARGV("-split", 1)//
MatCreate(Th, A, Pk);

fespace Wh(Th, Pk);                 // local finite element space
varf vPb([u, uB, p], [v, vB, q]) = int2d(Th)(grad(u)' * grad(v) + grad(uB)' * grad(vB) - div(u) * q - div(v) * p + 1e-10 * p * q) + on(1, 3, 5, u = 0, uB = 0) + on(2, u = y*(0.5-y), uB = 0);
real[int] rhs = vPb(0, Wh);

set(A, sparams = "-pc_type lu");
Wh<real> def(u);

A = vPb(Wh, Wh);
u[] = A^-1 * rhs;

macro def2(u)[u, u#B]// EOM
macro def1(u)u// EOM
plotMPI(Th, def2(u), [P2, P2], def2, real, cmm = "Global velocity");
plotMPI(Th, uC, P1, def1, real, cmm = "Global pressure");
