//  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]// 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[int] LL = [2,3, 2,1, 2,2];
mesh3 Th =  cube(10 * getARGV("-global", 5), getARGV("-global", 5), getARGV("-global", 5), [10 * x, y, z], label = LL);
Mat A;
macro ThRefinementFactor()getARGV("-split", 1)//
MatCreate(Th, A, Pk);

fespace Wh(Th, Pk);                 // local finite element space
real f = -9000.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)) = int3d(Th)(lambda * div(u) * div(v) + 2.0 * mu * (epsilon(u)' * epsilon(v))) + int3d(Th)(f * vC) + on(1, u = 0.0, uB = 0.0, uC = 0.0);

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

Wh<real> def(Rb)[6];
[Rb[0], RbB[0], RbC[0]] = [1, 0, 0];
[Rb[1], RbB[1], RbC[1]] = [0, 1, 0];
[Rb[2], RbB[2], RbC[2]] = [0, 0, 1];
[Rb[3], RbB[3], RbC[3]] = [y, -x, 0];
[Rb[4], RbB[4], RbC[4]] = [-z, 0, x];
[Rb[5], RbB[5], RbC[5]] = [0, z, -y];
set(A, sparams = "-pc_type gamg -ksp_type gmres -ksp_pc_side right -ksp_max_it 200 -pc_gamg_threshold 0.01", nearnullspace = Rb, bs = 3);
func real[int] funcRes(real[int]& inPETSc) {
    ChangeNumbering(A, u[], inPETSc, inverse = true, exchange = true);
    real[int] out(Wh.ndof);
    out = A * u[];
    real[int] rhs = vPb(0, Wh, tgv = -2);
    out -= rhs;
    real[int] outPETSc;
    ChangeNumbering(A, out, outPETSc);
    return outPETSc;
}
func int funcJ(real[int]& inPETSc) {
    A = vPb(Wh, Wh, tgv = -2);
    return 0;
}
real[int] xPETSc;
u[] = 0;
ChangeNumbering(A, u[], xPETSc);
SNESSolve(A, funcJ, funcRes, xPETSc, sparams = "-snes_monitor -snes_linesearch_monitor -snes_type newtonls");
ChangeNumbering(A, u[], xPETSc, inverse = true, exchange = true);
macro def1(u)u// EOM
plotMPI(Th, u, P1, def1, real, cmm = "Global solution");
real alpha = 1.0;
mesh3 ThMoved = movemesh3(Th, transfo = [x + alpha * u, y + alpha * uB, z + alpha * uC]);
u[] = mpirank;
plotMPI(ThMoved, u, P1, def1, real, cmm = "Global moved solution");
