// GenFem - A high-level finite element library
// 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 "interfaceXfemSolver.h"
#include "spaceReducer.h"
#include "genAlgorithms.h"
#include "genSField.h"
#include "genFunction.h"
#include "linearSystemCSR.h"
#include "linearSystemPETSc.h"
#include "linearSystemId.h"


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



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

void InterfaceXfemSolver::CheckProperties()
{
  for (int i=0; i<Data->Domains.size(); ++i)
  {
    if (Data->Domains[i]->type=="Interface")
    {
      InterfaceDomain* field = new InterfaceDomain(*(Data->Domains[i]));
      IDomains.push_back (field);
    }
  }
}

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


void InterfaceXfemSolver::CreateFunctionSpaces()
{
#ifdef VECT_SPACE
  if ( Data->dim==3 )
  {
    FSpace = genFSpace<genTensor1<double> >::Handle(new VectorLagrangeFSpace<double>());
  }
  if ( Data->dim==2 )
  {
    FSpace = genFSpace<genTensor1<double> >::Handle(new VectorLagrangeFSpace<double>(VectorLagrangeFSpace<double>::VECTOR_X, VectorLagrangeFSpace<double>::VECTOR_Y));
  }
#else
  FSpace = genFSpace<genTensor0<double> >::Handle(new ScalarLagrangeFSpace< genTensor0<double> >());
#endif
}

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

  for ( unsigned int i = 0; i < IDomains.size(); ++i)
  {
    SpaceReducer spaceReducer(IDomains[i]->group());
    spaceReducer.BuildLinearConstraints(FSpace, pAssembler);
//     spaceReducer.BuildLinearConstraintsGroups(FSpace, pAssembler);
    if(Data->Model->getNumRegions() && IDomains[i]->dim==1)
    {
      genFilterDofMultiEntities* filterDof = new genFilterDofMultiEntities(spaceReducer.FilteredEntities());
#ifdef VECT_SPACE
      FSpace = genFSpace<genTensor1<double> >::Handle(new FilteredFSpace<genTensor1<double> >(FSpace, filterDof));
#else
      FSpace = genFSpace<genTensor0<double> >::Handle(new FilteredFSpace<genTensor0<double> >(FSpace, filterDof));
#endif
    }
  }

  for ( unsigned int i = 0; i < Data->BCs.size(); i++ )
  {
    if ((Data->BCs[i]->kind=="Dirichlet")||(Data->BCs[i]->kind==""))
    {
#ifdef VECT_SPACE
      if (Data->BCs[i]->What=="Displacement")
      {
        std::cout <<  "Dirichlet BC Disp" << std::endl;
        genFilterDofComponent filter ( Data->BCs[i]->comp1 );
        FixNodalDofs ( FSpace, Data->BCs[i]->begin(), Data->BCs[i]->end(), *pAssembler, *(Data->BCs[i]->fscalar), filter );
      }
#else
     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());
        genSimpleFunctionConstant<genTensor0<double> > ori(1.0);
#ifdef Test // defined (or not) in spaceReduce.h
#else
        FixNodalDofsLC ( FSpace, Data->BCs[i]->begin(), Data->BCs[i]->end(), *pAssembler, *(Data->BCs[i]->fscalar), ori, genFilterDofTrivial());
#endif
      }
#endif
    }
  }

  // we number the dofs
  for ( unsigned int i = 0; i < IDomains.size(); ++i)
  {
    NumberDofs ( FSpace, IDomains[i]->begin(), IDomains[i]->end(), *pAssembler );
  }
}

void InterfaceXfemSolver::AssembleRHS()
{
  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==""))
    {
#ifdef VECT_SPACE
      if (Data->BCs[i]->What=="Force")
      {
        volume=0;
        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>(FSpace,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);
      }
#else
      if (Data->BCs[i]->What=="Flux")
      {
        volume=0;
        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 );
        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);
      }
#endif
    }
  }
}

