// main.cpp
// Darcy3d_CG_TetraP1_Proj2
// The main function for Darcy3d CG P1 on a tetrahedral mesh
// Project 2: A tetrahedral mesh generated by TetGen
// James Liu, Graham Harper, ColoState; 2014/07--2018/05

#include <ctime>
#include <iostream>

#include "LinSys.h"
#include "matrix.h"
#include "vector.h"

#include "Darcy3d_CG_TetraMesh.h"
#include "Darcy3d_SmplnPerm_TetraMesh.h"
#include "EqnBC_Poisson3d_Ex.h"
#include "GaussQuad.h"
#include "Interface_TetGen.h"
#include "IO_TetraMesh.h"
#include "TetraMesh.h"
#include "TetraMesh_ConvCP1toP0.h"


int main(int argc, const char * argv[])
{
  char filename1[64] = "Darcy3d_CG_TetraP1_NumerPres";
  char filename2[64] = "Darcy3d_CG_TetraP1_ProjPres";
  char varname1[64] = "NumerPres";
  char varname2[64] = "ProjPres";
  
  clock_t TimeStart = clock();
  
  GaussQuad GQTe, GQT;
  GQTe.setForTetrahedron(5);
  GQT.setForTriangle(3);
  
  // Generating a tetrahedral mesh by TetGen
  TetraMesh mesh;
  // TetGen parameters passed as a string
  Gen_TetraMesh_TetGen(mesh);
  std::cout << "Enriching the tetrahedral mesh ... " << "\n";
  mesh.enrich1();
  std::cout << "#elements= " << mesh.numberElements() << "  ";
  std::cout << "#nodes= " << mesh.numberNodes() << "\n\n";
  // mesh.save2file("TetraMesh.dat");

  // Sampling permeability
  Mat3 *PermK = new Mat3[mesh.numberElements()];
  Darcy3d_SmplnPerm_TetraMesh(PermK, fxnMatK, mesh, GQTe);
  std::cout << "PermK[0] = " << PermK[0] << "\n";

  // Assembling and solving
  // std::cout << "Assembling and solving...\n";
  int DOFs = mesh.numberNodes();
  std::cout << "DOFs = " << DOFs << "\n";
  SparseMatrix GlbMat;
  Vector DirichletVec(DOFs), NeumannVec(DOFs), GlbRHS(DOFs), sln(DOFs);
  Darcy3d_CG_TetraP1_AsmGlbMat(GlbMat, mesh, PermK);
  Darcy3d_CG_TetraP1_AsmGlbRHS(GlbRHS, fxnf, mesh, GQTe);
  Darcy3d_CG_TetraP1_AsmBndryConds(DirichletVec, NeumannVec,
                                   fxnpD, fxnuN, mesh, GQT);
  Darcy3d_CG_TetraP1_AdjLinSys(GlbMat, GlbRHS, DirichletVec, NeumannVec, mesh);
  int itr = 0;
  slvSpaSpdSysCG(sln, GlbMat, GlbRHS, itr, 5000, 1e-18, 1e-18);
  std::cout << "itr=" << itr << "\n";
  // GlbMat.save2file("GlbMat.dat");
  // GlbRHS.save2file("GlbRHS.dat");
  // DirichletVec.save2file("DirichletVec.dat");
  // sln.save2file("sln.dat");
  
  // Extracting info
  Vector NumerPresEmAvg(mesh.numberElements());
  TetraMesh_ConvCP1toP0(NumerPresEmAvg, sln, mesh);
  std::cout << "NumerPresEmAvg max = " << NumerPresEmAvg.l0norm() << "\n";
  // NumerPresEmAvg.save2file("NumerPresEmAvg.dat");
  
  // Projecting pressure for testing
  Vector IntpltPresNd(mesh.numberNodes());
  Vector ProjPresEmAvg(mesh.numberElements());
  Darcy3d_CG_TetraP1_IntpltPres(IntpltPresNd, fxnp, mesh);
  TetraMesh_ConvCP1toP0(ProjPresEmAvg, IntpltPresNd, mesh);
  std::cout << "ProjPresEmAvg max = " << ProjPresEmAvg.l0norm() << "\n";
  // ProjPresEmAvg.save2file("ProjPresEmAvg.dat");
  
  // Errors
  Vector ErrPresEm = NumerPresEmAvg - ProjPresEmAvg;
  std::cout << "Pressure elementwise error max = " << ErrPresEm.l0norm() << "\n";
  // ErrPresEm.save2file("ErrPresEm.dat");
  
  // Preparing for data files
  double *ScalarData = new double[mesh.numberElements()];
  //
  // Numerical pressure
  for (int i=1; i<=mesh.numberElements(); ++i)  ScalarData[i-1] = NumerPresEmAvg(i);
  Write_ScalarData_TetraMesh_silo(mesh, ScalarData, filename1, varname1);
  //
  // Projected pressure
  for (int i=1; i<=mesh.numberElements(); ++i)  ScalarData[i-1] = ProjPresEmAvg(i);
  Write_ScalarData_TetraMesh_silo(mesh, ScalarData, filename2, varname2);
  
  // Freeing...
  delete[] ScalarData;
  delete[] PermK;
  
  clock_t TimeEnd = clock();
  std::cout << "\n" << "Time taken = "
    << static_cast<double>(TimeEnd-TimeStart)/CLOCKS_PER_SEC << " seconds.\n";
  
  return(0);  // If successful
}
