// 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>.

#ifndef _NONLINEARXFEM_DOMAIN_H_
#define _NONLINEARXFEM_DOMAIN_H_

#include <string>
#include <vector>

#include "genSupports.h"
#include "genTerm.h"
#include "genTensors.h"
#include "genDataIO.h"
#include "genMatrix.h"
#include "savedGenTerm.h"

#include "savedGenTermManager.h"
#include "genAlgorithms.h"

#include "nonLinearDomain.h"
#include "elasticXfemDomain.h"

template <class T> class genInterfaceDomain : public NLDomain
{
public :
  typename genTerm<T,0>::Handle S; // Compliance
  typename genTerm<T,0>::Handle dS; // Tangent Compliance

  genTerm<genTensor1<double>,0>::Handle N; // Normal
  genTerm<genTensor1<double>,0>::Handle U; // Displacement

  int DomainInt, DomainExt;

  double eqfac;
  std::set<int> comps;

public :
  genInterfaceDomain(const genDomain &d, const double eqfac_=1.) : eqfac(eqfac_)
  {
    g = new genGroupOfElements(d.dim, d.tag);

    dim = d.dim;
    tag = d.tag;

    comps.insert(0);
    comps.insert(1);
    comps.insert(2);

//     History.defineMaxSteps(10);

    DomainInt = d.Integers.find("DomainInt")->second;
    DomainExt = d.Integers.find("DomainExt")->second;
    assert(DomainInt!=DomainExt);

    N = buildInterfaceBasis(begin(),end(),dim);

    genTerm<genTensor1<double>,0>::Handle InitDisp = genTerm<genTensor1<double>,0>::Handle(new ConstantField<genTensor1<double> >(genTensor1<double>(0.)));
    GaussQuadrature Integ_Val ( GaussQuadrature::Val );
    U = initGenTermInHistory<genTensor1<double> >("LagDisp", Integ_Val, InitDisp);

//- Spécifique à chaque ordre de Tenseur -------------------------------------------------------
//     double s = 0.;
//     if ( d.Scalars.find("k0") != d.Scalars.end() )
//       s = 1./d.Scalars.find("k0")->second;

//     S = genTerm<double,0>::Handle(new ConstantField<double>(s));
//     dS = genTerm<double,0>::Handle(new ConstantField<double>(s));
  }

  template <int n1, int n2> typename genTerm<genTensor0<double>,n1+n2>::Handle
  CrossLagTerm( const typename genTerm<genTensor1<double>,n1>::ConstHandle &FSpace1, const typename genTerm<genTensor1<double>,n2>::ConstHandle &FSpace2 );
  template <int n1, int n2> typename genTerm<genTensor0<double>, n1+n2>::Handle
  SelfLagTerm( const typename genTerm<genTensor1<double>,n1>::ConstHandle &SLag1, const typename genTerm<genTensor1<double>,n2>::ConstHandle &SLag2);

  virtual typename genTerm<T,0>::Handle Compliance(const genTerm<genTensor1<double>,0>::ConstHandle &deltaDisp, const genTerm<genTensor1<double>,0>::ConstHandle &Disp);
  virtual typename genTerm<T,0>::Handle Tangent(const genTerm<genTensor1<double>,0>::ConstHandle &deltaDisp, const genTerm<genTensor1<double>,0>::ConstHandle &Disp);
  template <int n> typename genTerm<genTensor1<double>,n>::Handle deltaDisplacement(const typename genTerm<genTensor1<double>,n>::ConstHandle &deltaForce);
  genTerm<genTensor1<double>,0>::Handle Displacement();
  template <int n> typename genTerm<genTensor0<double>,n>::Handle Internal(const typename genTerm<genTensor1<double>,n>::ConstHandle &FSpace);

  genInterfaceDomain() : NLDomain() {}
  virtual ~genInterfaceDomain() {}

  template <class Iterator> genTerm<genTensor1<double>,0>::Handle buildInterfaceBasis(Iterator itbegin, Iterator itend, int dimension); //WARNING a definir dans genAlgo de maniere + generale

public : //WARNING temporaire
  virtual void copyHistory() =0;
  virtual void saveHistory(const genTerm<genTensor1<double>,0>::ConstHandle &deltaDisp, const genTerm<genTensor1<double>,0>::ConstHandle &Disp,
                           const genTerm<genTensor1<double>,0>::ConstHandle &deltaForce, const genTerm<genTensor1<double>,0>::ConstHandle &Force) =0;
};

template<class T> template <int n1, int n2> typename genTerm<genTensor0<double>,n1+n2>::Handle
genInterfaceDomain<T>::CrossLagTerm( const typename genTerm<genTensor1<double>,n1>::ConstHandle &FSpace1, const typename genTerm<genTensor1<double>,n2>::ConstHandle &FSpace2)
{
  return typename genTerm<genTensor0<double>,n1+n2>::Handle (Compose<FullContrProdOp>(FSpace1,FSpace2)*eqfac);
}

template<class T> template <int n1, int n2> typename genTerm<genTensor0<double>,n1+n2>::Handle
genInterfaceDomain<T>::SelfLagTerm( const typename genTerm<genTensor1<double>,n1>::ConstHandle &FSpace1, const typename genTerm<genTensor1<double>,n2>::ConstHandle &FSpace2)
{
//   double eqfac2= eqfac*eqfac;
//   return typename genTerm<double,n1+n2>::Handle ( Compose<FullContrProdOp>(FSpace1,Compose<FullContrProdOp>(dS,FSpace2))*eqfac2);
  return typename genTerm<genTensor0<double>,n1+n2>::Handle ( Compose<FullContrProdOp>(FSpace1,Compose<FullContrProdOp>(dS,FSpace2))*eqfac);
}

