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

load "PETSc"
include "macro_ddm.idp"

real c = 6.25;
int N  = 80;

// domain: unit square
border aa(t=0,1) { x=t;   y=0;   };
border bb(t=0,1) { x=1;   y=t;   };
border cc(t=0,1) { x=1-t; y=1;   };
border dd(t=0,1) { x=0;   y=1-t; };

mesh M = buildmesh(aa(N)+bb(N)+cc(N)+dd(N));

load "Element_P3"
func Pk = P3;
Mat A;
macro MRefinementFactor()getARGV("-split", 1)//
MatCreate(M, A, Pk);
fespace Vh(M, Pk);
Vh u;
func BC = cos(pi*x)*cos(pi*y);
varf vInit(w, v) = on(aa, bb, cc, dd, w = BC);
varf vJ(w, v) = int2d(M)(dx(w)*dx(v) + dy(w)*dy(v) - c*exp(u)*w*v) + on(aa, bb, cc, dd, w = 0);
varf vRes(w, v) = int2d(M)(dx(u)*dx(v) + dy(u)*dy(v) - c*exp(u)*v) + on(aa, bb, cc, dd, w = u);
set(A, sparams = "-pc_type asm -sub_pc_type lu -ksp_converged_reason");
func real[int] funcRes(real[int]& inPETSc) {
    ChangeNumbering(A, u[], inPETSc, inverse = true, exchange = true);
    real[int] out(Vh.ndof);
    out = vRes(0, Vh, tgv = -2);
    real[int] outPETSc;
    ChangeNumbering(A, out, outPETSc);
    return outPETSc;
}
func int funcJ(real[int]& inPETSc) {
    ChangeNumbering(A, u[], inPETSc, inverse = true, exchange = true);
    A = vJ(Vh, Vh, tgv = -2);
    return 0;
}
real[int] bPETSc;
u[] = vInit(0, Vh, tgv = -2);
ChangeNumbering(A, u[], bPETSc);
real[int] xPETSc = bPETSc;
SNESSolve(A, funcJ, funcRes, bPETSc, xPETSc, sparams = "-snes_monitor -snes_linesearch_monitor -snes_type newtontr -snes_converged_reason -snes_view");
ChangeNumbering(A, u[], xPETSc, inverse = true, exchange = true);
macro def(u)u//
plotMPI(M, u, Pk, def, real, cmm = "Global solution");
