/****************************************************************************
*
*                                 M U E S L I   v 1.8
*
*
*     Copyright 2020 IMDEA Materials Institute, Getafe, Madrid, Spain
*     Contact: muesli.materials@imdea.org
*     Author: Ignacio Romero (ignacio.romero@imdea.org)
*
*     This file is part of MUESLI.
*
*     MUESLI is free software: you can redistribute it and/or modify
*     it under the terms of the GNU General Public License as published by
*     the Free Software Foundation, either version 3 of the License, or
*     (at your option) any later version.
*
*     MUESLI is distributed in the hope that it will be useful,
*     but WITHOUT ANY WARRANTY; without even the implied warranty of
*     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
*     GNU General Public License for more details.
*
*     You should have received a copy of the GNU General Public License
*     along with MUESLI.  If not, see <http://www.gnu.org/licenses/>.
*
****************************************************************************/


#include "johnsoncook.h"
#include "../Fcoupled/thermojc.h"
#include <string.h>
#include <cmath>

//Tolerances can be decreased in case convergence is affected
#define J2TOL1     1e-10
#define J2TOL2     1e-10
#define GAMMAITER1 10
#define GAMMAITER2 100
#define SQ23       0.816496580927726

using namespace std;
using namespace muesli;


johnsonCookMaterial::johnsonCookMaterial(const std::string& name,
                                         const materialProperties& cl)
:
finiteStrainMaterial(name, cl),
theThermoJCMaterial(0),
E(0.0), nu(0.0), lambda(0.0), mu(0.0), bulk(0.0),
cp(0.0), cs(0.0), _A(0.0), _B(0.0), _C(0.0), _M(0.0), _N(0.0),
_edot0(1.0), _curT(1.0), _modelRefTemp(1.0), _meltT(1.0), _maxDamage(0.99)
{
    muesli::assignValue(cl, "young",       E);
    muesli::assignValue(cl, "poisson",     nu);
    muesli::assignValue(cl, "lambda",      lambda);
    muesli::assignValue(cl, "mu",          mu);
    muesli::assignValue(cl, "a",           _A);
    muesli::assignValue(cl, "b",           _B);
    muesli::assignValue(cl, "c",           _C);
    muesli::assignValue(cl, "m",           _M);
    muesli::assignValue(cl, "n",           _N);
    muesli::assignValue(cl, "edot0",       _edot0);
    muesli::assignValue(cl, "temp",        _curT);
    muesli::assignValue(cl, "modelreftemp",_modelRefTemp);
    muesli::assignValue(cl, "melttemp",    _meltT);
    muesli::assignValue(cl, "maxDamage",   _maxDamage);
    if (cl.find("damagejc") != cl.end())
    {
        damageModelActivated = true;
    }
    else if (cl.find("damagejccustom") != cl.end())
    {
        damageModelActivated = true;
    }
    if ( cl.find("deletion") != cl.end() ) deletion = true;

    // E and nu have priority. If they are defined, define lambda and mu
    if (E*E > 0.0)
    {
        lambda = E*nu/(1.0-2.0*nu)/(1.0+nu);
        mu     = E/2.0/(1.0+nu);
    }
    else
    {
        nu = 0.5 * lambda / (lambda+mu);
        E  = mu*2.0*(1.0+nu);
    }
    
    // We set all the constants, so that later on all of them can be recovered fast
    bulk = lambda + 2.0/3.0 * mu;

    double rho = density();
    if (rho > 0.0)
    {
        cp = sqrt((lambda+2.0*mu)/rho);
        cs = sqrt(mu/rho);
    }
    
    theThermoJCMaterial = new thermoJCMaterial(name,cl);
}




// Alternative material creation method, for MP individual tests
johnsonCookMaterial::johnsonCookMaterial(const std::string& name, const double xE, const double xnu,
                                         const double rho, const double x_A, const double x_B,
                                         const double x_C, const double x_M, const double x_N,
                                         const double x_edot0, const double x_curT, const double x_modelReTtemp,
                                         const double x_meltT, const bool x_damageModelActivated,
                                         const double x_D1, const double x_D2, const double x_D3,
                                         const double x_D4, const double x_D5)