template<class T> typename genTerm<T,0>::Handle
genInterfaceDomain<T>::Compliance(const genTerm<genTensor1<double>,0>::ConstHandle &deltaDisp, const genTerm<genTensor1<double>,0>::ConstHandle &Disp)
{
  return typename genTerm<T,0>::Handle (S);
}

template<class T> typename genTerm<T,0>::Handle
genInterfaceDomain<T>::Tangent(const genTerm<genTensor1<double>,0>::ConstHandle &deltaDisp, const genTerm<genTensor1<double>,0>::ConstHandle &Disp)
{
  return typename genTerm<T,0>::Handle (dS);
}

template<class T> template <int n> typename genTerm<genTensor1<double>,n>::Handle
genInterfaceDomain<T>::deltaDisplacement(const typename genTerm<genTensor1<double>,n>::ConstHandle &deltaForce)
{
  return typename genTerm<genTensor1<double>,n>::Handle(Compose<FullContrProdOp>(dS,deltaForce));
}

template<class T> genTerm<genTensor1<double>,0>::Handle
genInterfaceDomain<T>::Displacement()
{
  return genTerm<genTensor1<double>,0>::Handle(U);
}

template<class T> template <int n> typename genTerm<genTensor0<double>,n>::Handle
genInterfaceDomain<T>::Internal(const typename genTerm<genTensor1<double>,n>::ConstHandle &FSpace)
{
  return typename genTerm<genTensor0<double>,n>::Handle(Compose<FullContrProdOp>(FSpace,U)*eqfac);
}

template<class T> template <class Iterator> genTerm<genTensor1<double>,0>::Handle
genInterfaceDomain<T>::buildInterfaceBasis(Iterator itbegin, Iterator itend, int dimension)
{
  std::shared_ptr<SavedGenTerm<genTensor1<double>,0> > N_(new SavedGenTerm<genTensor1<double>,0>);
//   std::shared_ptr<SavedGenTerm<genTensor2<double>,0> > T_(new SavedGenTerm<genTensor2<double>,0>);
  if (dimension==2) for (Iterator it = itbegin; it != itend; ++it)
  {
    MElement* e = *it;
    IntPt* GP;
    GaussQuadrature integrator(GaussQuadrature::Val);
    int npts = integrator.getIntPoints(e, &GP);
    std::vector<MVertex*> v;
    e->getVertices(v);

    //Vecteur normal à l'interface 2D
    genTensor1<double> tan1(v[1]->x()-v[0]->x(),v[1]->y()-v[0]->y(),v[1]->z()-v[0]->z());
    genTensor1<double> tan2(v[2]->x()-v[1]->x(),v[2]->y()-v[1]->y(),v[2]->z()-v[1]->z());
    genTensor1<double> normal = crossprod(tan1,tan2);
    normal.normalize();

    //Vecteurs tangents 1-2 à l'interface 2D
    tan1.normalize();
    tan2 = crossprod(normal,tan1);
    tan2.normalize();

    //Repere local
    genMatrix<double> axes(3,3);
    axes.set(0,0, tan1[0]); axes.set(0,1, tan2[0]); axes.set(0,2, normal[0]);
    axes.set(1,0, tan1[1]); axes.set(1,1, tan2[1]); axes.set(1,2, normal[1]);
    axes.set(2,0, tan1[2]); axes.set(2,1, tan2[2]); axes.set(2,2, normal[2]);
    genTensor2<double> axes_;
    axes_.setMat(axes);
    std::vector<genScalar<genTensor1<double> > > normal2(npts,normal); //WARNING cst per GP
//     std::vector<genScalar<genTensor2<double> > > axes2(npts,axes_); //WARNING cst per GP
    N_->set(e, npts, GP, normal2);
//     T_->set(e, npts, GP, axes2);
  }
  else if(dimension==1) for (Iterator it = itbegin; it != itend; ++it)
  {
    MElement* e = *it;
    IntPt* GP;
    GaussQuadrature integrator(GaussQuadrature::Val);
    int npts = integrator.getIntPoints(e, &GP);
    std::vector<MVertex*> v;
    e->getVertices(v);

    //Vecteur tangent à l'interface en 1D
    genTensor1<double> tan1(v[1]->x()-v[0]->x(),v[1]->y()-v[0]->y(),v[1]->z()-v[0]->z());
    tan1.normalize();

    //Utilisation de l'axe z pour construire la normale à l'interface
    genTensor1<double> tan2(0, 0, 1); // WARNING cela impose de travailler dans le plan xy
    genTensor1<double> normal = crossprod(tan1,tan2);
    normal.normalize();

    //Repere local
    genMatrix<double> axes(3,3);
    axes.set(0,0, tan1[0]); axes.set(0,1, tan2[0]); axes.set(0,2, normal[0]);
    axes.set(1,0, tan1[1]); axes.set(1,1, tan2[1]); axes.set(1,2, normal[1]);
    axes.set(2,0, tan1[2]); axes.set(2,1, tan2[2]); axes.set(2,2, normal[2]);
    genTensor2<double> axes_;
    axes_.setMat(axes);
    std::vector<genScalar<genTensor1<double> > > normal2(npts,normal); //WARNING cst per GP
//     std::vector<genScalar<genTensor2<double> > > axes2(npts,axes_); //WARNING cst per GP
    N_->set(e, npts, GP, normal2);
 //   T_->set(e, npts, GP, axes2);
  }
  N_->printfile("data_InterfaceNormal.txt");
//  T_->printfile("data_InterfaceBasis.txt");
  return N_;
}


