// Nonlinear_genTerm_xfem - A solver for nonlinear problems using X-FEM
// Copyright (C) 2010-2026 Eric Bechet, Frederic Duboeuf
//
// See the LICENSE file for license information and contributions.
// Please report all bugs and problems to <bechet@cadxfem.org> or <duboeuf@outlook.com>.

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

#include "nonLinearXfemSolver.h"
#include "SpaceReducerForStableMultipliers.h"
#include "spaceReducer.h"
#include "genAlgorithms.h"
#include "genSField.h"
// #include "linearSystemCSR.h"
// #include "linearSystemPETSc.h"
#include "linearSystemId.h"
// #include "linearSystemGmm.h"

// #include "Numeric.h"
// #include "FuncGradDisc.h"
// #include "FuncHeaviside.h"

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

NonLinearXfemSolver::~NonLinearXfemSolver()
{
}

void NonLinearXfemSolver::CheckProperties()
{
  NonLinearSolver::CheckProperties();

  double eqfac = Data->User.Scalars["ScaleElastic"];

  for (int i=0;i<Data->Domains.size();++i)
  {
    if (Data->Domains[i]->type=="PerfectInterface")
    {
//       std::vector<genSupport*> Domains;
//       Domains.push_back(EDomains[Data->Domains[i]->Integers["DomainInt"]]);
//       Domains.push_back(EDomains[Data->Domains[i]->Integers["DomainExt"]]);
// 
// //       BitriangularDamageInterfaceDomain* field = new BitriangularDamageInterfaceDomain(*(Data->Domains[i]), Domains, eqfac);
//       genPerfectInterfaceDomain* field = new genPerfectInterfaceDomain(*(Data->Domains[i]), Domains, eqfac);
      genPerfectInterfaceDomain* field = new genPerfectInterfaceDomain(*(Data->Domains[i]), eqfac);
      InterfaceDomains.push_back (field);
    }
    if (Data->Domains[i]->type=="BitriangularDamageInterface")
    {
//       std::vector<genSupport*> Domains;
//       Domains.push_back(EDomains[Data->Domains[i]->Integers["DomainInt"]]);
//       Domains.push_back(EDomains[Data->Domains[i]->Integers["DomainExt"]]);
//
// //       BitriangularDamageInterfaceDomain* field = new BitriangularDamageInterfaceDomain(*(Data->Domains[i]), Domains, eqfac);
//       genBitriangularDamageInterfaceDomain* field = new genBitriangularDamageInterfaceDomain(*(Data->Domains[i]), Domains, eqfac);
      genBitriangularDamageInterfaceDomain* field = new genBitriangularDamageInterfaceDomain(*(Data->Domains[i]), eqfac);
      InterfaceDomains.push_back (field);
    }
    else if (Data->Domains[i]->type=="PolynomialDamageInterface")
    {
// //      PolynomialDamageDomain* field = new PolynomialDamageDomain(*(Data->Domains[i]), DDomains, eqfac);
//      PolynomialDamageDomain* field = new PolynomialDamageDomain(*(Data->Domains[i]), eqfac);
//      InterfaceDomains.push_back (field);
    }
  }

  // Combines a boundary condition at a domain
  for (unsigned int j=0; j<Data->BCs.size(); ++j)
  {
    if (Data->BCs[j]->kind=="Lagrange")
    {
      if (Data->BCs[j]->What=="Displacement")
      {
        std::vector<genGroupOfElements*> contain;
        for (unsigned int i=0; i<EDomains.size(); ++i)
        {
          genFilterElementGroupVertex filterElt(EDomains[i]->group());
          genGroupOfElements* SF = new genGroupOfElements(Data->BCs[j]->group(), filterElt);
          contain.push_back(SF);
        }
        CombineSLMatDomain.push_back(contain);

        StableLagrangeMultiplierDomain* field = new StableLagrangeMultiplierDomain(Data->BCs[j], eqfac);
        StableLagMultDomains.push_back(field);
      }
    }
  }
}

void NonLinearXfemSolver::readInputFile ( const std::string &fileName )
{
  NonLinearSolver::readInputFile(fileName);
  printf("--> %d Interface Domains\n", (int)InterfaceDomains.size());
}



void NonLinearXfemSolver::BuildDirichlet(double facDirichlet)
{
  NonLinearSolver::BuildDirichlet(facDirichlet);

  int k=0;
  for ( unsigned int i = 0; i < Data->BCs.size(); ++i )
  {
    if (Data->BCs[i]->kind=="Lagrange")
    {
      if (Data->BCs[i]->What=="Displacement")
      {
        StableLagrangeMultiplierDomain* S = StableLagMultDomains[k];
        double fval[3];
        fval[0] = fval[1] = fval[2] = 0;
        fval[FinalBCs[i]->comp1] = (*FinalBCs[i]->fscalar)(0,0,0)*facDirichlet;
        S->f = new genSimpleFunctionConstant<genTensor1<double> >(genTensor1<double>(fval[0], fval[1], fval[2]));
        ++k;
      }
    }
//     if (Data->BCs[i]->kind=="Lagrange_GaetanMethod")
//     {
//       if (Data->BCs[i]->What=="Displacement")
//       {
// //        delete Data->BCs[i];
// //        Data->BCs[i] = new genBoundaryCondition(*FinalBCs[i]);
//         Data->BCs[i]->fscalar = new genSimpleFunctionOperatorProduct<double> (facDirichlet, FinalBCs[i]->fscalar );
//       }
//     }
  }
}

void NonLinearXfemSolver::BuildNeumann(double facNeumann)
{
  NonLinearSolver::BuildNeumann(facNeumann);
}


void NonLinearXfemSolver::BuildInitialFunctionSpaces()
{
  NonLinearSolver::BuildInitialFunctionSpaces();

  FSpaceStableLagMult.resize(StableLagMultDomains.size());
  for ( unsigned int i=0 ; i < FSpaceStableLagMult.size(); ++i)
  {
    VectorLagrangeFSpace<double>::Along vectLag[3];
    vectLag[0] = VectorLagrangeFSpace<double>::VECTOR_X;
    vectLag[1] = VectorLagrangeFSpace<double>::VECTOR_Y;
    vectLag[2] = VectorLagrangeFSpace<double>::VECTOR_Z;

    std::vector<int> comps(StableLagMultDomains[i]->comps.begin(), StableLagMultDomains[i]->comps.end());
    printf("StableLagMultDomains[%d] : %d components\n", i, (int)comps.size());

    switch (comps.size())
    {
      case 1 :
        FSpaceStableLagMult[i] = genFSpace<genTensor1<double> >::Handle(new VectorLagrangeFSpace<double>(vectLag[comps[0]]));
        break;
      case 2 :
        FSpaceStableLagMult[i] = genFSpace<genTensor1<double> >::Handle(new VectorLagrangeFSpace<double>(vectLag[comps[0]], vectLag[comps[1]]));
        break;
      case 3 :
        FSpaceStableLagMult[i] = genFSpace<genTensor1<double> >::Handle(new VectorLagrangeFSpace<double>());
        break;
      default :
        Msg::Error("StableLagrangeMultipliersSpace size error");
    }
  }

  FSpaceStableLagMultForInterface.resize(InterfaceDomains.size());
  for ( unsigned int i=0 ; i < FSpaceStableLagMultForInterface.size(); ++i)
  {
    switch (Data->dim)
    {
      case 1 :
        FSpaceStableLagMultForInterface[i] = genFSpace<genTensor1<double> >::Handle(new VectorLagrangeFSpace<double>(VectorLagrangeFSpace<double>::VECTOR_X));
        break;
      case 2 :
        FSpaceStableLagMultForInterface[i] = genFSpace<genTensor1<double> >::Handle(new VectorLagrangeFSpace<double>(VectorLagrangeFSpace<double>::VECTOR_X, VectorLagrangeFSpace<double>::VECTOR_Y));
        break;
      case 3 :
        FSpaceStableLagMultForInterface[i] = genFSpace<genTensor1<double> >::Handle(new VectorLagrangeFSpace<double>());
        break;
    }
  }


  for ( unsigned int i = 0; i < FSpaceStableLagMult.size(); ++i)
  {
#if 0
    std::vector<int> comps(StableLagMultDomains[i]->comps.begin(), StableLagMultDomains[i]->comps.end());
    SpaceReducerForStableMultipliers* spacereducer = new SpaceReducerForStableMultipliers(StableLagMultDomains[i]->group(), Data->tag+EDomains.size()+i);
    spacereducer->BuildLinearConstraints(comps, pAssembler);
#else
    SpaceReducer spaceReducer(StableLagMultDomains[i]->group());
    spaceReducer.BuildLinearConstraints(FSpaceStableLagMult[i], pAssembler);
#endif
  }

  for ( unsigned int i = 0; i < FSpaceStableLagMultForInterface.size(); ++i)
  {
#if 0
    std::vector<int> comps(InterfaceDomains[i]->comps.begin(), InterfaceDomains[i]->comps.end());
    SpaceReducerForStableMultipliers* spacereducer = new SpaceReducerForStableMultipliers(InterfaceDomains[i]->group(), Data->tag+EDomains.size()+StableLagMultDomains.size()+i);
    spacereducer->BuildLinearConstraints(comps, pAssembler);
#else
    SpaceReducer spaceReducer(InterfaceDomains[i]->group());
    spaceReducer.BuildLinearConstraints(FSpaceStableLagMultForInterface[i], pAssembler);
#endif
  }

//   for ( unsigned int i = 0; i < Data->BCs.size(); ++i)
//   {
//     if (Data->BCs[i]->kind=="Lagrange_GaetanMethod")
//     {
//        FSpaceLagMult = genFSpace<double>::Handle(new ScalarLagrangeFSpace<genTensor0<double> >()); // here we take parent vertex for dof...
// //      FSpaceLagMult = genFSpace<double>::Handle(new ScalarLagrangeFSpaceOfElement()); // here level set vertices ?
//     }
//   }

// we number the LagrangeMultipliers dofs
  for ( unsigned int i = 0; i < FSpaceStableLagMult.size(); ++i)
  {
    NumberDofs ( FSpaceStableLagMult[i], StableLagMultDomains[i]->begin(), StableLagMultDomains[i]->end(), *pAssembler );
  }
  for ( unsigned int i = 0; i < FSpaceStableLagMultForInterface.size(); ++i)
  {
    NumberDofs ( FSpaceStableLagMultForInterface[i], InterfaceDomains[i]->begin(), InterfaceDomains[i]->end(), *pAssembler );
  }
//   for (unsigned int i = 0; i<Data->BCs.size(); ++i)
//   {
//     if (Data->BCs[i]->kind=="Lagrange_GaetanMethod") // Gaetan method
//       NumberDofs ( FSpaceLagMult, Data->BCs[i]->begin(), Data->BCs[i]->end(), *pAssembler );
//   }

  DispJumps.resize(InterfaceDomains.size());
  EpsilonDispJumps.resize(InterfaceDomains.size());
  for ( unsigned int i = 0; i < InterfaceDomains.size(); ++i)
  {
    genTerm<genTensor1<double>,1>::Handle DispExt ( FSpaceDisps[InterfaceDomains[i]->DomainExt] );
    genTerm<genTensor1<double>,1>::Handle DispInt ( FSpaceDisps[InterfaceDomains[i]->DomainInt] );
    DispJumps[i] = DispExt + (DispInt*-1.);

    genTerm<genTensor2<double>,1>::Handle EpsilonDispExt ( Gradient( FSpaceDisps[InterfaceDomains[i]->DomainExt] ) );
    genTerm<genTensor2<double>,1>::Handle EpsilonDispInt ( Gradient( FSpaceDisps[InterfaceDomains[i]->DomainInt] ) );
    EpsilonDispJumps[i] = EpsilonDispExt + (EpsilonDispInt*-1.);
  }
}

void NonLinearXfemSolver::AssembleInitialRHS()
{
  NonLinearSolver::AssembleInitialRHS();
}

