// Membrane - A linear solver for membrane 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 "membraneSolver.h"
#include "genAlgorithms.h"
#include "genSField.h"
#include "membraneSpace.h"
#include "linearSystemCSR.h"
#include "linearSystemPETSc.h"
#include "linearSystemId.h"

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

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

void membraneSolver::CheckProperties()
{

  for (int i=0;i<Data->Domains.size();++i)
  {
    if(Data->Domains[i]->type!="IsotropicElastic") Msg::Error("Only IsotropicElastic for now !!");
/*
    if (Data->Domains[i]->type=="Elastic")
    {
      membraneDomain* field = new membraneDomain(*(Data->Domains[i]));
      EDomains.push_back (field);
    }
*/
//    else if (Data->Domains[i]->type=="IsotropicElastic")
//    {
      IsotropicMembraneDomain* field = new IsotropicMembraneDomain(*(Data->Domains[i]));
      EDomains.push_back (field);
//    }
/*    else if (Data->Domains[i]->type=="OrthotropicElastic")
    {
      OrthotropicMembraneDomain* field = new OrthotropicMembraneDomain(*(Data->Domains[i]));
      EDomains.push_back (field);
    }
*/
  }
}


void membraneSolver::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 membraneSolver::BuildFunctionSpaces()
{
  // Data->dim is unused but has to be 2 to create the groupOfElements GBDTAG
  FSpaceDisp = genFSpace<genTensor1<double> >::Handle(new membraneSpace()); // Cache it GBDTAG
 // FSpaceDisp= genTerm<genTensor1<double>,1>::Handle(new genCache<genTensor1<double>,1>(FSpaceDispNocache));


  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 )
  {
    std::cout << "Domain " << i << " : " << EDomains[i]->group()->size() << " elements"  << std::endl;
    NumberDofs ( FSpaceDisp, EDomains[i]->begin(), EDomains[i]->end(),*pAssembler );
  }
  genTerm<genTensor2<double>,1>::Handle G(Gradient(FSpaceDisp));
  genTerm<genTensor2<double>,1>::Handle CG(new genCache<genTensor2<double>,1>(G));
  EpsilonDisp = CG;
}

// Same as elastic_solver but not for Quadrature rule GBDTAG
void membraneSolver::AssembleRHS()
{
  GaussQuadrature Integ_Boundary ( 3 ); // CANNOT WORK FOR PRESSURE GBDTAG ??

  // 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 = FullContractedProductCompose(FSpaceDisp,Func);
        Assemble ( loadterm , Data->BCs[i]->begin(), Data->BCs[i]->end(), Integ_Boundary, *pAssembler );
      }
    }
  }
}

// Same as elastic_solver but != Quadrature Rule and domain (GenDom with getQuadrature and getEnergy ??)  GBDTAG
void membraneSolver::AssembleLHS()
{
  for (int i=0;i<EDomains.size();++i)
  {
    membraneDomain* E = EDomains[i];

    // 1st methode
    genTerm<genTensor2<double>,1>::Handle StressDisp = E->Stress<1>(FSpaceDisp);
    genTerm<genTensor2<double>,1>::Handle CS(new genCache<genTensor2<double>,1>(StressDisp));
    genTerm<genTensor0<double>,2>::Handle Total = E->Energy<1,1>(CS,EpsilonDisp);

    GaussQuadrature Integ_Bulk ( 3 );
    Assemble ( Total, E->begin(), E->end(), Integ_Bulk, *pAssembler );


    // 2nd methode
//     genTerm<genTensor2<double>,1>::Handle G = Gradient(FSpaceDisp);
//     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<double,2>::Handle Mech = FullContractedProductCompose( Transpose(CB), E->Stress<1>(CB) );
//
//     GaussQuadrature Integ_Bulk ( GaussQuadrature::GradGrad );
//     Assemble ( Mech, E->g()->begin(), E->g()->end(), Integ_Bulk, *pAssembler );


    // 3rd methode
//    genTerm<genTensor2<double>,1> &G = *Gradient(diffTerm<genTensor1<double>,1>::Handle(FSpaceDisp));
//    const genTerm<genTensor2<double>,1> &B = (G+Transpose(G))*0.5;
//    genCache<genTensor2<double>,1> CB(B); // use cache
//    genTerm<genTensor4<double>,0>::ConstHandle C = E->Tangent<1>(EpsilonDisp);

//     BilinearTermContractWithLaw<genTensor2<double>> K(B,C,B);
//     GaussQuadrature Integ_Bulk ( GaussQuadrature::GradGrad );
//     Assemble ( K ,E->_g->begin(),E->_g->end(),Integ_Bulk,*pAssembler );

//     genTerm<double,2>::Handle BtCB = FullContractedProductCompose(Transpose(CB),FullContractedProductCompose(C,CB));
//     GaussQuadrature Integ_Bulk ( GaussQuadrature::GradGrad );
//     Assemble ( BtCB ,E->_g->begin(),E->_g->end(),Integ_Bulk,*pAssembler );

//  -_Test_Manager_---------------------------------------
//      std::shared_ptr<SavedGenTerm<genTensor2<double>,0> > SaveStress(new SavedGenTerm<genTensor2<double>,0>());
//      genTerm<genTensor2<double>,0>::Handle SaveStress2(SaveStress);
//      std::shared_ptr<SavedGenTerm<genTensor2<double>,0> > SaveStress3 = std::dynamic_pointer_cast<SavedGenTerm<genTensor2<double>,0> >(SaveStress2);
//      SavedGenTerm<genTensor2<double>,0>* SaveStress4 = &*SaveStress;

//      std::shared_ptr<SavedGenTerm<genTensor4<double>,0> > saveC(new SavedGenTerm<genTensor4<double>,0>());
//      GaussQuadrature Integ_Val ( GaussQuadrature::Val );
//      SavePtGauss(*C,E->g()->begin(), E->g()->end(),Integ_Val,*saveC);
//      savedGenTermManager.addSavedGenTerm("C",saveC,0,0,i);
//      savedGenTermManager.print();
//  ----------------------------------------------------
  }
}