:
finiteStrainMaterial(name),
theThermoJCMaterial(0),
E(xE), nu(xnu), lambda(0.0), mu(0.0), bulk(0.0),
cp(0.0), cs(0.0), _A(x_A), _B(x_B), _C(x_C), _M(x_M), _N(x_N),
_edot0(x_edot0), _curT(x_curT), _modelRefTemp(x_modelReTtemp), _meltT(x_curT), _maxDamage(0.99)
{
    setDensity(rho);
    setReferenceTemperature(x_curT);

    // E and nu have priority. If they are defined, define lambda and mu
    if (E*E > 0.0)
    {
        lambda = E*nu/(1.0-2.0*nu)/(1.0+nu);
        mu     = E/2.0/(1.0+nu);
    }
    else
    {
        nu = 0.5 * lambda / (lambda+mu);
        E  = mu*2.0*(1.0+nu);
    }
    
    // We set all the constants, so that later on all of them can be recovered fast
    bulk = lambda + 2.0/3.0 * mu;
    
    if (rho > 0.0)
    {
        cp = sqrt((lambda+2.0*mu)/rho);
        cs = sqrt(mu/rho);
    }
    
    damageModelActivated = x_damageModelActivated;
    deletion = false;
    
    
    theThermoJCMaterial = new thermoJCMaterial(name, xE, xnu, rho, x_A,
                                               x_B, x_C, x_M, x_N, x_edot0,
                                               x_curT, x_modelReTtemp, x_meltT,
                                               x_damageModelActivated,
                                               x_D1, x_D2, x_D3, x_D4, x_D5);
}




johnsonCookMaterial::johnsonCookMaterial(const std::string& name, const double xE, const double xnu,
                                         const double rho, const double x_A, const double x_B,
                                         const double x_C, const double x_M, const double x_N,
                                         const double x_edot0, const double x_curT, const double x_modelReTtemp,
                                         const double x_meltT, const bool x_damageModelActivated,
                                         const double x_D1, const double x_D2, const double x_D3,
                                         const double x_D4, const double x_D5, const double x_D6)
:
finiteStrainMaterial(name),
theThermoJCMaterial(0),
E(xE), nu(xnu), lambda(0.0), mu(0.0), bulk(0.0),
cp(0.0), cs(0.0), _A(x_A), _B(x_B), _C(x_C), _M(x_M), _N(x_N),
_edot0(x_edot0), _curT(x_curT),_modelRefTemp(x_modelReTtemp), _meltT(x_curT), _maxDamage(0.99)
{
    setDensity(rho);
    setReferenceTemperature(x_curT);

    // E and nu have priority. If they are defined, define lambda and mu
    if (E*E > 0.0)
    {
        lambda = E*nu/(1.0-2.0*nu)/(1.0+nu);
        mu     = E/2.0/(1.0+nu);
    }
    else
    {
        nu = 0.5 * lambda / (lambda+mu);
        E  = mu*2.0*(1.0+nu);
    }
    
    // We set all the constants, so that later on all of them can be recovered fast
    bulk = lambda + 2.0/3.0 * mu;
    
    if (rho > 0.0)
    {
        cp = sqrt((lambda+2.0*mu)/rho);
        cs = sqrt(mu/rho);
    }
    
    damageModelActivated = x_damageModelActivated;
    deletion = false;
    
    
    theThermoJCMaterial = new thermoJCMaterial(name, xE, xnu, rho, x_A,
                                               x_B, x_C, x_M, x_N, x_edot0,
                                               x_curT, x_modelReTtemp, x_meltT,
                                               x_damageModelActivated,
                                               x_D1, x_D2, x_D3, x_D4, x_D5, x_D6);
}