class genPerfectInterfaceDomain : public genInterfaceDomain<genTensor0<double> >
{
public :
  genPerfectInterfaceDomain(const genDomain &d, const double eqfac_=1.) : genInterfaceDomain<genTensor0<double> >::genInterfaceDomain(d, eqfac_)
  {
//- Spécifique à chaque ordre de Tenseur - Definie dans le constructeur de base genInterfaceDomain ------------------------------------------------------
    double s = 0.;

    genTerm<genTensor0<double>,0>::Handle InitS = genTerm<genTensor0<double>,0>::Handle(new ConstantField<genTensor0<double> >(s));
    genTerm<genTensor0<double>,0>::Handle InitdS = genTerm<genTensor0<double>,0>::Handle(new ConstantField<genTensor0<double> >(s));
    GaussQuadrature Integ_Val ( GaussQuadrature::Val );
    S = initGenTermInHistory<genTensor0<double> >("Compliance", Integ_Val,InitS);
    dS = initGenTermInHistory<genTensor0<double> >("Tangent", Integ_Val,InitdS);
  }

  genPerfectInterfaceDomain () : genInterfaceDomain<genTensor0<double> >::genInterfaceDomain() {}
  virtual ~genPerfectInterfaceDomain() {};

public :
  virtual void copyHistory();
  virtual void saveHistory(const genTerm<genTensor1<double>,0>::ConstHandle &deltaDisp, const genTerm<genTensor1<double>,0>::ConstHandle &Disp,
                           const genTerm<genTensor1<double>,0>::ConstHandle &deltaForce, const genTerm<genTensor1<double>,0>::ConstHandle &Force);
};

inline void genPerfectInterfaceDomain::copyHistory()
{
  int step = genInterfaceDomain<genTensor0<double> >::History.getcurrentstep("Compliance");
  std::shared_ptr<SavedGenTerm<genTensor0<double>,0> > OldCompliance;
  genInterfaceDomain<genTensor0<double> >::History.getSavedGenTerm("Compliance", OldCompliance, step);
  genInterfaceDomain<genTensor0<double> >::History.addSavedGenTerm("Compliance", OldCompliance, step+1);

  std::shared_ptr<SavedGenTerm<genTensor0<double>,0> > OldTangent;
  genInterfaceDomain<genTensor0<double> >::History.getSavedGenTerm("Tangent", OldTangent, step);
  genInterfaceDomain<genTensor0<double> >::History.addSavedGenTerm("Tangent", OldTangent, step+1);

  std::shared_ptr<SavedGenTerm<genTensor1<double>,0> > OldDisp;
  genInterfaceDomain<genTensor0<double> >::History.getSavedGenTerm("LagDisp", OldDisp, step);
  genInterfaceDomain<genTensor0<double> >::History.addSavedGenTerm("LagDisp", OldDisp, step+1);
}

inline void genPerfectInterfaceDomain::saveHistory(const genTerm<genTensor1<double>,0>::ConstHandle &deltaDisp, const genTerm<genTensor1<double>,0>::ConstHandle &Disp,
                                            const genTerm<genTensor1<double>,0>::ConstHandle &deltaForce, const genTerm<genTensor1<double>,0>::ConstHandle &Force)
{
  //Save compliance
  genTerm<genTensor0<double>,0>::Handle newCompliance = genTerm<genTensor0<double>,0>::Handle(new ConstantField<genTensor0<double> >(0.));
  GaussQuadrature Integ_Val ( GaussQuadrature::Val );
  genInterfaceDomain<genTensor0<double> >::S = genInterfaceDomain<genTensor0<double> >::saveGenTermInHistory<genTensor0<double> >("Compliance", Integ_Val, newCompliance);

  //Save tangent
  genTerm<genTensor0<double>,0>::Handle newTangent = genTerm<genTensor0<double>,0>::Handle(new ConstantField<genTensor0<double> >(0.));
//   GaussQuadrature Integ_Val ( GaussQuadrature::Val );
  genInterfaceDomain<genTensor0<double> >::dS = genInterfaceDomain<genTensor0<double> >::saveGenTermInHistory<genTensor0<double> >("Tangent", Integ_Val, newTangent);

  //Save disp
  genTerm<genTensor1<double>,0>::Handle newDisp = genTerm<genTensor1<double>,0>::Handle(new ConstantField<genTensor1<double> >(genTensor1<double>()));
//   GaussQuadrature Integ_Val ( GaussQuadrature::Val );
  genInterfaceDomain<genTensor0<double> >::U = genInterfaceDomain<genTensor0<double> >::saveGenTermInHistory<genTensor1<double> >("LagDisp", Integ_Val, newDisp);
}



template <class T1,class T2> class genDamageInterfaceDomain : public genInterfaceDomain<T2>
{
public :
  typename genTerm<T1,0>::Handle D;

public :
  genDamageInterfaceDomain(const genDomain &d, const double eqfac_=1.) : genInterfaceDomain<T2>::genInterfaceDomain(d, eqfac_)
  {
    typename genTerm<T1,0>::Handle InitDamage = typename genTerm<T1,0>::Handle(new ConstantField<T1>(T1(0.)));
    GaussQuadrature Integ_Val ( GaussQuadrature::Val );
    D = genInterfaceDomain<T2>::template initGenTermInHistory<T1>("Damage", Integ_Val,InitDamage);
  }
  virtual typename genTerm<T2,0>::Handle Compliance(const genTerm<genTensor1<double>,0>::ConstHandle &deltaDisp, const genTerm<genTensor1<double>,0>::ConstHandle &Disp);
  virtual typename genTerm<T2,0>::Handle Tangent(const genTerm<genTensor1<double>,0>::ConstHandle &deltaDisp, const genTerm<genTensor1<double>,0>::ConstHandle &Disp);
  virtual typename genTerm<T1,0>::Handle deltaDamage(const genTerm<genTensor1<double>,0>::ConstHandle &deltaDisp, const genTerm<genTensor1<double>,0>::ConstHandle &Disp);
  virtual typename genTerm<T1,0>::Handle Damage();
  genDamageInterfaceDomain() : genInterfaceDomain<T2>::genInterfaceDomain() {}
  virtual ~genDamageInterfaceDomain() {}

public :
  virtual void copyHistory();
  virtual void saveHistory(const genTerm<genTensor1<double>,0>::ConstHandle &deltaDisp, const genTerm<genTensor1<double>,0>::ConstHandle &Disp,
                           const genTerm<genTensor1<double>,0>::ConstHandle &deltaForce, const genTerm<genTensor1<double>,0>::ConstHandle &Force);
};

