// main.cpp (MAC version)
// LinElas3d solved by WG(Q0^3,Q0^3;RT[0]^3,Q0) on a hexahedral mesh
// Proj1: A hexa. mesh generated by itself
// The single-matrix approach
// James Liu, 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 filename6[64] = "LinElas3d_WG_HexaQ03Q03RT03Q0_NumerStrs";
  char varname1[64] = "NumerDiv";
  char varname3x[64] = "NumerDspl_x";
  char varname3y[64] = "NumerDspl_y";
  char varname3z[64] = "NumerDspl_z";
  char varname6xx[64] = "NumerStress_xx";
  char varname6yy[64] = "NumerStress_yy";
  char varname6zz[64] = "NumerStress_zz";
  char varname6xy[64] = "NumerStress_xy";
  char varname6xz[64] = "NumerStress_xz";
  char varname6yz[64] = "NumerStress_yz";

  clock_t totalStart = clock();
  
  GaussQuad GQH, GQQ;
  GQH.setForBrick(3,3,3);
  GQQ.setForRectangle(3,3);

  // Setting up physical problem:
  // Physical parameters; Domain; ABC (Array for Boundary Conditions)
  // 1 Dirichlet, 2 Neumann, 3 Robin, 4,5,6: Symmetry in x,y,z respectively
  //
  // Ex.: Standard EqnBC_LinElas3d_Ex11sincoscos
  /*
  double lambda=1, mu=1;
  double  xa = 0,  xb = 1;
  double  yc = 0,  yd = 1;
  double  ze = 0,  zf = 1;
  Vector ABC(6);
  ABC(1) = 1;  ABC(2) = 1;
  ABC(3) = 1;  ABC(4) = 1;
  ABC(5) = 1;  ABC(6) = 1;
  */
  // Ex.: A variant of EqnBC_LinElas3d_Ex11sincoscos
  /*
  double lambda=1, mu=1;
  double  xa = 0,  xb = 1;
  double  yc = 0,  yd = 1;
  double  ze = 0,  zf = 1;
  Vector ABC(6);
  ABC(1) = 4;  ABC(2) = 1;
  ABC(3) = 5;  ABC(4) = 1;
  ABC(5) = 6;  ABC(6) = 1;
  */
  // Ex.: Mijuca_ComputMech_2004 p.479 Ex.6.7
  // double lambda=4.008376878002496e+05, mu=80.194000225029669;
  double E=240.56595979, nu=0.499899987;
  double lambda = E*nu/((1-2*nu)*(1+nu));
  double mu = E/(2*(1+nu));
  double  xa = 0,  xb = 1;
  double  yc = 0,  yd = 1;
  double  ze = 0,  zf = 1;
  Vector ABC(6);
  ABC(1) = 4;  // ABC(2) = 2;
  ABC(3) = 5;  // ABC(4) = 2;
  ABC(5) = 6;  ABC(6) = 2;

  // Generating a smooth hexahedral mesh on the unit cube
  int n = 8;
  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";
  // mesh.save2file("HexaMesh.dat");
  // std::cout << "Saving mesh data finished!\n";

  // JL20170427: The single-matrix approach
  // Assembling and solving
  int DOFs = 3*mesh.numberElements() + 3*mesh.numberFaces();
  std::cout << "DOFs= " << DOFs << "\n";
  
  clock_t assemblyStart = clock();
  SparseBlockMatrix GlbMat;
  Vector GlbVecSource, GlbVecDirichlet, GlbVecNeumann, GlbRHS, sln(DOFs), flag;
  LinElas3d_WG_HexaQ03Q03RT03Q0_AsmGlbMat1(GlbMat, lambda, mu, mesh, GQH, GQQ);
  LinElas3d_WG_HexaQ03Q03RT03Q0_AsmSource1(GlbVecSource, fxnf, mesh, GQH);
  LinElas3d_WG_HexaQ03Q03RT03Q0_AsmBndryConds1(GlbVecDirichlet, GlbVecNeumann, flag,
                                               fxnuD, fxntN, ABC, mesh, GQQ);
  GlbVecNeumann.save2file("GlbVecNeumann.dat");
  LinElas3d_WG_HexaQ03Q03RT03Q0_ModiLinSys1(GlbMat, GlbRHS, ABC, mesh, flag,                                    GlbVecSource, GlbVecDirichlet, GlbVecNeumann);
  clock_t assemblyEnd = clock();
  
  int itr = 0;
  clock_t solveStart = clock();
  slvSpaBlkSpdSysCG(sln, GlbMat, GlbRHS, itr, 5000, 1e-18, 1e-18);
  clock_t solveEnd = clock();
  std::cout << "itr=" << itr << "\n";
  // sln.save2file("sln.dat");
  
  // Postprocessing
  // Computing numerical displacement, stress,
  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);
  }
  // JL20171107: NOW FOR BRICK MESH ONLY
  FullMatrix NumerStrsEmAvg(mesh.numberElements(),6);
  for (int labele=1; labele<=mesh.numberElements(); labele++) {
    NumerStrsEmAvg(labele,1) = NumerStrsCofSym(labele, 1);        // sigmaxx
    NumerStrsEmAvg(labele,4) = NumerStrsCofSym(labele, 8);        // sigmayy
    NumerStrsEmAvg(labele,6) = NumerStrsCofSym(labele,15);        // sigmazz
    NumerStrsEmAvg(labele,2) = NumerStrsCofSym(labele, 2) * 0.5   // sigmaxy
                             + NumerStrsCofSym(labele, 7) * 0.5;
    NumerStrsEmAvg(labele,3) = NumerStrsCofSym(labele, 3) * 0.5   // sigmaxz
                             + NumerStrsCofSym(labele,13) * 0.5;
    NumerStrsEmAvg(labele,5) = NumerStrsCofSym(labele, 9) * 0.5   // sigmayz
                             + NumerStrsCofSym(labele,14) * 0.5;
  }

  // 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 numerical dilation
  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);
  //
  // For numerical stress (6 entries) (padded with numer.div. @beginning)
  for (int i=0; i<mesh.numberElements(); ++i)  ScalarData[i] = NumerDivCofQ0[i];
  Write_ScalarData_HexaMesh_silo(mesh, ScalarData, filename6, varname1);
  for (int i=0; i<mesh.numberElements(); ++i)  ScalarData[i] = NumerStrsEmAvg(i+1,1);
  Append_ScalarData_HexaMesh_silo(mesh, ScalarData, filename6, varname6xx);
  for (int i=0; i<mesh.numberElements(); ++i)  ScalarData[i] = NumerStrsEmAvg(i+1,2);
  Append_ScalarData_HexaMesh_silo(mesh, ScalarData, filename6, varname6xy);
  for (int i=0; i<mesh.numberElements(); ++i)  ScalarData[i] = NumerStrsEmAvg(i+1,3);
  Append_ScalarData_HexaMesh_silo(mesh, ScalarData, filename6, varname6xz);
  for (int i=0; i<mesh.numberElements(); ++i)  ScalarData[i] = NumerStrsEmAvg(i+1,4);
  Append_ScalarData_HexaMesh_silo(mesh, ScalarData, filename6, varname6yy);
  for (int i=0; i<mesh.numberElements(); ++i)  ScalarData[i] = NumerStrsEmAvg(i+1,5);
  Append_ScalarData_HexaMesh_silo(mesh, ScalarData, filename6, varname6yz);
  for (int i=0; i<mesh.numberElements(); ++i)  ScalarData[i] = NumerStrsEmAvg(i+1,6);
  Append_ScalarData_HexaMesh_silo(mesh, ScalarData, filename6, varname6zz);
  
  // Freeing...
  delete[] ProjDsplEm;
  
  // Finishing
  clock_t totalEnd = clock();
  /*
  std::cout << "Assembly time = "
    << static_cast<double>(assemblyEnd-assemblyStart)/CLOCKS_PER_SEC << " seconds \n";
  std::cout << "Solve time = "
    << static_cast<double>(solveEnd-solveStart)/CLOCKS_PER_SEC << " seconds \n";
  */
  std::cout << "Total time = "
    << static_cast<double>(totalEnd-totalStart)/CLOCKS_PER_SEC << " seconds \n";

  return(0);  // If successful
}
