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

load "PETSc-complex"
include "macro_ddm.idp"             // additional DDM functions

assert(mpisize == 1);
macro grad(u)[dx(u), dy(u)]// EOM   // two-dimensional gradient
func Pk = P1;                       // finite element space
func real wedge(real a, real b) {
    if(y < 0.4 + 0.1 * 0.75 * x)
        return 2;
    else if(y < 0.8 - 0.2 * 0.75 * x)
        return 1.5;
    else
        return 3;
}
real omega = 2 * pi * 5;
func f = 80 * 100 * exp(-20 * 100 * ((x-0.5)^2 + (y-0.25)^2));

mesh Th = square(getARGV("-global", 30), getARGV("-global", 30));
fespace Ph(Th, P0);
Ph val = wedge(x, y);
Ph k = omega / val;
varf vPb(u, v) = int2d(Th)(-k^2 * u * v + (grad(u)' * grad(v))) + int1d(Th, 2)(1i * k * u * v) + on(1, u = 0.0);

fespace Vh(Th, Pk);
Mat<complex> A(Vh.ndof);
A = vPb(Vh, Vh, tgv = -2);
int n = 10;
Vh<complex>[int]  sol(10);
complex[int] eigenvalues(10);
Vh<complex>[int] lsol(10);
Vh<complex>[int] rsol(10);
real[int] singularvalues(10);
int nconv = EPSSolve(A, values = eigenvalues, vectors = sol, sparams = "-eps_target 0 -eps_target_magnitude -eps_converged_reason -st_type sinvert -eps_nev " + n);
eigenvalues.resize(min(nconv, n));
sol.resize(min(nconv, n));
nconv = SVDSolve(A, values = singularvalues, lvectors = lsol, rvectors = rsol, sparams = "-svd_smallest -svd_converged_reason -svd_type lapack");
singularvalues.resize(min(nconv, n));
lsol.resize(min(nconv, n));
rsol.resize(min(nconv, n));
cout << eigenvalues << endl;
cout << singularvalues << endl;