template <class T1,class T2> typename genTerm<T2,0>::Handle
genDamageInterfaceDomain<T1,T2>::Compliance(const genTerm<genTensor1<double>,0>::ConstHandle &deltaDisp, const genTerm<genTensor1<double>,0>::ConstHandle &Disp)
{
//   genTerm<double,0>::Handle One = genTerm<double,0>::Handle(new ConstantField<double> (1.));
//   genTerm<double,0>::Handle State = One+(D*-1.);
//   return genTerm<double,0>::Handle (Divide(S,State));
  return typename genTerm<T2,0>::Handle (genInterfaceDomain<T2>::S);
}

template <class T1,class T2> typename genTerm<T2,0>::Handle
genDamageInterfaceDomain<T1,T2>::Tangent(const genTerm<genTensor1<double>,0>::ConstHandle &deltaDisp, const genTerm<genTensor1<double>,0>::ConstHandle &Disp)
{
// pour le moment valable uniquement pour la loi bitriangulaire
  return typename genTerm<T2,0>::Handle (genInterfaceDomain<T2>::dS);
}

template <class T1,class T2> typename genTerm<T1,0>::Handle
genDamageInterfaceDomain<T1,T2>::deltaDamage(const genTerm<genTensor1<double>,0>::ConstHandle &deltaDisp, const genTerm<genTensor1<double>,0>::ConstHandle &Disp)
{
  typename genTerm<T1,0>::Handle Zero = genTerm<genTensor0<double>,0>::Handle(new ConstantField<genTensor0<double> > (T1(0.)));
  return typename genTerm<T1,0>::Handle(Zero);
}

template <class T1,class T2> typename genTerm<T1,0>::Handle
genDamageInterfaceDomain<T1,T2>::Damage()
{
  return typename genTerm<T1,0>::Handle(D);
}

template <class T1,class T2> void
genDamageInterfaceDomain<T1,T2>::copyHistory()
{
  int step = genInterfaceDomain<T2>::History.getcurrentstep("Compliance");
  typename std::shared_ptr<SavedGenTerm<T2,0> > OldCompliance;
  genInterfaceDomain<T2>::History.getSavedGenTerm("Compliance", OldCompliance, step);
  genInterfaceDomain<T2>::History.addSavedGenTerm("Compliance", OldCompliance, step+1);

  typename std::shared_ptr<SavedGenTerm<T2,0> > OldTangent;
  genInterfaceDomain<T2>::History.getSavedGenTerm("Tangent", OldTangent, step);
  genInterfaceDomain<T2>::History.addSavedGenTerm("Tangent", OldTangent, step+1);

  typename std::shared_ptr<SavedGenTerm<genTensor1<double>,0> > OldDisp;
  genInterfaceDomain<T2>::History.getSavedGenTerm("LagDisp", OldDisp, step);
  genInterfaceDomain<T2>::History.addSavedGenTerm("LagDisp", OldDisp, step+1);

  typename std::shared_ptr<SavedGenTerm<T1,0> > OldDamage;
  genInterfaceDomain<T2>::History.getSavedGenTerm("Damage", OldDamage, step);
  genInterfaceDomain<T2>::History.addSavedGenTerm("Damage", OldDamage, step+1);
}

template <class T1,class T2> void
genDamageInterfaceDomain<T1,T2>::saveHistory(const genTerm<genTensor1<double>,0>::ConstHandle &deltaDisp, const genTerm<genTensor1<double>,0>::ConstHandle &Disp,
                                             const genTerm<genTensor1<double>,0>::ConstHandle &deltaForce, const genTerm<genTensor1<double>,0>::ConstHandle &Force)
{
  //Save damage
  typename genTerm<T1,0>::Handle oldDamage = Damage();
  typename genTerm<T1,0>::Handle dDamage = deltaDamage(deltaDisp,Disp);
  typename genTerm<T1,0>::Handle newDamage = oldDamage+dDamage;
  GaussQuadrature Integ_Val ( GaussQuadrature::Val );
  D = genInterfaceDomain<T2>::template saveGenTermInHistory<T1>("Damage", Integ_Val, newDamage);

  //Save compliance
  typename genTerm<T2,0>::Handle newCompliance = Compliance(deltaDisp,Disp);
  genInterfaceDomain<T2>::S = genInterfaceDomain<T2>::template saveGenTermInHistory<T2>("Compliance", Integ_Val, newCompliance);

  //Save tangent
  typename genTerm<T2,0>::Handle newTangent = Tangent(deltaDisp,Disp); // loi bitriangulaire
//   genTerm<genTensor1<double>,0>::Handle v = genTerm<genTensor1<double>,0>::Handle(new ConstantField<genTensor1<double> >(genTensor1<double>(0.,2.,0.)));
//   genTerm<double,0>::Handle newTangent = Compose<ContrProdOp>(Force,v); // loi racine carree
//   GaussQuadrature Integ_Val ( GaussQuadrature::Val );
  genInterfaceDomain<T2>::dS = genInterfaceDomain<T2>::template saveGenTermInHistory<T2>("Tangent", Integ_Val, newTangent);

  //Save disp
  genTerm<genTensor1<double>,0>::Handle oldDisp = genInterfaceDomain<T2>::Displacement();
  genTerm<genTensor1<double>,0>::Handle dDisp = genInterfaceDomain<T2>::template deltaDisplacement<0>(deltaForce);
  genTerm<genTensor1<double>,0>::Handle newDisp = oldDisp+dDisp;
//   GaussQuadrature Integ_Val ( GaussQuadrature::Val );
  genInterfaceDomain<T2>::U = genInterfaceDomain<T2>::template saveGenTermInHistory<genTensor1<double> >("LagDisp", Integ_Val, newDisp);
}


