// Nonlinear_genTerm - A solver for nonlinear problems using 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 _NONLINEAR_DOMAIN_H_
#define _NONLINEAR_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"


class NLDomain : public genSupport
{
public :
  SavedGenTermManager History;

public :
  NLDomain() : genSupport() {}
  virtual ~NLDomain() {};

protected :
  template <class T> typename genTerm<T,0>::Handle initGenTermInHistory(std::string name, QuadratureBase &integrator, const typename genTerm<T,0>::ConstHandle &gt);
  template <class T> typename genTerm<T,0>::Handle saveGenTermInHistory(std::string name, QuadratureBase &integrator, const typename genTerm<T,0>::ConstHandle &gt);
public :
  virtual void copyHistory() =0;
//   virtual void saveHistory(const genTerm<genTensor2<double> ,0>::ConstHandle &deltaStrain, const genTerm<genTensor2<double> ,0>::ConstHandle &Strain) =0;
//   virtual void saveHistory(const genTerm<genTensor1<double> ,0>::ConstHandle &deltaForce, const genTerm<genTensor1<double> ,0>::ConstHandle &Force) =0;
};


template <class T> typename genTerm<T,0>::Handle NLDomain::initGenTermInHistory(std::string name, QuadratureBase &integrator, const typename genTerm<T,0>::ConstHandle &gt)
{
  std::shared_ptr<SavedGenTerm<T,0> > InitGt(new SavedGenTerm<T,0>());
  SavePtGauss(*gt, begin(), end(), integrator, *InitGt);
  History.addSavedGenTerm(name, InitGt, 0);
  return InitGt;
}

template <class T> typename genTerm<T,0>::Handle NLDomain::saveGenTermInHistory(std::string name, QuadratureBase &integrator, const typename genTerm<T,0>::ConstHandle &gt)
{
  int step = History.getcurrentstep(name);
  std::shared_ptr<SavedGenTerm<T,0> > SaveGt(new SavedGenTerm<T,0>());
  SavePtGauss(*gt, begin(), end(), integrator, *SaveGt);
  History.addSavedGenTerm(name, SaveGt, step);

  std::ostringstream temp;
  temp << "data_" << name  << ".txt";
  SaveGt->printfile(temp.str());

  return SaveGt;
}
  


class NLElasticDomain : public NLDomain
{
public :
  genTerm<genTensor4<double> ,0>::Handle C; //elastic tensor
  genTerm<genTensor2<double>,0>::Handle S; //internal stress tensor

public :
  NLElasticDomain(const genDomain &d)
  {
    dim=d.dim;
    tag=d.tag;
    const genMatrix<double> & mC=d.Tensors.find("C")->second;
    genTensor4<double>  elast;
    FromVoigtNotation(mC,elast);
    C = genTerm<genTensor4<double> ,0>::Handle(new ConstantField<genTensor4<double> >(elast));

//    History.defineMaxSteps(2);

    genTerm<genTensor2<double>,0>::Handle InitStress = genTerm<genTensor2<double>,0>::Handle(new ConstantField<genTensor2<double> >(genTensor2<double>(0.)));
    GaussQuadrature Integ_Grad ( GaussQuadrature::Grad );
    S = initGenTermInHistory<genTensor2<double> >("Stress", Integ_Grad, InitStress);
  }
  template <int n> typename genTerm<genTensor2<double>,n>::Handle deltaStress(const typename genTerm<genTensor2<double>,n>::ConstHandle &deltaStrain);
  genTerm<genTensor2<double>,0>::Handle Stress();
  virtual genTerm<genTensor4<double> ,0>::Handle Tangent();
  template <int n> typename genTerm<genTensor0<double>,n>::Handle Force(const typename genTerm<genTensor2<double>,n>::ConstHandle &Strain);
  template <int n1, int n2> typename genTerm<genTensor0<double>, n1+n2>::Handle
  deltaEnergy( const typename genTerm<genTensor2<double>,n1>::ConstHandle &Strain1, const typename genTerm<genTensor2<double>,n2>::ConstHandle &Strain2 );
  NLElasticDomain() : NLDomain() {}
  virtual ~NLElasticDomain() {};
  
public :
  virtual void copyHistory();
  virtual void saveHistory(const genTerm<genTensor2<double> ,0>::ConstHandle &deltaStrain, const genTerm<genTensor2<double> ,0>::ConstHandle &Strain);
};