void NonLinearXfemSolver::AssembleInitialLHS()
{
  NonLinearSolver::AssembleInitialLHS();

  for (int i=0; i<EDomains.size(); ++i)
    for ( unsigned int j = 0; j < StableLagMultDomains.size(); ++j)
      if (CombineSLMatDomain[j][i]->size()!=0)
      {
        genGroupOfElements* SF = CombineSLMatDomain[j][i];
        StableLagrangeMultiplierDomain* S = StableLagMultDomains[j];
        genTerm<genTensor0<double>,2>::Handle CrossLagTerm1 = S->CrossLagTerm<1,1>(FSpaceDisps[i],FSpaceStableLagMult[j]);//pas de saut ici
        genTerm<genTensor0<double>,2>::Handle CrossLagTerm2 = S->CrossLagTerm<1,1>(FSpaceStableLagMult[j],FSpaceDisps[i]);
        genTerm<genTensor0<double>,2>::Handle SelfLagTerm = S->SelfLagTerm<1,1>(FSpaceStableLagMult[j],FSpaceStableLagMult[j]);

        GaussQuadrature Integ_LagrangeMult ( GaussQuadrature::ValVal );
        Assemble ( CrossLagTerm1, SF->begin(), SF->end(), Integ_LagrangeMult, *pAssembler );
        Assemble ( CrossLagTerm2, SF->begin(), SF->end(), Integ_LagrangeMult, *pAssembler );
        Assemble (  SelfLagTerm , SF->begin(), SF->end(), Integ_LagrangeMult, *pAssembler );
      }

  for ( unsigned int i = 0; i < InterfaceDomains.size(); ++i)
  {
    genInterfaceDomain<genTensor0<double> >* S = InterfaceDomains[i];
//     InterfaceDomain* S = InterfaceDomains[i];
    genTerm<genTensor0<double>,2>::Handle CrossLagTerm1 = S->CrossLagTerm<1,1>(DispJumps[i],FSpaceStableLagMultForInterface[i]);//WARNING probleme certainement dans la numerotation des dof associé -> combinaison linéaire ???
    genTerm<genTensor0<double>,2>::Handle CrossLagTerm2 = S->CrossLagTerm<1,1>(FSpaceStableLagMultForInterface[i],DispJumps[i]);
    genTerm<genTensor0<double>,2>::Handle SelfLagTerm = S->SelfLagTerm<1,1>(FSpaceStableLagMultForInterface[i],FSpaceStableLagMultForInterface[i]);

    GaussQuadrature Integ_LagrangeMult ( GaussQuadrature::ValVal );
    Assemble ( CrossLagTerm1*(-1), S->begin(), S->end(), Integ_LagrangeMult, *pAssembler );
    Assemble ( CrossLagTerm2*(-1), S->begin(), S->end(), Integ_LagrangeMult, *pAssembler );
    Assemble (  SelfLagTerm*(-1) , S->begin(), S->end(), Integ_LagrangeMult, *pAssembler );
  }

//   for (int i=0; i<EDomains.size(); ++i)
//     for (unsigned int j = 0; j<Data->BCs.size(); ++j)
//       if (Data->BCs[j]->kind=="Lagrange_GaetanMethod") // Gaetan method
//         if (CombineBCatDomain[j][i].size()!=0)
//         {
//           const genSimpleFunction<genTensor1<double> > &fvector = *(Data->BCs[j]->fvector);
//           const genSimpleFunction<double> &fscalar = *(Data->BCs[j]->fscalar);
//           genTensor1<double> v = fvector(0.,0.,0.);
//           double tau = fscalar(0.,0.,0.);
//           genTerm<genTensor1<double>,0>::Handle d(new ConstantField<genTensor1<double> >(v)); // genTerm<genTensor1<double>,0>::Handle d(new FunctionField<genTensor1<double> >(*(Data->BCs[i]->fvector)));
// 
//           genTerm<genTensor1<double>,1>::Handle G = Gradient(FSpaceLagMult);
//           genTerm<genTensor1<double>,1>::Handle GC(new genCache<genTensor1<double>,1>(G));
// 
//           genTerm<double,2>::Handle CrossLagTerm1 = Compose<TensProdOp>(Compose<FullContrProdOp>(FSpaceDisps[i],d),FSpaceLagMult);
//           genTerm<double,2>::Handle CrossLagTerm2 = Compose<TensProdOp>(FSpaceLagMult,Compose<FullContrProdOp>(FSpaceDisps[i],d));
//           genTerm<double,2>::Handle LapTerm = Compose<FullContrProdOp>(GC,GC)*tau;
// 
//           GaussQuadrature Integ_LagrangeMult ( GaussQuadrature::ValVal );
//           GaussQuadrature Integ_Laplace ( GaussQuadrature::GradGrad );
//           Assemble ( CrossLagTerm1, Data->BCs[j]->begin(), Data->BCs[j]->end(), Integ_LagrangeMult, *pAssembler );
//           Assemble ( CrossLagTerm2, Data->BCs[j]->begin(), Data->BCs[j]->end(), Integ_LagrangeMult, *pAssembler );
//           Assemble ( LapTerm, Data->BCs[j]->begin(), Data->BCs[j]->end(), Integ_Laplace, *pAssembler );
//         }

/*  NonLinearSolver::AssembleInitialLHS();

  for ( unsigned int i = 0; i < StableLagMultDomains.size(); ++i)
  {
    StableLagrangeMultiplierDomain* S = StableLagMultDomains[i];
    genTerm<double,2>::Handle CrossLagTerm1 = S->CrossLagTerm<1,1>(FSpaceDisp,FSpaceStableLagMult[i]);//pas de saut ici
    genTerm<double,2>::Handle CrossLagTerm2 = S->CrossLagTerm<1,1>(FSpaceStableLagMult[i],FSpaceDisp);
    genTerm<double,2>::Handle SelfLagTerm = S->SelfLagTerm<1,1>(FSpaceStableLagMult[i],FSpaceStableLagMult[i]);

    GaussQuadrature Integ_LagrangeMult ( GaussQuadrature::ValVal );
    Assemble ( CrossLagTerm1, S->begin(), S->end(), Integ_LagrangeMult, *pAssembler );
    Assemble ( CrossLagTerm2, S->begin(), S->end(), Integ_LagrangeMult, *pAssembler );
    Assemble (  SelfLagTerm , S->begin(), S->end(), Integ_LagrangeMult, *pAssembler );
  }
  for ( unsigned int i = 0; i < InterfaceDomains.size(); ++i)
  {
    InterfaceDomain* S = InterfaceDomains[i];
    genTerm<double,2>::Handle CrossLagTerm1 = S->CrossLagTerm<1,1>(DispJumps[i],FSpaceStableLagMultForInterface[i]);//WARNING probleme certainement dans la numerotation des dof associé -> combinaison linéaire ???
    genTerm<double,2>::Handle CrossLagTerm2 = S->CrossLagTerm<1,1>(FSpaceStableLagMultForInterface[i],DispJumps[i]);
    genTerm<double,2>::Handle SelfLagTerm = S->SelfLagTerm<1,1>(FSpaceStableLagMultForInterface[i],FSpaceStableLagMultForInterface[i]);

    GaussQuadrature Integ_LagrangeMult ( GaussQuadrature::ValVal );
    Assemble ( CrossLagTerm1, S->begin(), S->end(), Integ_LagrangeMult, *pAssembler );
    Assemble ( CrossLagTerm2, S->begin(), S->end(), Integ_LagrangeMult, *pAssembler );
    Assemble (  SelfLagTerm , S->begin(), S->end(), Integ_LagrangeMult, *pAssembler );
  }

  for (unsigned int i = 0; i<Data->BCs.size(); ++i)
  {
    if (Data->BCs[i]->kind=="Lagrange_GaetanMethod") // Gaetan method
    {

      const genSimpleFunction<genTensor1<double> > &fvector = *(Data->BCs[i]->fvector);
      const genSimpleFunction<double> &fscalar = *(Data->BCs[i]->fscalar);
      genTensor1<double> v = fvector(0.,0.,0.);
      double tau = fscalar(0.,0.,0.);
      genTerm<genTensor1<double>,0>::Handle d(new ConstantField<genTensor1<double> >(v)); // genTerm<genTensor1<double>,0>::Handle d(new FunctionField<genTensor1<double> >(*(Data->BCs[i]->fvector)));

      genTerm<genTensor1<double>,1>::Handle G = Gradient(FSpaceLagMult);
      genTerm<genTensor1<double>,1>::Handle GC(new genCache<genTensor1<double>,1>(G));

      genTerm<double,2>::Handle CrossLagTerm1 = Compose<TensProdOp>(Compose<FullContrProdOp>(FSpaceDisp,d),FSpaceLagMult);
      genTerm<double,2>::Handle CrossLagTerm2 = Compose<TensProdOp>(FSpaceLagMult,Compose<FullContrProdOp>(FSpaceDisp,d));
      genTerm<double,2>::Handle LapTerm = Compose<FullContrProdOp>(GC,GC)*tau;

      GaussQuadrature Integ_LagrangeMult ( GaussQuadrature::ValVal );
      GaussQuadrature Integ_Laplace ( GaussQuadrature::GradGrad );
      Assemble ( CrossLagTerm1, Data->BCs[i]->begin(), Data->BCs[i]->end(), Integ_LagrangeMult, *pAssembler );
      Assemble ( CrossLagTerm2, Data->BCs[i]->begin(), Data->BCs[i]->end(), Integ_LagrangeMult, *pAssembler );
      Assemble ( LapTerm, Data->BCs[i]->begin(), Data->BCs[i]->end(), Integ_Laplace, *pAssembler );
    }
  }*/
}



void NonLinearXfemSolver::BuildFunctionSpaces()
{
  NonLinearSolver::BuildFunctionSpaces();

  FSpaceStableLagMult.resize(StableLagMultDomains.size());
  for ( unsigned int i=0 ; i < FSpaceStableLagMult.size(); ++i)
  {
    VectorLagrangeFSpace<double>::Along vectLag[3];
    vectLag[0] = VectorLagrangeFSpace<double>::VECTOR_X;
    vectLag[1] = VectorLagrangeFSpace<double>::VECTOR_Y;
    vectLag[2] = VectorLagrangeFSpace<double>::VECTOR_Z;

    std::vector<int> comps(StableLagMultDomains[i]->comps.begin(), StableLagMultDomains[i]->comps.end());
    printf("StableLagMultDomains[%d] : %d components\n", i, (int)comps.size());

    switch (comps.size())
    {
      case 1 :
        FSpaceStableLagMult[i] = genFSpace<genTensor1<double> >::Handle(new VectorLagrangeFSpace<double>(vectLag[comps[0]]));
        break;
      case 2 :
        FSpaceStableLagMult[i] = genFSpace<genTensor1<double> >::Handle(new VectorLagrangeFSpace<double>(vectLag[comps[0]], vectLag[comps[1]]));
        break;
      case 3 :
        FSpaceStableLagMult[i] = genFSpace<genTensor1<double> >::Handle(new VectorLagrangeFSpace<double>());
        break;
      default :
        Msg::Error("StableLagrangeMultipliersSpace size error");
    }
  }

  FSpaceStableLagMultForInterface.resize(InterfaceDomains.size());
  for ( unsigned int i=0 ; i < FSpaceStableLagMultForInterface.size(); ++i)
  {
    switch (Data->dim)
    {
      case 1 :
        FSpaceStableLagMultForInterface[i] = genFSpace<genTensor1<double> >::Handle(new VectorLagrangeFSpace<double>(VectorLagrangeFSpace<double>::VECTOR_X));
        break;
      case 2 :
        FSpaceStableLagMultForInterface[i] = genFSpace<genTensor1<double> >::Handle(new VectorLagrangeFSpace<double>(VectorLagrangeFSpace<double>::VECTOR_X, VectorLagrangeFSpace<double>::VECTOR_Y));
        break;
      case 3 :
        FSpaceStableLagMultForInterface[i] = genFSpace<genTensor1<double> >::Handle(new VectorLagrangeFSpace<double>());
        break;
    }
  }

  for ( unsigned int i = 0; i < FSpaceStableLagMult.size(); ++i)
  {
#if 0
    std::vector<int> comps(StableLagMultDomains[i]->comps.begin(), StableLagMultDomains[i]->comps.end());
    SpaceReducerForStableMultipliers* spacereducer = new SpaceReducerForStableMultipliers(StableLagMultDomains[i]->group(), Data->tag+EDomains.size()+i);
    spacereducer->BuildLinearConstraints(comps, pAssembler);
#else
    SpaceReducer spaceReducer(StableLagMultDomains[i]->group());
    spaceReducer.BuildLinearConstraints(FSpaceStableLagMult[i], pAssembler);
#endif
  }

  for ( unsigned int i = 0; i < FSpaceStableLagMultForInterface.size(); ++i)
  {
#if 0
    std::vector<int> comps(InterfaceDomains[i]->comps.begin(), InterfaceDomains[i]->comps.end());
    SpaceReducerForStableMultipliers* spacereducer = new SpaceReducerForStableMultipliers(InterfaceDomains[i]->group(), Data->tag+EDomains.size()+StableLagMultDomains.size()+i);
    spacereducer->BuildLinearConstraints(comps, pAssembler);
#else
    SpaceReducer spaceReducer(InterfaceDomains[i]->group());
    spaceReducer.BuildLinearConstraints(FSpaceStableLagMultForInterface[i], pAssembler);
#endif
  }

//   for ( unsigned int i = 0; i < Data->BCs.size(); ++i)
//   {
//     if (Data->BCs[i]->kind=="Lagrange_GaetanMethod")
//     {
//        FSpaceLagMult = genFSpace<double>::Handle(new ScalarLagrangeFSpace<genTensor0<double> >()); // here we take parent vertex for dof...
// //      FSpaceLagMult = genFSpace<double>::Handle(new ScalarLagrangeFSpaceOfElement()); // here level set vertices ?
//     }
//   }

// we number the LagrangeMultipliers dofs
  for ( unsigned int i = 0; i < FSpaceStableLagMult.size(); ++i)
  {
    NumberDofs ( FSpaceStableLagMult[i], StableLagMultDomains[i]->begin(), StableLagMultDomains[i]->end(), *pAssembler );
  }
  for ( unsigned int i = 0; i < FSpaceStableLagMultForInterface.size(); ++i)
  {
    NumberDofs ( FSpaceStableLagMultForInterface[i], InterfaceDomains[i]->begin(), InterfaceDomains[i]->end(), *pAssembler );
  }
//   for (unsigned int i = 0; i<Data->BCs.size(); ++i)
//   {
//     if (Data->BCs[i]->kind=="Lagrange_GaetanMethod") // Gaetan method
//       NumberDofs ( FSpaceLagMult, Data->BCs[i]->begin(), Data->BCs[i]->end(), *pAssembler );
//   }


  DispJumps.resize(InterfaceDomains.size());
  EpsilonDispJumps.resize(InterfaceDomains.size());
  for ( unsigned int i = 0; i < InterfaceDomains.size(); ++i)
  {
    genTerm<genTensor1<double>,1>::Handle DispExt ( FSpaceDisps[InterfaceDomains[i]->DomainExt] );
    genTerm<genTensor1<double>,1>::Handle DispInt ( FSpaceDisps[InterfaceDomains[i]->DomainInt] );
    DispJumps[i] = DispExt + (DispInt*-1.);

    genTerm<genTensor2<double>,1>::Handle EpsilonDispExt ( Gradient( FSpaceDisps[InterfaceDomains[i]->DomainExt] ) );
    genTerm<genTensor2<double>,1>::Handle EpsilonDispInt ( Gradient( FSpaceDisps[InterfaceDomains[i]->DomainInt] ) );
    EpsilonDispJumps[i] = EpsilonDispExt + (EpsilonDispInt*-1.);
  }
}

