#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 Lagrange and Crouzeix-Raviart element
  //The result will be in a file name h_Vp_LA or h_Vp_CR
  //In this file on each ligne you have h vp_1 vp_2 vp_3 vp_4 vp_5 vp_6
  //choice(1) for Lagrange element
  //choice(0) for Crouzeix-Raviart element
  int choice(1);
  
  Int nbvp(6);   //number of eigenvalue computed
  Int nbSmax(8); //number of subdivision max
  Int nbN;       //number of _nnodes
  
  String name;
  FEType type;
  if(choice){
    name = "LA";
    type = _Lagrange;
  }
  else{
    name = "CR";
    type = _CrouzeixRaviart;
  }
  
  //Boucle on the number of subdivision
  std::ofstream myfile;
  myfile.precision(17);
  String fileName("h_Vp_"+name);
  myfile.open(fileName.c_str());
  for(int nbS=2; nbS<=nbSmax; nbS++){
    cout << "---------------------------" << endl;
    cout << "nbSubdiv = " << nbS << endl;
    nbN = pow(2,nbS)+1;
    
    //Mesh and domain
    Square sq(_origin=Point(0.,0.),_length=pi_,_nnodes=nbN,_domain_name="Omega",_side_names="Gamma");
    
    Mesh mesh(sq,_triangle,1,_structured);
    
    Domain Omega = mesh.domain("Omega");
    Domain Gamma = mesh.domain("Gamma");
    
    //Space and unknown
    Space V(Omega,type,1,"V");
    Unknown u(V, "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 << pi_*sqrt(2.)/pow(2,nbS);
    for(int k=1; k<=nbvp; k++){
        myfile << " " << eigs.value(k).real();
    }
    myfile << endl;
   }
  myfile.close();
  
  cout << "#### Fini ####" << endl;
  
  return 0;
}