class genBitriangularDamageInterfaceDomain : public genDamageInterfaceDomain<genTensor0<double>,genTensor0<double> >
{
public :
  //Données liés à la loi d'évolution bitriangulaire
  double Gic;
  double k0;
  double Y0t;

public :
  genBitriangularDamageInterfaceDomain(const genDomain &d, const double eqfac_=1.) : genDamageInterfaceDomain<genTensor0<double>,genTensor0<double> >::genDamageInterfaceDomain(d, eqfac_)
  {
    Gic = d.Scalars.find("Gic")->second;
    k0 = d.Scalars.find("k0")->second;
    Y0t = d.Scalars.find("Y0t")->second;

//- Spécifique à chaque ordre de Tenseur - Definie dans le constructeur de base genInterfaceDomain ------------------------------------------------------
    double s = 0.;
    if ( d.Scalars.find("k0") != d.Scalars.end() )
      s = 1./d.Scalars.find("k0")->second;

    genTerm<genTensor0<double>,0>::Handle InitS = genTerm<genTensor0<double>,0>::Handle(new ConstantField<genTensor0<double> >(s));
    genTerm<genTensor0<double>,0>::Handle InitdS = genTerm<genTensor0<double>,0>::Handle(new ConstantField<genTensor0<double> >(s));
    GaussQuadrature Integ_Val ( GaussQuadrature::Val );
    S = initGenTermInHistory<genTensor0<double> >("Compliance", Integ_Val,InitS);
    dS = initGenTermInHistory<genTensor0<double> >("Tangent", Integ_Val,InitdS);
  }

  virtual genTerm<genTensor0<double>,0>::Handle Tangent(const genTerm<genTensor1<double>,0>::ConstHandle &deltaForce, const genTerm<genTensor1<double>,0>::ConstHandle &Force);
  virtual genTerm<genTensor0<double>,0>::Handle deltaDamage(const genTerm<genTensor1<double>,0>::ConstHandle &deltaForce, const genTerm<genTensor1<double>,0>::ConstHandle &Force);
//  virtual genTerm<double,0>::Handle deltaDamage(const genTerm<genTensor2<double>,0>::ConstHandle &deltaStress, const genTerm<genTensor2<double>,0>::ConstHandle &Stress);
  genBitriangularDamageInterfaceDomain () : genDamageInterfaceDomain<genTensor0<double>,genTensor0<double> >::genDamageInterfaceDomain() {}
  virtual ~genBitriangularDamageInterfaceDomain() {};
};

genTerm<genTensor0<double>,0>::Handle inline genBitriangularDamageInterfaceDomain::Tangent(const genTerm<genTensor1<double>,0>::ConstHandle &deltaDisp, const genTerm<genTensor1<double>,0>::ConstHandle &Disp)
{
  std::shared_ptr<SavedGenTerm<genTensor0<double>,0> > dS_(new SavedGenTerm<genTensor0<double>,0>);
//   std::shared_ptr<SavedGenTerm<double,0> > dSign(new SavedGenTerm<double,0>);

  //Paramètres lois d'endommagement
  double sigmaMax = sqrt(2*Y0t*k0);
  double epsilonMax = sigmaMax/k0;
  double epsilonCri = 2*Gic/sigmaMax;
  double k1 = -pow(sigmaMax,2)*k0/(2*Gic*k0-pow(sigmaMax,2));

  std::cout<<"   __________________________________________"<<std::endl;
  std::cout<<"  |Gic "<<Gic<<" k0 "<<k0<<" Y0t "<<Y0t<<std::endl;
  std::cout<<"  |sigmaMax "<<sigmaMax<<" epsilonMax "<<epsilonMax<<" epsilonCri "<<epsilonCri<<" k1 "<<k1<<std::endl;
  std::cout<<"  |1/k0 "<<1./k0<<" 1/k1 "<<1./k1<<std::endl;
  double epsilonEq_Max = 0;
  bool printDamage = true;
  bool printFullDamage = true;

  //genTerm<genTensor1<double>,0>::ConstHandle TForce = Compose<ContrProdOp>(Compose<ContrProdOp>(Transpose(T),Force),T);
  double h = 1.;
  genTerm<genTensor2<double>,0>::ConstHandle Strain = (Apply<TensProdOp>(Disp,N)+Apply<TensProdOp>(N,Disp))*(0.5/h);

  for (std::set<MElement*>::const_iterator it = begin(); it != end(); ++it)
  {
    MElement* e = *it;
    IntPt* GP;
    GaussQuadrature integrator(GaussQuadrature::Val);
    int npts = integrator.getIntPoints(e, &GP);

    //Déformations sur l'élément
    std::vector<genScalar<genTensor2<double> > > s_vvals(npts);
    Strain->get(e, npts, GP, s_vvals);

    //Contraintes sur l'élément
    //std::vector<genScalar<genTensor1<double> > > f_vvals(npts);
    //TForce->get(e, npts, GP, f_vvals);

    //Endommagement (modification de k)
//     std::vector<genScalar<double> > d_vvals(npts);
//     D->get(e, npts, GP, d_vvals);

    // Compliance
    std::vector<genScalar<genTensor0<double> > > T_vvals(npts);
    S->get(e, npts, GP, T_vvals);

//     std::vector<genScalar<double> > Sign_vvals(npts);

    for (unsigned int i=0; i<npts; ++i)
    {
      //double ForceEq = f_vvals[i].norm()*eqfac;
      const genTensor2<double> &s = s_vvals[i]();
      double epsilonEq = sqrt( 0.5*( pow((s(0,0)-s(1,1)),2)
                                    +pow((s(1,1)-s(2,2)),2)
                                    +pow((s(2,2)-s(0,0)),2)
                                    +6*( pow(s(0,1),2)
                                        +pow(s(1,2),2)
                                        +pow(s(2,0),2))) ) * eqfac;

      if(epsilonEq > epsilonEq_Max) epsilonEq_Max = epsilonEq;

      if(epsilonEq<=epsilonMax)
      {
        T_vvals[i] = T_vvals[i];
      }
      else if(epsilonEq>epsilonMax)
      {
        if (printDamage) {std::cout<<"  |-*-*-*-*-*-*-*-*-*-*-DAMAGE-*-*-*-*-*-*-*-*-*- "<<std::endl; printDamage=false;}

        if (epsilonEq<epsilonCri)
        {
          T_vvals[i]() = 1./k1;
        }
        else
        {
          T_vvals[i]() = 1.E10;
          if (printFullDamage) {std::cout<<"  |******************-FULLDAMAGE-**************** "<<std::endl; printFullDamage=false;}
        }
      }
    }
    dS_->set(e, npts, GP, T_vvals);

//     dSign->set(e, npts, GP, Sign_vvals);

  }
  std::cout<<"  |epsilonEq_Max "<<epsilonEq_Max<<std::endl;

  return dS_;
}

