//ff-mpirun -np 8 elasticity-3d-deflateRBMs.edp -wg -ffddm_schwarz_method ras -ffddm_schwarz_coar_correction BNN -global 10

// If you have openmpi you may need to add the option --oversubscribe to allow more processes than the number of cores available on your computer

// for the make check:
// NBPROC 4
// PARAM -ffddm_schwarz_method ras -ffddm_schwarz_coarse_correction BNN -global 5

macro dimension 3// EOM            // 2D or 3D

include "ffddm.idp"

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 ThGlobal = cube(6 * getARGV("-global", 5), getARGV("-global", 5), getARGV("-global", 5), [6 * x, y, z], label = LL);

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;

macro Varf(varfName, meshName, PhName)
    varf varfName(def(u), def(v)) = int3d(meshName)(lambda * div(u) * div(v) + 2.0 * mu * (epsilon(u)' * epsilon(v))) + int3d(meshName)(f * vC) + on(dirichlet, u = 0.0, uB = 0.0, uC = 0.0); // EOM

vtgv = -2;
vtgvelim = -2;
vsym = 1;

ffddmbuild(E,ThGlobal,real,def,init,Pk,mpiCommWorld)

ffddmsetup(E,E,Varf,null)

real[int] rhs(1);

ffddmbuildrhs(E,Varf,rhs)

real[int] x0(rhs.n);
x0 = 0;

EVhi def(u);

if (mpirank == 0) cout << "RAS :" << endl; 

u[] = EfGMRES(x0, rhs, 1.e-6, 200, "right");

Ewritesummary

if (mpirank == 0) cout << "RAS + deflate rigid-body modes :" << endl; 

def(u) = [1, 0, 0]; ffddmdeflateglobalvector(E,u[]);
def(u) = [0, 1, 0]; ffddmdeflateglobalvector(E,u[]);
def(u) = [0, 0, 1]; ffddmdeflateglobalvector(E,u[]);
def(u) = [y, -x, 0]; ffddmdeflateglobalvector(E,u[]);
def(u) = [-z, 0, x]; ffddmdeflateglobalvector(E,u[]);
def(u) = [0, z, -y]; ffddmdeflateglobalvector(E,u[]);

ffddmbuildZTAZ(E,E,null)
Ecorr = ffddmcorrection;

u[] = EfGMRES(x0, rhs, 1.e-6, 200, "right");

Ewritesummary

ffddmplot(E,u, "Global solution");

EVhglob def(uglob);

EfromVhi(u[],EVhglob,uglob[])

real alpha = 20000.0;
EThglob = movemesh3(EThglob, transfo = [x + alpha * uglob, y + alpha * uglobB, z + alpha * uglobC]);
u[] = mpirank;

ffddmplot(E,u, "Deformed mesh");
