/****************************************************************************
*
*                                 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 "zerilliarmstrong.h"
#include "../Fcoupled/thermoza.h"
#include <string.h>
#include <cmath>
#include <stdio.h>
#include <iostream>

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




zerilliArmstrongMaterial :: zerilliArmstrongMaterial(const std::string& name,
                                                     const materialProperties& cl)
:
finiteStrainMaterial(name, cl),
theThermoZAMaterial(0),
E(0.0), nu(0.0), lambda(0.0), mu(0.0), bulk(0.0),
cp(0.0), cs(0.0), _C0 (0.0), _C1(0.0), _C2(0.0), _C3(0.0),
_C4(0.0), _C5(0.0), _n(0.0), _curT(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, "c0",           _C0);
    muesli::assignValue(cl, "c1",           _C1);
    muesli::assignValue(cl, "c2",           _C2);
    muesli::assignValue(cl, "c3",           _C3);
    muesli::assignValue(cl, "c4",           _C4);
    muesli::assignValue(cl, "c5",           _C5);
    muesli::assignValue(cl, "n",             _n);
    muesli::assignValue(cl, "temp",       _curT);
    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;
    
    if (rho > 0.0)
    {
        cp = sqrt((lambda+2.0*mu)/rho);
        cs = sqrt(mu/rho);
    }
    
    theThermoZAMaterial = new thermoZAMaterial(name,cl);
}




zerilliArmstrongMaterial :: zerilliArmstrongMaterial(const std::string& name, const double xE, const double xnu,
                                                     const double rho, const double x_C0, const double x_C1,
                                                     const double x_C2, const double x_C3, const double x_C4,
                                                     const double x_C5, const double x_n, const double x_curT,
                                                     const double x_edot0, 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),
theThermoZAMaterial(0),
E(xE), nu(xnu), lambda(0.0), mu(0.0), bulk(0.0),
cp(0.0), cs(0.0), _C0 (x_C0), _C1(x_C1), _C2(x_C2), _C3(x_C3),
_C4(x_C4), _C5(x_C5), _n(x_n), _curT(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;
        
        theThermoZAMaterial = new thermoZAMaterial(name, xE, xnu, rho,
                                                   x_C0, x_C1, x_C2, x_C3, x_C4, x_C5, x_n,
                                                   x_curT, x_edot0, x_modelReTtemp, x_meltT,
                                                   x_damageModelActivated,
                                                   x_D1, x_D2, x_D3, x_D4, x_D5);
}




zerilliArmstrongMaterial :: zerilliArmstrongMaterial(const std::string& name, const double xE, const double xnu,
                                                     const double rho, const double x_C0, const double x_C1,
                                                     const double x_C2, const double x_C3, const double x_C4,
                                                     const double x_C5, const double x_n, const double x_curT,
                                                     const double x_edot0, 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),
theThermoZAMaterial(0),
E(xE), nu(xnu), lambda(0.0), mu(0.0), bulk(0.0),
cp(0.0), cs(0.0), _C0 (x_C0), _C1(x_C1), _C2(x_C2), _C3(x_C3),
_C4(x_C4), _C5(x_C5), _n(x_n), _curT(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;

    
    theThermoZAMaterial = new thermoZAMaterial(name, xE, xnu, rho,
                                               x_C0, x_C1, x_C2, x_C3, x_C4, x_C5, x_n,
                                               x_curT, x_edot0, x_modelReTtemp, x_meltT,
                                               x_damageModelActivated,
                                               x_D1,x_D2, x_D3, x_D4, x_D5, x_D6);
}




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




bool zerilliArmstrongMaterial::check() const
{
    bool ok = theThermoZAMaterial->check();
    
    return ok;
}




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




double zerilliArmstrongMaterial::density() const
{
    return rho;
}




// This function is much faster than the one with string property names, because
// avoids string comparisons. It should be used.
double zerilliArmstrongMaterial::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 zerilliArmstrongMaterial :: print(std::ostream &of) const
{
    if (damageModelActivated)
    {
        of  << "\n Zerilli - Armstrong rate- and temperature (imposed)-dependent plasticity with Damage option ";
    }
    else
    {
        of  << "\n Zerilli - Armstrong 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                : " << rho;
    if (rho > 0.0)
    {
        of  << "\n   Wave velocities C_p    : " << cp;
        of  << "\n                   C_s    : " << cs;
    }
    if (damageModelActivated)
    {
    of  << "\n   The yield Kirchhoff stress is of the form:"
    << "\n    |tau| - (1-D) sqrt(2/3) (C1 + C2 eps^(1/2)) exp(-C3T + C4T ln(epsdot))+ C5eps^n+kl^(-1/2)+sigmaG";
    }
    else
    {
    of  << "\n   The yield Kirchhoff stress is of the form:"
    << "\n    |tau| - sqrt(2/3) (C1 + C2 eps^(1/2)) exp(-C3T + C4T ln(epsdot))+ C5eps^n+kl^(-1/2)+sigmaG";
    }
    
    of  << "\n   C0                      : " << _C0
    << "\n   C1                      : " << _C1
    << "\n   C2                      : " << _C2
    << "\n   C3                      : " << _C3
    << "\n   C4                      : " << _C4
    << "\n   C5                      : " << _C5
    << "\n   temp                    : " << _curT
    << "\n   n                       : " << _n;
    if (damageModelActivated)
    {
        of  <<"\n  maxDamage                  : " << _maxDamage;
    }
}




void zerilliArmstrongMaterial :: setRandom()
{
    theThermoZAMaterial->setRandom();
    
    this->E=theThermoZAMaterial->E;
    this->nu=theThermoZAMaterial->nu;
    this->_C0=theThermoZAMaterial->_C0;
    this->_C1=theThermoZAMaterial->_C1;
    this->_C2=theThermoZAMaterial->_C2;
    this->_C3=theThermoZAMaterial->_C3;
    this->_C4=theThermoZAMaterial->_C4;
    this->_C5=theThermoZAMaterial->_C5;
    this->_n=theThermoZAMaterial->_n;
    this->_curT=theThermoZAMaterial->_curT;
    
    setDensity(theThermoZAMaterial->material::density());
    setReferenceTemperature(theThermoZAMaterial->_curT);
    
    if (damageModelActivated)
    {
        this->_maxDamage = theThermoZAMaterial->_maxDamage;
    }

    this->lambda=theThermoZAMaterial->lambda;
    this->mu=theThermoZAMaterial->mu;
    this->cp=theThermoZAMaterial->cp;
    this->cs=theThermoZAMaterial->cs;
    this->bulk=theThermoZAMaterial->bulk;
}




bool zerilliArmstrongMaterial::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 zerilliArmstrongMaterial::waveVelocity() const
{
    return cp;
}




zerilliArmstrongMP::zerilliArmstrongMP(const zerilliArmstrongMaterial &m)
:
finiteStrainMP(m),
theElastoplasticMaterial(m),
theThermoZAMP(0)
{
    theThermoZAMP = new thermoZAMP(*m.theThermoZAMaterial);
}




// Brent method function
double zerilliArmstrongMP::brentroot(double a, double b, double Ga,
                                     double Gb, double eqpn, double ntbar,
                                     double mu, double C0,double C1,
                                     double C2, double C3,double C4,
                                     double C5, double Tc, double n, double dt)
{
    double root = theThermoZAMP->brentroot(a, b, Ga, Gb, eqpn, ntbar, mu,
                                           C0, C1, C2, C3, C4, C5, Tc, n, dt);
    
    return root;
}




void zerilliArmstrongMP::CauchyStress(istensor& sigma) const
{
    theThermoZAMP->mechCauchyStress(sigma);
}




void zerilliArmstrongMP::commitCurrentState()
{

    finiteStrainMP::commitCurrentState();
    
    theThermoZAMP->commitCurrentState();
}


void zerilliArmstrongMP::contractWithConvectedTangent(const ivector& V1, const ivector& V2, itensor& T) const
{
    itensor4 c;
    theThermoZAMP->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);
                }
}




void zerilliArmstrongMP :: convectedTangent(itensor4& C) const
{
    theThermoZAMP->mechConvectedTangent(C);
}




void zerilliArmstrongMP::convectedTangentTimesSymmetricTensor(const istensor& M,istensor& CM) const
{
    itensor4 C;
    theThermoZAMP->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 zerilliArmstrongMP::deviatoricEnergy() const
{
    return 0.0;
}




double zerilliArmstrongMP::dissipatedEnergy() const
{
    double dissEnergy = theThermoZAMP->dissipatedEnergy();
    
    return dissEnergy;
}




double zerilliArmstrongMP::effectiveStoredEnergy() const
{
    return theThermoZAMP->storedEnergy() + theThermoZAMP->dissipatedEnergy();
}




// Explicit radial return algorithm not tested within this code
// and not used in the "update state" function
void zerilliArmstrongMP::explicitRadialReturn(const ivector &taudev, double ep, double epdot)
{
    theThermoZAMP->explicitRadialReturn(taudev, ep, epdot);
}




double zerilliArmstrongMP::kineticPotential() const
{
    double kinPot = theThermoZAMP->kineticPotential();
    
    return kinPot;
}




double zerilliArmstrongMP::plasticSlip() const
{
    double iso_c = theThermoZAMP->plasticSlip();
    
    return iso_c;
}




// TaubarTR: trial deviatoric tau.
// This function finds dgamma root and equivalent plastic strain eqp in the current time step.
// Two-fold approach convergence algorithm, using in first place a NR scheme to find the root (dgamma)
// of G equation. In case this algorithm is not able (due to high derivative problems near the root
// associated to ZA equation), the function jumps to a Brent scheme, whose convergence is assured.
void zerilliArmstrongMP::plasticReturn(const ivector& taubarTR)
{
    theThermoZAMP->plasticReturn(taubarTR);
}




// Calculation of G function value given each iteration value of dgamma
void zerilliArmstrongMP::plasticReturnResidual(double mu, double C0, double C1, double C2, double C3,
                                               double C4, double C5, double Tc, double n, double eqpn,
                                               double tau, double dt, double dgamma, double& G)
{
    theThermoZAMP->plasticReturnResidual(mu, C0, C1, C2, C3, C4, C5, Tc, n, eqpn,
                                         tau, dt, dgamma, G);
}



// Calculation of G function derivative value given each iteration value of dgamma
void  zerilliArmstrongMP::plasticReturnTangent(double mu, double C0, double C1,double C2, double C3,
                                               double C4, double C5, double Tc, double n, double eqpn,
                                               double tau, double dt, double dgamma, double& DG)
{
    theThermoZAMP->plasticReturnTangent(mu, C0, C1, C2, C3, C4, C5, Tc, n, eqpn,
                                        tau, dt, dgamma, DG);
}




void zerilliArmstrongMP::resetCurrentState()
{
    finiteStrainMP::resetCurrentState();
    
    theThermoZAMP->resetCurrentState();
}




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




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




void zerilliArmstrongMP::spatialTangent(itensor4& Cs) const
{
    theThermoZAMP->mechSpatialTangent(Cs);
}




double zerilliArmstrongMP::storedEnergy() const
{
    double WpWe = theThermoZAMP->storedEnergy();

    return WpWe;
}




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




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




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




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




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