// main.cpp (MAC version)
// LinElas3d_WG_HexaQ03Q03RT03Q0_ProjSC
// The main function for LinElas3d: WG(Q0^3,Q0^3;RT[0]^3,Q0) on a hexahedral mesh
// The Schur complement approach
// A hexahedral mesh generated by itself
// James Liu, Graham Harper, ColoState; 2014/07--2018/05

#include <ctime>
#include <iostream>

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

#include "EqnBC_LinElas3d_Ex.h"
#include "GaussQuad.h"
#include "HexaMesh.h"
#include "HexaMesh_IntpltProj.h"
#include "IO_HexaMesh.h"
#include "LinElas3d_WG_HexaMesh.h"


int main(int argc, const char * argv[])
{
  char filename1[64] = "LinElas3d_WG_HexaQ03Q03RT03Q0_NumerDiv";
  char filename3[64] = "LinElas3d_WG_HexaQ03Q03RT03Q0_NumerDspl";
  char varname1[64] = "NumerDiv";
  char varname3x[64] = "NumerDspl_x";
  char varname3y[64] = "NumerDspl_y";
  char varname3z[64] = "NumerDspl_z";

  clock_t TimeStart = clock();
  
  double lambda=1, mu=1;
  
  GaussQuad GQH, GQQ;
  GQH.setForBrick(3,3,3);
  GQQ.setForRectangle(3,3);
  
  // Setting up a 3-dim domain: The unit cube
  double  xa = 0,  xb = 1;
  double  yc = 0,  yd = 1;
  double  ze = 0,  zf = 1;
  
  // Setting up ABC (Array for Boundary Conditions):
  // 1 Dirichlet, 2 Neumann, 3 Robin, 4,5,6: Symmetry in x,y,z respectively
  Vector ABC(6);
  // For EqnBC_LinElas3d_Ex11sincoscos
  ABC(1) = 1;  ABC(2) = 1;  ABC(3) = 1;
  ABC(4) = 1;  ABC(5) = 1;  ABC(6) = 1;

  // Generating a logically 3-dim rectangular hexahedral mesh on the unit cube
  int n = 16;
  int nx = n,  ny = n,  nz = n;
  double delta = 0;  // 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";

  // The Schur-complement approach: Assembling and Solving
  std::cout << "Assembling and Solving...\n";
  int DOFs = 3*mesh.numberElements() + 3*mesh.numberFaces();
  std::cout << "DOFs= " << DOFs << "\n";
  // std::cout << "DOFs_SC = " << 3*mesh.numberFaces() << "\n";
  Vector GlbVecE, GlbVecD, GlbVecN, GlbVecC, flag;
  BlockDiagMatrix GlbMatEE;
  SparseBlockMatrix GlbMatEC, GlbMatCE, GlbMatCC;
  LinElas3d_WG_HexaQ03Q03RT03Q0_AsmGlbMats(GlbMatEE, GlbMatEC, GlbMatCE, GlbMatCC,
                                           lambda, mu, mesh, GQH, GQQ);
  std::cout << "\n";
  LinElas3d_WG_HexaQ03Q03RT03Q0_AsmSource(GlbVecE, fxnf, mesh, GQH);
  LinElas3d_WG_HexaQ03Q03RT03Q0_AsmBndryConds(GlbVecD,GlbVecN,flag,
                                              fxnuD,fxntN,ABC,mesh,GQQ);
  LinElas3d_WG_HexaQ03Q03RT03Q0_ModiLinSys(GlbMatEC, GlbMatCE, GlbMatCC,
                                           GlbVecE, GlbVecC, mesh, GlbVecD, GlbVecN);
  Vector SlnEm(3*mesh.numberElements()), SlnFc(3*mesh.numberFaces());
  int itr=0;
  slvBDSchurCG(SlnEm, SlnFc, GlbMatEE, GlbMatEC, GlbMatCE, GlbMatCC,
               GlbVecE, GlbVecC, itr, 5000, 1e-18, 1e-18);
  std::cout << "itr= " << itr << "\n\n";

  // Concatenation
  Vector sln(DOFs);
  sln.setSubVector(1,SlnEm);
  sln.setSubVector(3*mesh.numberElements()+1,SlnFc);
  
  // Postprocessing
  // Computing numerical displacement, stress, dilation
  Vector NumerDsplEm(3*mesh.numberElements());
  Vector NumerDsplCofQ0x(mesh.numberElements());
  Vector NumerDsplCofQ0y(mesh.numberElements());
  Vector NumerDsplCofQ0z(mesh.numberElements());
  Vector NumerDivCofQ0(mesh.numberElements());
  FullMatrix NumerDsplCofQ03(mesh.numberElements(),3);
  FullMatrix NumerStrsCofSym(mesh.numberElements(),18);
  LinElas3d_WG_HexaQ03Q03RT03Q0_DsplStrsDiv(NumerDsplCofQ03, NumerStrsCofSym,
                                            NumerDivCofQ0,lambda,mu,mesh,sln,GQH,GQQ);
  for (int labele=1; labele<=mesh.numberElements(); labele++) {
    NumerDsplCofQ0x(labele) = NumerDsplCofQ03(labele,1);
    NumerDsplCofQ0y(labele) = NumerDsplCofQ03(labele,2);
    NumerDsplCofQ0z(labele) = NumerDsplCofQ03(labele,3);
  }
  
  // Projecting displacement and its divergence for testing
  PtVec3d *ProjDsplEm = new PtVec3d[mesh.numberElements()];
  Vector ProjDivEm(mesh.numberElements());
  HexaMesh_ProjSca_ElemQ0(ProjDivEm, fxndivu, mesh, GQH);
  HexaMesh_ProjVec_ElemQ03(ProjDsplEm, fxnu, mesh, GQH);
  
  // Computing errors
  double L2ErrDspl=0, L2ErrStrs=0, L2ErrDiv=0;
  LinElas3d_WG_HexaQ03Q03RT03Q0_Err(L2ErrDspl, L2ErrStrs, L2ErrDiv, fxnu, fxnsigma,fxndivu, mesh, NumerDsplCofQ03, NumerStrsCofSym, NumerDivCofQ0, GQH, GQQ);
  Vector ErrDiv = ProjDivEm - NumerDivCofQ0;
  std::cout << "L2ErrDspl= " << L2ErrDspl << "  ";
  std::cout << "L2ErrStrs= " << L2ErrStrs << "  ";
  std::cout << "L2ErrDiv= " << L2ErrDiv;
  std::cout << "\n\n";
  // ProjDivEm.save2file("ProjDivEm.dat");
  // NumerDivCofQ0.save2file("NumerDivCofQ0.dat");
  
  // Output & saving
  double *ScalarData = new double[mesh.numberElements()];
  // For the numerical diveregnce
  for (int i=0; i<mesh.numberElements(); ++i)  ScalarData[i] = NumerDivCofQ0[i];
  Write_ScalarData_HexaMesh_silo(mesh, ScalarData, filename1, varname1);
  //
  // For numerical displacement (3 components) (padded with numer.div. @beginning)
  for (int i=0; i<mesh.numberElements(); ++i)  ScalarData[i] = NumerDivCofQ0[i];
  Write_ScalarData_HexaMesh_silo(mesh, ScalarData, filename3, varname1);
  for (int i=0; i<mesh.numberElements(); ++i)  ScalarData[i] = NumerDsplCofQ0x[i];
  Append_ScalarData_HexaMesh_silo(mesh, ScalarData, filename3, varname3x);
  for (int i=0; i<mesh.numberElements(); ++i)  ScalarData[i] = NumerDsplCofQ0y[i];
  Append_ScalarData_HexaMesh_silo(mesh, ScalarData, filename3, varname3y);
  for (int i=0; i<mesh.numberElements(); ++i)  ScalarData[i] = NumerDsplCofQ0z[i];
  Append_ScalarData_HexaMesh_silo(mesh, ScalarData, filename3, varname3z);

  // Freeing...
  delete[] ProjDsplEm;
  delete[] ScalarData;
  
  // Finishing
  clock_t TimeEnd = clock();
  std::cout << "Time taken = "
    << static_cast<double>(TimeEnd-TimeStart)/CLOCKS_PER_SEC << " seconds \n";
  
  return(0);  // If successful
}