template <int n> typename genTerm<genTensor2<double>,n>::Handle NLElasticDomain::deltaStress(const typename genTerm<genTensor2<double>,n>::ConstHandle &deltaStrain)
{
  return typename genTerm<genTensor2<double>,n>::Handle(Compose<FullContrProdOp>(Tangent(),deltaStrain));
}

genTerm<genTensor2<double>,0>::Handle inline NLElasticDomain::Stress()
{
  return genTerm<genTensor2<double>,0>::Handle(S);
}

genTerm<genTensor4<double> ,0>::Handle inline NLElasticDomain::Tangent()
{
  return genTerm<genTensor4<double> ,0>::Handle(C);
}

template <int n> typename genTerm<genTensor0<double>,n>::Handle NLElasticDomain::Force(const typename genTerm<genTensor2<double>,n>::ConstHandle &Strain)
{
  return typename genTerm<genTensor0<double>,n>::Handle(Compose<FullContrProdOp>(Transpose(Strain),S));
}

template <int n1, int n2> typename genTerm<genTensor0<double>,n1+n2>::Handle
NLElasticDomain::deltaEnergy( const typename genTerm<genTensor2<double>,n1>::ConstHandle &deltaStrain1, const typename genTerm<genTensor2<double>,n2>::ConstHandle &deltaStrain2)
{
  return typename genTerm<genTensor0<double>,n1+n2>::Handle (Compose<FullContrProdOp>(Transpose(deltaStrain1),deltaStress<n2>(deltaStrain2)));
}

void inline NLElasticDomain::copyHistory()
{
  int step = History.getcurrentstep("Stress");
  std::shared_ptr<SavedGenTerm<genTensor2<double> ,0> > OldStress;
  History.getSavedGenTerm("Stress", OldStress, step);
  History.addSavedGenTerm("Stress", OldStress, step+1);
}

void inline NLElasticDomain::saveHistory(const genTerm<genTensor2<double> ,0>::ConstHandle &deltaStrain, const genTerm<genTensor2<double> ,0>::ConstHandle &Strain)
{
  //Save Stress
  genTerm<genTensor2<double>,0>::Handle newSress = Stress()+deltaStress<0>(deltaStrain);
  GaussQuadrature Integ_Grad ( GaussQuadrature::Grad );
  S = saveGenTermInHistory<genTensor2<double> >("Stress", Integ_Grad, newSress);
}
  


class IsotropicNLElasticDomain : public NLElasticDomain
{
public :
  double E, nu; // specific elastic datas (isotropic)

public :
  IsotropicNLElasticDomain(const genDomain &d)
  {
    dim=d.dim;
    tag=d.tag;
    const double & E_=d.Scalars.find("E")->second;
    const double & nu_=d.Scalars.find("nu")->second;
    E = E_; nu = nu_;
    update();

    genTerm<genTensor2<double>,0>::Handle InitStress = genTerm<genTensor2<double>,0>::Handle(new ConstantField<genTensor2<double> >(genTensor2<double>(0.)));
    GaussQuadrature Integ_Grad ( GaussQuadrature::Grad );
    S = initGenTermInHistory<genTensor2<double> >("Stress", Integ_Grad, InitStress);
  }
  IsotropicNLElasticDomain () : NLElasticDomain() {}
  ~IsotropicNLElasticDomain () {}
  void update ()
  {
    double FACT = E / (1 + nu);
    double C11 = FACT * (1 - nu) / (1 - 2 * nu);
    double C12 = FACT * nu / (1 - 2 * nu);
    double C44 = (C11 - C12) / 2;
//       FACT = E / (1 - nu * nu); // plane stress (plates)
//       C11  = FACT;
//       C12  = nu * FACT;
//       C44 = (1. - nu) * .5 * FACT;
    genMatrix<double> mC(6,6);
    for(int i = 0; i < 3; ++i)
    {
      mC(i,i) = C11;
      mC(i+3,i+3) = C44;
    }
    mC(1,0) = mC(0,1) = mC(2,0) = mC(0,2) = mC(1,2) = mC(2,1) = C12;

    genTensor4<double>  elast;
    FromVoigtNotation(mC,elast);

    C = genTerm<genTensor4<double> ,0>::Handle(new ConstantField<genTensor4<double> >(elast));
  }
};

