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

load "hpddm_substructuring"         // HPDDM plugin
macro dimension()2// EOM            // 2D or 3D
include "macro_ddm_substructuring.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
macro BC(u, val)u = val, u#B = val// EOM // Dirichlet boundary conditions
func Pk = [P2, P2];                 // finite element space

int labNeumann = getARGV("-labNeumann", 2);                     // interface between subdomains
int s = getARGV("-split", 1);                                   // refinement factor

mpiComm comm;
int p = getARGV("-hpddm_level_2_p", 1);
bool excluded = splitComm(mpiCommWorld, p, comm, topology = getARGV("-hpddm_level_2_topology", 0), exclude = (usedARGV("-hpddm_level_2_exclude") != -1));

if(verbosity > 0 && mpirank == 0) {
    cout << " --- " << mpirank << "/" << mpisize;
    cout << " - elasticity-2d-substructuring.edp - input parameters: refinement factor = " << s << endl;
}

mesh Th = square(1, 1);
fespace Wh(Th, Pk);           // local finite element space
int[int][int] intersection;   // local-to-neighbors renumbering
int[int] interfaceNb;         //
{
    int[int] l = [2, 1, 2, 2];
    Th = square(2 * getARGV("-global", 10), getARGV("-global", 10), [2 * x, y], label = l);    // global mesh
    buildSubstructuring(Th, interfaceNb, 10, 1, labNeumann, s, intersection, Pk, BC, comm, excluded)
}

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);
matrix<real> Mat = vPb(Wh, Wh, sym = 1); // local operator
real[int] rhs = vPb(0, Wh);                  // local right-hand side

bdd A(Mat, intersection, communicator = comm);
// this can be replaced by feti
bool adaptive = isSetOption("geneo_nu") || isSetOption("geneo_threshold");

pair ret;
if(mpisize == 1)
    renumber(A, Mat, interfaceNb, effort = rhs);
else {
    if(excluded)
        AttachCoarseOperator(mpiCommWorld, A);
    else {
        Wh[int] def(Rb)(0);
        real[int] float(Wh.ndof);
        varf floatingPb(def(u), def(v)) = on(1, BC(u, 1.0));
        float = floatingPb(0, Wh);
        if(float.max < 0.9 && !adaptive) {
            Rb.resize(3);
            [Rb[0], RbB[0]] = [1, 0];
            [Rb[1], RbB[1]] = [0, 1];
            [Rb[2], RbB[2]] = [y, -x];
        }
        if(getARGV("-hpddm_substructuring_scaling", 0) != 2)
            renumber(A, Mat, interfaceNb, R = Rb, effort = rhs);
        else {
            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));
            Wh def(muFunc) = mu;
            renumber(A, Mat, interfaceNb, R = Rb, effort = rhs, rho = muFunc[]);
        }
        AttachCoarseOperator(mpiCommWorld, A, R = Rb, ret = ret);
    }
}

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

if(mpisize == mpiSize(comm))
    u[] = A^-1 * rhs;
else
    DDM(A, rhs, u[], excluded = excluded);

if(!excluded)
    OriginalNumbering(A, u[], interfaceNb);
plotMPI(Th, def(u), Pk, def, real, cmm = "Global solution");