genTerm<genTensor0<double>,0>::Handle inline genBitriangularDamageInterfaceDomain::deltaDamage(const genTerm<genTensor1<double>,0>::ConstHandle &deltaDisp, const genTerm<genTensor1<double>,0>::ConstHandle &Disp)
{
  std::shared_ptr<SavedGenTerm<genTensor0<double>,0> > dD(new SavedGenTerm<genTensor0<double>,0>);

  //Paramètres lois d'endommagement
  double sigmaMax = sqrt(2*Y0t*k0);
  double epsilonMax = sigmaMax/k0;
  double epsilonCri = 2*Gic/sigmaMax;
  double k1 = pow(sigmaMax,2)*k0/(2*Gic*k0-pow(sigmaMax,2));

  //genTerm<genTensor1<double>,0>::ConstHandle TForce = Compose<ContrProdOp>(Compose<ContrProdOp>(Transpose(T),Force),T);
  double h = 1.;
  genTerm<genTensor2<double>,0>::ConstHandle Strain = Apply<TensProdOp>(Disp,N)*(1./h);

  for (std::set<MElement*>::const_iterator it = begin(); it != end(); ++it)
  {
    MElement* e = *it;
    IntPt* GP;
    GaussQuadrature integrator(GaussQuadrature::Val);
    int npts = integrator.getIntPoints(e, &GP);

    //Déformations sur l'élément
    std::vector<genScalar<genTensor2<double> > > s_vvals(npts);
    Strain->get(e, npts, GP, s_vvals);

    //Endommagement (modification de k)
    std::vector<genScalar<genTensor0<double> > > d_vvals(npts);
    D->get(e, npts, GP, d_vvals);
    for (unsigned int i=0; i<npts; ++i)
    {
      //double ForceEq = f_vvals[i].norm()*eqfac;
      const genTensor2<double> &s = s_vvals[i]();
      double epsilonEq = sqrt( 0.5*( pow((s(0,0)-s(1,1)),2)
                                    +pow((s(1,1)-s(2,2)),2)
                                    +pow((s(2,2)-s(0,0)),2)
                                    +6*( pow(s(0,1),2)
                                        +pow(s(1,2),2)
                                        +pow(s(2,0),2))) ) * eqfac;

      if(epsilonEq<=epsilonMax)
      {
        d_vvals[i]() = 0.;
      }
      else if(epsilonEq>epsilonMax)
      {
        if (epsilonEq<epsilonCri)
        {
          double newD = (1+k1/k0)*(1-sigmaMax/(k0*epsilonEq));
          d_vvals[i]() = newD - d_vvals[i]();
        }
        else
        {
          d_vvals[i]() = 0.99999;
//           std::cout<<"WARNING full Damage in element "<< e->getNum()
//                    <<" for the GP "<< i
//                    <<" with the value (epsilonEq-epsilonCri)/epsilonCri="<< ((epsilonEq-epsilonCri)/epsilonCri*100) <<"%"
//                    <<std::endl;
        }
      }
    }
    dD->set(e, npts, GP, d_vvals);
  }
  return dD;
}

