// thermal - A linear solver for elastic problems using FEM
// 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>.
//
// Initial design: Frederic Duboeuf (rev.2285)


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

#include "thermalSolver.h"
#include "genAlgorithms.h"
#include "genSField.h"
#include "linearSystemCSR.h"
#include "linearSystemPETSc.h"
#include "linearSystemId.h"
#include "linearSystemFull.h"

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


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

void ThermalSolver::CheckProperties()
{
  for (int i=0; i<Data->Domains.size(); ++i)
  {
    if (Data->Domains[i]->type=="Thermal")
    {
      ThermalDomain* field = new ThermalDomain(*(Data->Domains[i]));
      TDomains.push_back (field);
    }
    else if (Data->Domains[i]->type=="IsotropicThermal")
    {
      IsotropicThermalDomain* field = new IsotropicThermalDomain(*(Data->Domains[i]));
      TDomains.push_back (field);
    }
  }
}


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




void ThermalSolver::CreateFunctionSpaces()
{
  std::cout << "\nCreateFunctionSpaces :" << std::endl;

  FSpace = genFSpace<genTensor0<double> >::Handle(new ScalarLagrangeFSpace<genTensor0<double> >());

  std::cout << "Thermal Field : iField=" << FSpace->getIncidentSpaceTag() << std::endl;
}


void ThermalSolver::BuildFunctionSpaces()
{
  std::cout << "\nBuildFunctionSpaces :" << std::endl;
  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=="Temperature")
      {
        std::cout <<  "Dirichlet BC Temp" << std::endl;
        FixNodalDofs ( FSpace, Data->BCs[i]->begin(), Data->BCs[i]->end(), *pAssembler, *(Data->BCs[i]->fscalar), genFilterDofTrivial() );
      }
    }
  }

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



void ThermalSolver::AssembleRHS()
{
  std::cout << "\nAssembleRHS :" << std::endl;
  GaussQuadrature Integ_Boundary ( GaussQuadrature::Val );
  double volume=0;
  genTerm<genTensor0<double>,0>::Handle Volume(new ConstantField<genTensor0<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=="Flux")
      {
        std::cout << "Neumann BC Flux" << std::endl;
        genTerm<genTensor0<double>,0>::Handle Func(new FunctionField<genTensor0<double> >(*(Data->BCs[i]->fscalar)));
        genTerm<genTensor0<double>,1>::Handle loadterm = Compose<FullContrProdOp>(FSpace,Func);
        Assemble ( loadterm, Data->BCs[i]->begin(), Data->BCs[i]->end(), Integ_Boundary, *pAssembler );
        volume=0;
        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 ThermalSolver::AssembleLHS()
{
  std::cout << "\nAssembleLHS :" << std::endl;
  double volume=0;
  genTerm<genTensor0<double>,0>::Handle Volume(new ConstantField<genTensor0<double> >(1.0));
  for (int i=0;i<TDomains.size();++i)
  {
    ThermalDomain* T = TDomains[i];

    genTerm<genTensor1<double>,1>::Handle G = Gradient(FSpace);
    genTerm<genTensor1<double>,1>::Handle St = Compose<FullContrProdOp>(T->C,G);
    genTerm<genTensor0<double>,2>::Handle Energy = Compose<FullContrProdOp>(G,St);

    GaussQuadrature Integ_Bulk ( GaussQuadrature::GradGrad );
    Assemble ( Energy, T->begin(), T->end(), Integ_Bulk, *pAssembler );
    volume=0;
    Assemble ( Volume,T->begin(), T->end(), Integ_Bulk, volume);
    printf("Support(%dD,%d) = %f [L^%d]\n", T->dim, T->tag, volume, T->dim);
  }
}

void ThermalSolver::AssembleResidual()
{
  linearSystemId<double>* buf= new linearSystemId<double>;
  dofManager<double>* pAssembler2 = new dofManager<double> ( buf );

  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=="Temperature")
      {
        std::cout <<  "Dirichlet BC Temp" << std::endl;
        FixNodalDofs ( FSpace, Data->BCs[i]->begin(), Data->BCs[i]->end(), *pAssembler2, *(Data->BCs[i]->fscalar), genFilterDofTrivial() );
      }
    }
  }

  // we number the dofs : when a dof is numbered, it cannot be numbered
  // again with another number.
  for ( unsigned int i = 0; i < TDomains.size(); ++i )
  {
    NumberDofs ( FSpace, TDomains[i]->begin(), TDomains[i]->end(),*pAssembler2 );
  }

  for (int i=0;i<TDomains.size();++i)
  {
    ThermalDomain*  T= TDomains[i];

    genTerm<genTensor1<double>,1>::Handle G = Gradient(FSpace);

    diffTerm<genTensor0<double> ,0>::Handle Field ( new genSField<genTensor0<double> >(pAssembler, FSpace) );
    genTerm<genTensor1<double> ,0>::Handle GG = Gradient(Field);
    genTerm<genTensor1<double>,0>::Handle St = Compose<FullContrProdOp>(T->C,GG);

    genTerm<genTensor0<double>,1>::Handle Energy = Compose<FullContrProdOp>(G,St);

    GaussQuadrature Integ_Bulk ( GaussQuadrature::GradGrad );
    Assemble ( Energy, T->begin(), T->end(), Integ_Bulk, *pAssembler2 );
  }

  GaussQuadrature Integ_Boundary ( GaussQuadrature::Val );
  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=="Flux")
      {
        std::cout << "Neumann BC Flux" << std::endl;
        genTerm<genTensor0<double>,0>::Handle Func(new FunctionField<genTensor0<double> >(*(Data->BCs[i]->fscalar)));
        genTerm<genTensor0<double>,1>::Handle loadterm = Compose<FullContrProdOp>(FSpace,Func)*-1.;
        
        Assemble ( loadterm, Data->BCs[i]->begin(), Data->BCs[i]->end(), Integ_Boundary, *pAssembler );
      }
    }
  }