void NonLinearXfemSolver::AssembleRHS()
{
  NonLinearSolver::AssembleRHS();

  GaussQuadrature Integ_Boundary ( GaussQuadrature::Val );
  for (int i=0; i<EDomains.size(); ++i)
    for ( unsigned int j = 0; j < StableLagMultDomains.size(); ++j)
      if (CombineSLMatDomain[j][i]->size()!=0)
      {
        genGroupOfElements* SF = CombineSLMatDomain[j][i];
        StableLagrangeMultiplierDomain* S = StableLagMultDomains[j];

    //     std::shared_ptr<SavedGenTerm<genTensor1<double>,0> > OldL;
    //     savedGenTermManager.getSavedGenTerm("Lagrangian",OldL,_step,0,i);
    //     genTerm<double,1>::Handle CrossLagTerm1 = S->CrossLagTerm<1,0>(FSpaceDisps[i],OldL);
    //     Assemble ( CrossLagTerm1*(-1), SF->begin(), SF->end(), Integ_Boundary, *pAssembler );

        //Specifique aux conditions limites
    //     genTerm<genTensor1<double>,0>::Handle Func(new FunctionField<genTensor1<double> >(*(S->f));
    //     genTerm<double,1>::Handle loadterm = Compose<FullContrProdOp>(FSpaceStableLagMult[j],Func);
        genTerm<genTensor0<double>,1>::Handle loadterm = S->Lterm<1>(FSpaceStableLagMult[j]);
        Assemble ( loadterm, SF->begin(), SF->end(), Integ_Boundary, *pAssembler );

    //     std::shared_ptr<SavedGenTerm<genTensor1<double>,0> > OldD;
    //     savedGenTermManager.getSavedGenTerm("StableLagMultDisplacement",OldD,_step,0,i);
    //     genTerm<double,1>::Handle CrossLagTerm2 = S->CrossLagTerm<1,0>(FSpaceStableLagMult[j],OldD);
    //     Assemble ( CrossLagTerm2*(-1), SF->begin(), SF->end(), Integ_Boundary, *pAssembler );

        //Specifique a l'interface quand 1/k!=0
    //    std::shared_ptr<SavedGenTerm<genTensor1<double>,0> > OldU;
    //    savedGenTermManager.getSavedGenTerm("InterfaceLagrangian",OldU,_step,0,i);
    //    genTerm<double,1>::Handle SelfLagTerm = S->Internal(FSpaceStableLagMultForInterface[i]);
    //    Assemble ( SelfLagTerm*(-1), SF->begin(), S-F>end(), Integ_Boundary, *pAssembler );

        for (int k=0; k<(int)S->comps.size(); ++k) std::cout << "Lagrange BC Disp" << std::endl;
      }

  for ( unsigned int i = 0; i < InterfaceDomains.size(); ++i)
  {
    genInterfaceDomain<genTensor0<double> >* I = InterfaceDomains[i];
//     InterfaceDomain* I = InterfaceDomains[i];

    std::shared_ptr<SavedGenTerm<genTensor1<double>,0> > OldL;
    savedGenTermManager.getSavedGenTerm("InterfaceLagrangian",OldL,_step,0,i);
    genTerm<genTensor0<double>,1>::Handle CrossLagTerm1 = I->CrossLagTerm<1,0>(DispJumps[i],OldL);
    Assemble ( CrossLagTerm1, I->begin(), I->end(), Integ_Boundary, *pAssembler );

    //Specifique aux conditions limites
//     genTerm<genTensor1<double>,0>::Handle Func(new FunctionField<genTensor1<double> >(*(I->f));
//     genTerm<double,1>::Handle loadterm = Compose<FullContrProdOp>(FSpaceStableLagMultForInterface[i],Func);
//    genTerm<double,1>::Handle loadterm = I->Lterm<1>(FSpaceStableLagMultForInterface[i]);
//    Assemble ( loadterm, I->begin(), I->end(), Integ_Boundary, *pAssembler );

    std::shared_ptr<SavedGenTerm<genTensor1<double>,0> > OldD;
    savedGenTermManager.getSavedGenTerm("InterfaceDisplacement",OldD,_step,0,i);
    genTerm<genTensor0<double>,1>::Handle CrossLagTerm2 = I->CrossLagTerm<1,0>(FSpaceStableLagMultForInterface[i],OldD);
    Assemble ( CrossLagTerm2, I->begin(), I->end(), Integ_Boundary, *pAssembler );

    //Specifique a l'interface quand 1/k!=0
    //std::shared_ptr<SavedGenTerm<genTensor1<double>,0> > OldU;
    //savedGenTermManager.getSavedGenTerm("InterfaceLagrangian",OldU,_step,0,i);
    //genTerm<double,1>::Handle SelfLagTerm = Compose<FullContrProdOp>(FSpaceStableLagMultForInterface[i],OldU);
    genTerm<genTensor0<double>,1>::Handle SelfLagTerm = I->Internal<1>(FSpaceStableLagMultForInterface[i]);
    Assemble ( SelfLagTerm, I->begin(), I->end(), Integ_Boundary, *pAssembler );

//     std::shared_ptr<SavedGenTerm<genTensor1<double>,0> > OldL;
//     savedGenTermManager.getSavedGenTerm("InterfaceLagrangian",OldL,_step,0,i);
//     genTerm<double,1>::Handle SelfLagTerm = I->SelfLagTerm<1,0>(FSpaceStableLagMultForInterface[i],OldL);
//     Assemble ( SelfLagTerm, I->begin(), I->end(), Integ_Boundary, *pAssembler );
  }

//   for (unsigned int i=0; i<Data->BCs.size(); ++i)
//   {
//     if (Data->BCs[i]->kind=="Lagrange_GaetanMethod")
//     {
//       std::cout << "Dirichlet BC Disp" << std::endl;
//       genTerm<double,0>::Handle Func(new FunctionField<double>(*(Data->BCs[i]->fscalar)));
//       genTerm<double,1>::Handle loadterm = Compose<FullContrProdOp>(FSpaceLagMult,Func);
//       Assemble ( loadterm , Data->BCs[i]->begin(), Data->BCs[i]->end(), Integ_Boundary, *pAssembler );
//     }
//   }

  
  /*NonLinearSolver::AssembleRHS();

  GaussQuadrature Integ_Boundary ( GaussQuadrature::Val );
  for ( unsigned int i = 0; i < StableLagMultDomains.size(); ++i)
  {
    StableLagrangeMultiplierDomain* S = StableLagMultDomains[i];
    
//     std::shared_ptr<SavedGenTerm<genTensor1<double>,0> > OldL;
//     savedGenTermManager.getSavedGenTerm("Lagrangian",OldL,_step,0,i);
//     genTerm<double,1>::Handle CrossLagTerm1 = S->CrossLagTerm<1,0>(FSpaceDisp,OldL);
//     Assemble ( CrossLagTerm1*(-1), S->begin(), S->end(), Integ_Boundary, *pAssembler );
    
    //Specifique aux conditions limites
//     genTerm<genTensor1<double>,0>::Handle Func(new FunctionField<genTensor1<double> >(*(S->f));
//     genTerm<double,1>::Handle loadterm = Compose<FullContrProdOp>(FSpaceStableLagMult[i],Func);
    genTerm<double,1>::Handle loadterm = S->Lterm<1>(FSpaceStableLagMult[i]);
    Assemble ( loadterm, S->begin(), S->end(), Integ_Boundary, *pAssembler );

//     std::shared_ptr<SavedGenTerm<genTensor1<double>,0> > OldD;
//     savedGenTermManager.getSavedGenTerm("StableLagMultDisplacement",OldD,_step,0,i);
//     genTerm<double,1>::Handle CrossLagTerm2 = S->CrossLagTerm<1,0>(FSpaceStableLagMult[i],OldD);
//     Assemble ( CrossLagTerm2*(-1), S->begin(), S->end(), Integ_Boundary, *pAssembler );

    //Specifique a l'interface quand 1/k!=0
//    std::shared_ptr<SavedGenTerm<genTensor1<double>,0> > OldU;
//    savedGenTermManager.getSavedGenTerm("InterfaceLagrangian",OldU,_step,0,i);
//    genTerm<double,1>::Handle SelfLagTerm = S->Internal(FSpaceStableLagMultForInterface[i]);
//    Assemble ( SelfLagTerm*(-1), S->begin(), S->end(), Integ_Boundary, *pAssembler );
    
    for (int j=0; j<3; ++j) if (S->comp[j]) std::cout << "Lagrange BC Disp" << std::endl;
  }

  for ( unsigned int i = 0; i < InterfaceDomains.size(); ++i)
  {
    InterfaceDomain* S = InterfaceDomains[i];

    std::shared_ptr<SavedGenTerm<genTensor1<double>,0> > OldL;
    savedGenTermManager.getSavedGenTerm("InterfaceLagrangian",OldL,_step,0,i);
    genTerm<double,1>::Handle CrossLagTerm1 = S->CrossLagTerm<1,0>(FSpaceDisp,OldL);
    Assemble ( CrossLagTerm1*(-1), S->begin(), S->end(), Integ_Boundary, *pAssembler );

    //Specifique aux conditions limites
//     genTerm<genTensor1<double>,0>::Handle Func(new FunctionField<genTensor1<double> >(*(S->f));
//     genTerm<double,1>::Handle loadterm = Compose<FullContrProdOp>(FSpaceStableLagMultForInterface[i],Func);
//    genTerm<double,1>::Handle loadterm = S->Lterm<1>(FSpaceStableLagMultForInterface[i]);
//    Assemble ( loadterm, S->begin(), S->end(), Integ_Boundary, *pAssembler );

    std::shared_ptr<SavedGenTerm<genTensor1<double>,0> > OldD;
    savedGenTermManager.getSavedGenTerm("InterfaceDisplacement",OldD,_step,0,i);
    genTerm<double,1>::Handle CrossLagTerm2 = S->CrossLagTerm<1,0>(FSpaceStableLagMultForInterface[i],OldD);
    Assemble ( CrossLagTerm2*(-1), S->begin(), S->end(), Integ_Boundary, *pAssembler );

    //Specifique a l'interface quand 1/k!=0
    std::shared_ptr<SavedGenTerm<genTensor1<double>,0> > OldU;
    savedGenTermManager.getSavedGenTerm("InterfaceLagrangian",OldU,_step,0,i);
    genTerm<double,1>::Handle SelfLagTerm = S->Internal<1>(FSpaceStableLagMultForInterface[i]);
    Assemble ( SelfLagTerm*(-1), S->begin(), S->end(), Integ_Boundary, *pAssembler );
  }

  for (unsigned int i=0; i<Data->BCs.size(); ++i)
  {
    if (Data->BCs[i]->kind=="Lagrange_GaetanMethod")
    {
      std::cout << "Dirichlet BC Disp" << std::endl;
      genTerm<double,0>::Handle Func(new FunctionField<double>(*(Data->BCs[i]->fscalar)));
      genTerm<double,1>::Handle loadterm = Compose<FullContrProdOp>(FSpaceLagMult,Func);
      Assemble ( loadterm , Data->BCs[i]->begin(), Data->BCs[i]->end(), Integ_Boundary, *pAssembler );
    }
  }*/
}