void InterfaceXfemSolver::AssembleLHS()
{
  GaussQuadrature Integ_Bulk ( GaussQuadrature::GradGrad );
  GaussQuadrature Integ_Boundary ( GaussQuadrature::Val );
  double volume=0;
  genTerm<genTensor0<double>,0>::Handle Volume(new ConstantField<genTensor0<double> >(1.0));
  for ( unsigned int i = 0; i < IDomains.size(); ++i)
  {
    InterfaceDomain* I = IDomains[i];
#ifdef VECT_SPACE
    genTerm<genTensor2<double>,1>::Handle G = Gradient(FSpace);
    genTerm<genTensor2<double>,1>::Handle B = (G+Transpose(G))*0.5;
    genTerm<genTensor2<double>,1>::Handle CB ( new genCache<genTensor2<double> ,1>(B) ); // use cache

    genTerm<genTensor2<double>,1>::Handle St = CB*I->k;

    genTerm<genTensor0<double>,2>::Handle Energy = Compose<FullContrProdOp>( Transpose(CB), St );

    Assemble ( Energy, I->begin(), I->end(), Integ_Bulk, *pAssembler );
#else
    genTerm<genTensor1<double>,1>::Handle Gs = Gradient(FSpace);

    genTerm<genTensor1<double>,1>::Handle Sts = Gs*I->k;

    genTerm<genTensor0<double>,2>::Handle Energys = Compose<FullContrProdOp>( Gs, Sts );

    Assemble ( Energys, I->begin(), I->end(), Integ_Bulk, *pAssembler );
#endif
    Assemble ( Volume, I->begin(), I->end(), Integ_Boundary, volume);
    printf("Support(%dD,%d) = %f [L^%d]\n", I->dim, I->tag, volume, I->dim);
  }
}

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

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

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

class CFF : public genSimpleFunction<genTensor0<double> >
{
  typedef typename TensorialTraits<genTensor0<double> >::Scalar Scalar;
  typedef typename TensorialTraits<genTensor0<double> >::ScalarType ScalarType;
  typedef typename TensorialTraits<genTensor0<double> >::ValType ValType;
  typedef typename TensorialTraits<genTensor0<double> >::GradType GradType;
  typedef typename TensorialTraits<genTensor0<double> >::HessType HessType;
  virtual ValType operator()(double x, double y, double z) const { return 3*y;/*+z;*/ };
  virtual GradType grad(double x, double y, double z) const {return GradType(0.0,3.0,0.0);}
  virtual HessType hess(double x, double y, double z) const {return HessType(0.0);}
};




void InterfaceXfemSolver::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();
  CFF f;
  diffTerm<genTensor0<double>,0>::Handle Func(new genFunction<genTensor0<double> >(f));
#ifdef Test // defined (or not) in spaceReduce.h
  for ( unsigned int i = 0; i < IDomains.size(); ++i)
  {
    InterfaceDomain* I = IDomains[i];
    GaussQuadrature Integ_Bulk ( GaussQuadrature::ValVal);
    genTerm<genTensor0<double>,2>::Handle L2proj = Compose<FullContrProdOp>( FSpace, FSpace );
    Assemble ( L2proj, I->begin(), I->end(), Integ_Bulk, *pAssembler );
    genTerm<genTensor0<double>,1>::Handle loadterm = Compose<FullContrProdOp>(FSpace,Func);
    Assemble ( loadterm, I->begin(), I->end(), Integ_Bulk, *pAssembler );
  }
  if ( Data->solvertype == 1 )
  {
    exportKb();
  }
#else
  BuildLinearSystem();
#endif


  pAssembler->systemSolve();
  printf ( "-- done solving! \n" );
  int ndofs=pAssembler->sizeOfR();
  double energ=0;
  double energex=0;
  double sqsol=0;
  double sqsolex=0;
  double volume=0;
