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

load "hpddm"                        // HPDDM plugin
macro dimension()3// EOM            // 2D or 3D
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
real Sqrt = sqrt(2.0);
macro epsilon(u)[dx(u), dy(u#B), dz(u#C), (dz(u#B) + dy(u#C)) / Sqrt, (dz(u) + dx(u#C)) / Sqrt, (dy(u) + dx(u#B)) / Sqrt]// EOM
macro div(u)(dx(u) + dy(u#B) + dz(u#C))// EOM
func Pk = [P1, P1, P1];             // finite element space

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

int[int] LL = [2,3, 2,1, 2,2];
mesh3 Th = cube(1, 1, 1, [x, y, z]);
fespace Wh(Th, Pk);           // local finite element space
int[int][int] intersection;   // local-to-neighbors renumbering
real[int] D;                  // partition of unity
{
    Th = cube(6 * getARGV("-global", 5), getARGV("-global", 5), getARGV("-global", 5), [6 * x, y, z], label = LL);
    build(Th, s, intersection, D, Pk, mpiCommWorld)
}

real f = -9000.0;
real strain = 100.0;
real Young = 2.0e11;
real poisson = 0.35;
real tmp = 1.0 + poisson;
real mu = Young  / (2.0 * tmp);
real lambda = Young * poisson / (tmp * (1.0 - 2.0 * poisson));
int dirichlet = 1;
varf vPb(def(u), def(v)) = int3d(Th)(lambda * div(u) * div(v) + 2.0 * mu * (epsilon(u)' * epsilon(v))) + int3d(Th)(f * vC) + on(dirichlet, u = 0.0, uB = 0.0, uC = 0.0);
matrix<real> Mat = vPb(Wh, Wh, tgv = -1);
real[int] rhs = vPb(0, Wh, tgv = -1);

schwarz A(Mat, intersection, D);
set(A, sparams = "-hpddm_schwarz_coarse_correction balanced -hpddm_geneo_nu 10");

int solver = getOption("schwarz_method");
if(mpisize > 1 && isSetOption("schwarz_coarse_correction")) { // two-level Schwarz methods
    matrix<real> noPen = vPb(Wh, Wh, sym = 1);
    AttachCoarseOperator(mpiCommWorld, A, A = noPen);
}

Wh<real> def(u);    // local solution

u[] = A^-1 * rhs;

real[int] err(u[].n);
err = A * u[];      // global matrix-vector product
err -= rhs;

mpiBarrier(mpiCommWorld);

macro def1(u)u// EOM
plotMPI(Th, u, P1, def1, real, cmm = "Global solution");
real alpha = 20000.0;
mesh3 ThMoved = movemesh3(Th, transfo = [x + alpha * u, y + alpha * uB, z + alpha * uC]);
u[] = err;
plotMPI(Th, u, P1, def1, real, cmm = "Global residual");
u[] = mpirank;
plotMPI(ThMoved, u, P1, def1, real, cmm = "Global moved solution");
