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

load "Element_Mixte3d"
load "PETSc"
macro dimension()3//
include "macro_ddm.idp"

include "cube.idp"
int[int] LL = [1,2, 1,2, 1,1];
mesh3 Th = cube(getARGV("-global", 10), getARGV("-global", 10), getARGV("-global", 10), [x, y, z], label = LL); // global mesh
mesh3 ThGlobal = Th;

Mat A;
int[int] n2o;
macro ThN2O()n2o//
DmeshCreate(Th);
func Pk = [RT03d, P0];
fespace Wh(Th, Pk);
{
macro def(u)[u, u#B, u#C, u#D]//
macro init(u)[u, u, u, u]//
macro ThPostProcessD(D) {
    Wh def(u), def(v);
    varf onG(def(u), def(v)) = on(-111111, u = 10 + x - z, uB = 100 + y - x, uC = 1000 + z - y);
    v[] = D;
    u[] = onG(0, Wh);
    for [j, dj : u[]] dj = abs(dj) > 1e-2 ? 0.0 : 1.0;
    def(u) = [u, uB, uC, vD];
    D = u[];
}// EOM
MatCreate(Th, A, Pk);
}

macro div(u1, u2, u3)(dx(u1) + dy(u2) + dz(u3))//
varf vMixedLaplace([u1, u2, u3, p], [v1, v2, v3, q]) = int3d(Th)(u1 * v1 + u2 * v2 + u3 * v3 + p*div(v1, v2, v3) + div(u1, u2, u3)*q)
                                                     - int3d(Th)(q)
                                                     + int2d(Th, 1)(v1*N.x + v2*N.y + v3*N.z)
                                                     + on(2, u1 = 2.0, u2 = 3.0, u3 = 4.0);

real[int] rhs = vMixedLaplace(0, Wh, tgv = -1);
A = vMixedLaplace(Wh, Wh, tgv = -1);
set(A, sparams ="-pc_type lu");
Wh [u1, u2, u3, u4];
u1[] = A^-1 * rhs;
fespace WhGlobal(ThGlobal, Pk);
int[int] rest = restrict(Wh, WhGlobal, n2o);
WhGlobal [uG1, uG2, uG3, uG4];
real[int] reduce(uG1[].n);
{
    real[int] tmp;
    ChangeNumbering(A, u1[], tmp);
    ChangeNumbering(A, u1[], tmp, inverse = true, exchange = false);
}
for[i, v : rest] reduce[v] = u1[][i];
mpiReduce(reduce, uG1[], processor(0), mpiSUM);
if(mpirank == 0) {
    medit("Global solution", ThGlobal, uG4);
}
[u1, u2, u3, u4] = [0.0, 0.0, 0.0, 1.0];
u1[] += 1.0;
set(A, sparams = "-pc_type fieldsplit -ksp_monitor -ksp_view " +
    "-ksp_type fgmres " +
    "-pc_fieldsplit_type schur " +
    "-pc_fieldsplit_schur_fact_type full " +
    "-fieldsplit_0_ksp_type preonly " +
    "-fieldsplit_1_ksp_type gmres " +
    "-fieldsplit_1_ksp_converged_reason -fieldsplit_1_ksp_rtol 1e-4 " +
    "-pc_fieldsplit_schur_precondition selfp " +
    "-fieldsplit_1_pc_type hypre", fields = u1[]);
u1[] = A^-1 * rhs;
{
    real[int] tmp;
    ChangeNumbering(A, u1[], tmp);
    ChangeNumbering(A, u1[], tmp, inverse = true, exchange = false);
}
for[i, v : rest] reduce[v] = u1[][i];
mpiReduce(reduce, uG1[], processor(0), mpiSUM);
if(mpirank == 0) {
    medit("Global solution", ThGlobal, uG4);
}
