// This codes shows how to solve an eigenvalue problems
// in different coordinate systems. The example here is
// the Laplace equation on a toroidal surface embedded
// in R3. The cartesian coordinates on the torus read
//     x = (r1+r2*cos(u))*cos(v)
//     y = (r1+r2*cos(u))*sin(v)
//     z = r2*sin(v)
// with r1>r2 the major and minor radii and u,v \in [0,2\pi]
// The toroidal surface thus maps to a square domain of
// side 2\pi and periodic boundary conditions in both
// coordinates.
//
// Laplace-Beltrami on a torus, using SLEPc.
//
// Is the eigenspectrum known analytically?
//
// Usage :
// mpirun -np 4 FreeFem++-mpi -wg laplace-2d-torus-SLEPc.edp \
//        -split 1 -npts 400 -nev 15 -sigma 0.0
//
// Authors: Julien Garaud <julien.garaud@gmail.com>
//          Pierre Jolivet <pierre.jolivet@enseeiht.fr>

/***************************************/
/*         Geometry parameters         */
/***************************************/
int[int] Labels = [1,2,3,4]; // labels : bottom, right, top, left sides
int[int] labPeriodic = [Labels[0],Labels[2],Labels[1],Labels[3]];

/**************************************/
/*     Load PETSc & SLEPc macros      */
/**************************************/
load "PETSc"                             // PETSc plugin
macro dimension(      )2// End Of Macro  // 2D or 3D
include "macro_ddm.idp"                  // Additional DDM functions
macro def(i)i// EOM
macro init(i)i// EOM

macro Pk() P1, periodic=[[Labels[0],x],[Labels[2],x],[Labels[1],y],[Labels[3],y]]//EOM

/***************************************/
/*    Options for distributed solver   */
/***************************************/
int s      = getARGV("-split", 1)   ; // Refinement factor
//
int Npts   = getARGV("-npts" , 400) ; // Number of points on the perimeter
//
int nEV    = getARGV("-nev"  ,  5)  ; // Number of eigenvalues
real sigma = getARGV("-sigma", 0.0) ; // Shift
//
real r1 = getARGV("-r1",2.0); // Large radius of the torus
real r2 = getARGV("-r2",1.0); // Small radius of the torus
if (r1<=r2){if(!mpirank) cout << "Error r1<r2 is a wrong dimensions of the torus"<<endl; exit(1);}
/***************************************/
/*     Verbosity and passed options    */
/***************************************/
if(verbosity > 0 && mpirank == 0) {
    cout << "********************************************"  << endl
        << " --- " << mpirank << "/" << mpisize
        << "- laplace-2d-torus-SLEPc.edp "        << endl
        << "********************************************"  << endl
        << "- input parameters: "                          << endl
        << "  refinement factor  = " << s                  << endl
        << "********************************************"  << endl
        << "  nb of pts on perimeter = " << Npts           << endl
        << "********************************************"  << endl
        << "  nb of eigenvalues  = " << nEV                << endl
        << "  value of the shift = " << sigma              << endl
        << "********************************************"  << endl
        << "  Large radius of the torus = " << r1          << endl
        << "  Small radius of the torus = " << r2          << endl
        << "********************************************"  << endl
        << endl;
}
/***************************************/
/*   ###############################   */
/***************************************/
mesh Th = square(1, 1);      // Local mesh
int[int][int] intersection;   // Local-to-neighbors renumbering
real[int] D;                  // Partition of unity

/***************************************/
/*         Finite Element space        */
/***************************************/
// Definition of the finite element space on the domain Th
// P1 are the first order Lagrange elements
fespace Vh(Th, Pk);      // local finite element space

