// elastic_genTerm - A linear solver for elastic problems using FEM & new genTensors
// Copyright (C) 2010-2026 Eric Bechet
//
// See the LICENSE file for license information and contributions.
// Please report all bugs and problems to <bechet@cadxfem.org>.

#include <cstring>
#include <cstdio>
#include <fstream>

#include "elasticSolver_ngt.h"
#include "genAlgorithms.h"
#include "genSField.h"
#include "linearSystemCSR.h"
#include "linearSystemPETSc.h"
#include "linearSystemId.h"

#if defined(HAVE_POST)
#include "PView.h"
#include "PViewData.h"
#endif


ElasticSolver_ngt::~ElasticSolver_ngt()
{
  if (pAssembler) delete pAssembler;
  for (int i=0; i<EDomains.size(); ++i) delete EDomains[i];
}

void ElasticSolver_ngt::CheckProperties()
{
  for (int i=0; i<Data->Domains.size(); ++i)
  {
    if (Data->Domains[i]->type=="Elastic")
    {
      ElasticDomain_ngt* field = new ElasticDomain_ngt(*(Data->Domains[i]));
      EDomains.push_back (field);
    }
    else if (Data->Domains[i]->type=="IsotropicElastic")
    {
      IsotropicElasticDomain_ngt* field = new IsotropicElasticDomain_ngt(*(Data->Domains[i]));
      EDomains.push_back (field);
    }
    else if (Data->Domains[i]->type=="OrthotropicElastic")
    {
      OrthotropicElasticDomain_ngt* field = new OrthotropicElasticDomain_ngt(*(Data->Domains[i]));
      EDomains.push_back (field);
    }
  }
}


void ElasticSolver_ngt::readInputFile ( const std::string &fileName )
{
  modelname = fileName;
  std::cout<<"model name : " << modelname << std::endl ;
  Data->readInputFile(fileName);
  CheckProperties();
  printf("--> %d Domains\n", (int)EDomains.size());
  printf("--> %d BCs \n", (int)Data->BCs.size());
}




void ElasticSolver_ngt::CreateFunctionSpaces()
{
  if ( Data->dim==3 )
    FSpaceDisp = genFSpace<genTensor<1,double> >::Handle(new VectorLagrangeFSpaceNG<double>());
//     FSpaceDisp = genFSpace<genTensor<1,double> >::Handle(new VectorLagrangeFSpaceOriented<double>());
  if ( Data->dim==2 )
    FSpaceDisp = genFSpace<genTensor<1,double> >::Handle(new VectorLagrangeFSpaceNG<double>(VectorLagrangeFSpaceNG<double>::VECTOR_X, VectorLagrangeFSpaceNG<double>::VECTOR_Y));

  std::cout << "Displacement Field : iField=" << FSpaceDisp->getIncidentSpaceTag() << std::endl;

  genTerm<genTensor<2,double> ,1>::Handle G(Gradient(FSpaceDisp)*0.5);
  EpsilonDisp = (G+Transpose(G));
}


void ElasticSolver_ngt::BuildFunctionSpaces()
{
  std::cout << *(DofIdManager::getInstance());

  for ( unsigned int i = 0; i < Data->BCs.size(); i++ )
  {
    if ((Data->BCs[i]->kind=="Dirichlet")||(Data->BCs[i]->kind==""))
    {
      if (Data->BCs[i]->What=="Displacement")
      {
        std::cout <<  "Dirichlet BC Disp" << std::endl;
        genFilterDofComponent filter ( Data->BCs[i]->comp1 );
        FixNodalDofs ( FSpaceDisp, Data->BCs[i]->begin(), Data->BCs[i]->end(), *pAssembler, *(Data->BCs[i]->fscalar), filter );
      }
    }
  }

  // we number the dofs : when a dof is numbered, it cannot be numbered
  // again with another number.
  for ( unsigned int i = 0; i < EDomains.size(); ++i )
  {
    ElasticDomain_ngt* E = EDomains[i];
    std::cout << "Domain " << i << " : dim=" << E->dim << ", tag=" << E->tag << ", " << E->group()->size() << " elements" << std::endl;
    NumberDofs ( FSpaceDisp, E->begin(), E->end(), *pAssembler );
  }
}