/*
class genBitriangularDamageInterfaceDomain : public genDamageInterfaceDomain<genTensor2<double>,genTensor2<double> >
{
public :
  //Données liés à la loi d'évolution bitriangulaire
  std::vector<double> Gic;
  std::vector<double> k0;
  std::vector<double> Y0t;

public :
  genBitriangularDamageInterfaceDomain(const genDomain &d, const double eqfac_=1.) : genDamageInterfaceDomain<double,double>(d, eqfac_)
  {
    Gic = d.Vectors.find("Gic")->second;
    k0 = d.Vectors.find("k0")->second;
    Y0t = d.Vectors.find("Y0t")->second;

//- Spécifique à chaque ordre de Tenseur - Definie dans le constructeur de base genInterfaceDomain ------------------------------------------------------
    genMatrix<double> m(3,3);
    for (int i = 0; i < 3; i++)
        m(i, i) = k0(i);
    genTensor2<double> s(m);
    genTerm<genTensor2<double>,0>::Handle InitS = genTerm<genTensor2<double>,0>::Handle(new ConstantField<genTensor2<double> >(s));
    genTerm<genTensor2<double>,0>::Handle InitdS = genTerm<genTensor2<double>,0>::Handle(new ConstantField<genTensor2<double> >(s));
    GaussQuadrature Integ_Val ( GaussQuadrature::Val );
    S = initGenTermInHistory<genTensor2<double> >("Compliance", Integ_Val,InitS);
    dS = initGenTermInHistory<genTensor2<double> >("Tangent", Integ_Val,InitdS);
  }

  virtual genTerm<genTensor2<double>,0>::Handle Tangent(const genTerm<genTensor1<double>,0>::ConstHandle &deltaForce, const genTerm<genTensor1<double>,0>::ConstHandle &Force);
  virtual genTerm<genTensor2<double>,0>::Handle deltaDamage(const genTerm<genTensor1<double>,0>::ConstHandle &deltaForce, const genTerm<genTensor1<double>,0>::ConstHandle &Force);
//  virtual genTerm<double,0>::Handle deltaDamage(const genTerm<genTensor2<double>,0>::ConstHandle &deltaStress, const genTerm<genTensor2<double>,0>::ConstHandle &Stress);
  genBitriangularDamageInterfaceDomain () : genDamageInterfaceDomain<double,double>() {}
  virtual ~genBitriangularDamageInterfaceDomain() {};
};

genTerm<genTensor2<double>,0>::Handle inline genBitriangularDamageInterfaceDomain::Tangent(const genTerm<genTensor1<double>,0>::ConstHandle &deltaDisp, const genTerm<genTensor1<double>,0>::ConstHandle &Disp)
{
  std::shared_ptr<SavedGenTerm<double,0> > dS_(new SavedGenTerm<double,0>);
//   std::shared_ptr<SavedGenTerm<double,0> > dSign(new SavedGenTerm<double,0>);

  //Paramètres lois d'endommagement
  double sigmaMax = sqrt(2*Y0t*k0);
  double epsilonMax = sigmaMax/k0;
  double epsilonCri = 2*Gic/sigmaMax;
  double k1 = -pow(sigmaMax,2)*k0/(2*Gic*k0-pow(sigmaMax,2));

  std::cout<<"   __________________________________________"<<std::endl;
  std::cout<<"  |Gic "<<Gic<<" k0 "<<k0<<" Y0t "<<Y0t<<std::endl;
  std::cout<<"  |sigmaMax "<<sigmaMax<<" epsilonMax "<<epsilonMax<<" epsilonCri "<<epsilonCri<<" k1 "<<k1<<std::endl;
  std::cout<<"  |1/k0 "<<1./k0<<" 1/k1 "<<1./k1<<std::endl;
  double epsilonEq_Max = 0;

  //genTerm<genTensor1<double>,0>::ConstHandle TForce = Compose<ContrProdOp>(Compose<ContrProdOp>(Transpose(T),Force),T);
  double h = 1.;
  genTerm<genTensor2<double>,0>::ConstHandle Strain = (Apply<TensProdOp>(Disp,N)+Apply<TensProdOp>(N,Disp))*(0.5/h);

  for (std::set<MElement*>::const_iterator it = begin(); it != end(); ++it)
  {
    MElement* e = *it;
    IntPt* GP;
    GaussQuadrature integrator(GaussQuadrature::Val);
    int npts = integrator.getIntPoints(e, &GP);

    //Déformations sur l'élément
    std::vector<genScalar<genTensor2<double> > > s_vvals(npts);
    Strain->get(e, npts, GP, s_vvals);

    //Contraintes sur l'élément
    //std::vector<genScalar<genTensor1<double> > > f_vvals(npts);
    //TForce->get(e, npts, GP, f_vvals);

    //Endommagement (modification de k)
//     std::vector<genScalar<double> > d_vvals(npts);
//     D->get(e, npts, GP, d_vvals);

    // Compliance
    std::vector<genScalar<double> > T_vvals(npts);
    S->get(e, npts, GP, T_vvals);

//     std::vector<genScalar<double> > Sign_vvals(npts);

    for (unsigned int i=0; i<npts; ++i)
    {
      //double ForceEq = f_vvals[i].norm()*eqfac;

      const genTensor2<double> &s = s_vvals[i]();
      double epsilonEq = sqrt( 0.5*( pow((s(0,0)-s(1,1)),2)
                                    +pow((s(1,1)-s(2,2)),2)
                                    +pow((s(2,2)-s(0,0)),2)
                                    +6*( pow(s(0,1),2)
                                        +pow(s(1,2),2)
                                        +pow(s(2,0),2))) ) * eqfac;


      if(epsilonEq > epsilonEq_Max) epsilonEq_Max = epsilonEq;

      if(epsilonEq<=epsilonMax)
      {
//         Sign_vvals[i] = 1.;
        T_vvals[i] = T_vvals[i];
      }
      else if(epsilonEq>epsilonMax)
      {
        if (it == begin()) std::cout<<"  |-*-*-*-*-*-*-*-*-*-*-DAMAGE-*-*-*-*-*-*-*-*-*- "<<std::endl;
//         Sign_vvals[i] = -1.;
        if (epsilonEq<epsilonCri)
        {
          T_vvals[i] = 1./k1;
        }
        else
        {
          T_vvals[i] = 1.E10;
//           std::cout<<"  |WARNING full Damage in element "<< e->getNum()
//                    <<"  | for the GP "<< i
//                    <<"  | with the value (epsilonEq-epsilonCri)/epsilonCri="<< ((epsilonEq-epsilonCri)/epsilonCri*100) <<"%"
//                    <<std::endl;
        }
      }
    }
    dS_->set(e, npts, GP, T_vvals);

//     dSign->set(e, npts, GP, Sign_vvals);

  }
  std::cout<<"  |epsilonEq_Max "<<epsilonEq_Max<<std::endl;

  return dS_;
}

genTerm<genTensor2<double>,0>::Handle inline genBitriangularDamageInterfaceDomain::deltaDamage(const genTerm<genTensor1<double>,0>::ConstHandle &deltaDisp, const genTerm<genTensor1<double>,0>::ConstHandle &Disp)
{
  std::shared_ptr<SavedGenTerm<double,0> > dD(new SavedGenTerm<double,0>);

  //Paramètres lois d'endommagement
  double sigmaMax = sqrt(2*Y0t*k0);
  double epsilonMax = sigmaMax/k0;
  double epsilonCri = 2*Gic/sigmaMax;
  double k1 = pow(sigmaMax,2)*k0/(2*Gic*k0-pow(sigmaMax,2));

  //genTerm<genTensor1<double>,0>::ConstHandle TForce = Compose<ContrProdOp>(Compose<ContrProdOp>(Transpose(T),Force),T);
  double h = 1.;
  genTerm<genTensor2<double>,0>::ConstHandle Strain = Apply<TensProdOp>(Disp,N)*(1./h);

  for (std::set<MElement*>::const_iterator it = begin(); it != end(); ++it)
  {
    MElement* e = *it;
    IntPt* GP;
    GaussQuadrature integrator(GaussQuadrature::Val);
    int npts = integrator.getIntPoints(e, &GP);

    //Déformations sur l'élément
    std::vector<genScalar<genTensor2<double> > > s_vvals(npts);
    Strain->get(e, npts, GP, s_vvals);

    //Contraintes sur l'élément
    //std::vector<genScalar<genTensor1<double> > > f_vvals(npts);
    //TForce->get(e, npts, GP, f_vvals);

    //Endommagement (modification de k)
    std::vector<genScalar<double> > d_vvals(npts);
    D->get(e, npts, GP, d_vvals);
    for (unsigned int i=0; i<npts; ++i)
    {
      //double ForceEq = f_vvals[i].norm()*eqfac;
      const genTensor2<double> &s = s_vvals[i]();
      double epsilonEq = sqrt( 0.5*( pow((s(0,0)-s(1,1)),2)
                                    +pow((s(1,1)-s(2,2)),2)
                                    +pow((s(2,2)-s(0,0)),2)
                                    +6*( pow(s(0,1),2)
                                        +pow(s(1,2),2)
                                        +pow(s(2,0),2))) ) * eqfac;

      if(epsilonEq<=epsilonMax)
      {
        d_vvals[i] = 0;
      }
      else if(epsilonEq>epsilonMax)
      {
        //double epsilonEq = epsilonMax+fabs(sigmaMax-ForceEq)/k1;
        if (epsilonEq<epsilonCri)
        {
          double newD = (1+k1/k0)*(1-sigmaMax/(k0*epsilonEq));
          d_vvals[i] = newD - d_vvals[i];
        }
        else
        {
          d_vvals[i] = 0.99999;
          std::cout<<"WARNING full Damage in element "<< e->getNum()
                   <<" for the GP "<< i
                   <<" with the value (epsilonEq-epsilonCri)/epsilonCri="<< ((epsilonEq-epsilonCri)/epsilonCri*100) <<"%"
                   <<std::endl;
        }
      }
    }
    dD->set(e, npts, GP, d_vvals);
  }
  return dD;
}
*/


