// JL20170131: Need to check Max.flux.discrepancy
// Darcy3d_WG_TetraP0P0RT0_Proj2
// Main function for Darcy3d WG(P0,P0;RT0) on a tetrahedral mesh
// Project 2: 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_SmplnPerm_TetraMesh.h"
#include "Darcy3d_WG_TetraMesh.h"
#include "EqnBC_Poisson3d_Ex.h"
#include "GaussQuad.h"
#include "Interface_TetGen.h"
#include "IO_TetraMesh.h"
#include "mat3.h"
#include "TetraMesh.h"


int main(int argc, const char * argv[])
{
  char filename1[64] = "Darcy3d_WG_TetraP0P0RT0_NumerPres";
  char filename2[64] = "Darcy3d_WG_TetraP0P0RT0_ProjPres";
  char varname1[64] = "NumerPres";
  char varname2[64] = "ProjPres";
  
  clock_t TimeStart = clock();
  
  GaussQuad GQTe, GQT;
  GQTe.setForTetrahedron(5);
  GQT.setForTriangle(3);
  
  // Generating tetrahedral mesh
  TetraMesh mesh;
  // TetGen parameters passed as a string
  Gen_TetraMesh_TetGen(mesh);
  mesh.enrich3();
  std::cout << "#elements = " << mesh.numberElements() << "  ";
  std::cout << "#faces = " << mesh.numberFaces() << "\n";
  std::cout << "\n";

  // Sampling the permeability
  Mat3 *PermK = new Mat3[mesh.numberElements()];
  Darcy3d_SmplnPerm_TetraMesh(PermK, fxnMatK, mesh, GQTe);
  // std::cout << "PermK[1] = " << PermK[1] << "\n";
  
  // Assembling and solving: The single-matrix approach
  std::cout << "Assembling and solving...\n";
  int DOFs = mesh.numberElements() + mesh.numberFaces();
  std::cout << "DOFs=" << DOFs << "\n";
  SparseMatrix GlbMat;
  Vector GlbVecSource, GlbVecDirichlet, GlbVecNeumann, GlbRHS;
  Darcy3d_WG_TetraP0P0RT0_AsmGlbMat1(GlbMat, mesh, PermK, GQTe);
  Darcy3d_WG_TetraP0P0RT0_AsmSource1(GlbVecSource, fxnf, mesh, GQTe);
  Darcy3d_WG_TetraP0P0RT0_AsmBndryConds1(GlbVecDirichlet, GlbVecNeumann,
                                         fxnpD, fxnuN, mesh, GQT);
  Darcy3d_WG_TetraP0P0RT0_ModiLinSys1(GlbMat, GlbRHS, GlbVecSource,
                                      GlbVecDirichlet, GlbVecNeumann, mesh);
  int itr = 0,  MaxItr = 5000;
  double thres = 1e-18,  tol = 1e-18;
  Vector sln(DOFs);  // The size (DOFs) must be specified
  slvSpaSpdSysCG(sln, GlbMat, GlbRHS, itr, MaxItr, thres, tol);
  std::cout << "itr=" << itr << "\n";
  std::cout << "\n";
  
  // Post-processing for numerical pressure, velocity, flux
  Vector NumerPresEm;
  FullMatrix NumerVelCofRT0, NumerFlux;
  Darcy3d_WG_TetraP0P0RT0_PresVelFlux(NumerPresEm, NumerVelCofRT0, NumerFlux,
                                      mesh, PermK, sln, GQTe, GQT);
  
  // Projecting pressure for testing
  Vector ProjPresEm, ProjPresFc;
  Darcy3d_WG_TetraP0P0RT0_ProjPres(ProjPresEm, ProjPresFc, fxnp, fxnpD,
                                   mesh, GQTe, GQT);
  
  std::cout << "NumerPresEm max = " << NumerPresEm.l0norm() << "\n";
  std::cout << "ProjPresEm  max = " << ProjPresEm.l0norm() << "\n";
  // NumerPresEm.save2file("NumerPresEm.dat");
  // ProjPresEm.save2file("ProjPresEm.dat");
  
  // Computing errors
  Vector ErrPresEm = NumerPresEm - ProjPresEm;
  std::cout << "Error max in pres. = " << ErrPresEm.l0norm() << "\n";
  //
  double L2ErrPres, L2ErrVel;
  Darcy3d_WG_TetraP0P0RT0_Err(L2ErrPres, L2ErrVel, mesh,
                              NumerPresEm, NumerVelCofRT0, fxnp, fxnu);
  std::cout << "L2ErrPres= " << L2ErrPres << "  ";
  std::cout << "L2ErrVel= " << L2ErrVel << "\n";
  std::cout << "\n";

  // Output & saving
  double *ScalarData = new double[mesh.numberElements()];
  // Numerical pressure
  for (int i=1; i<=mesh.numberElements(); ++i)  ScalarData[i-1] = NumerPresEm(i);
  Write_ScalarData_TetraMesh_silo(mesh, ScalarData, filename1, varname1);
  // Projected pressure
  for (int i=1; i<=mesh.numberElements(); ++i)  ScalarData[i-1] = ProjPresEm(i);
  Write_ScalarData_TetraMesh_silo(mesh, ScalarData, filename2, varname2);
  
  // Freeing...
  delete[] ScalarData;
  delete[] PermK;
  
  // Finishing
  clock_t TimeEnd = clock();
  std::cout << "Time taken = "
    << static_cast<double>(TimeEnd-TimeStart)/CLOCKS_PER_SEC << " seconds.\n";

  return(0);  // If successful
}