//   GaussQuadrature Integ_Bulk ( GaussQuadrature::ValVal );
  GaussQuadrature Integ_Bulk ( GaussQuadrature::GradGrad );
  GaussQuadrature Integ_Boundary ( GaussQuadrature::Val );
  GaussQuadrature Integ_Bulk2 ( GaussQuadrature::ValVal );
  genTerm<genTensor0<double>,0>::Handle Volume(new ConstantField<genTensor0<double> >(1.0));

  for ( unsigned int i = 0; i < IDomains.size(); i++ )
  {
    InterfaceDomain* I = IDomains[i];
#ifdef VECT_SPACE
    diffTerm<genTensor1<double>,0>::Handle Field ( new genSField<genTensor1<double> >(pAssembler, FSpace) );
    genTerm<genTensor2<double>,0>::Handle G = Gradient(Field);
    genTerm<genTensor2<double>,0>::Handle B = (G+Transpose(G))*0.5;
    genTerm<genTensor2<double>,0>::Handle CB ( new genCache<genTensor2<double>,0>(B) );
    genTerm<genTensor2<double>,0>::Handle St = CB*I->k;
    genTerm<genTensor0<double>,0>::Handle Energy = Compose<FullContrProdOp>( Transpose(CB), St );
#else
    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 St = G*I->k;
    genTerm<genTensor0<double>,0>::Handle Energy = Compose<FullContrProdOp>(G,St);
    genTerm<genTensor0<double>,0>::Handle Ival = Compose<FullContrProdOp>(Field,Field);

    genTerm<genTensor1<double>,0>::Handle Gex = Gradient(Func);
    genTerm<genTensor1<double>,0>::Handle Stex = Gex*I->k;
    genTerm<genTensor0<double>,0>::Handle Energyex = Compose<FullContrProdOp>(Gex,Stex);
    genTerm<genTensor0<double>,0>::Handle Ivalex = Compose<FullContrProdOp>(Func,Func);
#endif
    Assemble ( Energy, I->begin(), I->end(), Integ_Bulk, energ);
    Assemble ( Volume, I->begin(), I->end(), Integ_Boundary, volume);
#ifdef Test // defined (or not) in spaceReduce.h
    Assemble ( Energyex, I->begin(), I->end(), Integ_Bulk, energex);
    Assemble ( Ival, I->begin(), I->end(), Integ_Bulk2, sqsol);
    Assemble ( Ivalex, I->begin(), I->end(), Integ_Bulk2, sqsolex);
    printf("sqsol ex=%.15f\n", sqsolex);
    printf("elastic energy ex=%.15f\n", energex);
    printf("sqsol =%.15f\n", sqsol);
#endif
  }
  printf("elastic energy=%.15f\n", energ);
  printf("Total volume = %.15f [L^%d]\n", volume, Data->dim);
}


void InterfaceXfemSolver::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* InterfaceXfemSolver::buildView(const std::string &postFileName)
{
#ifdef VECT_SPACE
  genTerm<genTensor1<double>,0>::ConstHandle Field(new genSField<genTensor1<double> >( pAssembler, FSpace ));
#else
  genTerm<genTensor0<double>,0>::ConstHandle Field(new genSField<genTensor0<double> >( pAssembler, FSpace ));
#endif
  return buildViewNodal(Field, IDomains, postFileName);
}

PView* InterfaceXfemSolver::buildViewG(const std::string &postFileName)
{
#ifdef VECT_SPACE
  diffTerm<genTensor1<double>,0>::ConstHandle Field(new genSField<genTensor1<double> >( pAssembler, FSpace ));
  genTerm<genTensor2<double>,0>::ConstHandle GField=Gradient(Field);
#else
  diffTerm<genTensor0<double>,0>::ConstHandle Field(new genSField<genTensor0<double> >( pAssembler, FSpace ));
  genTerm<genTensor1<double>,0>::ConstHandle GField=Gradient(Field);
#endif
  return buildViewElement(GField, IDomains, postFileName);
}