/*
void membraneSolver::AssembleResidual()
{

//   linearSystemCSRGmm<double>* buf= new linearSystemCSRGmm<double>;
  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=="Displacement")
      {
        genFilterDofComponent filter ( Data->BCs[i]->comp1 );
        FixNodalDofs ( FSpaceDisp, Data->BCs[i]->begin(), Data->BCs[i]->end(), *pAssembler2, *(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 )
  {
    NumberDofs ( FSpaceDisp, EDomains[i]->begin(), EDomains[i]->end(),*pAssembler2 );
  }

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

    // 1st methode
    genTerm<genTensor2<double>,1>::Handle G = Gradient(FSpaceDisp);
    genTerm<genTensor2<double>,1>::Handle B = (G+Transpose(G))*0.5;

    diffTerm<genTensor1<double>,0>::Handle Field ( new genSField<genTensor1<double>>(pAssembler, FSpaceDisp) );
    genTerm<genTensor2<double>,0>::Handle GG = Gradient(Field);
    genTerm<genTensor2<double>,0>::Handle BB = (GG+Transpose(GG))*0.5;
    genTerm<double,1>::Handle Mech = FullContractedProductCompose( Transpose(B), E->Stress<0>(BB) );

    GaussQuadrature Integ_Bulk ( GaussQuadrature::GradGrad );
    Assemble ( Mech, E->begin(), E->end(), Integ_Bulk, *pAssembler2 );

    // 2nd methode
//    diffTerm<genTensor1<double>,0>::Handle Field ( new genSField<genTensor1<double>>(pAssembler, FSpaceDisp) );
//     genTerm<genTensor2<double>,0>::Handle GG = Gradient(Field);
//     genTerm<genTensor2<double>,0>::Handle BB = (GG+Transpose(GG))*0.5;
//    GaussQuadrature Integ_Bulk( GaussQuadrature::GradGrad );
//     Assemble(*E->Energy<1,0>(EpsilonDisp,BB), E->g()->begin(), E->g()->end(), Integ_Bulk, *pAssembler2 );

    // 3rd methode;
//    genTerm<genTensor4<double>,0>::Handle C = E->Tangent<0>(CB);
//    diffTerm<genTensor1<double>,0>::Handle Field ( new genSField<genTensor1<double>>(pAssembler, FSpaceDisp) );
//    genTerm<genTensor2<double>,1>::Handle G = Gradient(FSpaceDisp);
//    genTerm<genTensor2<double>,1>::Handle B = (G+Transpose(G))*0.5;
//    genTerm<genTensor2<double>,0>::Handle GG = Gradient(Field);
//    genTerm<genTensor2<double>,0>::Handle BB = (GG+Transpose(GG))*0.5;
//    genTerm<double,1>::Handle BtCB = FullContractedProductCompose(Transpose(B),FullContractedProductCompose(C,BB));
//    Assemble ( *BtHB ,E->_g->begin(),E->_g->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=="Force")
      {
        genTerm<genTensor1<double>,0>::Handle Func(new FunctionField<genTensor1<double>>(*(Data->BCs[i]->fvector)));
        genTerm<double,1>::Handle loadterm = FullContractedProductCompose(FSpaceDisp,Func)*-1;
        Assemble ( loadterm , Data->BCs[i]->begin(), Data->BCs[i]->end(), Integ_Boundary, *pAssembler2 );
      }
    }
  }

//   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<genTensor1<double>> saveResidual(FSpaceDisp);
  for ( unsigned int i = 0; i < EDomains.size(); ++i )
  {
    SaveUnknown(*pAssembler2,EDomains[i]->begin(), EDomains[i]->end(),saveResidual);
  }
  saveResidual.printfile("saveResidual.txt");


  delete buf;
  delete pAssembler2;

}
*/


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

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

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




