// C++ Interface: material law
//
// Description: general class for nonlocal porous material law
//
// Author:  <V.D. Nguyen>, (C) 2018
//
// Copyright: See COPYING file that comes with this distribution
//
//

#include "mlawNonLocalPorous.h"
#include "STensorOperations.h"

mlawNonLocalPorosity::mlawNonLocalPorosity(const int num,const double E,const double nu, const double rho, 
                           const double fVinitial, const J2IsotropicHardening &j2IH,
                           const double tol, const bool matrixbyPerturbation, const double pert) : 
                                               materialLaw(num,true),
                                               _E(E), _nu(nu), _rho(rho),
                                               _lambda((E*nu)/(1.+nu)/(1.-2.*nu)),
                                               _mu(0.5*E/(1.+nu)),_K(E/3./(1.-2.*nu)),
                                               _K3(3.*_K), _mu3(3.*_mu),
                                               _mu2(2.*_mu),
                                               _orderlogexp(1), _tol(tol),
                                               _perturbationfactor(pert),
                                               _tangentByPerturbation(matrixbyPerturbation ),
                                               _fVinitial(fVinitial), _randMaxbound(1.), _randMinbound(1.),
                                               _I4(1.,1.), _I(1.),_withSubstepping(false),_maxAttemptSubstepping(1)
{
  // Internal laws
  _j2IH  = j2IH.clone();
  _viscoplastlaw = NULL;
  
  _cLLaw = new ZeroCLengthLaw(num);
  _gdnLawContainer.clear();
  _coalescenceLaw = new NoCoalescenceLaw(num);

  // Fill elastic tensor
  Cel*=0.;
  Cel(0,0,0,0) = _lambda + _mu2;
  Cel(1,1,0,0) = _lambda;
  Cel(2,2,0,0) = _lambda;
  Cel(0,0,1,1) = _lambda;
  Cel(1,1,1,1) = _lambda + _mu2;
  Cel(2,2,1,1) = _lambda;
  Cel(0,0,2,2) = _lambda;
  Cel(1,1,2,2) = _lambda;
  Cel(2,2,2,2) = _lambda + _mu2;

  Cel(1,0,1,0) = _mu;
  Cel(2,0,2,0) = _mu;
  Cel(0,1,0,1) = _mu;
  Cel(2,1,2,1) = _mu;
  Cel(0,2,0,2) = _mu;
  Cel(1,2,1,2) = _mu;

  Cel(0,1,1,0) = _mu;
  Cel(0,2,2,0) = _mu;
  Cel(1,0,0,1) = _mu;
  Cel(1,2,2,1) = _mu;
  Cel(2,0,0,2) = _mu;
  Cel(2,1,1,2) = _mu;
    
  // Initialise tools
  _Idev = _I4;
  STensor3 mIon3(-1./3);
  STensor43 mIIon3;
  tensprod(_I,mIon3, mIIon3);
  _Idev += mIIon3;
}


mlawNonLocalPorosity::mlawNonLocalPorosity(const int num,const double E,const double nu, const double rho, 
                           const double fVinitial, const J2IsotropicHardening &j2IH, const CLengthLaw &cLLaw,
                           const double tol, const bool matrixbyPerturbation, const double pert) : 
                                               materialLaw(num,true),
                                               _E(E), _nu(nu), _rho(rho),
                                               _lambda((E*nu)/(1.+nu)/(1.-2.*nu)),
                                               _mu(0.5*E/(1.+nu)),_K(E/3./(1.-2.*nu)),
                                               _K3(3.*_K), _mu3(3.*_mu),
                                               _mu2(2.*_mu),
                                               _orderlogexp(1), _tol(tol),
                                               _perturbationfactor(pert),
                                               _tangentByPerturbation(matrixbyPerturbation ),
                                               _fVinitial(fVinitial), _randMaxbound(1.), _randMinbound(1.),
                                               _I4(1.,1.), _I(1.),_withSubstepping(false),_maxAttemptSubstepping(1)
{
  // Internal laws
  _j2IH  = j2IH.clone();
  _viscoplastlaw = NULL;
  
  _cLLaw = cLLaw.clone();
  _gdnLawContainer.clear();
  _coalescenceLaw = new NoCoalescenceLaw(num);

  // Fill elastic tensor
  Cel=0.;
  Cel(0,0,0,0) = _lambda + _mu2;
  Cel(1,1,0,0) = _lambda;
  Cel(2,2,0,0) = _lambda;
  Cel(0,0,1,1) = _lambda;
  Cel(1,1,1,1) = _lambda + _mu2;
  Cel(2,2,1,1) = _lambda;
  Cel(0,0,2,2) = _lambda;
  Cel(1,1,2,2) = _lambda;
  Cel(2,2,2,2) = _lambda + _mu2;

  Cel(1,0,1,0) = _mu;
  Cel(2,0,2,0) = _mu;
  Cel(0,1,0,1) = _mu;
  Cel(2,1,2,1) = _mu;
  Cel(0,2,0,2) = _mu;
  Cel(1,2,1,2) = _mu;

  Cel(0,1,1,0) = _mu;
  Cel(0,2,2,0) = _mu;
  Cel(1,0,0,1) = _mu;
  Cel(1,2,2,1) = _mu;
  Cel(2,0,0,2) = _mu;
  Cel(2,1,1,2) = _mu;
    
  // Initialise tools
  _Idev = _I4;
  STensor3 mIon3(-1./3);
  STensor43 mIIon3;
  tensprod(_I,mIon3, mIIon3);
  _Idev += mIIon3;
}