class OrthotropicNLElasticDomain : public NLElasticDomain
{
public :
  std::vector<double> E, nu, G; // E1,E2,E3, nu12,nu13,nu23, G12,G13,G23 (orthotropic_3D)

public :
  OrthotropicNLElasticDomain(const genDomain &d)
  {
    dim=d.dim;
    tag=d.tag;
    const std::vector<double> & E_=d.Vectors.find("E")->second;
    const std::vector<double> & nu_=d.Vectors.find("nu")->second;
    const std::vector<double> & G_=d.Vectors.find("G")->second;
    E = E_; nu = nu_; G = G_;
    update();

    genTerm<genTensor2<double>,0>::Handle InitStress = genTerm<genTensor2<double>,0>::Handle(new ConstantField<genTensor2<double> >(genTensor2<double>(0.)));
    GaussQuadrature Integ_Grad ( GaussQuadrature::Grad );
    S = initGenTermInHistory<genTensor2<double> >("Stress", Integ_Grad, InitStress);
  }
  OrthotropicNLElasticDomain () :  NLElasticDomain() {}
  ~OrthotropicNLElasticDomain () {}
  void update ()
  {
    double S11 = 1. / E[0];
    double S22 = 1. / E[1];
    double S33 = 1. / E[2];
    double S12 = - nu[0] / E[0];
    double S13 = - nu[1] / E[1];
    double S23 = - nu[2] / E[2];

    double dS = S11*S22*S33 - S11*S23*S23 - S22*S13*S13 - S33*S12*S12 + 2*S12*S23*S13;

    genMatrix<double> mC(6,6);
    mC(0,0) = (S22*S33 - S23*S23)/dS;
    mC(1,1) = (S33*S11 - S13*S13)/dS;
    mC(2,2) = (S11*S22 - S12*S12)/dS;
    mC(0,1) = mC(1,0) = (S13*S23 - S12*S33)/dS;
    mC(0,2) = mC(2,0) = (S12*S23 - S13*S22)/dS;
    mC(1,2) = mC(2,1) = (S12*S13 - S23*S11)/dS;
    mC(3,3) = G[2];
    mC(4,4) = G[1];
    mC(5,5) = G[0];

    genTensor4<double>  elast;
    FromVoigtNotation(mC,elast);

    C = genTerm<genTensor4<double> ,0>::Handle(new ConstantField<genTensor4<double> >(elast));
  }
};



class DamageDomain : public NLElasticDomain
{
public :
  genTerm<genTensor0<double>,0>::Handle D;

public :
  DamageDomain(const genDomain &d) : NLElasticDomain(d)
  {
    genTerm<genTensor0<double>,0>::Handle InitDamage = genTerm<genTensor0<double>,0>::Handle(new ConstantField<genTensor0<double> >(0.));
    GaussQuadrature Integ_Val ( GaussQuadrature::Val );
    D = initGenTermInHistory<genTensor0<double> >("Damage", Integ_Val,InitDamage);
  }
  genTerm<genTensor4<double> ,0>::Handle Tangent();
  virtual genTerm<genTensor0<double>,0>::Handle deltaDamage(const genTerm<genTensor2<double> ,0>::ConstHandle &deltaStrain, const genTerm<genTensor2<double> ,0>::ConstHandle &Strain);
  genTerm<genTensor0<double>,0>::Handle Damage();
  DamageDomain () : NLElasticDomain() {}
  ~DamageDomain () {}

  virtual void copyHistory();
  virtual void saveHistory(const genTerm<genTensor2<double> ,0>::ConstHandle &deltaStrain, const genTerm<genTensor2<double> ,0>::ConstHandle &Strain);
};

genTerm<genTensor4<double> ,0>::Handle inline DamageDomain::Tangent()
{
  genTerm<genTensor0<double>,0>::Handle One = genTerm<genTensor0<double>,0>::Handle(new ConstantField<genTensor0<double> > (1.));
  genTerm<genTensor0<double>,0>::Handle State = One+(D*-1.);
  return genTerm<genTensor4<double> ,0>::Handle (Apply<TensProdOp>(C,State));
}