void NonLinearXfemSolver::AssembleLHS()
{
  NonLinearSolver::AssembleLHS();

  for (int i=0; i<EDomains.size(); ++i)
    for ( unsigned int j = 0; j < StableLagMultDomains.size(); ++j)
      if (CombineSLMatDomain[j][i]->size()!=0)
      {
        genGroupOfElements* SF = CombineSLMatDomain[j][i];
        StableLagrangeMultiplierDomain* S = StableLagMultDomains[j];
        genTerm<genTensor0<double>,2>::Handle CrossLagTerm1 = S->CrossLagTerm<1,1>(FSpaceDisps[i],FSpaceStableLagMult[j]);
        genTerm<genTensor0<double>,2>::Handle CrossLagTerm2 = S->CrossLagTerm<1,1>(FSpaceStableLagMult[j],FSpaceDisps[i]);
        genTerm<genTensor0<double>,2>::Handle SelfLagTerm = S->SelfLagTerm<1,1>(FSpaceStableLagMult[j],FSpaceStableLagMult[j]);

        GaussQuadrature Integ_LagrangeMult ( GaussQuadrature::ValVal );
        Assemble ( CrossLagTerm1, SF->begin(), SF->end(), Integ_LagrangeMult, *pAssembler );
        Assemble ( CrossLagTerm2, SF->begin(), SF->end(), Integ_LagrangeMult, *pAssembler );
        Assemble (  SelfLagTerm , SF->begin(), SF->end(), Integ_LagrangeMult, *pAssembler );
      }

  for ( unsigned int i = 0; i < InterfaceDomains.size(); ++i)
  {
    genInterfaceDomain<genTensor0<double> >* S = InterfaceDomains[i];
//     InterfaceDomain* S = InterfaceDomains[i];
    genTerm<genTensor0<double>,2>::Handle CrossLagTerm1 = S->CrossLagTerm<1,1>(DispJumps[i],FSpaceStableLagMultForInterface[i]);
    genTerm<genTensor0<double>,2>::Handle CrossLagTerm2 = S->CrossLagTerm<1,1>(FSpaceStableLagMultForInterface[i],DispJumps[i]);
    genTerm<genTensor0<double>,2>::Handle SelfLagTerm = S->SelfLagTerm<1,1>(FSpaceStableLagMultForInterface[i],FSpaceStableLagMultForInterface[i]);

    GaussQuadrature Integ_LagrangeMult ( GaussQuadrature::ValVal );
    Assemble ( CrossLagTerm1*(-1), S->begin(), S->end(), Integ_LagrangeMult, *pAssembler );
    Assemble ( CrossLagTerm2*(-1), S->begin(), S->end(), Integ_LagrangeMult, *pAssembler );
    Assemble (  SelfLagTerm*(-1) , S->begin(), S->end(), Integ_LagrangeMult, *pAssembler );
  }

//   for (int i=0; i<EDomains.size(); ++i)
//     for (unsigned int j = 0; j<Data->BCs.size(); ++j)
//       if (Data->BCs[j]->kind=="Lagrange_GaetanMethod") // Gaetan method
//         if (CombineBCatDomain[j][i].size()!=0)
//         {
//           const genSimpleFunction<genTensor1<double> > &fvector = *(Data->BCs[j]->fvector);
//           const genSimpleFunction<double> &fscalar = *(Data->BCs[j]->fscalar);
//           genTensor1<double> v = fvector(0.,0.,0.);
//           double tau = fscalar(0.,0.,0.);
//           genTerm<genTensor1<double>,0>::Handle d(new ConstantField<genTensor1<double> >(v)); // genTerm<genTensor1<double>,0>::Handle d(new FunctionField<genTensor1<double> >(*(Data->BCs[i]->fvector)));
// 
//           genTerm<genTensor1<double>,1>::Handle G = Gradient(FSpaceLagMult);
//           genTerm<genTensor1<double>,1>::Handle GC(new genCache<genTensor1<double>,1>(G));
// 
//           genTerm<double,2>::Handle CrossLagTerm1 = Compose<TensProdOp>(Compose<FullContrProdOp>(FSpaceDisps[i],d),FSpaceLagMult);
//           genTerm<double,2>::Handle CrossLagTerm2 = Compose<TensProdOp>(FSpaceLagMult,Compose<FullContrProdOp>(FSpaceDisps[i],d));
//           genTerm<double,2>::Handle LapTerm = Compose<FullContrProdOp>(GC,GC)*tau;
// 
//           GaussQuadrature Integ_LagrangeMult ( GaussQuadrature::ValVal );
//           GaussQuadrature Integ_Laplace ( GaussQuadrature::GradGrad );
//           Assemble ( CrossLagTerm1, Data->BCs[j]->begin(), Data->BCs[j]->end(), Integ_LagrangeMult, *pAssembler );
//           Assemble ( CrossLagTerm2, Data->BCs[j]->begin(), Data->BCs[j]->end(), Integ_LagrangeMult, *pAssembler );
//           Assemble ( LapTerm, Data->BCs[j]->begin(), Data->BCs[j]->end(), Integ_Laplace, *pAssembler );
//         }
  
  /*NonLinearSolver::AssembleLHS();

  for ( unsigned int i = 0; i < StableLagMultDomains.size(); ++i)
  {
    StableLagrangeMultiplierDomain* S = StableLagMultDomains[i];
    genTerm<double,2>::Handle CrossLagTerm1 = S->CrossLagTerm<1,1>(FSpaceDisp,FSpaceStableLagMult[i]);
    genTerm<double,2>::Handle CrossLagTerm2 = S->CrossLagTerm<1,1>(FSpaceStableLagMult[i],FSpaceDisp);
    genTerm<double,2>::Handle SelfLagTerm = S->SelfLagTerm<1,1>(FSpaceStableLagMult[i],FSpaceStableLagMult[i]);

    GaussQuadrature Integ_LagrangeMult ( GaussQuadrature::ValVal );
    Assemble ( CrossLagTerm1, S->begin(), S->end(), Integ_LagrangeMult, *pAssembler );
    Assemble ( CrossLagTerm2, S->begin(), S->end(), Integ_LagrangeMult, *pAssembler );
    Assemble (  SelfLagTerm , S->begin(), S->end(), Integ_LagrangeMult, *pAssembler );
  }

  for ( unsigned int i = 0; i < InterfaceDomains.size(); ++i)
  {
    InterfaceDomain* S = InterfaceDomains[i];
    genTerm<double,2>::Handle CrossLagTerm1 = S->CrossLagTerm<1,1>(FSpaceDisp,FSpaceStableLagMultForInterface[i]);
    genTerm<double,2>::Handle CrossLagTerm2 = S->CrossLagTerm<1,1>(FSpaceStableLagMultForInterface[i],FSpaceDisp);
    genTerm<double,2>::Handle SelfLagTerm = S->SelfLagTerm<1,1>(FSpaceStableLagMultForInterface[i],FSpaceStableLagMultForInterface[i]);

    GaussQuadrature Integ_LagrangeMult ( GaussQuadrature::ValVal );
    Assemble ( CrossLagTerm1, S->begin(), S->end(), Integ_LagrangeMult, *pAssembler );
    Assemble ( CrossLagTerm2, S->begin(), S->end(), Integ_LagrangeMult, *pAssembler );
    Assemble (  SelfLagTerm , S->begin(), S->end(), Integ_LagrangeMult, *pAssembler );
  }

  for (unsigned int i = 0; i<Data->BCs.size(); ++i)
  {
    if (Data->BCs[i]->kind=="Lagrange_GaetanMethod") // Gaetan method
    {

      const genSimpleFunction<genTensor1<double> > &fvector = *(Data->BCs[i]->fvector);
      const genSimpleFunction<double> &fscalar = *(Data->BCs[i]->fscalar);
      genTensor1<double> v = fvector(0.,0.,0.);
      double tau = fscalar(0.,0.,0.);
      genTerm<genTensor1<double>,0>::Handle d(new ConstantField<genTensor1<double> >(v)); // genTerm<genTensor1<double>,0>::Handle d(new FunctionField<genTensor1<double> >(*(Data->BCs[i]->fvector)));

      genTerm<genTensor1<double>,1>::Handle G = Gradient(FSpaceLagMult);
      genTerm<genTensor1<double>,1>::Handle GC(new genCache<genTensor1<double>,1>(G));

      genTerm<double,2>::Handle CrossLagTerm1 = Compose<TensProdOp>(Compose<FullContrProdOp>(FSpaceDisp,d),FSpaceLagMult);
      genTerm<double,2>::Handle CrossLagTerm2 = Compose<TensProdOp>(FSpaceLagMult,Compose<FullContrProdOp>(FSpaceDisp,d));
      genTerm<double,2>::Handle LapTerm = Compose<FullContrProdOp>(GC,GC)*tau;

      GaussQuadrature Integ_LagrangeMult ( GaussQuadrature::ValVal );
      GaussQuadrature Integ_Laplace ( GaussQuadrature::GradGrad );
      Assemble ( CrossLagTerm1, Data->BCs[i]->begin(), Data->BCs[i]->end(), Integ_LagrangeMult, *pAssembler );
      Assemble ( CrossLagTerm2, Data->BCs[i]->begin(), Data->BCs[i]->end(), Integ_LagrangeMult, *pAssembler );
      Assemble ( LapTerm, Data->BCs[i]->begin(), Data->BCs[i]->end(), Integ_Laplace, *pAssembler );
    }
  }*/
}