mlawNonLocalPorosity::mlawNonLocalPorosity(const mlawNonLocalPorosity &source) : 
                                                         materialLaw(source),
                                                         _E(source._E), _nu(source._nu), _rho(source._rho), 
                                                         _lambda(source._lambda),
                                                         _mu(source._mu), _K(source._K),
                                                         _K3(source._K3), _mu3(source._mu3),
                                                         _mu2(source._mu2),
                                                         _orderlogexp(source._orderlogexp), _tol(source._tol),
                                                         _perturbationfactor(source._perturbationfactor),
                                                         _tangentByPerturbation(source._tangentByPerturbation),
                                                         _fVinitial(source._fVinitial), _randMaxbound(source._randMaxbound),
                                                         _randMinbound(source._randMinbound),
                                                         Cel(source.Cel), _I4(source._I4), _I(source._I), _Idev(source._Idev),
                                                         _withSubstepping(source._withSubstepping),
                                                         _maxAttemptSubstepping(source._maxAttemptSubstepping)


{
  // Internal laws
    // Isotropic hardening
  _j2IH = NULL;
  if(source._j2IH != NULL) _j2IH=source._j2IH->clone();
    // Viscoplastic law
  _viscoplastlaw = NULL;
  if(source._viscoplastlaw != NULL) _viscoplastlaw = source._viscoplastlaw->clone();
    // Non-local lenght
  _cLLaw = NULL;
  if(source._cLLaw != NULL) _cLLaw=source._cLLaw->clone();
    // Nucleation law vector
  _gdnLawContainer.resize(source._gdnLawContainer.size());
  for(int i=0; i < source._gdnLawContainer.size(); i++)
  {
    if(source._gdnLawContainer[i] != NULL){_gdnLawContainer[i] = source._gdnLawContainer[i]->clone();};
  }
    // Coalescence law
  _coalescenceLaw = NULL;
  if(source._coalescenceLaw != NULL) _coalescenceLaw=source._coalescenceLaw->clone();
}

mlawNonLocalPorosity::~mlawNonLocalPorosity() 
{ 
  if(_j2IH!=NULL) delete _j2IH; _j2IH = NULL;
  
  if (_viscoplastlaw != NULL) delete _viscoplastlaw; _viscoplastlaw = NULL;   
  
  if(_cLLaw!= NULL) delete _cLLaw; _cLLaw = NULL;
  
  for(int i=0; i < _gdnLawContainer.size(); i++)
  {
    if(_gdnLawContainer[i] != NULL){delete _gdnLawContainer[i]; _gdnLawContainer[i] = NULL;};
  }
  
  if(_coalescenceLaw != NULL) delete _coalescenceLaw;
}


/* Option settings */
void mlawNonLocalPorosity::setOrderForLogExp(const int newOrder)
{
  _orderlogexp = newOrder;
  Msg::Info("mlawNonLocalPorosity :: order = %d is used to approximate log and exp");
};


void mlawNonLocalPorosity::setViscosityLaw(const viscosityLaw& vico)
{
  if (_viscoplastlaw != NULL) delete _viscoplastlaw;
  _viscoplastlaw = vico.clone();
};

void mlawNonLocalPorosity::setNucleationLaw(const NucleationLaw& added_gdnLaw)
{   
  _gdnLawContainer.push_back(added_gdnLaw.clone());
};

void mlawNonLocalPorosity::setScatterredInitialPorosity(double f0min, double f0max)
{
  Msg::Info("porous law :: Add scatter in initial porosity");
  _randMinbound = f0min;
  if(f0min < 0.){Msg::Fatal("mlawNonLocalPorosity :: Wrong scatter parameter fV0_min");}
  _randMaxbound = f0max;
  if(f0max < f0min){Msg::Fatal("mlawNonLocalPorosity :: Wrong scatter parameter fV0_max");}
};

