// main.cpp
// LinElas3d_CG_TetraP13_Proj1
// The main function for sloving lin.elas.eqn. by CGP 1^3 on a tetrahedral mesh
// Project 1: A tetrahedral 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 "IO_TetraMesh.h"
#include "LinElas3d_CG_TetraP13.h"
#include "TetraMesh.h"


int main(int argc, const char * argv[])
{
  char filename1[64] = "LinElas3d_CG_TetraP13_NumerDspl1";
  char filename2[64] = "LinElas3d_CG_TetraP13_ProjDspl1";
  char filename3[64] = "LinElas3d_CG_TetraP13_ProjDiv";
  char varname1[64] = "NumerDspl1";
  char varname2[64] = "ProjDspl1";
  char varname3[64] = "ProjDiv";

  clock_t TimeStart = clock();

  GaussQuad GQTe, GQT;
  GQTe.setForTetrahedron(5);
  GQT.setForTriangle(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;
  
  // Generating a uniform tetrahedral mesh
  std::cout << "Generating a uniform tetrahedral mesh ... " << "\n";
  int n = 8;
  double  nx = n,  ny = n,  nz = n;
  TetraMesh mesh(xa,xb,nx,yc,yd,ny,ze,zf,nz);
  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");

  // Setting up the problem
  double lambda = 1,  mu = 1;
  
  // Assembling and solving
  std::cout << "Assembling and solving...\n";
  int DOFs = 3*mesh.numberNodes();
  std::cout << "DOFs=" << DOFs << "\n";
  SparseBlockMatrix GlbMat;
  Vector GlbVecSource(DOFs), GlbVecDirichlet(DOFs), GlbVecNeumann(DOFs);
  LinElas3d_CG_TetraP13_AsmGlbMat(GlbMat, lambda, mu, mesh);
  LinElas3d_CG_TetraP13_AsmSource(GlbVecSource, fxnf, mesh, GQTe);
  GlbVecSource.save2file("GlbVecSource.dat");
  LinElas3d_CG_TetraP13_AsmBndryConds(GlbVecDirichlet, GlbVecNeumann,
                                      fxnuD, mesh, GQT);
  Vector GlbRHS(DOFs);
  LinElas3d_CG_TetraP13_ModiLinSys(GlbMat, GlbRHS, mesh, GlbVecSource,
                                   GlbVecDirichlet, GlbVecNeumann);
  // Solving by Conjugate Gradient
  int itr = 0,  MaxItr = 5000;
  double thres = 1e-18,  tol = 1e-18;
  Vector sln(DOFs);
  slvSpaBlkSpdSysCG(sln, GlbMat, GlbRHS, itr, MaxItr, thres, tol);
  // BlockDiagMatrix B = diagonal(GlbMat);
  // slvSpaBlkSpdSysPCG(sln, GlbMat, GlbRHS, B, itr, MaxItr, thres, tol);
  // GlbMat.save2file("GlbMat.dat");
  // GlbVecSource.save2file("GlbVecSource.dat");
  // sln.save2file("sln.dat");
  std::cout << "#iterations=" << itr << "\n";
  std::cout << "sln.max.=" << sln.l0norm() << "\n\n";

  // Post-processing
  PtVec3d *NumerDspl = new PtVec3d[mesh.numberElements()];
  double *NumerDiv = new double[mesh.numberElements()];
  Mat3 *NumerStrs = new Mat3[mesh.numberElements()];
  LinElas3d_CG_TetraP13_DsplDivStrs(NumerDspl,NumerDiv,NumerStrs,lambda,mu,mesh,sln);
  
  // Projected displacement & divergence for verification
  PtVec3d *ProjDspl = new PtVec3d[mesh.numberElements()];
  Vector ProjDiv(mesh.numberElements());
  /*
  LinElas3d_CG_TetraP13_ProjDsplDivStrs(ProjDspl, ProjDiv, lambda, mu,
                                        fxnu, fxnsigma, fxndivu, mesh, GQTe);
  */
  
  // Computing errors
  double L2ErrDspl = 0,  L2ErrStrs = 0;
  LinElas3d_CG_TetraP13_Err(L2ErrDspl, L2ErrStrs, lambda, mu,
                            fxnu, fxnsigma, sln, mesh, GQTe);
  std::cout << "L2ErrDspl=" << L2ErrDspl << "  ";
  std::cout << "L2ErrStrs=" << L2ErrStrs << "\n";

  // Preparing for data files
  double *ScalarData = new double[mesh.numberElements()];
  //
  // Numerical displacement: 1st component
  for (int i=1; i<=mesh.numberElements(); ++i) {
    ScalarData[i-1] = NumerDspl[i-1].xCrd();
  }
  Write_ScalarData_TetraMesh_silo(mesh, ScalarData, filename1, varname1);
  //
  // Projected displacement: 1st component
  for (int i=1; i<=mesh.numberElements(); ++i) {
    ScalarData[i-1] = ProjDspl[i-1].xCrd();
  }
  Write_ScalarData_TetraMesh_silo(mesh, ScalarData, filename2, varname2);
  
  // Freeing...
  delete[] NumerDspl, NumerDiv, NumerStrs;
  delete[] ProjDspl;
  delete[] ScalarData;
  
  // Finishing ...
  clock_t TimeEnd = clock();
  std::cout << "\n" << "Time taken = "
    << static_cast<double>(TimeEnd-TimeStart)/CLOCKS_PER_SEC << " seconds.\n";

  return(0);  // If successful
}