void NonLinearXfemSolver::InitField()
{
  //--------------------------------------------------------------
  //_Init_GaussPoint_Values_--------------------------------------
  //--------------------------------------------------------------


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

    genTensor1<double> v(0.);
    ConstantField<genTensor1<double> > vZero ( v );
    GaussQuadrature Integ_Val ( GaussQuadrature::Val );


    // Init Displacement
    std::shared_ptr<SavedGenTerm<genTensor1<double>,0> > InitDIncr(new SavedGenTerm<genTensor1<double>,0>());
    std::shared_ptr<SavedGenTerm<genTensor1<double>,0> > InitD(new SavedGenTerm<genTensor1<double>,0>());

    SavePtGauss(vZero, D->begin(), D->end(), Integ_Val, *InitDIncr);
    SavePtGauss(vZero, D->begin(), D->end(), Integ_Val, *InitD);

    savedGenTermManager.addSavedGenTerm("DisplacementIncr",InitDIncr,0,0,i);
    savedGenTermManager.addSavedGenTerm("Displacement",InitD,0,0,i);

    printData(InitDIncr,"DisplacementIncr",i);
    printData(InitD,"Displacement",i);


    // Init StableLagMultDomains Displacement
    for ( unsigned int j = 0; j < StableLagMultDomains.size(); ++j )
      if (CombineSLMatDomain[j][i]->size()!=0)
      {
        genGroupOfElements* SF = CombineSLMatDomain[j][i];
        StableLagrangeMultiplierDomain* S = StableLagMultDomains[j];

        std::shared_ptr<SavedGenTerm<genTensor1<double>,0> > InitSLMDIncr(new SavedGenTerm<genTensor1<double>,0>());
        std::shared_ptr<SavedGenTerm<genTensor1<double>,0> > InitSLMD(new SavedGenTerm<genTensor1<double>,0>());

        SavePtGauss(vZero, SF->begin(), SF->end(), Integ_Val, *InitSLMDIncr);
        SavePtGauss(vZero, SF->begin(), SF->end(), Integ_Val, *InitSLMD);

        savedGenTermManager.addSavedGenTerm("StableLagMultDisplacementIncr",InitSLMDIncr,0,0,j);
        savedGenTermManager.addSavedGenTerm("StableLagMultDisplacement",InitSLMD,0,0,j);

        printData(InitSLMDIncr,"StableLagMultDisplacementIncr",j);
        printData(InitSLMD,"StableLagMultDisplacement",j);
      }


    // Init Lagrangian
    for ( unsigned int j = 0; j < StableLagMultDomains.size(); ++j )
      if (CombineSLMatDomain[j][i]->size()!=0)
      {
        StableLagrangeMultiplierDomain* I = StableLagMultDomains[j];

        std::shared_ptr<SavedGenTerm<genTensor1<double>,0> > InitLIncr(new SavedGenTerm<genTensor1<double>,0>());
        std::shared_ptr<SavedGenTerm<genTensor1<double>,0> > InitL(new SavedGenTerm<genTensor1<double>,0>());

        SavePtGauss(vZero, I->begin(), I->end(), Integ_Val, *InitLIncr);
        SavePtGauss(vZero, I->begin(), I->end(), Integ_Val, *InitL);

        savedGenTermManager.addSavedGenTerm("LagrangianIncr",InitLIncr,0,0,j);
        savedGenTermManager.addSavedGenTerm("Lagrangian",InitL,0,0,j);

        printData(InitLIncr,"LagrangianIncr",j);
        printData(InitL,"Lagrangian",j);
      }
  }

  for ( unsigned int i = 0; i < InterfaceDomains.size(); ++i )
  {
    genInterfaceDomain<genTensor0<double> >* I = InterfaceDomains[i];
//     InterfaceDomain* I = InterfaceDomains[i];

    genTensor1<double> v(0.);
    ConstantField<genTensor1<double> > vZero ( v );
    GaussQuadrature Integ_Val ( GaussQuadrature::Val );


    // Init InterfaceDomains Displacement
    std::shared_ptr<SavedGenTerm<genTensor1<double>,0> > InitDIncr(new SavedGenTerm<genTensor1<double>,0>());
    std::shared_ptr<SavedGenTerm<genTensor1<double>,0> > InitD(new SavedGenTerm<genTensor1<double>,0>());

    SavePtGauss(vZero, I->begin(), I->end(), Integ_Val, *InitDIncr);
    SavePtGauss(vZero, I->begin(), I->end(), Integ_Val, *InitD);

    savedGenTermManager.addSavedGenTerm("InterfaceDisplacementIncr",InitDIncr,0,0,i);
    savedGenTermManager.addSavedGenTerm("InterfaceDisplacement",InitD,0,0,i);

    printData(InitDIncr,"InterfaceDisplacementIncr",i);
    printData(InitD,"InterfaceDisplacement",i);


    // Init Interface Lagrangian
    std::shared_ptr<SavedGenTerm<genTensor1<double>,0> > InitLIncr(new SavedGenTerm<genTensor1<double>,0>());
    std::shared_ptr<SavedGenTerm<genTensor1<double>,0> > InitL(new SavedGenTerm<genTensor1<double>,0>());

    SavePtGauss(vZero, I->begin(), I->end(), Integ_Val, *InitLIncr);
    SavePtGauss(vZero, I->begin(), I->end(), Integ_Val, *InitL);

    savedGenTermManager.addSavedGenTerm("InterfaceLagrangianIncr",InitLIncr,0,0,i);
    savedGenTermManager.addSavedGenTerm("InterfaceLagrangian",InitL,0,0,i);

    printData(InitLIncr,"InterfaceLagrangianIncr",i);
    printData(InitL,"InterfaceLagrangian",i);
  }

  NonLinearSolver::InitField();

  //--------------------------------------------------------------
  //_Init_Unknown_Values_-----------------------------------------
  //--------------------------------------------------------------

  // Initialise StableLagMult, StableLagMult increment and StableLagMult increment on a step
  for ( unsigned int i = 0; i < StableLagMultDomains.size(); ++i )
  {
    std::shared_ptr<SavedUnknownTerm<genTensor1<double> > > InitSLagIncr(new SavedUnknownTerm<genTensor1<double> >(FSpaceStableLagMult[i]));
    std::shared_ptr<SavedUnknownTerm<genTensor1<double> > > InitSLagIncrS(new SavedUnknownTerm<genTensor1<double> >(FSpaceStableLagMult[i]));
    std::shared_ptr<SavedUnknownTerm<genTensor1<double> > > InitSLag(new SavedUnknownTerm<genTensor1<double> >(FSpaceStableLagMult[i]));

    SaveUnknown(*pAssembler,StableLagMultDomains[i]->begin(), StableLagMultDomains[i]->end(),*InitSLagIncr);
//     genTerm<double,1>::Handle SLagIncr = (InitSLagIncr)*0.;
//     SaveModifUnknown(*SLagIncr,StableLagMultDomains[i]->begin(), StableLagMultDomains[i]->end(),*InitSLagIncr);
    *InitSLagIncrS = *InitSLagIncr;
    *InitSLag = *InitSLagIncr;

    savedGenTermManager.addSavedGenTerm("SLagIncr",InitSLagIncr,0,0,i);
    savedGenTermManager.addSavedGenTerm("SLagIncrStep",InitSLagIncrS,0,0,i);
    savedGenTermManager.addSavedGenTerm("SLag",InitSLag,0,0,i);

    printData(InitSLagIncr,"SLagIncr",i);
    printData(InitSLagIncrS,"SLagIncrStep",i);
    printData(InitSLag,"SLag",i);
  }

  for ( unsigned int i = 0; i < InterfaceDomains.size(); ++i )
  {
    std::shared_ptr<SavedUnknownTerm<genTensor1<double> > > InitSLag(new SavedUnknownTerm<genTensor1<double> >(FSpaceStableLagMultForInterface[i]));
    std::shared_ptr<SavedUnknownTerm<genTensor1<double> > > InitSLagIncr(new SavedUnknownTerm<genTensor1<double> >(FSpaceStableLagMultForInterface[i]));
    std::shared_ptr<SavedUnknownTerm<genTensor1<double> > > InitSLagIncrS(new SavedUnknownTerm<genTensor1<double> >(FSpaceStableLagMultForInterface[i]));

    SaveUnknown(*pAssembler,InterfaceDomains[i]->begin(), InterfaceDomains[i]->end(),*InitSLag);
    *InitSLagIncr = *InitSLag;
    genTerm<genTensor0<double>,1>::Handle SLagIncr = InitSLagIncr*0.;
    SaveModifUnknown(*SLagIncr,InterfaceDomains[i]->begin(), InterfaceDomains[i]->end(),*InitSLagIncr);
    *InitSLagIncrS = *InitSLagIncr;

    savedGenTermManager.addSavedGenTerm("ILag",InitSLag,0,0,i);
    savedGenTermManager.addSavedGenTerm("ILagIncr",InitSLagIncr,0,0,i);
    savedGenTermManager.addSavedGenTerm("ILagIncrStep",InitSLagIncrS,0,0,i);

    printData(InitSLag,"ILag",i);
    printData(InitSLagIncr,"ILagIncr",i);
    printData(InitSLagIncrS,"ILagIncrStep",i);
  }
}

void NonLinearXfemSolver::CopyField()
{
  //--------------------------------------------------------------
  //_Copy_GaussPoint_Values_--------------------------------------
  //--------------------------------------------------------------

  for ( unsigned int i = 0; i < EDomains.size(); ++i )
  {
    
    // Copy Displacement
    std::shared_ptr<SavedGenTerm<genTensor1<double>,0> > CopyDIncr;
    std::shared_ptr<SavedGenTerm<genTensor1<double>,0> > CopyD;

    if(i==0) savedGenTermManager.getSavedGenTerm("DisplacementIncr",CopyDIncr,_step-1,0,i);
    else savedGenTermManager.getSavedGenTerm("DisplacementIncr",CopyDIncr,_step,1,i);
    if(i==0) savedGenTermManager.getSavedGenTerm("Displacement",CopyD,_step-1,0,i);
    else savedGenTermManager.getSavedGenTerm("Displacement",CopyD,_step,1,i);

    savedGenTermManager.addSavedGenTerm("DisplacementIncr",CopyDIncr,_step,0,i);
    savedGenTermManager.addSavedGenTerm("Displacement",CopyD,_step,0,i);

    
    // Copy StableLagMultDomains Displacement
    for ( unsigned int j = 0; j < StableLagMultDomains.size(); ++j )
      if (CombineSLMatDomain[j][i]->size()!=0)
      {
        genGroupOfElements* SF = CombineSLMatDomain[j][i];
        StableLagrangeMultiplierDomain* S = StableLagMultDomains[j];

        std::shared_ptr<SavedGenTerm<genTensor1<double>,0> > CopySLMDIncr;
        std::shared_ptr<SavedGenTerm<genTensor1<double>,0> > CopySLMD;

        if(j==0) savedGenTermManager.getSavedGenTerm("StableLagMultDisplacementIncr",CopySLMDIncr,_step-1,0,j);
        else savedGenTermManager.getSavedGenTerm("StableLagMultDisplacementIncr",CopySLMDIncr,_step,1,j);
        if(j==0) savedGenTermManager.getSavedGenTerm("StableLagMultDisplacement",CopySLMD,_step-1,0,j);
        else savedGenTermManager.getSavedGenTerm("StableLagMultDisplacement",CopySLMD,_step,1,j);

        savedGenTermManager.addSavedGenTerm("StableLagMultDisplacementIncr",CopySLMDIncr,_step,0,j);
        savedGenTermManager.addSavedGenTerm("StableLagMultDisplacement",CopySLMD,_step,0,j);
      }

    // Copy Lagrangian
    for ( unsigned int j = 0; j < StableLagMultDomains.size(); ++j )
      if (CombineSLMatDomain[j][i]->size()!=0)
      {
        genGroupOfElements* SF = CombineSLMatDomain[j][i];
        StableLagrangeMultiplierDomain* S = StableLagMultDomains[j];

        std::shared_ptr<SavedGenTerm<genTensor1<double>,0> > CopyLIncr;
        std::shared_ptr<SavedGenTerm<genTensor1<double>,0> > CopyL;

        if(j==0) savedGenTermManager.getSavedGenTerm("LagrangianIncr",CopyLIncr,_step-1,0,j);
        else savedGenTermManager.getSavedGenTerm("LagrangianIncr",CopyLIncr,_step,1,j);
        if(j==0) savedGenTermManager.getSavedGenTerm("Lagrangian",CopyL,_step-1,0,j);
        else savedGenTermManager.getSavedGenTerm("Lagrangian",CopyL,_step,1,j);

        savedGenTermManager.addSavedGenTerm("LagrangianIncr",CopyLIncr,_step,0,j);
        savedGenTermManager.addSavedGenTerm("Lagrangian",CopyL,_step,0,j);
      }
  }

  // Copy InterfaceDomains Displacement
  for ( unsigned int i = 0; i < InterfaceDomains.size(); ++i )
  {
    genInterfaceDomain<genTensor0<double> >* I = InterfaceDomains[i];
//     InterfaceDomain* I = InterfaceDomains[i];

    std::shared_ptr<SavedGenTerm<genTensor1<double>,0> > CopyIDIncr;
    std::shared_ptr<SavedGenTerm<genTensor1<double>,0> > CopyID;
    
    if(i==0) savedGenTermManager.getSavedGenTerm("InterfaceDisplacementIncr",CopyIDIncr,_step-1,0,i);
    else savedGenTermManager.getSavedGenTerm("InterfaceDisplacementIncr",CopyIDIncr,_step,1,i);
    if(i==0) savedGenTermManager.getSavedGenTerm("InterfaceDisplacement",CopyID,_step-1,0,i);
    else savedGenTermManager.getSavedGenTerm("InterfaceDisplacement",CopyID,_step,1,i);

    savedGenTermManager.addSavedGenTerm("InterfaceDisplacementIncr",CopyIDIncr,_step,0,i);
    savedGenTermManager.addSavedGenTerm("InterfaceDisplacement",CopyID,_step,0,i);
  }
  
  // Copy Interface Lagrangian
  for ( unsigned int i = 0; i < InterfaceDomains.size(); ++i )
  {
    genInterfaceDomain<genTensor0<double> >* I = InterfaceDomains[i];
//     InterfaceDomain* I = InterfaceDomains[i];

    std::shared_ptr<SavedGenTerm<genTensor1<double>,0> > CopyLIncr;
    std::shared_ptr<SavedGenTerm<genTensor1<double>,0> > CopyL;

    if(i==0) savedGenTermManager.getSavedGenTerm("InterfaceLagrangianIncr",CopyLIncr,_step-1,0,i);
    else savedGenTermManager.getSavedGenTerm("InterfaceLagrangianIncr",CopyLIncr,_step,1,i);
    if(i==0) savedGenTermManager.getSavedGenTerm("InterfaceLagrangian",CopyL,_step-1,0,i);
    else savedGenTermManager.getSavedGenTerm("InterfaceLagrangian",CopyL,_step,1,i);

    savedGenTermManager.addSavedGenTerm("InterfaceLagrangianIncr",CopyLIncr,_step,0,i);
    savedGenTermManager.addSavedGenTerm("InterfaceLagrangian",CopyL,_step,0,i);
  }


  
  NonLinearSolver::CopyField();

  //--------------------------------------------------------------
  //_Copy_Unknown_Values_-----------------------------------------
  //--------------------------------------------------------------
  
  // Copy StableLagMult, StableLagMult increment and StableLagMult increment on a step
  for ( unsigned int i = 0; i < StableLagMultDomains.size(); ++i )
  {
    std::shared_ptr<SavedUnknownTerm<genTensor1<double> > > CopySLagIncr;
    std::shared_ptr<SavedUnknownTerm<genTensor1<double> > > CopySLagIncrS(new SavedUnknownTerm<genTensor1<double> >(FSpaceStableLagMult[i]));
    std::shared_ptr<SavedUnknownTerm<genTensor1<double> > > CopySLag;

    if(i==0) savedGenTermManager.getSavedGenTerm("SLagIncr",CopySLagIncr,_step-1,0,i);
    else savedGenTermManager.getSavedGenTerm("SLagIncr",CopySLagIncr,_step,1,i);
    *CopySLagIncrS = *CopySLagIncr;
    genTerm<genTensor0<double>,1>::Handle SLagIncrS = CopySLagIncrS*0.;
    SaveModifUnknown(*SLagIncrS,StableLagMultDomains[i]->begin(), StableLagMultDomains[i]->end(),*CopySLagIncrS);
    if(i==0) savedGenTermManager.getSavedGenTerm("SLag",CopySLag,_step-1,0,i);
    else savedGenTermManager.getSavedGenTerm("SLag",CopySLag,_step,1,i);

    savedGenTermManager.addSavedGenTerm("SLagIncr",CopySLagIncr,_step,0,i);
    savedGenTermManager.addSavedGenTerm("SLagIncrStep",CopySLagIncrS,_step,0,i);
    savedGenTermManager.addSavedGenTerm("SLag",CopySLag,_step,0,i);
  }

  for ( unsigned int i = 0; i < InterfaceDomains.size(); ++i )
  {
    genInterfaceDomain<genTensor0<double> >* I = InterfaceDomains[i];
//     InterfaceDomain* I = InterfaceDomains[i];
    std::shared_ptr<SavedUnknownTerm<genTensor1<double> > > CopySLagIncr;
    std::shared_ptr<SavedUnknownTerm<genTensor1<double> > > CopySLagIncrS(new SavedUnknownTerm<genTensor1<double> >(FSpaceStableLagMultForInterface[i]));
    std::shared_ptr<SavedUnknownTerm<genTensor1<double> > > CopySLag;

    if(i==0) savedGenTermManager.getSavedGenTerm("ILagIncr",CopySLagIncr,_step-1,0,i);
    else savedGenTermManager.getSavedGenTerm("ILagIncr",CopySLagIncr,_step,1,i);
    *CopySLagIncrS = *CopySLagIncr;
    genTerm<genTensor0<double>,1>::Handle SLagIncrS = CopySLagIncrS*0.;
    SaveModifUnknown(*SLagIncrS,I->begin(), I->end(),*CopySLagIncrS);
    if(i==0) savedGenTermManager.getSavedGenTerm("ILag",CopySLag,_step-1,0,i);
    else savedGenTermManager.getSavedGenTerm("ILag",CopySLag,_step,1,i);

    savedGenTermManager.addSavedGenTerm("ILagIncr",CopySLagIncr,_step,0,i);
    savedGenTermManager.addSavedGenTerm("ILagIncrStep",CopySLagIncrS,_step,0,i);
    savedGenTermManager.addSavedGenTerm("ILag",CopySLag,_step,0,i);

    I->copyHistory();
  }

}