double johnsonCookMaterial::characteristicStiffness() const
{
    return E;
}




bool johnsonCookMaterial::check() const
{
    bool ok = theThermoJCMaterial->check();
    
    return ok;
}




muesli::finiteStrainMP* johnsonCookMaterial::createMaterialPoint() const
{
    muesli::finiteStrainMP *mp = new johnsonCookMP(*this);
    
    return mp;
}




// This function is much faster than the one with string property names, since it
// avoids string comparisons. It should be used.
double johnsonCookMaterial::getProperty(const propertyName p) const
{
    double ret=0.0;
    
    switch (p)
    {
        case PR_LAMBDA:     ret = lambda;   break;
        case PR_MU:         ret = mu;       break;
        case PR_YOUNG:      ret = E;        break;
        case PR_POISSON:    ret = nu;       break;
        case PR_BULK:       ret = bulk;     break;
        case PR_CP:         ret = cp;       break;
        case PR_CS:         ret = cs;       break;
            
        default:
            std::cout << "\n Error in elasticIsotropicMaterial. Property not defined";
    }
    return ret;
}




void johnsonCookMaterial::print(std::ostream &of) const
{
    if (damageModelActivated)
    {
        of  << "\n Johnson - Cook rate- and temperature (imposed)-dependent plasticity with Damage option ";
    }
    else
    {
        of  << "\n Johnson - Cook rate- and temperature (imposed)-dependent  plasticity ";
    }
    of << "\n   Young modulus:  E      : " << E
    << "\n   Poisson ratio:  nu     : " << nu
    << "\n   Lame constants: Lambda : " << lambda
    << "\n                   Mu     : " << mu
    << "\n   Bulk modulus:   K      : " << bulk
    << "\n   Density                : " << density()
    << "\n   Wave velocities C_p    : " << cp
    << "\n                   C_s    : " << cs;

    if (damageModelActivated)
    {
        of  << "\n   The yield Kirchhoff stress is of the form:"
        << "\n    |tau| - (1-D) sqrt(2/3) (A+B e^N) (1+C ln (edot/edot0))(1-theta^M)"
        << "\n    theta = (T - T0)/(Tm - T0)";
    }
    else
    {
        of  << "\n   The yield Kirchhoff stress is of the form:"
        << "\n    |tau| - sqrt(2/3) (A+B e^N) (1+C ln (edot/edot0))(1-theta^M)"
        << "\n    theta = (T - T0)/(Tm - T0)";
    }
    
    of  << "\n   A                      : " << _A
    << "\n   B                      : " << _B
    << "\n   C                      : " << _C
    << "\n   N                      : " << _N
    << "\n   dot{eps}_0             : " << _edot0
    << "\n   T                      : " << _curT
    << "\n   T0                     : " << _modelRefTemp
    << "\n   Tm                     : " << _meltT
    << "\n   M                      : " << _M;
    if (damageModelActivated)
    {
        of  <<"\n  maxDamage                  : " << _maxDamage;
    }
}




void johnsonCookMaterial::setRandom()
{
    theThermoJCMaterial->setRandom();
    
    this->E=theThermoJCMaterial->E;
    this->nu=theThermoJCMaterial->nu;
    this->_A=theThermoJCMaterial->_A;
    this->_B=theThermoJCMaterial->_B;
    this->_C=theThermoJCMaterial->_C;
    this->_M=theThermoJCMaterial->_M;
    this->_N=theThermoJCMaterial->_N;
    this->_edot0=theThermoJCMaterial->_edot0;
    this->_meltT=theThermoJCMaterial->_meltT;
    this->_modelRefTemp=theThermoJCMaterial->_modelRefTemp;
    this->_curT=theThermoJCMaterial->_curT;

    setDensity(theThermoJCMaterial->material::density());
    setReferenceTemperature(theThermoJCMaterial->_curT);
    
    if (damageModelActivated)
    {
        this->_maxDamage = theThermoJCMaterial->_maxDamage;
    }
    
    this->lambda=theThermoJCMaterial->lambda;
    this->mu=theThermoJCMaterial->mu;
    this->cp=theThermoJCMaterial->cp;
    this->cs=theThermoJCMaterial->cs;
    this->bulk=theThermoJCMaterial->bulk;
}




