//  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]// EOM          // vector field definition
macro init(i)[i, i]// EOM           // vector field initialization
real Sqrt = sqrt(2.0);
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]; // finite element space

int[int] l = [2, 1, 2, 2];
mesh Th = square(2 * getARGV("-global", 10), getARGV("-global", 10), [2 * x, y], label = l); // global mesh
Mat A;
macro ThRefinementFactor()getARGV("-split", 1)//
MatCreate(Th, A, Pk);

real f = -90000.0;
real strain = 100.0;
real Young = 1.0e8;
real poisson = 0.45;
real tmp = 1.0 + poisson;
real mu = Young  / (2.0 * tmp);
real lambda = Young * poisson / (tmp * (1.0 - 2.0 * poisson));
varf vPb(def(u), def(v)) = int2d(Th)(lambda * div(u) * div(v) + 2.0 * mu * (epsilon(u)' * epsilon(v))) + int2d(Th)(f * vB) + on(1, u = 0.0, uB = 0.0);
fespace Wh(Th, Pk);                 // local finite element space
matrix Loc = vPb(Wh, Wh);
real[int] rhs = vPb(0, Wh);

set(A, sparams = "-ksp_view", bs = 2);
Wh<real> def(u);                    // local solution

A = Loc;
u[] = A^-1 * rhs;
if(mpirank == 0)
    cout << "KSPGetIterationNumber = " << KSPGetIterationNumber(A) << endl;

real[int] err = A * u[];            // global matrix-vector product
exchange(A, rhs, scaled = true);
err -= rhs;

plotMPI(Th, def(u), Pk, def, real, cmm = "Global solution");
u[] = err;
macro def1(u)u// EOM
plotMPI(Th, u, P2, def1, real, cmm = "Global residual");

Wh<real> def(Rb)[3];
[Rb[0], RbB[0]] = [1, 0];
[Rb[1], RbB[1]] = [0, 1];
[Rb[2], RbB[2]] = [y, -x];
set(A, sparams = "-pc_type gamg -ksp_type gmres -ksp_max_it 200", nearnullspace = Rb);
u[] = 0.0;
real[int] residuals(200);
KSPSetResidualHistory(A, residuals);
u[] = A^-1 * rhs;
if(mpirank == 0)
    cout << "KSPResidualHistory = " << residuals(0:KSPGetIterationNumber(A)) << endl;
plotMPI(Th, def(u), Pk, def, real, cmm = "Global solution");
real alpha = 1.0e+1;
mesh ThMoved = movemesh(Th, [x + alpha * u, y + alpha * uB]);
u[] = mpirank;
plotMPI(ThMoved, u, P2, def1, real, cmm = "Global moved solution");
if(HasType("PC", "hpddm")) {
    u[] = 0;
    set(A, sparams = "-pc_type hpddm -pc_hpddm_levels_1_eps_nev 20 -pc_hpddm_levels_1_sub_pc_type lu -pc_hpddm_define_subdomains -pc_hpddm_has_neumann -pc_hpddm_coarse_p 2 -ksp_monitor");
    u[] = A^-1 * rhs;
}
if(HasType("PC", "fancy_bug")) {
    exit(1);
}
if(HasType("KSP", "fancy_bug_bis")) {
    exit(2);
}
