// elastic_genTerm_dfloat - A solver for non linear elastic 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>.

#ifndef _ELASTIC_DOMAIN_DFLOAT_H_
#define _ELASTIC_DOMAIN_DFLOAT_H_

#include "mechanicsDomain.h"


class ElasticDomainDfloat : public MechanicsDomain
{
public :
  genTerm<genTensor4<double >,0>::ConstHandle C; //elastic tensor
  genTensor4<double> elast;
//  here , C isnt constant anymore. C=C_0( 1+k*||grad u||^2) to test.
public :
  ElasticDomainDfloat (const genDomain &s) : MechanicsDomain(s)
  {
    const genMatrix<double> & mC = s.Tensors.find("C")->second;
    FromVoigtNotation(mC,elast);
    C = genTerm<genTensor4<double>,0>::Handle(new ConstantField<genTensor4<double> >(elast));
  }
  ElasticDomainDfloat() : MechanicsDomain() {}
  virtual ~ElasticDomainDfloat() {}

  Macro_Stress
};

template <class scalar, int n> typename genTerm<genTensor2<scalar>,n>::ConstHandle ElasticDomainDfloat::UserStress(const typename genTerm<genTensor2<scalar>,n>::ConstHandle &Eps)
{
  typename genTerm<genTensor4<scalar>,0>::ConstHandle scalarC=Build<genTensor4<scalar> >(C);
  return Compose<FullContrProdOp>(Eps,scalarC);
}


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

public :
  IsotropicElasticDomainDfloat(const genDomain &s) : ElasticDomainDfloat(s)
  {
    const double & E_=s.Scalars.find("E")->second;
    const double & nu_=s.Scalars.find("nu")->second;
    E = E_; nu = nu_;
    update();
  }
  IsotropicElasticDomainDfloat () : ElasticDomainDfloat() {}
  virtual ~IsotropicElasticDomainDfloat() {}
  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;
    FromVoigtNotation(mC,elast);
    C = genTerm<genTensor4<double>,0>::Handle(new ConstantField<genTensor4<double> >(elast));
  }
};

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

public :
  OrthotropicElasticDomainDfloat(const genDomain &s) : ElasticDomainDfloat(s)
  {
    const std::vector<double> & E_=s.Vectors.find("E")->second;
    const std::vector<double> & nu_=s.Vectors.find("nu")->second;
    const std::vector<double> & G_=s.Vectors.find("G")->second;
    E = E_; nu = nu_; G = G_;
    update();
  }
  OrthotropicElasticDomainDfloat () :  ElasticDomainDfloat() {}
  ~OrthotropicElasticDomainDfloat () {}
  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];

    FromVoigtNotation(mC,elast);

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





class IsotropicPerfectPlasticDomainDfloat : public IsotropicElasticDomainDfloat
{
public :
  double YieldStress;
  double Hardening;
  genTerm<genTensor2<double>,0>::ConstHandle NomEps;

public :
  IsotropicPerfectPlasticDomainDfloat(const genDomain &s) : IsotropicElasticDomainDfloat(s)
  {
    YieldStress=s.Scalars.find("YieldStress")->second;
    Hardening=s.Scalars.find("Hardening")->second;
    NomEps = genTerm<genTensor2<double>,0>::ConstHandle(new ConstantField<genTensor2<double> >(genTensor2<double>(0.)));
    UpdateHistory(NomEps);
  }
  IsotropicPerfectPlasticDomainDfloat() : IsotropicElasticDomainDfloat() {}
  virtual ~IsotropicPerfectPlasticDomainDfloat() {}

  Macro_Stress
  Macro_History
};


template<class T1, class T2=T1 > class PlasticStressOp
{
  public :
};

template<class scalar, int N> class PlasticStressOp<genTensor2<scalar,N> > // computes the stress in function of epsilon
{
public :
  typedef genTensor2<scalar,N> ValType;
protected :
  double YieldStress;
  double Hardening;
  genTensor4<scalar,N> elast;
public :
  PlasticStressOp(const double &YieldStress_, const double &Hardening_, const genTensor4<double,N> &elast_) : YieldStress(YieldStress_), Hardening(Hardening_)
  {
    for (int i=0;i<N;++i)
      for (int j=0;j<N;++j)
        for (int k=0;k<N;++k)
          for (int l=0;l<N;++l)
            elast(i,j,k,l)=scalar(elast_(i,j,k,l));
  };
  void operator ()(const genTensor2<scalar,N> &strain, const genTensor2<scalar,N> &nominalstrain, ValType &stress) const
  {
    ValType nstress=elast*nominalstrain;
    stress=elast*strain;
    scalar VM(0.0);
    scalar VM2=nstress*nstress*1.5;
    if (VM2!=0.0) VM=sqrt(VM2);
    if (VM>scalar(YieldStress))
    {
      scalar coeff=(scalar(YieldStress)+(VM-scalar(YieldStress))*Hardening)/VM;
      stress*=coeff;
    }
  }
};


template <class scalar, int n> typename genTerm<genTensor2<scalar>,n>::ConstHandle IsotropicPerfectPlasticDomainDfloat::UserStress(const typename genTerm<genTensor2<scalar>,n>::ConstHandle &Eps)
{
  PlasticStressOp<genTensor2<scalar> > Pop(YieldStress,Hardening,elast);
  typename genTerm<genTensor2<scalar>,0>::ConstHandle scalarNomEps=Build<genTensor2<scalar> >(NomEps);
  return Apply(Eps,scalarNomEps,Pop);
}


void inline IsotropicPerfectPlasticDomainDfloat::CopyHistory()
{
  int step = History.getcurrentstep("NomEps");
  std::shared_ptr<SavedGenTerm<genTensor2<double>,0> > CurrentNomEps;
  History.getSavedGenTerm("NomEps", CurrentNomEps, step);
  History.addSavedGenTerm("NomEps", CurrentNomEps, step+1);
}

void inline IsotropicPerfectPlasticDomainDfloat::UpdateHistory(const genTerm<genTensor2<double>,0>::ConstHandle &Eps)
{
  genTerm<genTensor2<double>,0>::ConstHandle newNomEps = NomEps+Eps;
  GaussQuadrature Integ_Grad ( GaussQuadrature::Grad );
  NomEps = saveGenTermInHistory<genTensor2<double> >("NomEps", Integ_Grad, newNomEps);
}

#endif// _ELASTIC_DOMAIN_DFLOAT_H_