void membraneSolver::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" );
  BuildFunctionSpaces();
  BuildLinearSystem();

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

//  -_Test_SaveDof_---------------------------------------
    SavedDiffTerm<genTensor1<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<genTensor1<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<genTensor1<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<genTensor1<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<genTensor1<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;
  GaussQuadrature Integ_Bulk ( 3 );
  for ( unsigned int i = 0; i < EDomains.size(); i++ )
  {
    membraneDomain* E = EDomains[i];

    diffTerm<genTensor1<double>,0>::Handle Field ( new genSField<genTensor1<double> >(pAssembler, FSpaceDisp) );
    genTerm<genTensor2<double>,0>::Handle G = Gradient(Field);
    genTerm<genTensor2<double>,0>::Handle CG(new genCache<genTensor2<double>,0>(G));
    genTerm<genTensor2<double>,0>::Handle StressDisp = E->Stress<0>(Field);
    genTerm<genTensor2<double>,0>::Handle CS(new genCache<genTensor2<double>,0>(StressDisp));

    genTerm<genTensor0<double>,0>::Handle Energy = E->Energy<0,0>(CS,CG);

 //   GaussQuadrature Integ_Bulk( GaussQuadrature::GradGrad );
    Assemble ( Energy, E->begin(), E->end(), Integ_Bulk, energ);

//     genTerm<genTensor4<double>,0>::ConstHandle C = E->C;
//     genTerm<double,0>::Handle BtCB = FullContractedProductCompose(Transpose(CB),FullContractedProductCompose(C,CB));
//     Assemble ( *BtCB, E->begin(), E->end(), Integ_Bulk, energ);
  }
  printf ( "elastic energy=%f\n", energ/2. );


//   for ( unsigned int i = 0; i < EDomains.size(); i++ )
//   {
//     diffTerm<genTensor1<double>,0>::Handle Field (new genSField<genTensor1<double>>(pAssembler, FSpaceDisp) );
//     SavedgenTerm<genTensor4<double>,0>* H;
//     _savedgenTermManager.getSavedgenTerm("H",H,0,0,i);
//
//     const genTerm<genTensor2<double>,0> &G=Gradient(Field);
//     const genTerm<genTensor2<double>,0> &Strain=(G+Transpose(G))*0.5;
//     const genTerm<genTensor2<double>,0> &Stress=FullContractedProductCompose(*H,Strain);
//
//     SavedgenTerm<genTensor2<double>,0>* saveStress;
//     _savedgenTermManager.addSavedgenTerm("Stress",saveStress,0,0,i);
//
//     GaussQuadrature Integ_Val ( GaussQuadrature::Grad );
//     SavePtGauss(Stress,EDomains[i]->_g->begin(), EDomains[i]->_g->end(),Integ_Val,*saveStress);
//
//     //_savedgenTermManager.print();
//     //saveStress->print("saveStress");
//   }
//   printf ( "--SaveField--\n" );
}


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

PView* membraneSolver::buildElasticEnergyView ( const std::string &postFileName )
{
  diffTerm<genTensor1<double>,0>::Handle Field(new genSField<genTensor1<double> >( pAssembler, FSpaceDisp ));
  genTerm<genTensor2<double>,0>::Handle G = Gradient(Field);
  genTerm<genTensor2<double>,0>::Handle CG(new genCache<genTensor2<double>,0>(G));

  std::vector<genTerm<genTensor0<double>,0>::ConstHandle> Energy;
  for ( unsigned int i = 0; i < EDomains.size(); i++ )
  {
    membraneDomain* E = EDomains[i];
    genTerm<genTensor2<double>,0>::Handle StressDisp = E->Stress<0>(Field);
    genTerm<genTensor2<double>,0>::Handle CS(new genCache<genTensor2<double>,0>(StressDisp));
    Energy.push_back(E->Energy<0,0>(CS,CG));
  }

  return buildViewElement(Energy, EDomains, postFileName);
}