void mlawNonLocalPorosity::setCoalescenceLaw(const CoalescenceLaw& added_coalsLaw)
{
  if (_coalescenceLaw != NULL) delete _coalescenceLaw;
  _coalescenceLaw = added_coalsLaw.clone();
};

double mlawNonLocalPorosity::soundSpeed() const
{
  double factornu = (1.-_nu)/((1.+_nu)*(1.-2.*_nu));
  return sqrt(_E*factornu/_rho);
}

  
double mlawNonLocalPorosity::frand(const double a,const double b){
  return (rand()/(double)RAND_MAX) * (b-a) + a;
}

void mlawNonLocalPorosity::setSubStepping(const bool fl, const int maxNumStep){
  _withSubstepping = fl;
  _maxAttemptSubstepping = maxNumStep;
  if (fl){
    Msg::Info("with substeppint, try maximal %d steps", _maxAttemptSubstepping);
  }
};

double mlawNonLocalPorosity::deformationEnergy(const double K, const double G, const STensor3 &C) const
{
  double Jac= sqrt((C(0,0) * (C(1,1) * C(2,2) - C(1,2) * C(2,1)) -
          C(0,1) * (C(1,0) * C(2,2) - C(1,2) * C(2,0)) +
          C(0,2) * (C(1,0) * C(2,1) - C(1,1) * C(2,0))));
  double lnJ = log(Jac);
  static STensor3 logCdev;
  STensorOperation::logSTensor3(C,_orderlogexp,logCdev);
  double trace = logCdev.trace();
  logCdev(0,0)-=trace/3.;
  logCdev(1,1)-=trace/3.;
  logCdev(2,2)-=trace/3.;

  double Psy = K*0.5*lnJ*lnJ+G*0.25*dot(logCdev,logCdev);
  return Psy;
}

void mlawNonLocalPorosity::elasticPredictor(const double K, const double G, const STensor3& F1, const STensor3& Fppr,
                                  STensor3& kCor,  STensor3& kcorDev, double& ppr,
                                  STensor43 &dKCordEe, STensor3& Fepr, STensor3& Cepr, STensor3& Eepr, 
                                  STensor43 &Le, STensor63 &dLe, const bool stiff
                                  ) const{
  static STensor3 invFp0;
  STensorOperation::inverseSTensor3(Fppr,invFp0);
  STensorOperation::multSTensor3(F1,invFp0,Fepr);
  STensorOperation::multSTensor3FirstTranspose(Fepr,Fepr,Cepr);
  STensorOperation::logSTensor3(Cepr,_orderlogexp,Eepr,&Le,&dLe); 
  Eepr *= 0.5;
  
  corKirchhoffStress(K,G,kcorDev, ppr, Eepr);
  kCor = kcorDev;
  kCor(0,0) += ppr;
  kCor(1,1) += ppr;
  kCor(2,2) += ppr;
  
  if (stiff){
    for (int i=0; i< 3; i++){
      for (int j=0; j<3; j++){
        for (int k=0; k<3; k++){
          for (int l=0; l<3; l++){
             dKCordEe(i,j,k,l) = K*_I(i,j)*_I(k,l) + 2.*G*_Idev(i,j,k,l);
          }
        }
      }
    }
  }
};

void mlawNonLocalPorosity::corKirchhoffStress(const double K, const double G, STensor3 &KirDev_, double &p, const STensor3 &E_) const
{
  /* Compute Kirchoff stress for natural strain tensor */
  // Pressure part
  double trace = E_.trace();
  double val = (2.*G)*trace/3.;
  p = K*trace;
  // Deviatoric part
  KirDev_ = E_;
  KirDev_*= (2.*G);
  KirDev_(0,0)-=val;
  KirDev_(1,1)-=val;
  KirDev_(2,2)-=val;
};