void NonLinearXfemSolver::SaveField()
{
  // Change InterfaceLagIncr and InterfaceDispIncr
 /* for ( unsigned int i = 0; i < InterfaceDomains.size(); ++i )
  {
    InterfaceDomain* I = InterfaceDomains[i];
    std::shared_ptr<SavedUnknownTerm<genTensor1<double> > > SaveSLagIncr(new SavedUnknownTerm<genTensor1<double> >(FSpaceStableLagMultForInterface[i]));

    std::shared_ptr<SavedUnknownTerm<genTensor1<double> > > SaveSLagIncr2(new SavedUnknownTerm<genTensor1<double> >(FSpaceStableLagMultForInterface[i]));
    SaveUnknown(*pAssembler,I->begin(), I->end(),*SaveSLagIncr2);
    *SaveSLagIncr = *SaveSLagIncr2;
    SaveModifUnknown(*(Compose<TensProdOp>(SaveSLagIncr2,I->Sign)+SaveSLagIncr2*-1.),I->begin(), I->end(),*SaveSLagIncr);
    SavedUnknownAddToSolution(*SaveSLagIncr,*pAssembler);


    std::shared_ptr<SavedUnknownTerm<genTensor1<double> > > SaveDExtIncr(new SavedUnknownTerm<genTensor1<double> >(FSpaceDisps[InterfaceDomains[i]->DomainExt]));
    std::shared_ptr<SavedUnknownTerm<genTensor1<double> > > SaveDIntIncr(new SavedUnknownTerm<genTensor1<double> >(FSpaceDisps[InterfaceDomains[i]->DomainInt]));

    std::shared_ptr<SavedUnknownTerm<genTensor1<double> > > SaveDExtIncr2(new SavedUnknownTerm<genTensor1<double> >(FSpaceDisps[InterfaceDomains[i]->DomainExt]));
    std::shared_ptr<SavedUnknownTerm<genTensor1<double> > > SaveDIntIncr2(new SavedUnknownTerm<genTensor1<double> >(FSpaceDisps[InterfaceDomains[i]->DomainInt]));
    SaveUnknown(*pAssembler,I->begin(), I->end(),*SaveDExtIncr2);
    SaveUnknown(*pAssembler,I->begin(), I->end(),*SaveDIntIncr2);
    *SaveDExtIncr = *SaveDExtIncr2;
    *SaveDIntIncr = *SaveDIntIncr2;
    SaveModifUnknown(*(Compose<TensProdOp>(SaveDExtIncr2,I->Sign)+SaveDExtIncr2*-1.),I->begin(), I->end(),*SaveDExtIncr);
    SaveModifUnknown(*(Compose<TensProdOp>(SaveDIntIncr2,I->Sign)+SaveDIntIncr2*-1.),I->begin(), I->end(),*SaveDIntIncr);
    SavedUnknownAddToSolution(*SaveDExtIncr,*pAssembler);
    SavedUnknownAddToSolution(*SaveDIntIncr,*pAssembler);
  }*/


  //--------------------------------------------------------------
  //_Save_GaussPoint_Values_--------------------------------------
  //--------------------------------------------------------------


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

    // Save Displacement
    std::shared_ptr<SavedGenTerm<genTensor1<double>,0> > SaveDIncr(new SavedGenTerm<genTensor1<double>,0>());
    std::shared_ptr<SavedGenTerm<genTensor1<double>,0> > SaveD(new SavedGenTerm<genTensor1<double>,0>());

    diffTerm<genTensor1<double>,0>::Handle DIncr ( new genSField<genTensor1<double> >(pAssembler, FSpaceDisps[i]) );
    GaussQuadrature Integ_Val ( GaussQuadrature::Val );
    SavePtGauss(*DIncr, D->begin(), D->end(), Integ_Val, *SaveDIncr);

    std::shared_ptr<SavedGenTerm<genTensor1<double>,0> > OldD;
    savedGenTermManager.getSavedGenTerm("Displacement",OldD,_step,0,i);
    genTerm<genTensor1<double>,0>::Handle NewD = OldD + OldD;
    SavePtGauss(*NewD, D->begin(), D->end(), Integ_Val, *SaveD);

    savedGenTermManager.addSavedGenTerm("DisplacementIncr",SaveDIncr,_step,0,i);
    savedGenTermManager.addSavedGenTerm("Displacement",SaveD,_step,0,i);

    printData(SaveDIncr,"DisplacementIncr",i);
    printData(SaveD,"Displacement",i);


    // Save StableLagMultDomains Displacement
    for ( unsigned int j = 0; j < StableLagMultDomains.size(); ++j )
      if (CombineSLMatDomain[j][i]->size()!=0)
      {
        genGroupOfElements* SF = CombineSLMatDomain[j][i];
        StableLagrangeMultiplierDomain* S = StableLagMultDomains[j];

        std::shared_ptr<SavedGenTerm<genTensor1<double>,0> > SaveSLMDIncr(new SavedGenTerm<genTensor1<double>,0>());
        std::shared_ptr<SavedGenTerm<genTensor1<double>,0> > SaveSLMD(new SavedGenTerm<genTensor1<double>,0>());

        diffTerm<genTensor1<double>,0>::Handle DIncr ( new genSField<genTensor1<double> >(pAssembler, FSpaceDisps[i]) );
        GaussQuadrature Integ_Val ( GaussQuadrature::Val );
        SavePtGauss(*DIncr, SF->begin(), SF->end(), Integ_Val, *SaveSLMDIncr);

        std::shared_ptr<SavedGenTerm<genTensor1<double>,0> > OldD;
        savedGenTermManager.getSavedGenTerm("StableLagMultDisplacement",OldD,_step,0,j);
        genTerm<genTensor1<double>,0>::Handle NewD = OldD + DIncr;
        SavePtGauss(*NewD, SF->begin(), SF->end(), Integ_Val, *SaveSLMD);

        savedGenTermManager.addSavedGenTerm("StableLagMultDisplacementIncr",SaveSLMDIncr,_step,0,j);
        savedGenTermManager.addSavedGenTerm("StableLagMultDisplacement",SaveSLMD,_step,0,j);

        printData(SaveSLMDIncr,"StableLagMultDisplacementIncr",j);
        printData(SaveSLMD,"StableLagMultDisplacement",j);
      }

    // Save Lagrangian
    for ( unsigned int j = 0; j < StableLagMultDomains.size(); ++j )
      if (CombineSLMatDomain[j][i]->size()!=0)
      {
        genGroupOfElements* SF = CombineSLMatDomain[j][i];
        StableLagrangeMultiplierDomain* S = StableLagMultDomains[j];

        std::shared_ptr<SavedGenTerm<genTensor1<double>,0> > SaveLIncr(new SavedGenTerm<genTensor1<double>,0>());
        std::shared_ptr<SavedGenTerm<genTensor1<double>,0> > SaveL(new SavedGenTerm<genTensor1<double>,0>());

        diffTerm<genTensor1<double>,0>::Handle LIncr ( new genSField<genTensor1<double> >(pAssembler, FSpaceStableLagMult[j]) );
        GaussQuadrature Integ_Val ( GaussQuadrature::Val );
        SavePtGauss(*LIncr, SF->begin(), SF->end(), Integ_Val, *SaveLIncr);
        std::shared_ptr<SavedGenTerm<genTensor1<double>,0> > OldL;
        savedGenTermManager.getSavedGenTerm("Lagrangian",OldL,_step,0,j);
        genTerm<genTensor1<double>,0>::Handle NewL = OldL + LIncr;
        SavePtGauss(*NewL, SF->begin(), SF->end(), Integ_Val, *SaveL);

        savedGenTermManager.addSavedGenTerm("LagrangianIncr",SaveLIncr,_step,0,j);
        savedGenTermManager.addSavedGenTerm("Lagrangian",SaveL,_step,0,j);

        printData(SaveLIncr,"LagrangianIncr",j);
        printData(SaveL,"Lagrangian",j);
      }
  }

  // Save InterfaceDomains Displacement
  for ( unsigned int i = 0; i < InterfaceDomains.size(); ++i )
  {
    genInterfaceDomain<genTensor0<double> >* I = InterfaceDomains[i];
//     InterfaceDomain* I = InterfaceDomains[i];

    std::shared_ptr<SavedGenTerm<genTensor1<double>,0> > SaveIDIncr(new SavedGenTerm<genTensor1<double>,0>());
    std::shared_ptr<SavedGenTerm<genTensor1<double>,0> > SaveID(new SavedGenTerm<genTensor1<double>,0>());

    genTerm<genTensor1<double>,0>::Handle IDExtIncr ( new genSField<genTensor1<double> >(pAssembler, FSpaceDisps[InterfaceDomains[i]->DomainExt]) );
    genTerm<genTensor1<double>,0>::Handle IDIntIncr ( new genSField<genTensor1<double> >(pAssembler, FSpaceDisps[InterfaceDomains[i]->DomainInt]) );
    genTerm<genTensor1<double>,0>::Handle IDIncr = IDExtIncr + (IDIntIncr*-1.);
//     diffTerm<genTensor1<double>,0>::Handle IDIncr ( new genSField<genTensor1<double> >(pAssembler, DispJumps[i]) );
    GaussQuadrature Integ_Val ( GaussQuadrature::Val );
    SavePtGauss(*IDIncr, I->begin(), I->end(), Integ_Val, *SaveIDIncr);
//     SavePtGauss(*Compose<TensProdOp>(DIncr,I->Sign), I->begin(), I->end(), Integ_Val, *SaveIDIncr);

    std::shared_ptr<SavedGenTerm<genTensor1<double>,0> > OldID;
    savedGenTermManager.getSavedGenTerm("InterfaceDisplacement",OldID,_step,0,i);
    genTerm<genTensor1<double>,0>::Handle NewID = OldID + IDIncr;
    SavePtGauss(*NewID, I->begin(), I->end(), Integ_Val, *SaveID);

    savedGenTermManager.addSavedGenTerm("InterfaceDisplacementIncr",SaveIDIncr,_step,0,i);
    savedGenTermManager.addSavedGenTerm("InterfaceDisplacement",SaveID,_step,0,i);

    printData(SaveIDIncr,"InterfaceDisplacementIncr",i);
    printData(SaveID,"InterfaceDisplacement",i);

    
    // Save Interface Lagrangian
    std::shared_ptr<SavedGenTerm<genTensor1<double>,0> > SaveLIncr(new SavedGenTerm<genTensor1<double>,0>());
    std::shared_ptr<SavedGenTerm<genTensor1<double>,0> > SaveL(new SavedGenTerm<genTensor1<double>,0>());

    diffTerm<genTensor1<double>,0>::Handle LIncr ( new genSField<genTensor1<double> >(pAssembler, FSpaceStableLagMultForInterface[i]) );
    SavePtGauss(*LIncr, I->begin(), I->end(), Integ_Val, *SaveLIncr);
//     SavePtGauss(*Compose<TensProdOp>(LIncr,I->Sign), I->begin(), I->end(), Integ_Val, *SaveLIncr);
    
    std::shared_ptr<SavedGenTerm<genTensor1<double>,0> > OldL;
    savedGenTermManager.getSavedGenTerm("InterfaceLagrangian",OldL,_step,0,i);
    genTerm<genTensor1<double>,0>::Handle NewL = OldL + LIncr;
    SavePtGauss(*NewL, I->begin(), I->end(), Integ_Val, *SaveL);

    savedGenTermManager.addSavedGenTerm("InterfaceLagrangianIncr",SaveLIncr,_step,0,i);
    savedGenTermManager.addSavedGenTerm("InterfaceLagrangian",SaveL,_step,0,i);

    printData(SaveLIncr,"InterfaceLagrangianIncr",i);
    printData(SaveL,"InterfaceLagrangian",i);
    
    I->saveHistory(SaveIDIncr,SaveID,SaveLIncr,SaveL);
  }

  
  NonLinearSolver::SaveField();

  
  //--------------------------------------------------------------
  //_Save_Unknown_Values_-----------------------------------------
  //--------------------------------------------------------------

  // Save StableLagMult, StableLagMult increment and StableLagMult increment on a step
  for ( unsigned int i = 0; i < StableLagMultDomains.size(); ++i )
  {
    std::shared_ptr<SavedUnknownTerm<genTensor1<double> > > SaveSLagIncr(new SavedUnknownTerm<genTensor1<double> >(FSpaceStableLagMult[i]));
    std::shared_ptr<SavedUnknownTerm<genTensor1<double> > > SaveSLagIncrS(new SavedUnknownTerm<genTensor1<double> >(FSpaceStableLagMult[i]));
    std::shared_ptr<SavedUnknownTerm<genTensor1<double> > > SaveSLag(new SavedUnknownTerm<genTensor1<double> >(FSpaceStableLagMult[i]));

    SaveUnknown(*pAssembler,StableLagMultDomains[i]->begin(), StableLagMultDomains[i]->end(),*SaveSLagIncr);
    std::shared_ptr<SavedUnknownTerm<genTensor1<double> > > OldSLagIncrS;
    savedGenTermManager.getSavedGenTerm("SLagIncrStep",OldSLagIncrS,_step,0,i);
    *SaveSLagIncrS = *OldSLagIncrS;
    genTerm<genTensor0<double>,1>::Handle SLagIncrS = OldSLagIncrS+SaveSLagIncr;
    SaveModifUnknown(*SLagIncrS,StableLagMultDomains[i]->begin(), StableLagMultDomains[i]->end(),*SaveSLagIncrS);
    std::shared_ptr<SavedUnknownTerm<genTensor1<double> > > OldSLag;
    savedGenTermManager.getSavedGenTerm("SLag",OldSLag,_step,0,i);
    SavedUnknownAddToSolution(*OldSLag,*pAssembler);
    SaveUnknown(*pAssembler, StableLagMultDomains[i]->begin(), StableLagMultDomains[i]->end(), *SaveSLag);

    savedGenTermManager.addSavedGenTerm("SLagIncr",SaveSLagIncr,_step,0,i);
    savedGenTermManager.addSavedGenTerm("SLagIncrStep",SaveSLagIncrS,_step,0,i);
    savedGenTermManager.addSavedGenTerm("SLag",SaveSLag,_step,0,i);

    printData(SaveSLagIncr,"SLagIncr",i);
    printData(SaveSLagIncrS,"SLagIncrStep",i);
    printData(SaveSLag,"SLag",i);
  }

  for ( unsigned int i = 0; i < InterfaceDomains.size(); ++i )
  {
    genInterfaceDomain<genTensor0<double> >* I = InterfaceDomains[i];
//     InterfaceDomain* I = InterfaceDomains[i];
    std::shared_ptr<SavedUnknownTerm<genTensor1<double> > > SaveSLagIncr(new SavedUnknownTerm<genTensor1<double> >(FSpaceStableLagMultForInterface[i]));
    std::shared_ptr<SavedUnknownTerm<genTensor1<double> > > SaveSLagIncrS(new SavedUnknownTerm<genTensor1<double> >(FSpaceStableLagMultForInterface[i]));
    std::shared_ptr<SavedUnknownTerm<genTensor1<double> > > SaveSLag(new SavedUnknownTerm<genTensor1<double> >(FSpaceStableLagMultForInterface[i]));

    SaveUnknown(*pAssembler,I->begin(), I->end(),*SaveSLagIncr);
//     std::shared_ptr<SavedUnknownTerm<genTensor1<double> > > SaveSLagIncr2(new SavedUnknownTerm<genTensor1<double> >(FSpaceStableLagMultForInterface[i]));
//     SaveUnknown(*pAssembler,I->begin(), I->end(),*SaveSLagIncr2);
//     *SaveSLagIncr = *SaveSLagIncr2;
//     SaveModifUnknown(*Compose<TensProdOp>(SaveSLagIncr2,I->Sign),I->begin(), I->end(),*SaveSLagIncr);

    std::shared_ptr<SavedUnknownTerm<genTensor1<double> > > OldSLagIncrS;
    savedGenTermManager.getSavedGenTerm("ILagIncrStep",OldSLagIncrS,_step,0,i);
    *SaveSLagIncrS = *OldSLagIncrS;
    genTerm<genTensor0<double>,1>::Handle SLagIncrS = OldSLagIncrS+SaveSLagIncr;
    SaveModifUnknown(*SLagIncrS,I->begin(), I->end(),*SaveSLagIncrS);
    std::shared_ptr<SavedUnknownTerm<genTensor1<double> > > OldSLag;
    savedGenTermManager.getSavedGenTerm("ILag",OldSLag,_step,0,i);
    SavedUnknownAddToSolution(*OldSLag,*pAssembler);
    SaveUnknown(*pAssembler, I->begin(), I->end(), *SaveSLag);

    savedGenTermManager.addSavedGenTerm("ILagIncr",SaveSLagIncr,_step,0,i);
    savedGenTermManager.addSavedGenTerm("ILagIncrStep",SaveSLagIncrS,_step,0,i);
    savedGenTermManager.addSavedGenTerm("ILag",SaveSLag,_step,0,i);

    printData(SaveSLagIncr,"ILagIncr",i);
    printData(SaveSLagIncrS,"ILagIncrStep",i);
    printData(SaveSLag,"ILag",i);
  }

