// Piezo - A linear solver for piezo problem 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>.

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

#include "piezoSolver.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


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

void PiezoSolver::CheckProperties()
{
  for (int i=0;i<Data->Domains.size();++i)
  {
    if (Data->Domains[i]->type=="Piezoelectric")
    {
      PiezoDomain* field = new PiezoDomain(*(Data->Domains[i]));
      PDomains.push_back (field);
    }
  }
}


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




void PiezoSolver::BuildFunctionSpaces()
{
  if ( Data->dim==3 )
    FSpaceDisp = genFSpace<genTensor1<double> >::Handle(new VectorLagrangeFSpace<double>());
  if ( Data->dim==2 )
    FSpaceDisp = genFSpace<genTensor1<double> >::Handle(new VectorLagrangeFSpace<double>(VectorLagrangeFSpace<double>::VECTOR_Y, VectorLagrangeFSpace<double>::VECTOR_Z));

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

  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 );
      }
      if (Data->BCs[i]->What=="Potential")
      {
        std::cout <<  "Dirichlet BC Pot " << std::endl;
        FixNodalDofs ( FSpacePot, 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 < PDomains.size(); ++i )
  {
    NumberDofs ( FSpaceDisp, PDomains[i]->begin(), PDomains[i]->end(),*pAssembler );
    NumberDofs ( FSpacePot, PDomains[i]->begin(), PDomains[i]->end(),*pAssembler );
  }
  genTerm<genTensor2<double>,1>::Handle G(Gradient(FSpaceDisp)*0.5);
  genTerm<genTensor2<double>,1>::Handle GC(new genCache<genTensor2<double>,1>(G));
  Eps= (GC+Transpose(GC));
  EField=Gradient(FSpacePot)*-1.0;
}




void PiezoSolver::AssembleRHS()
{
  GaussQuadrature Integ_Boundary ( GaussQuadrature::Val );

  // 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")
      {
        std::cout << "Neumann BC Force" << std::endl;
        genTerm<genTensor1<double>,0>::Handle Func(new FunctionField<genTensor1<double> >(*(Data->BCs[i]->fvector)));
        genTerm<genTensor0<double>,1>::Handle loadterm(Compose<FullContrProdOp>(FSpaceDisp,Func));
        Assemble ( loadterm , Data->BCs[i]->begin(), Data->BCs[i]->end(), Integ_Boundary, *pAssembler );
      }
      if (Data->BCs[i]->What=="Charge")
      {
        std::cout << "Neumann BC Charge" << std::endl;
        genTerm<genTensor0<double>,0>::Handle Func(new FunctionField<genTensor0<double> >(*(Data->BCs[i]->fscalar)));
        genTerm<genTensor0<double>,1>::Handle loadterm(Compose<FullContrProdOp>(FSpacePot,Func));
        Assemble ( loadterm , Data->BCs[i]->begin(), Data->BCs[i]->end(), Integ_Boundary, *pAssembler );
      }
    }
  }
}

void PiezoSolver::AssembleLHS()
{
  for (int i=0;i<PDomains.size();++i)
  {
    PiezoDomain* P = PDomains[i];

//     genTerm<double,2>::Handle BtCB=Compose<FullContrProdOp>(Transpose(Eps),Compose<FullContrProdOp>(P->C,Eps));
//     genTerm<double,2>::Handle BtEP=Compose<FullContrProdOp>(Transpose(Eps),Compose<FullContrProdOp>(Transpose(P->E),EField));
//     genTerm<double,2>::Handle PtEB=Compose<FullContrProdOp>(Transpose(EField),Compose<FullContrProdOp>(P->E,Eps));
//     genTerm<double,2>::Handle PtKP=Compose<FullContrProdOp>(Transpose(EField),Compose<FullContrProdOp>(P->K,EField));
//
//     genTerm<double,2>::Handle Mech=Compose<FullContrProdOp>(Transpose(Eps),P->Stress<1>(Eps,EField));
//     genTerm<double,2>::Handle Elec=Compose<FullContrProdOp>(Transpose(EField),P->Edisp<1>(Eps,EField));
//
//     genTerm<double,2>::Handle Total=(BtCB+BtEP*-1.0)+(PtEB*-1.0+PtKP*-1);
//     genTerm<double,2>::Handle Total=Mech+Elec*-1.0;
    genTerm<genTensor0<double>,2>::Handle Total=P->Enthalpy<1,1>(Eps,EField,Eps,EField);
    // works now !!!

    GaussQuadrature Integ_Bulk ( GaussQuadrature::GradGrad );
/*     Assemble ( BtCB ,_Domains[i]->_g->begin(),_Domains[i]->_g->end(),Integ_Bulk,*_pAssembler );
    Assemble ( BtEP ,_Domains[i]->_g->begin(),_Domains[i]->_g->end(),Integ_Bulk,*_pAssembler );
    Assemble ( PtEB ,_Domains[i]->_g->begin(),_Domains[i]->_g->end(),Integ_Bulk,*_pAssembler );
    Assemble ( PtKP ,_Domains[i]->_g->begin(),_Domains[i]->_g->end(),Integ_Bulk,*_pAssembler );
*/
/*    Assemble ( Mech ,_Domains[i]->_g->begin(),_Domains[i]->_g->end(),Integ_Bulk,*_pAssembler );
    Assemble ( Elec ,_Domains[i]->_g->begin(),_Domains[i]->_g->end(),Integ_Bulk,*_pAssembler );
*/
    Assemble ( Total ,P->begin(),P->end(),Integ_Bulk,*pAssembler );
  }
}

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

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

  if ( Data->solvertype == 1 )
    exportKb();
}

void PiezoSolver::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;
  }

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



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

  pAssembler->systemSolve();

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

}


void PiezoSolver::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 );
      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* PiezoSolver::buildDisplacementView ( const std::string &postFileName )
{
  genTerm<genTensor1<double>,0>::ConstHandle Field(new genSField<genTensor1<double> >( pAssembler, FSpaceDisp ));
  return buildViewNodal(Field,PDomains,postFileName);
}

PView* PiezoSolver::buildPotentialView ( const std::string &postFileName )
{
  genTerm<genTensor0<double>,0>::ConstHandle Field(new genSField<genTensor0<double> >( pAssembler, FSpacePot ));
  return buildViewNodal(Field,PDomains,postFileName);
}

PView* PiezoSolver::buildElasticEnergyView ( const std::string &postFileName )
{
  diffTerm<genTensor1<double>,0>::Handle Field(new genSField<genTensor1<double> >( pAssembler, FSpaceDisp ));
  diffTerm<genTensor0<double>,0>::Handle Fieldp(new genSField<genTensor0<double> >( pAssembler, FSpacePot ));
  genTerm<genTensor2<double>,0>::Handle G(Gradient(Field)*0.5);
  genTerm<genTensor2<double>,0>::Handle CB( (G+Transpose(G)) );
  genTerm<genTensor1<double>,0>::Handle CP=Gradient(Fieldp)*-1.0;

  std::vector<genTerm<genTensor0<double>,0>::ConstHandle> Energy;
  for ( unsigned int i = 0; i < PDomains.size(); i++ )
  {
    PiezoDomain* P = PDomains[i];
    Energy.push_back(P->Energy<0,0>(CB,CP,CB,CP));
  }

  return buildViewElement(Energy, PDomains, postFileName);
}