#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 matMass_LA or matMass_CR
  //In this file on each ligne you have i j m_{i,j}
  //choice(1) for Lagrange element
  //choice(0) for Crouzeix-Raviart element
  int choice(1);
  
  String name;
  FEType type;
  if(choice){
    name = "LA";
    type = _Lagrange;
  }
  else{
    name = "CR";
    type = _CrouzeixRaviart;
  }
  
  //Mesh and domain
  Square sq(_origin=Point(0.,0.),_length=pi_,_nnodes=5,_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");
  
  B.saveToFile("matMass_"+name,_coo);
  
  cout << "#### Fini ####" << endl;
  
  return 0;
}