void ElasticSolver_ngt::AssembleRHS()
{
  GaussQuadrature Integ_Boundary ( GaussQuadrature::Val );
  double volume=0;
  genTerm<genTensor<0,double>,0>::Handle Volume(new ConstantField<genTensor<0,double> >(1.0));

  // build the external force vector
  for ( unsigned int i = 0; i < Data->BCs.size(); i++ )
  {
    if ((Data->BCs[i]->kind=="Neumann")||(Data->BCs[i]->kind==""))
    {
      if (Data->BCs[i]->What=="Force")
      {
        volume=0;
        std::cout << "Neumann BC Force" << std::endl;
        genTerm<genTensor<1,double> ,0>::Handle Func(new FunctionField<genTensor<1,double> >(*(Data->BCs[i]->fvector)));
        genTerm<genTensor<0,double>,1>::Handle loadterm = Compose<FullContrProdOp>(FSpaceDisp,Func);
        Assemble ( loadterm , Data->BCs[i]->begin(), Data->BCs[i]->end(), Integ_Boundary, *pAssembler );
        Assemble ( Volume, Data->BCs[i]->begin(), Data->BCs[i]->end(), Integ_Boundary, volume);
        printf("Support(%dD,%d) = %f [L^%d]\n", Data->BCs[i]->dim, Data->BCs[i]->tag, volume, Data->BCs[i]->dim);
      }
    }
  }
}


void ElasticSolver_ngt::AssembleLHS()
{
  double volume=0;
  genTerm<genTensor<0,double>,0>::Handle Volume(new ConstantField<genTensor<0,double> >(1.0));
  for (int i=0;i<EDomains.size();++i)
  {
    ElasticDomain_ngt* E = EDomains[i];

    // 1st methode
    genTerm<genTensor<0,double>,2>::Handle Total = E->Energy<1,1>(EpsilonDisp,EpsilonDisp);
    
    GaussQuadrature Integ_Bulk ( GaussQuadrature::GradGrad );
    Assemble ( Total, E->begin(), E->end(), Integ_Bulk, *pAssembler );
    volume=0;
    Assemble ( Volume,E->begin(), E->end(), Integ_Bulk,volume);
    printf("Support(%dD,%d) = %f [L^%d]\n", E->dim, E->tag, volume, E->dim);
  }
}



void ElasticSolver_ngt::BuildLinearSystem()
{
  AssembleRHS();
  AssembleLHS();

  printf ( "nDofs=%d\n",pAssembler->sizeOfR() );
  printf ( "-- done assembling!\n" );

  if ( Data->solvertype == 1 )
  {
    exportKbCompressed();
//     exportKbToMatrixMarketFormat();
  }
}




void ElasticSolver_ngt::solve()
{
  linearSystem<double>* lsys=NULL;
  if ( Data->solvertype == 2 )
  {
#if defined(HAVE_TAUCS)
    lsys = new linearSystemCSRTaucs<double>;
#else
    printf ( "Taucs is not installed : Gmm is chosen to solve\n" );
    linearSystemCSRGmm<double>* buf= new linearSystemCSRGmm<double>;
    buf->setNoisy ( 2 );
    lsys=buf;
    Data->solvertype = 1;
#endif
  }
  else if ( Data->solvertype == 3 )
  {
#if defined(HAVE_PETSC)
    lsys = new linearSystemPETSc<double>;
#else
    printf ( "Petsc is not installed : Gmm is chosen to solve\n" );
    linearSystemCSRGmm<double>* buf= new linearSystemCSRGmm<double>;
    buf->setNoisy ( 2 );
    lsys=buf;
    Data->solvertype = 1;
#endif
  }
  else if ( Data->solvertype == 1 )
  {
    linearSystemCSRGmm<double>* buf= new linearSystemCSRGmm<double>;
    buf->setNoisy ( 2 );
    lsys=buf;
    Data->solvertype = 1;
  }


  if ( pAssembler ) delete pAssembler;
  pAssembler = new dofManager<double> ( lsys );

  printf ( "-- start solving\n" );
  CreateFunctionSpaces();
  BuildFunctionSpaces();
  BuildLinearSystem();

  pAssembler->systemSolve();
//  AssembleResidual();
  printf ( "-- done solving! \n" );

//  -_Test_SaveDof_---------------------------------------
//     SavedDiffTerm<genTensor<1,double> > saveDof(FSpaceDisp);
//     for ( unsigned int i = 0; i < EDomains.size(); ++i )
//     {
//       SaveDof(*pAssembler,EDomains[i]->begin(), EDomains[i]->end(),saveDof);
//     }
//     saveDof.printfile("saveDof");

//  -_Test_SaveUnknown_---------------------------------------
//     SavedUnknownTerm<genTensor<1,double> > saveUnknown(FSpaceDisp);
//     for ( unsigned int i = 0; i < EDomains.size(); ++i )
//     {
//       SaveUnknown(*pAssembler,EDomains[i]->_g->begin(), EDomains[i]->_g->end(),saveUnknown);
//     }
//     saveUnknown.printfile("saveUnknown");

//  -_Test_SaveModifUnknown_-----------------------------------
//     SavedUnknownTerm<genTensor<1,double> > saveUnknown2(saveUnknown);
//     for ( unsigned int i = 0; i < EDomains.size(); ++i )
//     {
//       const genTerm<double,1>& saveUnknown_2 = saveUnknown+saveUnknown;
//       SaveModifUnknown(saveUnknown_2,EDomains[i]->_g->begin(), EDomains[i]->_g->end(),saveUnknown2);
//     }
//     saveUnknown2.printfile("saveUnknown*2");

//  -_Test_SaveSolution_---------------------------------------
//     std::string linSystName("A");
//     linearSystem<double>* linearSyst = pAssembler->getLinearSystem(linSystName);
//     linearSyst->zeroSolution();
//     SavedUnknownAddToSolution(saveUnknown2,*pAssembler);
//
//     SavedUnknownTerm<genTensor<1,double> > saveSolution(FSpaceDisp);
//     for ( unsigned int i = 0; i < EDomains.size(); ++i )
//     {
//       SaveUnknown(*pAssembler,EDomains[i]->_g->begin(), EDomains[i]->_g->end(),saveSolution);
//     }
//     saveSolution.printfile("saveSolution");
//
//     SaveddiffTerm<genTensor<1,double> > saveDof2(FSpaceDisp);
//     for ( unsigned int i = 0; i < EDomains.size(); ++i )
//     {
//       SaveDof(*pAssembler,EDomains[i]->_g->begin(), EDomains[i]->_g->end(),saveDof2);
//     }
//     saveDof2.printfile("saveDof2");
//  ----------------------------------------------------

  double energ=0;
  double volume=0;
  GaussQuadrature Integ_Bulk ( GaussQuadrature::GradGrad );
  genTerm<genTensor<0,double>,0>::Handle Volume(new ConstantField<genTensor<0,double> >(1.0));

  for ( unsigned int i = 0; i < EDomains.size(); i++ )
  {
    ElasticDomain_ngt* E = EDomains[i];

    diffTerm<genTensor<1,double> ,0>::Handle Field ( new genSField<genTensor<1,double> >(pAssembler, FSpaceDisp) );
    genTerm<genTensor<2,double> ,0>::Handle G = Gradient(Field);
    genTerm<genTensor<2,double> ,0>::Handle B = (G+Transpose(G))*0.5;
    genTerm<genTensor<0,double> ,0>::Handle Total (Compose<FullContrProdOp>(Field,Field));

    genTerm<genTensor<0,double>,0>::Handle Energy = E->Energy<0,0>(B,B);
    Assemble ( Energy, E->begin(), E->end(), Integ_Bulk, energ);
    Assemble ( Volume, E->begin(), E->end(), Integ_Bulk, volume);
  }
  printf("elastic energy=%f\n", energ/2.);
  printf("Total volume = %f [L^%d]\n", volume, Data->dim);

}