void mlawNonLocalPorosity::tangentComputationLocal(STensor43& dStressDF,
                              const bool plastic, 
                              const STensor3& F,  // current F
                              const STensor3& corKir,  // cor Kir
                              const STensor3& S, //  second PK
                              const STensor3& Fepr, const STensor3& Fppr, // predictor
                              const STensor43& Lpr,
                              const STensor3& Fe, const STensor3& Fp, // corrector
                              const STensor43& L, const STensor63& dL, // corrector value
                              const STensor43& DcorKirDEepr,
                              const STensor43& dFpDEepr, 
                              const STensor43& EprToF, const STensor3& Fpinv
                              ) const{
  // P = Fe. (Kcor:Le) . invFp
  // need to compute dFedF, dKcordF, dinvFpdF, dLedF
  
  static STensor43 DcorKirDF, dFpdF;
  STensorOperation::multSTensor43(DcorKirDEepr,EprToF,DcorKirDF);
  if (plastic){
    STensorOperation::multSTensor43(dFpDEepr,EprToF,dFpdF);
  }
  else{
    STensorOperation::zero(dFpdF);
  }

  // done DcorKirDF, DFpDF

  static STensor43 DinvFpdF;
  STensorOperation::zero(DinvFpdF);
  if (plastic){
    for (int i=0; i<3; i++){
      for (int j=0; j<3; j++){
        for (int p=0; p<3; p++){
          for (int q=0; q<3; q++){
            for (int k=0; k<3; k++){
              for (int l=0; l<3; l++){
                DinvFpdF(i,j,k,l) -= Fpinv(i,p)*dFpdF(p,q,k,l)*Fpinv(q,j);
              }
            }
          }
        }
      }
    }
  }

  static STensor43 dFedF;
  STensorOperation::zero(dFedF);
  for (int i=0; i<3; i++){
    for (int j=0; j<3; j++){
      for (int k=0; k<3; k++){
        dFedF(i,j,i,k) += Fpinv(k,j);
        if (plastic){
          for (int l=0; l<3; l++){
            for (int p=0; p<3; p++){
              dFedF(i,j,k,l) += F(i,p)*DinvFpdF(p,j,k,l);
            }
          }
        }
      }
    }
  }
  static STensor63 DLDF;
  STensorOperation::zero(DLDF);
  for (int i=0; i<3; i++){
    for (int j=0; j<3; j++){
      for (int k=0; k<3; k++){
        for (int l=0; l<3; l++){
          for (int r=0; r<3; r++){
            for (int s=0; s<3; s++){
              for (int a=0; a<3; a++){
                for (int p=0; p<3; p++){
                  for (int q=0; q<3; q++){
                    DLDF(i,j,k,l,p,q) += dL(i,j,k,l,r,s)*2.*Fe(a,r)*dFedF(a,s,p,q);
                  }
                }
              }
            }
          }
        }
      }
    }
  }
  //

  static STensor43 DSDF; // S = corKir:L
  STensorOperation::zero(DSDF);
  for (int i=0; i<3; i++){
    for (int j=0; j<3; j++){
      for (int r=0; r<3; r++){
        for (int s=0; s<3; s++){
          for (int k=0; k<3; k++){
            for (int l=0; l<3; l++){
              DSDF(i,j,k,l) += DcorKirDF(r,s,k,l)*L(r,s,i,j) + corKir(r,s)*DLDF(r,s,i,j,k,l);
            }
          }
        }
      }
    }
  }


  // compute mechanical tengent
  STensorOperation::zero(dStressDF);
  for (int i=0; i<3; i++){
    for (int j=0; j<3; j++){
      for (int k=0; k<3; k++){
        for (int l=0; l<3; l++){
          for (int p=0; p<3; p++){
            for (int q=0; q<3; q++){
              dStressDF(i,j,p,q) += (dFedF(i,k,p,q)*S(k,l)*Fpinv(j,l) + Fe(i,k)*DSDF(k,l,p,q)*Fpinv(j,l)+Fe(i,k)*S(k,l)*DinvFpdF(j,l,p,q));
            }
          }
        }
      }
    }
  }                              
}

