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

load "hpddm"                        // HPDDM plugin
load "medit"
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
/*# DiffMacros #*/
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
/*# DiffMacrosEnd #*/
func Pk = [P1, P1, P1];             // finite element space

/*# DDMoptions #*/
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
/*# SchwarzMethod #*/
int[int][int] intersection;   // local-to-neighbors renumbering
real[int] D;                  // partition of unity
{
    mesh3 ThBorder;
    Th = cube(10 * getARGV("-global", 5), getARGV("-global", 5), getARGV("-global", 5), [10 * x, y, z], label = LL);      // global mesh
    buildOverlap(Th, ThBorder, 10, s, 1, intersection, D, Pk, mpiCommWorld, false)
}

real f = -9000.0;
real strain = 100.0;
real Young = 2.0e11; // steel
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));
real[int] rhs;                              // local right-hand side
matrix<real> Mat;                           // local operator
varf vPb(def(u), def(v)) = intN(Th)(lambda * div(u) * div(v) + 2.0 * mu * (epsilon(u)' * epsilon(v))) + on(1, u = 0.0, uB = 0.0, uC = 0.0);
Mat = vPb(Wh, Wh, sym = 1, tgv = -2);

schwarz A(Mat, intersection, D);
set(A, sparams = "-hpddm_schwarz_method ras -hpddm_schwarz_coarse_correction balanced -hpddm_variant right -hpddm_verbosity 1 -hpddm_geneo_nu 20");

{
    if(mpisize > 1 &&
       isSetOption("schwarz_coarse_correction")) { // two-level Schwarz methods
        varf vPbNoPen(def(u), def(v)) = intN(Th)(lambda * div(u) * div(v) + 2.0 * mu * (epsilon(u)' * epsilon(v))) + on(1, u = 0.0, uB = 0.0, uC = 0.0);
        matrix<real> noPen = vPbNoPen(Wh, Wh, sym = 1, tgv = -2);
        AttachCoarseOperator(mpiCommWorld, A, A = noPen);
    }
}

real alpha = 5.0e3;
Wh def(u), def(d);
real[int] uBlock(u[].n * 2);
real[int] rhsBlock(u[].n * 2);

varf vPbC(def(u), def(v)) = intN(Th)(f * vC) + on(1, u = 0.0, uB = 0.0, uC = 0.0);
rhsBlock(0 * u[].n:1 * u[].n - 1) = vPbC(0, Wh, tgv = -2);

varf vPbB(def(u), def(v)) = intN(Th)(f * vB) + on(1, u = 0.0, uB = 0.0, uC = 0.0);
rhsBlock(1 * u[].n:2 * u[].n - 1) = vPbB(0, Wh, tgv = -2);

uBlock = A^-1 * rhsBlock;

u[] = uBlock(0 * u[].n:1 * u[].n - 1);
d[] = uBlock(1 * u[].n:2 * u[].n - 1);
mesh3 ThMoved1 = movemesh3(Th, transfo = [x + alpha * u, y + alpha * uB, z + alpha * uC]);
mesh3 ThMoved2 = movemesh3(Th, transfo = [x + alpha * d, y + alpha * dB, z + alpha * dC]);
macro def1(u)u// EOM
plotMPI(ThMoved1, u, P1, def1, real, cmm = "Global moved solution");
plotMPI(ThMoved2, u, P1, def1, real, cmm = "Global moved solution");