bool johnsonCookMaterial::test(std::ostream &of)
{
    bool isok=true;
    setRandom();
    
    muesli::finiteStrainMP* p = this->createMaterialPoint();
    p->setRandom();
    
    isok = p->testImplementation(of, false);
    delete p;
    return isok;
}




double johnsonCookMaterial::waveVelocity() const
{
    return cp;
}




johnsonCookMP::johnsonCookMP(const johnsonCookMaterial &m)
:
finiteStrainMP(m),
theElastoplasticMaterial(m),
theThermoJCMP(0)
{
    theThermoJCMP = new thermoJCMP(*m.theThermoJCMaterial);
}




// Brent method function, as backup for Newton-Raphson, in case it fails
double johnsonCookMP::brentroot(double a, double b, double Ga,
                                double Gb, double eqpn, double ntbar,
                                double mu, double A,double B,
                                double C, double N,double edot0,
                                double dt, double tempterm)
{
    double root = theThermoJCMP->brentroot(a, b, Ga, Gb, eqpn, ntbar, mu,
                                           A, B, C, N, edot0, dt, tempterm);
    
    return root;
}




void johnsonCookMP::CauchyStress(istensor& sigma) const
{
    theThermoJCMP->mechCauchyStress(sigma);
}




void johnsonCookMP::commitCurrentState()
{
    finiteStrainMP::commitCurrentState();
    
    theThermoJCMP->commitCurrentState();
}




void johnsonCookMP::contractWithConvectedTangent(const ivector& V1, const ivector& V2, itensor& T) const
{
    itensor4 c;
    theThermoJCMP->mechConvectedTangent(c);
    
    T.setZero();
    for (unsigned i=0; i<3; i++)
        for (unsigned j=0; j<3; j++)
            for (unsigned k=0; k<3; k++)
                for (unsigned l=0; l<3; l++)
                {
                    T(i,k) += c(i,j,k,l)*V1(j)*V2(l);
                }
}




// Calls generic FiniteStrain numerical convected tangent
void johnsonCookMP::convectedTangent(itensor4& C) const
{
    theThermoJCMP->mechConvectedTangent(C);
}




void johnsonCookMP::convectedTangentTimesSymmetricTensor(const istensor& M,istensor& CM) const
{
    itensor4 C;
    theThermoJCMP->mechConvectedTangent(C);
    
    CM.setZero();
    for (unsigned i=0; i<3; i++)
        for (unsigned j=0; j<3; j++)
            for (unsigned k=0; k<3; k++)
                for (unsigned l=0; l<3; l++)
                {
                    CM(i,j) += C(i,j,k,l)*M(k,l);
                }
}




double johnsonCookMP::deviatoricEnergy() const
{
    return 0.0;
}




double johnsonCookMP::dissipatedEnergy() const
{
    double dissEnergy = theThermoJCMP->dissipatedEnergy();
    
    return dissEnergy;
}




double johnsonCookMP::dissipatedEnergyDTheta() const
{
    double dissEnergyDTheta = theThermoJCMP->dissipatedEnergyDtheta();
    
    return dissEnergyDTheta;
}




double johnsonCookMP::effectiveStoredEnergy() const
{
    return theThermoJCMP->storedEnergy() + theThermoJCMP->dissipatedEnergy();
}




void johnsonCookMP::explicitRadialReturn(const ivector &taudev, double ep, double epdot)
{
    theThermoJCMP->explicitRadialReturn(taudev, ep, epdot);
}




double johnsonCookMP::kineticPotential() const
{
    double kinPot = theThermoJCMP->kineticPotential();
    
    return kinPot;
}




