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

load "mshmet"
load "mmg"
load "medit"
load "PETSc"
macro dimension()3// EOM            // 2D or 3D
include "macro_ddm.idp"             // additional DDM functions

int n = 12;
int[int] L=[1,1,1,1], LU=[0,2], LD=[1,1], LR=[0,1];
mesh3 Th3 = buildlayers(square(n, n, region=0, label=L), n, zbound=[0,1], labelmid=LD, labelup=LU, labeldown=LR);
Th3 = trunc(Th3, x < 0.5 || y < 0.5 || z < 0.5, label=3);

macro grad(u)[dx(u), dy(u), dz(u)]// EOM    // three-dimensional gradient

varf vPoisson(u, v) = int3d(Th3)(grad(u)' * grad(v)) + int3d(Th3)(v) + on(1,2,3, u=0);

mesh3 ThBackup = Th3;
Mat A;
int[int] n2o;
NewMacro Th3N2O() n2o EndMacro
DmeshCreate(Th3);
MatCreate(Th3, A, P1);
fespace Vh(Th3, P1);
Vh u;
real err = 1.0e-2;

int iMax = 4;
int[int] rt(2);
rt = [1,2];
for(int i = 0; i < iMax; ++i) {
    A = vPoisson(Vh, Vh, tgv = -2);
    set(A, sparams = "-ksp_monitor -pc_type hypre");
    real[int] rhs = vPoisson(0, Vh, tgv = -1);
    u[] = A^-1 * rhs;
    int[int] fforder = [1];
    savevtk("laplace-adapt-3d.vtu", Th3, u, bin = 1, order = fforder, append = i ? true : false);
    if(i < iMax - 1) {
        fespace VhBackup(ThBackup, P1);
        VhBackup h, hReduced;
        u[] .*= A.D;
        int[int] rest = restrict(Vh, VhBackup, n2o);
        for[i, v : rest] h[][v] = u[][i];
        mpiAllReduce(h[], hReduced[], mpiCommWorld, mpiSUM);
        if(mpirank == 0) {
            real[int] met = mshmet(ThBackup, hReduced, aniso = 0, hmin = 1.0e-3, hmax = 2.0e-1, err = err);
            plot(hReduced);
            Th3 = mmg3d(ThBackup, metric = met, hgrad = 3.0, requiredTriangle = rt);
            medit("U-adapt-iso-" + i, Th3, hReduced);
        }
        broadcast(processor(0), Th3);
        ThBackup = Th3;
        Mat Adapt;
        MatCreate(Th3, Adapt, P1);
        A = Adapt;
        u = 0.0;
        err *= 0.5;
    }
}