void mlawNonLocalPorosity::tangentComputation(STensor43& dStressDF, STensor3& dStressDtildefV, 
                              STensor3& DfVdF, double& dfVDNonLocalPorosity,
                              const bool plastic, 
                              const STensor3& F,  // current F
                              const STensor3& corKir,  // cor Kir
                              const STensor3& S, //  second PK
                              const STensor3& Fepr, const STensor3& Fppr, // predictor
                              const STensor43& Lpr,
                              const STensor3& Fe, const STensor3& Fp, // corrector
                              const STensor43& L, const STensor63& dL, // corrector value
                              const STensor43& DcorKirDEepr, const STensor3& DcorKirDtildefV,  // small strain tangents
                              const STensor43& dFpDEepr, const STensor3& DFpDtildefV, 
                              const STensor3& DfVDEpr, const double& DfVDtildefV,
                              const STensor43& EprToF, const STensor3& Fpinv
                              ) const{
  // P = Fe. (Kcor:Le) . invFp
  // need to compute dFedF, dKcordF, dinvFpdF, dLedF
  
  STensorOperation::multSTensor3STensor43(DfVDEpr,EprToF,DfVdF);
  dfVDNonLocalPorosity = DfVDtildefV;

  static STensor43 DcorKirDF, dFpdF;
  STensorOperation::multSTensor43(DcorKirDEepr,EprToF,DcorKirDF);
  if (plastic){
    STensorOperation::multSTensor43(dFpDEepr,EprToF,dFpdF);
  }
  else{
    STensorOperation::zero(dFpdF);
  }

  // done DcorKirDF, DcorKirDtildefV, DFpDF, DFpDtildefV

  static STensor43 DinvFpdF;
  static STensor3 DinvFpDtildefV;
  STensorOperation::zero(DinvFpdF);
  STensorOperation::zero(DinvFpDtildefV);
  if (plastic){
    for (int i=0; i<3; i++){
      for (int j=0; j<3; j++){
        for (int p=0; p<3; p++){
          for (int q=0; q<3; q++){
            DinvFpDtildefV(i,j) -= Fpinv(i,p)*DFpDtildefV(p,q)*Fpinv(q,j);
            for (int k=0; k<3; k++){
              for (int l=0; l<3; l++){
                DinvFpdF(i,j,k,l) -= Fpinv(i,p)*dFpdF(p,q,k,l)*Fpinv(q,j);
              }
            }
          }
        }
      }
    }
  }

  static STensor43 dFedF;
  static STensor3 dFeDtildefV;
  STensorOperation::zero(dFedF);
  STensorOperation::zero(dFeDtildefV);
  for (int i=0; i<3; i++){
    for (int j=0; j<3; j++){
      for (int k=0; k<3; k++){
        dFedF(i,j,i,k) += Fpinv(k,j);
        if (plastic){
          dFeDtildefV(i,j) += F(i,k)*DinvFpDtildefV(k,j);
          for (int l=0; l<3; l++){
            for (int p=0; p<3; p++){
              dFedF(i,j,k,l) += F(i,p)*DinvFpdF(p,j,k,l);
            }
          }
        }
      }
    }
  }
  static STensor63 DLDF;
  static STensor43 DLDtildefV;
  STensorOperation::zero(DLDF);
  STensorOperation::zero(DLDtildefV);
  for (int i=0; i<3; i++){
    for (int j=0; j<3; j++){
      for (int k=0; k<3; k++){
        for (int l=0; l<3; l++){
          for (int r=0; r<3; r++){
            for (int s=0; s<3; s++){
              for (int a=0; a<3; a++){
                DLDtildefV(i,j,k,l) += dL(i,j,k,l,r,s)*2.*Fe(a,r)*dFeDtildefV(a,s);
                for (int p=0; p<3; p++){
                  for (int q=0; q<3; q++){
                    DLDF(i,j,k,l,p,q) += dL(i,j,k,l,r,s)*2.*Fe(a,r)*dFedF(a,s,p,q);
                  }
                }
              }
            }
          }
        }
      }
    }
  }
  //

  static STensor43 DSDF; // S = corKir:L
  static STensor3 DSDtildefV;
  STensorOperation::zero(DSDF);
  STensorOperation::zero(DSDtildefV);
  for (int i=0; i<3; i++){
    for (int j=0; j<3; j++){
      for (int r=0; r<3; r++){
        for (int s=0; s<3; s++){
          DSDtildefV(i,j) += DcorKirDtildefV(r,s)*L(r,s,i,j) + corKir(r,s)*DLDtildefV(r,s,i,j);
          for (int k=0; k<3; k++){
            for (int l=0; l<3; l++){
              DSDF(i,j,k,l) += DcorKirDF(r,s,k,l)*L(r,s,i,j) + corKir(r,s)*DLDF(r,s,i,j,k,l);
            }
          }
        }
      }
    }
  }


  // compute mechanical tengent
  STensorOperation::zero(dStressDF);
  STensorOperation::zero(dStressDtildefV);
  for (int i=0; i<3; i++){
    for (int j=0; j<3; j++){
      for (int k=0; k<3; k++){
        for (int l=0; l<3; l++){
          dStressDtildefV(i,j) += (dFeDtildefV(i,k)*S(k,l)*Fpinv(j,l) + Fe(i,k)*DSDtildefV(k,l)*Fpinv(j,l)+Fe(i,k)*S(k,l)*DinvFpDtildefV(j,l));
          for (int p=0; p<3; p++){
            for (int q=0; q<3; q++){
              dStressDF(i,j,p,q) += (dFedF(i,k,p,q)*S(k,l)*Fpinv(j,l) + Fe(i,k)*DSDF(k,l,p,q)*Fpinv(j,l)+Fe(i,k)*S(k,l)*DinvFpdF(j,l,p,q));
            }
          }
        }
      }
    }
  }
};