//   FILE* f = fopen ( "residual.txt", "w" );
//   double valeur;
//   std::string sysname = "A";
//   for ( int i = 0 ; i < pAssembler2->sizeOfR() ; i ++ )
//   {
//     pAssembler2->getLinearSystem ( sysname )->getFromRightHandSide ( i,valeur );
//     fprintf ( f,"%+e\n",valeur ) ;
//   }
//   fclose ( f );

  pAssembler2->systemSolve();

  SavedUnknownTerm<genTensor0<double> > saveResidual(FSpace);
  for ( unsigned int i = 0; i < TDomains.size(); ++i )
  {
    SaveUnknown(*pAssembler2,TDomains[i]->begin(), TDomains[i]->end(),saveResidual);
  }
  saveResidual.printfile("saveResidual.txt");


  delete buf;
  delete pAssembler2;
}

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

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

  if ( Data->solvertype == 1 )
  {
    exportKb();
  }
  printf ( "-- done BuildLinearSystem!\n" );
}



void ThermalSolver::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;
  }
  else
  {
    printf ( "Default linearSystemFull is chosen to solve\n" );
    linearSystemFull<double>* buf= new linearSystemFull<double>;
    lsys=buf;
    Data->solvertype = 0;
  }


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

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

  pAssembler->systemSolve();

  printf ( "-- done solving! \n" );

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

  for ( unsigned int i = 0; i < TDomains.size(); i++ )
  {
    ThermalDomain* T = TDomains[i];

    diffTerm<genTensor0<double>,0>::Handle Field( new genSField<genTensor0<double> >(pAssembler, FSpace) );
    genTerm<genTensor1<double>,0>::Handle G = Gradient(Field);
    genTerm<genTensor1<double>,0>::Handle CG ( new genCache<genTensor1<double>,0>(G) ); // use cache
    genTerm<genTensor1<double>,0>::Handle St = Compose<FullContrProdOp>(T->C,CG);
    genTerm<genTensor0<double>,0>::Handle Energy = Compose<FullContrProdOp>(CG,St);

    Assemble ( Energy, T->begin(), T->end(), Integ_Bulk, energ);
    Assemble ( Volume, T->begin(), T->end(), Integ_Bulk, volume);
  }
  printf("Total volume = %f [L^%d]\n", volume, Data->dim);
  printf("Thermal energy = %f\n", energ);
}



void ThermalSolver::exportKb()
{
  bool exportSystem = false;

  std::map<std::string,int>::iterator it;
  it = Data->User.Integers.find("ExportSystem");
  if ( it!=Data->User.Integers.end() )
    if (it->second==1)
      exportSystem = true;

  if (!exportSystem) return;

  FILE* f;

  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);
      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);
}


PView* ThermalSolver::buildTemperatureView(const std::string &postFileName)
{
  genTerm<genTensor0<double>,0>::ConstHandle Field(new genSField<genTensor0<double> >( pAssembler, FSpace ));

  return buildViewNodal(Field,TDomains,postFileName);
}

PView* ThermalSolver::buildThermalStrainView(const std::string &postFileName)
{
  diffTerm<genTensor0<double>,0>::Handle Field(new genSField<genTensor0<double> >( pAssembler, FSpace ));
  genTerm<genTensor1<double>,0>::ConstHandle G = Gradient(Field);

  return buildViewElementNode(G, TDomains, postFileName);
}

PView* ThermalSolver::buildThermalStressView(const std::string &postFileName)
{
  diffTerm<genTensor0<double>,0>::Handle Field(new genSField<genTensor0<double> >( pAssembler, FSpace ));
  genTerm<genTensor1<double>,0>::Handle G = Gradient(Field);

  std::vector<genTerm<genTensor1<double>,0>::ConstHandle> St;
  for( unsigned int i = 0; i < TDomains.size(); i++ )
  {
    ThermalDomain* T = TDomains[i];
    St.push_back(Compose<FullContrProdOp>(T->C,G));
  }

  return buildViewElementNode(St, TDomains, postFileName);
}

PView* ThermalSolver::buildThermalEnergyView(const std::string &postFileName)
{
  diffTerm<genTensor0<double>,0>::Handle Field(new genSField<genTensor0<double> >( pAssembler, FSpace ));
  genTerm<genTensor1<double>,0>::Handle G = Gradient(Field);
  genTerm<genTensor1<double>,0>::Handle CG(new genCache<genTensor1<double>,0>(G));

  std::vector<genTerm<genTensor0<double>,0>::ConstHandle> Energy;
  for ( unsigned int i = 0; i < TDomains.size(); i++ )
  {
    ThermalDomain* T = TDomains[i];
    genTerm<genTensor1<double>,0>::Handle St = Compose<FullContrProdOp>(T->C,CG);
    Energy.push_back(Compose<FullContrProdOp>(CG,St));
  }

  return buildViewElement(Energy, TDomains, postFileName);
}