void ElasticSolver_ngt::exportKb()
{
  FILE* f = fopen("K.txt", "w");
  double valeur;
  std::string sysname = "A";
  for (int i=0; i < pAssembler->sizeOfR(); ++i)
  {
    for (int j=0; j < pAssembler->sizeOfR(); ++j)
    {
      pAssembler->getLinearSystem(sysname)->getFromMatrix(i,j,valeur);
      if (valeur)
      fprintf(f, "%+e ", valeur);
    }
    fprintf(f ,"\n");
  }
  fclose(f);

  f = fopen("b.txt", "w");
  for ( int i=0 ; i < pAssembler->sizeOfR(); ++i)
  {
    pAssembler->getLinearSystem(sysname)->getFromRightHandSide(i,valeur);
    fprintf(f, "%+e\n", valeur) ;
  }
  fclose(f);
}

void ElasticSolver_ngt::exportKbCompressed()
{
  FILE* f = fopen("K.txt", "w");
  double valeur;
  std::string sysname = "A";
  for (int i=0; i < pAssembler->sizeOfR(); ++i)
  {
    bool first=true;
    int nz=0;
    for (int j=0; j < pAssembler->sizeOfR(); ++j)
    {
      pAssembler->getLinearSystem(sysname)->getFromMatrix(i,j,valeur);
      if (valeur!=0.)
      {
        if (nz>3) fprintf(f, "Z%d ", nz);
        else
          for (int k=0;k<nz;++k) fprintf(f, "%d ", 0);
        nz=0;
        fprintf(f, "%+e ", valeur);
      }
      else
        nz++;
    }
    if (nz>3) fprintf(f, "Z%d ", nz);
     else
      for (int k=0;k<nz;++k) fprintf(f, "%d ", 0);
    fprintf(f ,"\n");
  }
  fclose(f);

  f = fopen("b.txt", "w");
  for ( int i=0 ; i < pAssembler->sizeOfR(); ++i)
  {
    pAssembler->getLinearSystem(sysname)->getFromRightHandSide(i,valeur);
    fprintf(f, "%+e\n", valeur) ;
  }
  fclose(f);
}


