// main.cpp for Darcy3d_WG_HexaQ0Q0RT0_Proj1 (MAC version)
// The main function for Darcy3d: WG(Q0,Q0;RT[0]) on a hexahedral mesh
// The single-matrix approach
// A hexahedral mesh is generated by perturbing a brick mesh
// 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_HexaMesh.h"
#include "Darcy3d_WG_HexaMesh.h"
#include "EqnBC_Poisson3d_Ex.h"
#include "GaussQuad.h"
#include "HexaMesh.h"
#include "IO_HexaMesh.h"
#include "mat3.h"


int main(int argc, const char * argv[])
{
  char filename1a[64] = "Darcy3d_WG_HexaQ0Q0RT0_NumerPres";
  char filename2a[64] = "Darcy3d_WG_HexaQ0Q0RT0_ProjPres";
  char varname1a[64] = "NumerPres";
  char varname2a[64] = "ProjPres";
  
  clock_t TimeStart = clock();
  
  GaussQuad GQH, GQQ;
  GQH.setForBrick(3,3,3);
  GQQ.setForRectangle(3,3);

  // Setting up a 3-dim brick domain: The unit cube
  double  xa = 0,  xb = 1;
  double  yc = 0,  yd = 1;
  double  ze = 0,  zf = 1;

  // Generating a logically rectangular hexahedral mesh on the brick domain
  int n = 16;
  int nx = n,  ny = n,  nz = n;
  double delta = 1;  // Perturbation magnitude
  HexaMesh mesh(xa, xb, yc, yd, ze, zf, nx, ny, nz, delta);
/*
  int m = 3,  nz = 8;
  HexaMesh mesh(xa, xb, yc, yd, ze, zf, m, nz);
*/
  mesh.enrich();
  std::cout << "#elements= " << mesh.numberElements() << "  ";
  std::cout << "#faces= " << mesh.numberFaces() << "\n";
  // mesh.save2file("HexaMesh.dat");
  // std::cout << "Saving mesh data finished!\n";
  std::cout << "\n";

  // Setting up a 3-dim domain: A cylindrical sector domain
  /*
   double  ra = 1,  rb = 2;
   double  tc = 0,  td = PI/2.0;
   double  zb = 0,  zt = 1;
   // Generating a logically rectangular hexahedral mesh for the cylindrical sector domain
   int n = 32;
   int nr = n,  nt = n,  nz = n;
   HexaMesh mesh(ra, rb, tc, td, zb, zt, nr, nt, nz);
   mesh.enrich();
   std::cout << "#elements= " << mesh.numberElements() << "  ";
   std::cout << "#faces= " << mesh.numberFaces() << "\n";
   // mesh.save2file("HexaMesh.dat");
  */

  // Sampling permeability
  std::cout << "Sampling permeability...\n";
  Mat3 *PermK = new Mat3[mesh.numberElements()];
  Darcy3d_SmplnPerm_HexaMesh(PermK, fxnMatK, mesh, GQH);
  // std::cout << "PermK = " << PermK[mesh.numberElements()-1] << "\n";
  std::cout << "\n";
  
  // JL20161102: The single-matrix approach
  // Assembling and solving
  std::cout << "Assembling and solving...\n";
  int DOFs = mesh.numberElements() + mesh.numberFaces();
  std::cout << "DOFs= " << DOFs << "\n";
  SparseMatrix GlbMat;
  Vector GlbVecSource, GlbVecDirichlet, GlbVecNeumann, GlbRHS, sln(DOFs);
  Darcy3d_WG_HexaQ0Q0RT0_AsmGlbMat1(GlbMat, mesh, PermK, GQH, GQQ);
  Darcy3d_WG_HexaQ0Q0RT0_AsmSource1(GlbVecSource, fxnf, mesh, GQH);
  Darcy3d_WG_HexaQ0Q0RT0_AsmBndryConds1(GlbVecDirichlet, GlbVecNeumann,
                                        fxnpD, fxnuN, mesh, GQQ);
  Darcy3d_WG_HexaQ0Q0RT0_ModiLinSys1(GlbMat, GlbRHS, GlbVecSource,
                                     GlbVecDirichlet, GlbVecNeumann, mesh);
  int itr = 0;
  slvSpaSpdSysCG(sln, GlbMat, GlbRHS, itr, 5000, 1e-18, 1e-18);
  std::cout << "itr=" << itr << "\n";
  std::cout << "\n";
  
  // Postprocessing
  // Computing numerical pressure, velocity, flux
  Vector NumerPresEm(mesh.numberElements());
  FullMatrix NumerVelCofRT0(mesh.numberElements(),6);
  FullMatrix NumerFlux(mesh.numberElements(),6);
  Darcy3d_WG_HexaQ0Q0RT0_PresVelFlux(NumerPresEm, NumerVelCofRT0, NumerFlux,
                                     mesh, PermK, sln, GQH, GQQ);
  
  // Projecting pressure for testing
  Vector ProjPresEm(mesh.numberElements());
  Vector ProjPresFc(mesh.numberFaces());
  Darcy3d_WG_HexaQ0Q0RT0_ProjPres(ProjPresEm, ProjPresFc,
                                  fxnp, fxnpD, mesh, GQH, GQQ);

  //
  Darcy3d_WG_HexaQ0Q0RT0_VeriNmlFluxCnty(mesh, NumerFlux);
  std::cout << "NumerPresEm max = " << NumerPresEm.l0norm() << "\n";
  std::cout << "ProjPresEm  max = " << ProjPresEm.l0norm()  << "\n";

  // Computing errors
  Vector ErrPresEm = NumerPresEm - ProjPresEm;
  std::cout << "Err. press. max = " << ErrPresEm.l0norm() << "\n";
  //
  double L2ErrPres, L2ErrVel, L2ErrFlux;
  Darcy3d_WG_HexaQ0Q0RT0_Err(L2ErrPres, L2ErrVel, L2ErrFlux, mesh,
                             NumerPresEm, NumerVelCofRT0, fxnp, fxnu);
  std::cout << "L2ErrPres= " << L2ErrPres << "  ";
  std::cout << "L2ErrVel= " << L2ErrVel << "  ";
  std::cout << "L2ErrFlux= " << L2ErrFlux << "\n";
  std::cout << "\n";

  // Output & saving
  double *ScalarData = new double[mesh.numberElements()];
  // For the numerical pressure
  for (int i=0; i<mesh.numberElements(); ++i)  ScalarData[i] = NumerPresEm[i];
  Write_ScalarData_HexaMesh_silo(mesh, ScalarData, filename1a, varname1a);
  // For the projected pressure
  for (int i=0; i<mesh.numberElements(); ++i)  ScalarData[i] = ProjPresEm[i];
  Write_ScalarData_HexaMesh_silo(mesh, ScalarData, filename2a, varname2a);
  
  // 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
}