/*genTerm<double,0>::Handle inline genBitriangularDamageInterfaceDomain::deltaDamage(const genTerm<genTensor2<double>,0>::ConstHandle &deltaStress, const genTerm<genTensor2<double>,0>::ConstHandle &Stress)
{
  std::shared_ptr<SavedGenTerm<double,0> > dD(new SavedGenTerm<double,0>);

  //Paramètres lois d'endommagement
  double sigmaMax = sqrt(2*Y0t*k0);
  double epsilonMax = sqrt(2*Y0t/k0);//=sigmaMax/k0;
  double k1 = pow(sigmaMax,2)*k0/(2*Gic*k0-pow(sigmaMax,2));

  genTerm<genTensor2<double>,0>::ConstHandle StressB = Compose<ContrProdOp>(Compose<ContrProdOp>(Transpose(T),Stress),T);

  for (std::set<MElement*>::const_iterator it = begin(); it != end(); ++it)
  {
    MElement* e = *it;
    IntPt* GP;
    GaussQuadrature integrator(GaussQuadrature::Val);
    int npts = integrator.getIntPoints(e, &GP);

    //Contraintes sur l'élément
    std::vector<genScalar<genTensor2<double> > > s_vvals(npts);
    StressB->get(e, npts, GP, s_vvals);

    //Endommagement (modification de k)
    std::vector<genScalar<double> > d_vvals(npts);
    D->get(e, npts, GP, d_vvals);
    for (unsigned int i=0; i<npts; ++i)
    {
      double sigmaEq = sqrt( 0.5*( pow((s_vvals[i](0,0)-s_vvals[i](1,1)),2)
                                       +pow((s_vvals[i](1,1)-s_vvals[i](2,2)),2)
                                       +pow((s_vvals[i](2,2)-s_vvals[i](0,0)),2)
                                       +6*( pow(s_vvals[i](0,1),2)
                                           +pow(s_vvals[i](1,2),2)
                                           +pow(s_vvals[i](2,0),2))) );
      if(d_vvals[i]==0 && sigmaEq<=sigmaMax)
      {
        d_vvals[i] = 0;
      }
      else if(sigmaEq>epsilonMax)
      {
        double epsilonEq=(sigmaMax-sigmaEq)/k1+sigmaMax/k0;
        double newD=(1+k1/k0)*(1-sigmaMax/(k0*epsilonEq));
        d_vvals[i] = newD - d_vvals[i];
      }
    }
    dD->set(e, npts, GP, d_vvals);
  }
  return dD;
}
*/


#endif// _NONLINEARXFEM_DOMAIN_H_