#include "xlife++.h"
#include <iostream>
#include <fstream>
using namespace xlifepp;
using namespace std;

int main(int argc, char** argv)
{
  init(_lang=en); //mandatory initialization of xlifepp
  
  //Change the value of choice to choose between uniform and Gauss-Lobatto interpolation
  //The result will be in a file name degVp_UI or deg_Vp_GL
  //In this file on each ligne you have degree vp_1 vp_2 vp_3 vp_4 vp_5 vp_6
  //choice(1) for uniform interpolation
  //choice(0) for Gauss-Lobatto interpolation
  int choice(1);
  
  Int nbvp(6);  //number of eigenvalue computed
  Int dMin(4);  //degree min
  Int dMax(20); //degree max
  
  String name;
  FESubType type;
  if(choice){
    name = "UI";
    type = _standard;
  }
  else{
    name = "GL";
    type = _GaussLobatto;
  }
  
  //Mesh and domain
  Square sq(_origin=Point(0.,0.),_length=pi_,_nnodes=2,_domain_name="Omega",_side_names="Gamma");
  
  Mesh mesh(sq,_quadrangle,1,_structured);
  
  Domain Omega = mesh.domain("Omega");
  Domain Gamma = mesh.domain("Gamma");
  
  //Boucle on the degree
  std::ofstream myfile;
  myfile.precision(17);
  String fileName("degVp_"+name);
  myfile.open(fileName.c_str());
  for(int d=dMin; d<=dMax; d++){
    cout << "------------------------------" << endl;
    cout << "degree = " << d << endl;
    
    //Space and unknown
    Space Vd(Omega,interpolation(Lagrange,type,d,H1),"Vd");
    Unknown u(Vd, "u");
    TestFunction v(u, "v");
    
    //Bilinear forme
    BilinearForm auv = intg(Omega,grad(u)|grad(v));
    BilinearForm buv = intg(Omega,u*v);
    EssentialConditions ecs = ((u|Gamma) = 0);
    
    //TermMatrix and Solver
    cpuTime();
    TermMatrix A(auv, ecs, ReductionMethod(_pseudoReduction,10.), "A");
    TermMatrix B(buv, ecs, ReductionMethod(_pseudoReduction,0.1), "B");
    cpuTime("Time TermMatrix");
    
    EigenElements eigs = eigenSolve(A,B,_nev=nbvp,_sigma=1.);
    cpuTime("Time eigenSolve");
    
    myfile << d;
    for(int k=1; k<=nbvp; k++){
        myfile << " " << eigs.value(k).real();
    }
    myfile << endl;
   }
  myfile.close();
  
  cout << ">>>> Fini <<<<" << endl;
  
  return 0;
}