// JL20180518: TO BE REVISED FOR USING Array for Boundary Conditions (ABC)
// main.cpp for Darcy3d_WG_HexaQ0Q0RT0_ProjSC (MAC version)
// The main function for Darcy3d: WG(Q0,Q0;RT[0]) on a hexahedral mesh
// Project SC: Schur complement
// 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;

  // JL20180518: TO BE REVISED FOR USING Array for Boundary Conditions (ABC)
  // Setting up ABC (Array for Boundary Conditions)
  Vector ABC(6);  // 6 boundary pieces
  ABC(1) = 1;  ABC(2) = 1;  ABC(3) = 1;
  ABC(4) = 1;  ABC(5) = 1;  ABC(6) = 1;

  // Generating a smooth hexahedral mesh on a 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);
  //
  mesh.enrich();
  std::cout << "#elements= " << mesh.numberElements() << "  ";
  std::cout << "#faces= " << mesh.numberFaces() << "\n";
  // mesh.save2file("HexaMesh.dat");
  // std::cout << "Saving mesh data finished!\n";
  
  //
  int DOFs = mesh.numberElements() + mesh.numberFaces();
  std::cout << "DOFs= " << DOFs << "\n";
  std::cout << "\n";
  
  // Sampling permeability
  // std::cout << "Sampling permeability ... ";
  Mat3 *PermK = new Mat3[mesh.numberElements()];
  Darcy3d_SmplnPerm_HexaMesh(PermK, fxnMatK, mesh, GQH);
  // std::cout << "PermK = " << PermK[mesh.numberElements()-1] << "\n";
  // std::cout << "done! \n\n";

  // The Schur-complement approach
  // Assembling and solving
  std::cout << "Schur-complement approach: Assembling and solving ... \n";
  BlockDiagMatrix GlbMatEE;
  SparseBlockMatrix GlbMatEF, GlbMatFE, GlbMatFF;
  Vector GlbVecE, GlbVecF, GlbVecD, GlbVecN;
  Darcy3d_WG_HexaQ0Q0RT0_AsmGlbMats(GlbMatEE, GlbMatEF, GlbMatFE, GlbMatFF,
                                    mesh, PermK, GQH, GQQ);
  std::cout << "\n";
  Darcy3d_WG_HexaQ0Q0RT0_AsmSource(GlbVecE, fxnf, mesh, GQH);
  Darcy3d_WG_HexaQ0Q0RT0_AsmBndryConds(GlbVecD,GlbVecN,fxnpD,fxnuN, // ABC,
                                       mesh,GQQ);
  Darcy3d_WG_HexaQ0Q0RT0_ModiLinSys(GlbMatEF, GlbMatFE, GlbMatFF,
                                    GlbVecE, GlbVecF, GlbVecD, GlbVecN, mesh);
  Vector SlnEm(mesh.numberElements());
  Vector SlnFc(mesh.numberFaces());
  int itr = 0;
  slvBDSchurCG(SlnEm, SlnFc, GlbMatEE, GlbMatEF, GlbMatFE, GlbMatFF,
               GlbVecE, GlbVecF, itr, 5000, 1e-18, 1e-18);
  std::cout << "itr= " << itr << "\n\n";

  // JL20161102: Concatenation
  Vector sln(DOFs);
  sln.setSubVector(1,SlnEm);
  sln.setSubVector(mesh.numberElements()+1,SlnFc);
  
  // Postprocessing
  // Extracting elementwise numerical pressure averages
  Vector NumerPresEm(mesh.numberElements());
  FullMatrix NmlFlux(mesh.numberElements(),6);
  FullMatrix NumerVelCofRT0(mesh.numberElements(),6);
  Darcy3d_WG_HexaQ0Q0RT0_PresVelFlux(NumerPresEm, NumerVelCofRT0, NmlFlux,
                                     mesh, PermK, sln, GQH, GQQ);
  Darcy3d_WG_HexaQ0Q0RT0_VeriNmlFluxCnty(mesh, NmlFlux);
  
  // Projecting pressure for testing
  Vector ProjPresEm(mesh.numberElements());
  Vector ProjPresFc(mesh.numberFaces());
  Darcy3d_WG_HexaQ0Q0RT0_ProjPres(ProjPresEm, ProjPresFc,
                                  fxnp, fxnpD, mesh, GQH, GQQ);

  std::cout << "Numer.Pres.Em. max = " << NumerPresEm.l0norm() << "\n";
  std::cout << " Proj.Pres.Em. max = " << ProjPresEm.l0norm()  << "\n";
  
  // Computing errors
  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=1; i<=mesh.numberElements(); ++i)  ScalarData[i-1] = NumerPresEm(i);
  Write_ScalarData_HexaMesh_silo(mesh, ScalarData, filename1a, varname1a);
  //
  // For the projected pressure
  for (int i=1; i<=mesh.numberElements(); ++i)  ScalarData[i-1] = 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
}