//   savedGenTermManager.print();;
}


bool NonLinearXfemSolver::ResidualForceCriterion(const double &normFext, double &normResidual)
{
  linearSystemId<double>* buf= new linearSystemId<double>;
  dofManager<double>* pAssembler2 = new dofManager<double> ( buf );

  
  //---------------------------------------------------------------------------------------------------
  // BuildFunctionSpaces ------------------------------------------------------------------------------
  //---------------------------------------------------------------------------------------------------


  for ( unsigned int i=0; i<EDomains.size(); ++i )
  {
    // Displacement - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    for ( unsigned int j = 0; j < Data->BCs.size(); ++j )
    {
      if (CombineBCatDomain[j][i]->size()!=0)
      if ((Data->BCs[j]->kind=="Dirichlet")||(Data->BCs[j]->kind==""))
      {
        if (Data->BCs[j]->What=="Displacement")
        {
          genGroupOfElements* EF = CombineBCatDomain[j][i];
          genFilterDofComponent filterDof ( Data->BCs[j]->comp1 );
          FixNodalDofs ( FSpaceDisps[i], EF->begin(), EF->end(), *pAssembler2, *(Data->BCs[j]->fscalar), filterDof );
        }
      }
    }
  }
  
  for ( unsigned int i=0; i<EDomains.size(); ++i )
  {
    NumberDofs ( FSpaceDisps[i], EDomains[i]->begin(), EDomains[i]->end(),*pAssembler2 );
  }

  // Lagrange Multiplier - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  for ( unsigned int i = 0; i < FSpaceStableLagMult.size(); ++i)
  {
#if 0
    std::vector<int> comps(StableLagMultDomains[i]->comps.begin(), StableLagMultDomains[i]->comps.end());
    SpaceReducerForStableMultipliers* spacereducer = new SpaceReducerForStableMultipliers(StableLagMultDomains[i]->group(), Data->tag+EDomains.size()+i);
    spacereducer->BuildLinearConstraints(comps, pAssembler2);
#else
    SpaceReducer spaceReducer(StableLagMultDomains[i]->group());
    spaceReducer.BuildLinearConstraints(FSpaceStableLagMult[i], pAssembler2);
#endif
  }

  for ( unsigned int i = 0; i < FSpaceStableLagMultForInterface.size(); ++i)
  {
#if 0
    std::vector<int> comps(InterfaceDomains[i]->comps.begin(), InterfaceDomains[i]->comps.end());
    SpaceReducerForStableMultipliers* spacereducer = new SpaceReducerForStableMultipliers(InterfaceDomains[i]->group(), Data->tag+EDomains.size()+StableLagMultDomains.size()+i);
    spacereducer->BuildLinearConstraints(comps, pAssembler2);
#else
    SpaceReducer spaceReducer(InterfaceDomains[i]->group());
    spaceReducer.BuildLinearConstraints(FSpaceStableLagMultForInterface[i], pAssembler2);
#endif
  }

  for ( unsigned int i = 0; i < FSpaceStableLagMult.size(); ++i)
  {
    NumberDofs ( FSpaceStableLagMult[i], StableLagMultDomains[i]->begin(), StableLagMultDomains[i]->end(), *pAssembler2 );
  }
  for ( unsigned int i = 0; i < FSpaceStableLagMultForInterface.size(); ++i)
  {
    NumberDofs ( FSpaceStableLagMultForInterface[i], InterfaceDomains[i]->begin(), InterfaceDomains[i]->end(), *pAssembler2 );
  }
//   for (unsigned int i = 0; i<Data->BCs.size(); ++i)
//   {
//     if (Data->BCs[i]->kind=="Lagrange_GaetanMethod") // Gaetan method
//       NumberDofs ( FSpaceLagMult, Data->BCs[i]->begin(), Data->BCs[i]->end(), *pAssembler2 );
//   }


  //---------------------------------------------------------------------------------------------------
  // Assemble -----------------------------------------------------------------------------------------
  //---------------------------------------------------------------------------------------------------

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

    // build the internal force vector
    std::shared_ptr<SavedGenTerm<genTensor2<double>,0> > lastG;
    savedGenTermManager.getSavedGenTerm("Grad",lastG,_step,0,i);
    genTerm<genTensor2<double>,0>::Handle lastG2 = lastG;
    genTerm<genTensor2<double>,1>::Handle G = Gradient(FSpaceDisps[i]);
    #ifndef LARGE_STRAIN
    genTerm<genTensor2<double>,1>::Handle Bl = ( G+Transpose(G) )*0.5;
    #else
    genTerm<genTensor2<double>,1>::Handle Bl = ( G+Transpose(G) + Compose<ContrProdOp>(Transpose(lastG2),G) + Compose<ContrProdOp>(Transpose(G),lastG2) )*0.5;
    #endif
    GaussQuadrature Integ_Bulk ( GaussQuadrature::GradGrad );
    
    genTerm<genTensor0<double>,1>::Handle Fint = E->Force<1>(Bl);

    Assemble ( Fint ,E->begin(),E->end(),Integ_Bulk,*pAssembler2 );


    // build the external force vector
    GaussQuadrature Integ_Boundary ( GaussQuadrature::Val );
    for ( unsigned int j = 0; j < Data->BCs.size(); ++j )
    {
      if (CombineBCatDomain[j][i]->size()!=0)
      if ((Data->BCs[j]->kind=="Neumann")||(Data->BCs[j]->kind==""))
      {
        if (Data->BCs[j]->What=="Force")
        {
          genTerm<genTensor1<double>,0>::Handle Func(new FunctionField<genTensor1<double> >(*(Data->BCs[j]->fvector)));
          genTerm<genTensor0<double>,1>::Handle loadterm = Compose<FullContrProdOp>(FSpaceDisps[i],Func)*-1;
          genGroupOfElements* EF = CombineBCatDomain[j][i];
          Assemble ( loadterm , EF->begin(), EF->end(), Integ_Boundary, *pAssembler2 );
        }
      }
    }


    for ( unsigned int j = 0; j < StableLagMultDomains.size(); ++j)
      if (CombineSLMatDomain[j][i]->size()!=0)
      {
        genGroupOfElements* SF = CombineSLMatDomain[j][i];
        StableLagrangeMultiplierDomain* S = StableLagMultDomains[j];

        std::shared_ptr<SavedGenTerm<genTensor1<double>,0> > OldL;
        savedGenTermManager.getSavedGenTerm("Lagrangian",OldL,_step,0,j);
        genTerm<genTensor0<double>,1>::Handle CrossLagTerm1 = S->CrossLagTerm<1,0>(FSpaceDisps[i],OldL);
        Assemble ( CrossLagTerm1, SF->begin(), SF->end(), Integ_Boundary, *pAssembler2 );

        //Specifique aux conditions limites
//      genTerm<genTensor1<double>,0>::Handle Func(new FunctionField<genTensor1<double> >(*(S->f));
//      genTerm<double,1>::Handle loadterm = Compose<FullContrProdOp>(FSpaceStableLagMult[i],Func);
//      genTerm<double,1>::Handle loadterm = S->Lterm<1>(FSpaceStableLagMult[i]);
//      Assemble ( loadterm*(-1), S->begin(), S->end(), Integ_Boundary, *pAssembler2 );
// 
//      std::shared_ptr<SavedGenTerm<genTensor1<double>,0> > OldD;
//      savedGenTermManager.getSavedGenTerm("StableLagMultDisplacement",OldD,_step,0,i);
//      genTerm<double,1>::Handle CrossLagTerm2 = S->CrossLagTerm<1,0>(FSpaceStableLagMult[i],OldD);
//      Assemble ( CrossLagTerm2, S->begin(), S->end(), Integ_Boundary, *pAssembler2 );
// 
//      //Specifique a l'interface quand 1/k!=0
//      std::shared_ptr<SavedGenTerm<genTensor1<double>,0> > OldU;
//      savedGenTermManager.getSavedGenTerm("InterfaceLagrangian",OldU,_step,0,i);
//      genTerm<double,1>::Handle SelfLagTerm = S->Internal(FSpaceStableLagMultForInterface[i]);
//      Assemble ( SelfLagTerm, S->begin(), S->end(), Integ_Boundary, *pAssembler2 );
      }
  }

  for ( unsigned int i = 0; i < InterfaceDomains.size(); ++i)
  {
    genInterfaceDomain<genTensor0<double> >* I = InterfaceDomains[i];
//     InterfaceDomain* I = InterfaceDomains[i];

    GaussQuadrature Integ_Boundary ( GaussQuadrature::Val );
    std::shared_ptr<SavedGenTerm<genTensor1<double>,0> > OldL;
    savedGenTermManager.getSavedGenTerm("InterfaceLagrangian",OldL,_step,0,i);
    genTerm<genTensor0<double>,1>::Handle CrossLagTerm1 = I->CrossLagTerm<1,0>(DispJumps[i],OldL);
    Assemble ( CrossLagTerm1*-1, I->begin(), I->end(), Integ_Boundary, *pAssembler2 );

    //Specifique aux conditions limites
//     genTerm<genTensor1<double>,0>::Handle Func(new FunctionField<genTensor1<double> >(*(I->f));
//     genTerm<double,1>::Handle loadterm = Compose<FullContrProdOp>(FSpaceStableLagMultForInterface[i],Func);
//    genTerm<double,1>::Handle loadterm = I->Lterm<1>(FSpaceStableLagMultForInterface[i]);
//    Assemble ( loadterm*(-1), I->begin(), I->end(), Integ_Boundary, *pAssembler2 );

//     std::shared_ptr<SavedGenTerm<genTensor1<double>,0> > OldD;
//     savedGenTermManager.getSavedGenTerm("InterfaceDisplacement",OldD,_step,0,i);
//     genTerm<double,1>::Handle CrossLagTerm2 = I->CrossLagTerm<1,0>(FSpaceStableLagMultForInterface[i],OldD);
//     Assemble ( CrossLagTerm2, I->begin(), I->end(), Integ_Boundary, *pAssembler2 );

    //Specifique a l'interface quand 1/k!=0
//      std::shared_ptr<SavedGenTerm<genTensor1<double>,0> > OldU;
//      savedGenTermManager.getSavedGenTerm("InterfaceLagrangian",OldU,_step,0,i);
//      genTerm<double,1>::Handle SelfLagTerm = I->Internal<1>(FSpaceStableLagMultForInterface[i]);
//      Assemble ( SelfLagTerm, I->begin(), I->end(), Integ_Boundary, *pAssembler2 );
  }