double johnsonCookMP::plasticSlip() const
{
    double iso_c = theThermoJCMP->plasticSlip();
    
    return iso_c;
}




void johnsonCookMP::plasticReturn(const ivector& taubarTR)
{
    theThermoJCMP->plasticReturn(taubarTR);
}




// Calculation of G function value given each iteration value of dgamma
void johnsonCookMP::plasticReturnResidual(double mu, double A, double B, double C,
                                          double N, double edot0, double eqpn, double tempterm,
                                          double tau, double dt, double dgamma, double& G)
{
    theThermoJCMP->plasticReturnResidual(mu, A, B, C, N, edot0, eqpn, tempterm,
                                         tau, dt, dgamma, G);
}




// Calculation of G function derivative value given each iteration value of dgamma
void  johnsonCookMP::plasticReturnTangent(double mu, double A, double B, double C,
                                          double N, double edot0, double eqpn, double tempterm,
                                          double tau, double dt, double dgamma, double& DG)
{
    theThermoJCMP->plasticReturnTangent(mu, A, B, C, N, edot0, eqpn, tempterm,
                                        tau, dt, dgamma, DG);
}




void johnsonCookMP::resetCurrentState()
{
    finiteStrainMP::resetCurrentState();
    
    theThermoJCMP->resetCurrentState();
}




void johnsonCookMP::setConvergedState(const double theTime, const itensor& F,
                                      const double iso, const ivector& kine,
                                      const istensor& be)
{
    tn     = theTime;
    Fn     = F;
    Jn     = Fn.determinant();
    
    theThermoJCMP->setConvergedState(theTime, F, iso, kine, be);
}




void johnsonCookMP::setRandom()
{
    theThermoJCMP->setRandom();
    
    Fc = Fn = theThermoJCMP->getDefGrad();
    Jc = Fc.determinant();
    if (theElastoplasticMaterial.damageModelActivated)
    {
        this->Dc=theThermoJCMP->getDamage();
    }
}




void johnsonCookMP::spatialTangent(itensor4& Cs) const
{
    theThermoJCMP->mechSpatialTangent(Cs);
}




double johnsonCookMP::storedEnergy() const
{
    double WpWe = theThermoJCMP->storedEnergy();

    return WpWe;
}




void johnsonCookMP::updateCurrentState(const double theTime, const istensor& C)
{
    ivector zero (0.0, 0.0, 0.0);
    theThermoJCMP->updateCurrentState(theTime, C, zero, theElastoplasticMaterial._curT);
    
    if (theElastoplasticMaterial.damageModelActivated)
    {
        Dc = theThermoJCMP->Dc;
        fullyDamaged = theThermoJCMP->fullyDamaged;
    }
}




void johnsonCookMP::updateCurrentState(const double theTime, const itensor& F)
{
    tc = theTime;
    Fc = F;
    Jc = F.determinant();
    
    ivector zero (0.0, 0.0, 0.0);
    theThermoJCMP->updateCurrentState(theTime, F, zero, theElastoplasticMaterial._curT);
    
    if (theElastoplasticMaterial.damageModelActivated)
    {
        Dc = theThermoJCMP->Dc;
        fullyDamaged = theThermoJCMP->fullyDamaged;
    }
}




double johnsonCookMP::volumetricEnergy() const
{
    return 0.0;
}




// Yield function in principal Kirchhoff space
double johnsonCookMP::yieldFunction(const ivector& tau,
                                    const double&  eps,
                                    const double&  epsdot) const
{
    double yield = theThermoJCMP->yieldFunction(tau, eps, epsdot);
    
    return yield;
}




// Yield function in principal Kirchhoff space, when Damage model is in use
double johnsonCookMP::yieldFunctionDMG(const ivector& tau,
                                       const double&  eps,
                                       const double&  epsdot,
                                       const double& D) const
{
    double yield = theThermoJCMP->yieldFunctionDMG(tau, eps, epsdot, D);
    
    return yield;
}