/***************************************/
/*            ##############           */
/***************************************/
{ // Construction of the rectangular domain
    int Thetapts = int(0.25*Npts); // pts on the x-axis sides
    int Phipts   = int(0.25*Npts); // pts on the y-axis sides

    Th = square(Thetapts,Phipts,[2.0*x*pi,2.0*pi*y],label=Labels);
    // .....
    buildPeriodic(Th,       // The local mesh
            s,        // Refinement factor
            intersection, // local-to-neighbors renumbering
            D,        // partition of unity
            Pk,           // FE-space
            mpiCommWorld, // Communicator
            labPeriodic   // Array of labels for periodic boundaries
         );
}
/***********************************************************************/
/*         Coordinate dependant differential operators                 */
/***********************************************************************/
/*              Parametric coordinates in the (u,v)-plane              */
/*                                                                     */
/*                  u  --> x  in [0,2*pi]                              */
/*                  v  --> y  in [0,2*pi]                              */
/*                                                                     */
/*          Jacobian determinant on the bi-periodic square             */
/*        The det(J) = r2*(r1+r2*cos(u)) --> r2(r1+r2*cos(x))          */
/*                                                                     */
macro Jac()( r2*(r1+r2*cos(x)) )            // End Of Macro /*         */
     /*  The Jacobian                                                  */
/*                                                                     */
/*       The gradiant operator on the induced toroidal surface         */
/*                                                                     */
/*       grad =   1/r2*d/du              ->  1/r2*d/dx                 */
/*                1/(r1+r2*cos(u))*d/dv  ->  1/(r1+r2*cos(x)*d/dy      */
/*                                                                     */
macro Grad(u) [dx(u)/r2,dy(u)/(r1+r2*cos(x))]      // End Of Macro /*  */
     /*  The Gradient operator                                         */
macro Lap(u,v) ( Grad(u)'*Grad(v)) //')      // End Of Macro /*        */
     /*  The Laplace-Beltrami operator on the induced surface          */ 
/*                                                                     */
/***********************************************************************/
/***************************************/
/*          Problem parameters         */
/***************************************/

/***************************************/
/*          Problem definition         */
/***************************************/
varf vA(uh,vh)= int2d(Th)        //  Definion of  the problem
  (Jac*(Lap(uh,vh)))//  Bilinear form
  ;

varf vB(uh,vh)= int2d(Th)  //  Definion of  the problem
  ( Jac*uh*vh )           //  Bilinear form
  ;

matrix<real> A = vA(Vh,Vh);
matrix<real> B = vB(Vh,Vh);
/***************************************/
/*      Build distributed matrices     */
/***************************************/
Mat DistA(A, intersection, D);
Mat DistB(DistA, B);
/***************************************/
/*          Problem resolution         */
/***************************************/
real[int]     EigenVAL(0);       // array to store eigenvalues
Vh<real>[int] def(EigenVEC)(1); // array to store eigenvectors

string ssparams =            // Parameters for the distributed EigenValue solver
  " -eps_nev " + nEV       + // Number of eigenvalues
  " -eps_type krylovschur" +
  " -eps_target "+ sigma   + // Shift value
  " -st_type sinvert "     +
  " -st_pc_type cholesky "       +
  " -st_matstructure same "       +
  " -eps_view"             +
  " -eps_gen_hermitian"      // The problem is symmetric
  ;

int k = EPSSolve
(DistA,              // matrix OP = A − sigma*B
 DistB,              //
 vectors = EigenVEC, // Array to store the FEM-EigenFunctions
 values  = EigenVAL, // Array to store the EigenValues
 sparams = ssparams  // Parameters for the distributed EigenValue solver
);

k=min(k,nEV); // some time the number of converged eigen value
              // can be greater than nev;
/***************************************/
/*           View the solution         */
/***************************************/
Vh<real> Temp;

for(int i=0;i<k;i++){
    if(!mpirank) cout << " Eigenvalue #"+i+" = "+EigenVAL[i]<<endl;
    Temp = EigenVEC[i];
    plotMPI(Th,     // The local mesh
            Temp,   // The local solution
            Pk,     // Local FE-space
            def,    // Macro for field definition
            real,   // Type
            cmm = "Psi("+i+")  EV = "+EigenVAL[i]
           )
}
