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

load "PETSc"                        // PETSc plugin
macro dimension()3// EOM            // 2D or 3D
include "macro_ddm.idp"             // additional DDM functions

macro def(i)[i, i#B, i#C, i#D]//    // vector field definition
macro init(i)[i, i, i, i]// EOM     // vector field initialization
macro grad(u)[dx(u), dy(u), dz(u)]//// three-dimensional gradient
real Sqrt = sqrt(2.);
macro div(u)(dx(u) + dy(u#B) + dz(u#C))// EOM
func Pk = [P2, P2, P2, P1];         // finite element space

string solver;
if(!HasType("MATSOLVER", "mumps") && !HasType("MATSOLVER", "superlu"))
    exit(0);
else
    solver = (HasType("MATSOLVER", "mumps") ? "mumps" : "superlu");

mesh3 Th;
{
    mesh ThGlobal2d = square(getARGV("-global", 12), getARGV("-global", 12), [x, y]);    // global mesh
    ThGlobal2d = trunc(ThGlobal2d, (x <= 0.5) || (y <= 0.5), label = 5);
    ThGlobal2d = trunc(ThGlobal2d, (y >= 0.25) || (x >= 0.25), label = 5);
    mesh Th2d = movemesh(ThGlobal2d, [-x, y]);
    ThGlobal2d = ThGlobal2d + Th2d;
    Th = buildlayers(ThGlobal2d, getARGV("-global", 12) / 2, zbound = [0, 0.4]);
}
Mat A;
macro ThRefinementFactor()getARGV("-split", 1)//
MatCreate(Th, A, Pk);

fespace Wh(Th, Pk);                 // local finite element space
varf vPb([u, uB, uC, p], [v, vB, vC, q]) = int3d(Th)(grad(u)' * grad(v) + grad(uB)' * grad(vB) + grad(uC)' * grad(vC) - div(u) * q - div(v) * p + 1e-10 * p * q) + on(0, 1, 3, 5, u = 0, uB = 0, uC = 0) + on(2, u = 1000*y*(0.5-y)*z*(0.4-z), uB = 0, uC = 0);
real[int] rhs = vPb(0, Wh, tgv = -1);
Wh<real> def(u) = [1.0, 1.0, 1.0, 2.0];
string[int] names(2);
names[0] = "velocity";
names[1] = "pressure";

A = vPb(Wh, Wh, tgv = -1);
set(A, sparams = "-ksp_type fgmres -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type lower -pc_fieldsplit_detect_saddle_point -pc_fieldsplit_schur_precondition selfp -fieldsplit_velocity_sub_pc_type lu " + " -fieldsplit_pressure_sub_pc_type lu -fieldsplit_velocity_sub_pc_factor_mat_solver_type " + solver + " -fieldsplit_pressure_sub_pc_factor_mat_solver_type " + solver + " -ksp_monitor -ksp_view " + " -fieldsplit_velocity_ksp_type gmres -fieldsplit_velocity_ksp_max_it 5 -fieldsplit_pressure_ksp_type gmres -fieldsplit_pressure_ksp_max_it 5 -ksp_rtol 1e-6", fields = u[], names = names);
u[] = 0.0;
u[] = A^-1 * rhs;
macro def1(u)u//
plotMPI(Th, u, P2, def1, real, cmm = "Global solution");