void ElasticSolver_ngt::exportKbToMatrixMarketFormat()
{
  double valeur;
  std::string sysname = "A";
  long seek;
  int numVal;

  int size_i = pAssembler->sizeOfR();
  int size_j = pAssembler->sizeOfR();

  FILE* f = fopen("K.mtx", "w");
  fprintf(f, "%%%%MatrixMarket matrix coordinate real general\n");
  fprintf(f, "%%\n");
  fprintf(f, "%% MM File Characteristics :\n");
  fprintf(f, "%% http://people.sc.fsu.edu/~jburkardt/data/mm/mm.html\n");

  seek = ftell(f);
  fprintf(f, "%d %d %10d\n", size_i, size_j, 0);

  numVal=0;
  for (int j=0; j < size_j; ++j)
  {
    for (int i=0; i < size_i; ++i)
    {
      pAssembler->getLinearSystem(sysname)->getFromMatrix(i,j,valeur);
      if (valeur>1.e-100)
      {
        fprintf(f, "%d %d %+1.13e\n", i+1 ,j+1, valeur);
        ++numVal;
      }
    }
  }
  fseek(f, seek, SEEK_SET);
  fprintf(f, "%d %d %10d\n", size_i, size_j, numVal);
  fseek(f, 0, SEEK_END);
  fclose(f);


  f = fopen("b.mtx", "w");

  fprintf(f, "%%%%MatrixMarket matrix coordinate real general\n");
  fprintf(f, "%%\n");
  fprintf(f, "%% MM File Characteristics :\n");
  fprintf(f, "%% http://people.sc.fsu.edu/~jburkardt/data/mm/mm.html\n");


  seek = ftell(f);
  fprintf(f, "%d 2 %10d\n", size_i, 0);

  numVal=0;
  for ( int i=0 ; i < size_i; ++i)
  {
    pAssembler->getLinearSystem(sysname)->getFromRightHandSide(i,valeur);
    if (valeur>1.e-100)
      {
        fprintf(f, "%d 1 %+1.13e\n", i+1, valeur);
        ++numVal;
      }
  }
  fseek(f, seek, SEEK_SET);
  fprintf(f, "%d 2 %10d\n", size_i, numVal);
  fseek(f, 0, SEEK_END);
  fclose(f);
}


PView* ElasticSolver_ngt::buildDisplacementView(const std::string &postFileName)
{
  genTerm<genTensor<1,double>,0>::ConstHandle Field(new genSField<genTensor<1,double> >( pAssembler, FSpaceDisp ));
  return buildViewNodal(Field,EDomains,postFileName);
}

PView* ElasticSolver_ngt::buildElasticEnergyView(const std::string &postFileName)
{
  diffTerm<genTensor<1,double> ,0>::Handle Field(new genSField<genTensor<1,double> >( pAssembler, FSpaceDisp ));
  genTerm<genTensor<2,double> ,0>::Handle G = Gradient(Field);
  genTerm<genTensor<2,double> ,0>::Handle B = (G+Transpose(G))*0.5;

  std::vector<genTerm<genTensor<0,double>,0>::ConstHandle> Energy;
  for ( unsigned int i = 0; i < EDomains.size(); i++ )
  {
    ElasticDomain_ngt* E = EDomains[i];
    Energy.push_back(E->Energy<0,0>(B,B));
  }
  return buildViewElement(Energy, EDomains, postFileName);
}

PView* ElasticSolver_ngt::buildStressView(const std::string &postFileName)
{
  diffTerm<genTensor<1,double>,0>::Handle Field(new genSField<genTensor<1,double> >( pAssembler, FSpaceDisp ));
  genTerm<genTensor<2,double>,0>::Handle G = Gradient(Field);
  genTerm<genTensor<2,double>,0>::Handle B = (G+Transpose(G))*0.5;

  std::vector<genTerm<genTensor<2,double>,0>::ConstHandle> Stress;
  for( unsigned int i = 0; i < EDomains.size(); i++ )
  {
    ElasticDomain_ngt* E = EDomains[i];
    Stress.push_back(E->Stress<0>(B));
  }

  return buildViewElementNode(Stress, EDomains, postFileName);
//   return buildViewElement(Stress, EDomains, postFileName);
}

PView* ElasticSolver_ngt::buildStrainView(const std::string &postFileName)
{
  diffTerm<genTensor<1,double>,0>::Handle Field(new genSField<genTensor<1,double> >( pAssembler, FSpaceDisp ));
  genTerm<genTensor<2,double>,0>::Handle G = Gradient(Field);
  genTerm<genTensor<2,double>,0>::ConstHandle Strain = (G+Transpose(G))*0.5;

  return buildViewElementNode(Strain, EDomains, postFileName);
//   return buildViewElement(*Strain, EDomains, postFileName);
}