//   for (unsigned int i=0; i<Data->BCs.size(); ++i)
//   {
//     if (Data->BCs[i]->kind=="Lagrange_GaetanMethod")
//     {
//       std::cout << "Dirichlet BC Disp" << std::endl;
//       GaussQuadrature Integ_Boundary ( GaussQuadrature::Val );
//       genTerm<double,0>::Handle Func(new FunctionField<double>(*(Data->BCs[i]->fscalar)));
//       genTerm<double,1>::Handle loadterm = Compose<FullContrProdOp>(FSpaceLagMult,Func);
//       Assemble ( loadterm*(-1) , Data->BCs[i]->begin(), Data->BCs[i]->end(), Integ_Boundary, *pAssembler2 );
//     }
//   }



  //---------------------------------------------------------------------------------------------------
  // Solve --------------------------------------------------------------------------------------------
  //---------------------------------------------------------------------------------------------------
  pAssembler2->systemSolve();

  
  //---------------------------------------------------------------------------------------------------
  // _Critere_Residual_ --------------------------------------------------------------------------------------------
  //---------------------------------------------------------------------------------------------------
  double criterion = Data->User.Scalars["tolerance"]*normFext;
  bool convergeance = true;
  double max_normResidual = 0;
  for ( unsigned int i=0; i<EDomains.size(); ++i )
  {
    genSField<genTensor1<double> > Residual(pAssembler2,FSpaceDisps[i]);
    bool conv = ValidateCriterion(Residual, *EDomains[i], criterion, normResidual);
//  bool convergeance = ValidateCriterion(pAssembler2, criterion, normResidual);

    if(normResidual>max_normResidual) max_normResidual = normResidual;
    if(!conv) convergeance = conv;
  }
  normResidual = max_normResidual;

//   delete pAssembler2;
//   delete buf;

  return convergeance;
}

void NonLinearXfemSolver::postprocessing(int step)
{
  NonLinearSolver::postprocessing(step);

  std::ostringstream temp;
  temp << "interface displacement " << step ;
  PView* pv = buildInterfaceDisplacementView( temp.str() );
  temp.str(""); temp << "intdisp_" << step << ".msh";
  pv->getData()->writeMSH(temp.str(), 3);
  delete pv;
  
  temp.str("");
  temp << "interface lagrangian " << step ;
  pv = buildInterfaceLagrangianView( temp.str() );
  temp.str(""); temp << "intlag_" << step << ".msh";
  pv->getData()->writeMSH(temp.str(), 3);
  delete pv;
}

PView* NonLinearXfemSolver::buildInterfaceDisplacementView ( const std::string &postFileName )
{
  std::vector<genTerm<genTensor1<double>,0>::ConstHandle> InterfaceDisplacement;
  for ( unsigned int i = 0; i < InterfaceDomains.size(); ++i )
  {
    genTerm<genTensor1<double>,0>::Handle DExtIncr ( new genSField<genTensor1<double> >(pAssembler, FSpaceDisps[InterfaceDomains[i]->DomainExt]) );
    genTerm<genTensor1<double>,0>::Handle DIntIncr ( new genSField<genTensor1<double> >(pAssembler, FSpaceDisps[InterfaceDomains[i]->DomainInt]) );
    genTerm<genTensor1<double>,0>::Handle IntDisp = DExtIncr + (DIntIncr*-1.);
    InterfaceDisplacement.push_back(IntDisp);
  }

  return buildViewNodal(InterfaceDisplacement,InterfaceDomains,postFileName);
}

PView* NonLinearXfemSolver::buildInterfaceLagrangianView ( const std::string &postFileName )
{
  std::vector<genTerm<genTensor1<double>,0>::ConstHandle> InterfaceLagrangian;
  for ( unsigned int i = 0; i < InterfaceDomains.size(); ++i )
  {
    diffTerm<genTensor1<double>,0>::Handle Field(new genSField<genTensor1<double> >( pAssembler, FSpaceStableLagMultForInterface[i] ));
    genTerm<genTensor1<double>,0>::Handle IntLag = Field;
    InterfaceLagrangian.push_back(IntLag);
  }

  return buildViewNodal(InterfaceLagrangian,InterfaceDomains,postFileName);
}

/*PView* NonLinearXfemSolver::buildStableLagrangeMultiplierView ( const std::string &postFileName )
{
  std::vector<genTerm<genTensor1<double>,0>::Handle> Fields;
  for ( unsigned int i = 0; i < StableLagMultDomains.size(); ++i )
  {
    Fields.push_back(new genSField<genTensor1<double> >( pAssembler, FSpaceStableLagMult[i] ));
  }

  return buildViewNodal(Fields,StableLagMultDomains,postFileName);
}*/

/*PView* NonLinearXfemSolver::buildStableLagrangeMultiplierView ( const std::string &postFileName )
{
  std::cout <<  "build Stable Lagrange Multiplier View"<< std::endl;

  // store the StableLagMult values <vertex,value>
  std::map<int, std::vector<double> > data;
  PView* pv = NULL;

  // for new elements
  std::vector<MVertex*> vecVertices;
  std::vector <MElement*> vecElements;
  std::map<int, std::vector<MElement*> > Cut_Elements;

  // copy of _pModel
  GModel* pModelPost = new GModel ( *_pModel );
  GModel::setCurrent ( pModelPost );

  // tolerance for discontinuity point calcul
  double eps = 1e-6;

  // entity associate to cut elements
  std::map<MElement*, GEntity*> ElementCutEntity;
  typedef std::set<GFace*, GEntityLessThan>::iterator fiter;
  for ( fiter it = pModelPost->firstFace(); it != pModelPost->lastFace(); ++it )
    for ( unsigned int i = 0; i < ( *it )->polygons.size(); i++ )
      ElementCutEntity[ ( *it )->polygons[i]] = ( *it );

  // REGION
  typedef std::set<GRegion*, GEntityLessThan>::iterator riter;
  for ( riter it = pModelPost->firstRegion(); it != pModelPost->lastRegion(); ++it )
    for ( unsigned int i = 0; i < ( *it )->polyhedra.size(); i++ )
      ElementCutEntity[ ( *it )->polyhedra[i]] = ( *it );


  for ( unsigned int i = 0; i < _StableLagrangeMultipliersSpace.size(); ++i )
  {
    SolverField<genTensor1<double> > Field(_pAssembler, _StableLagrangeMultipliersSpace[i]);
    groupOfElements::elementContainer::const_iterator it;

    for ( it = _StableLagrangeMultiplierFields[i]._g->begin(); it != _StableLagrangeMultiplierFields[i]._g->end(); ++it )
    {
      MElement* e=*it;
      for ( int j = 0; j < e->getNumVertices(); ++j )
      {
        MVertex* v = e->getVertex(j);
        SPoint3 p ( v->x(), v->y(), v->z() );
        double pnt[3];
        pnt[0]=p.x();
        pnt[1]=p.y();
        pnt[2]=p.z();
        double uvw[3];
        e->xyz2uvw ( pnt,uvw );
        genTensor1<double> val;
        Field.f( e,uvw[0],uvw[1],uvw[2],val );
        std::vector<double> vec ( 3 );
        vec[0]=val(0);
        vec[1]=val(1);
        vec[2]=val(2);
        data[ e->getVertex(j)->getNum() ]=vec;
        // printf("Vertex index: %d, force[ %lf %lf %lf ]\n", e->getVertex ( j )->getNum(), vec[0], vec[1], vec[2]);
        // La valeur pour un vertex dans deux elt distincts est identique
      }
    }

    if (pv == NULL)
      pv = new PView (postFileName, "NodeData", _pModel, data, i);
    else
      pv->addStep(_pModel, data, i, 3);

    int temp=0;
    std::vector<double> vec(3,0);
    std::map<int, std::vector<double> >::iterator itdata;
    for( itdata = data.begin(); itdata != data.end(); ++itdata)
    {
      vec[0] += (*itdata).second[0];
      vec[1] += (*itdata).second[1];
      vec[2] += (*itdata).second[2];
      temp++;
    }

    printf("nb de lag: %d\n", temp);
    double eqfac = _elasticFields[0]->scale();
    vec[0] *= eqfac/(double)temp;
    vec[1] *= eqfac/(double)temp;
    vec[2] *= eqfac/(double)temp;

    printf("LS n°%d : { %lf %lf %lf }\n ", i, vec[0], vec[1], vec[2] );
  }

  std::map < int, std::map<int, std::string > > physical; // empty

  // store new elements and vertices
  pModelPost->storeChain ( vecVertices,_dim, Cut_Elements, physical );

  return pv;
}*/


/*PView* NonLinearXfemSolver::buildDamageView ( const std::string &postFileName )
{

  for ( unsigned int i = 0; i < InterfaceDomains.size(); ++i )
  {
    InterfaceDomain* I = InterfaceDomains[i];
    
  }
  genSField<genTensor1<double> > Field ( pAssembler, FSpaceDisp );
  return buildViewNodal(Field,EDomains,postFileName);
}*/