genTerm<genTensor0<double>,0>::Handle inline DamageDomain::deltaDamage(const genTerm<genTensor2<double> ,0>::ConstHandle &deltaStrain, const genTerm<genTensor2<double> ,0>::ConstHandle &Strain)
{
  genTerm<genTensor0<double>,0>::Handle Zero = genTerm<genTensor0<double>,0>::Handle(new ConstantField<genTensor0<double> > (0.));
  return genTerm<genTensor0<double>,0>::Handle(Zero);
}

genTerm<genTensor0<double>,0>::Handle inline DamageDomain::Damage()
{
  return genTerm<genTensor0<double>,0>::Handle(D);
}

void inline DamageDomain::copyHistory()
{
  int step = History.getcurrentstep("Stress");
  std::shared_ptr<SavedGenTerm<genTensor2<double> ,0> > OldStress;
  History.getSavedGenTerm("Stress", OldStress, step);
  History.addSavedGenTerm("Stress", OldStress, step+1);

  std::shared_ptr<SavedGenTerm<genTensor0<double>,0> > OldDamage;
  History.getSavedGenTerm("Damage", OldDamage, step);
  History.addSavedGenTerm("Damage", OldDamage, step+1);
}

void inline DamageDomain::saveHistory(const genTerm<genTensor2<double> ,0>::ConstHandle &deltaStrain, const genTerm<genTensor2<double> ,0>::ConstHandle &Strain)
{
  //Save damage
  genTerm<genTensor0<double>,0>::Handle oldDamage = Damage();
  genTerm<genTensor0<double>,0>::Handle dDamage = deltaDamage(deltaStrain,Strain);
  genTerm<genTensor0<double>,0>::Handle newDamage = oldDamage+dDamage;
  GaussQuadrature Integ_Val ( GaussQuadrature::Val );
  D = saveGenTermInHistory<genTensor0<double> >("Damage", Integ_Val, newDamage);
  
  //Save Stress
  genTerm<genTensor2<double> ,0>::Handle oldSress = Stress();
  genTerm<genTensor2<double> ,0>::Handle dSress = deltaStress<0>(deltaStrain);
  genTerm<genTensor2<double> ,0>::Handle newSress = oldSress+dSress;
  GaussQuadrature Integ_Grad ( GaussQuadrature::Grad );
  S = saveGenTermInHistory<genTensor2<double> >("Stress", Integ_Grad, newSress);
}



class BitriangularDamageDomain : public DamageDomain
{
public :
  //Données liés à la loi d'évolution bitriangulaire
  double Gic;
  double k0;
  double Y0t;

public :
  BitriangularDamageDomain(const genDomain &d) : DamageDomain(d)
  {
    Gic=d.Scalars.find("Gic")->second;
    k0=d.Scalars.find("k0")->second;
    Y0t=d.Scalars.find("Y0t")->second;
  }
  virtual genTerm<genTensor0<double>,0>::Handle deltaDamage(const genTerm<genTensor2<double> ,0>::ConstHandle &deltaStrain, const genTerm<genTensor2<double> ,0>::ConstHandle &Strain);
  BitriangularDamageDomain () : DamageDomain() {}
  virtual ~BitriangularDamageDomain() {};
};

genTerm<genTensor0<double>,0>::Handle inline BitriangularDamageDomain::deltaDamage(const genTerm<genTensor2<double> ,0>::ConstHandle &deltaStrain, const genTerm<genTensor2<double> ,0>::ConstHandle &Strain)
{
  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 = sqrt(2*Y0t/k0);//=sigmaMax/k0;
  double k1 = pow(sigmaMax,2)*k0/(2*Gic*k0-pow(sigmaMax,2));

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

    //Deformation 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)
    {
      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))) );
      if(epsilonEq<=epsilonMax)
      {
        d_vvals[i] = genTensor0<double>(0.);
      }
      else if(epsilonEq>epsilonMax)
      {
        double newD=(1+k1/k0)*(1-sigmaMax/(k0*epsilonEq));
        d_vvals[i] = genTensor0<double>(newD) - d_vvals[i];
      }
    }
    dD->set(e, npts, GP, d_vvals);
  }
  return dD;
}

#endif// _NONLINEAR_DOMAIN_H_