/****************************************************************************
*
*                                 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 "thermoza.h"
#include "muesli/Damage/damagemodellib.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;


thermoZAMaterial::thermoZAMaterial(const std::string& name,
                                   const materialProperties& cl)
:
thermofiniteStrainMaterial(name, cl, "thermoZA"),
theDamageModel(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), taylorQuinney(0.9), _maxDamage(0.99),
damageModelActivated(false), deletion(false)
{
    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, "taylor_quinney", taylorQuinney);
    muesli::assignValue(cl, "maxDamage", _maxDamage);
    if (cl.find("damagejc") != cl.end())
    {
        damageModelActivated = true;
        theDamageModel = new JCDamageModel(name, cl);
    }
    else if (cl.find("damagejccustom") != cl.end())
    {
        damageModelActivated = true;
        theDamageModel = new JCCustomDamageModel(name, cl);
    }
    else if (cl.find("damagejccustom2") != cl.end())
    {
        damageModelActivated = true;
        theDamageModel = new JCCustom2DamageModel(name, cl);
    }
    
    if ( cl.find("deletion") != cl.end() ) deletion = true;
    
    muesli::assignValue(cl, "heat_supply", _heatSupply);
    muesli::assignValue(cl, "thermal_expansion", _thermalExpansion);
    muesli::assignValue(cl, "heat_capacity", _thermalCapacity);
    muesli::assignValue(cl, "conductivity", _thermalConductivity);
    
    
    // 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 = material::density();
    if (rho > 0.0)
    {
        cp = sqrt((lambda+2.0*mu)/rho);
        cs = sqrt(mu/rho);
    }
}




// Alternative material creation method, for MP individual tests
thermoZAMaterial::thermoZAMaterial(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_modelRefTemp, 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)
:
thermofiniteStrainMaterial(name),
theDamageModel(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), _maxDamage(0.99),
damageModelActivated(x_damageModelActivated), deletion(false)
{
    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);
    }
    
    // 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);
    }
    
    if (damageModelActivated) theDamageModel = new JCDamageModel(name, x_edot0,
                                                                 x_modelRefTemp, x_meltT,
                                                                 x_D1, x_D2, x_D3,
                                                                 x_D4, x_D5);
}




thermoZAMaterial::thermoZAMaterial(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_modelRefTemp, 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)
:
thermofiniteStrainMaterial(name),
theDamageModel(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), _maxDamage(0.99),
damageModelActivated(x_damageModelActivated), deletion(false)
{
    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);
    }
    
    // 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);
    }
    
    if (damageModelActivated) theDamageModel = new JCCustomDamageModel(name, x_edot0,
                                                                       x_modelRefTemp, x_meltT,
                                                                       x_D1, x_D2, x_D3,
                                                                       x_D4, x_D5, x_D6);
}




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




bool thermoZAMaterial::check() const
{
    if (_thermalConductivity <= 0.0 && mu > 0.0 && lambda+2.0*mu > 0.0) return true;
    else return false;
}




muesli::thermofiniteStrainMP* thermoZAMaterial::createMaterialPoint() const
{
    thermoZAMP *mp = new thermoZAMP(*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 thermoZAMaterial::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;
        case PR_CONDUCTIVITY:   ret = _thermalConductivity; break;
        case PR_THERMAL_EXP:    ret = _thermalExpansion; break;
        case PR_THERMAL_CAP:    ret = _thermalCapacity; break;
        case PR_HEAT_SUPPLY:    ret = _heatSupply; break;
            
        default:
            std::cout << "\n Error in Thermo-ElasticPlasticMaterial. Property not defined";
    }
    return ret;
}




void thermoZAMaterial::print(std::ostream &of) const
{
    if (damageModelActivated)
    {
        of  << "\n Zerilli - Armstrong rate- and temperature (variable)-dependent plasticity with Damage option ";
    }
    else
    {
        of  << "\n Zerilli - Armstrong rate- and temperature (variable)-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) (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;
    }
    
    of  << "\n   Conductivity : " << _thermalConductivity;
    of  << "\n   Heat supply : " << _heatSupply;
    of  << "\n   Thermal expansion coef. : " << _thermalExpansion;
    of  << "\n   Thermal capacity : " << _thermalCapacity;
    of  << "\n";
}




void thermoZAMaterial::setRandom()
{
    material::setRandom();

    E      = muesli::randomUniform(1.0, 10.0);
    nu     = muesli::randomUniform(0.05, 0.45);
    _C0    = muesli::randomUniform(1.0, 10.0);
    _C1    = muesli::randomUniform(1.0, 10.0);
    if(_C1>5.5)
    {
        _C2 = 0.0;
        _C5 = muesli::randomUniform(1.0, 10.0);
    }
    else
    {
        _C2 = muesli::randomUniform(1.0, 10.0);
        _C1 = 0.0;
        _C5 = 0.0;
    }
    _C3    = muesli::randomUniform(1.0, 10.0);
    _C4    = muesli::randomUniform(1.0, 10.0);
    _n     = muesli::randomUniform(1.0, 10.0);
    _curT  = muesli::randomUniform(400.0, 500.0);

    if (damageModelActivated)
    {
        _maxDamage = muesli::randomUniform(0.0, 1.0);
        theDamageModel->setRandom();
    }

    double rho = density();
    lambda = E*nu/(1.0-2.0*nu)/(1.0+nu);
    mu     = E/2.0/(1.0+nu);
    cp     = sqrt((lambda+2.0*mu)/rho);
    cs     = sqrt(2.0*mu/rho);
    bulk   = lambda + 2.0/3.0 * mu;
    _thermalConductivity = muesli::randomUniform(1.0, 10.0);
    _thermalExpansion    = 1e-6 * muesli::randomUniform(1.0, 10.0);
    _thermalCapacity     = muesli::randomUniform(10.0, 100.0);
}




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




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




thermoZAMP::thermoZAMP(const thermoZAMaterial &m)
:
thermofiniteStrainMP(m),
thethermoFiniteStrainMaterial(m),
tn(0.0), tc(0.0),
Jn(1.0), Jc(1.0),
Dn(0.0), Dc(0.0),
fullyDamaged(false),

iso_n(0.0), iso_c(0.0),
epdot_n(0.0), epdot_c(0.0),
dgamma(0.0)
{
    Fn = Fc = itensor::identity();
    be_n = be_c = istensor::identity();
    tau.setZero();
    lambda2 = lambda2TR = ivector(1.0, 1.0, 1.0);
    tau.setZero();
    nubarTR.setZero();
    nn[0] = ivector(1.0, 0.0, 0.0);
    nn[1] = ivector(0.0, 1.0, 0.0);
    nn[2] = ivector(0.0, 0.0, 1.0);
    
    currT = m._curT;
    temp_c = currT;
    temp_n = temp_c;
    GradT_n.setZero();
    GradT_c.setZero();
}




// Brent method function, as backup for Newton-Raphson, in case it fails
double thermoZAMP::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 s(0.0);
    double d(0.0);
    size_t count=1;
    
    if (fabs(Ga) < fabs(Gb))
    {
        double ch=a;
        a=b; b=ch;
        double Gch=Ga;
        Ga=Gb; Gb=Gch;
    }
    
    double c = a;
    double Gs = Gb;
    double Gc = Ga;
    int flag=1;
    
    
    while( !(Gs == 0.0 ) && !((fabs(Ga) < J2TOL2) && (fabs(Gb) < J2TOL2))  && count < GAMMAITER2)
    {
        if((Ga != Gc) && (Gb != Gc))
        {
            s=((a*Gb*Gc/((Ga-Gb)*(Ga-Gc)))+(b*Ga*Gc/((Gb-Ga)*(Gb-Gc)))+(c*Ga*Gb/((Gc-Ga)*(Gc-Gb))));
        }
        else
        {
            s=b-(Gb*((b-a)/(Gb-Ga)));
        }
        
        if(!(((3.0*a+b)/4.0)<s<b) || ((flag==1) && (fabs(s-b)>=fabs((b-c)/2.0))) || ((flag==0) && (fabs(s-b)>= fabs((c-d)/2.0))) || ((flag==1) && (fabs(b-c)< J2TOL2)) || ((flag==0) && (fabs(c-d)< J2TOL2)))
        {
            s=(a+b)/2.0;
            flag=1;
        }
        else
        {
            flag=0;
        }
        
        plasticReturnResidual(mu, C0, C1, C2, C3, C4, C5, Tc, n, eqpn, ntbar, dt, s, Gs);
        d=c;
        c=b;
        Gc=Gb;
        
        if(Ga*Gs < 0.0)
        {
            b=s;
            Gb=Gs;
        }
        else
        {
            a=s;
            Ga=Gs;
        }
        
        if (fabs(Ga) < fabs(Gb))
        {
            double cp=a;
            a=b; b=cp;
            double Gcp=Ga;
            Ga=Gb; Gb=Gcp;
        }
        
        count++;
        
    }
    
    if (fabs(Ga) > J2TOL2 || fabs(Gb) > J2TOL2)
    {
        return b;
    }
    return s;
}




void thermoZAMP::CauchyStress(istensor& sigma) const
{
    double Jinv = 1/Jc;
    istensor S;
    secondPiolaKirchhoffStress(S);
    sigma = Jinv * istensor::FSFt(Fc, S);
}




void thermoZAMP::commitCurrentState()
{
    be_n    = be_c;
    iso_n   = iso_c;
    epdot_n = epdot_c;
    tn = tc;
    Fn = Fc;
    Jn = Jc;
    if (thethermoFiniteStrainMaterial.damageModelActivated) Dn = Dc;
    
    time_n  = time_c;
    GradT_n = GradT_c;
    temp_n  = temp_c;
}




// Calls generic numerical convected tangent
void thermoZAMP::convectedTangent(itensor4& ctg) const
{
    this->mechConvectedTangent(ctg);
    
    istensor C = istensor::tensorTransposedTimesTensor(Fc);
    istensor Cinv = C.inverse();

    const double alpha = thethermoFiniteStrainMaterial._thermalExpansion;
    const double bulk  = thethermoFiniteStrainMaterial.getProperty(muesli::PR_BULK);
    const double dtheta= temp_c - thethermoFiniteStrainMaterial.referenceTemperature();
    const double f     = 6.0 * alpha * bulk * dtheta;
    
    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++)
                {
                    ctg(i,j,k,l) += f * 0.5 * ( Cinv(i,k)* Cinv(l,j) + Cinv(i,l) * Cinv(k,j) );
                }
}




void thermoZAMP::convectedTangentTimesSymmetricTensor(const istensor& M,istensor& CM) const
{
    itensor4 C;
    convectedTangent(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);
                }
}




void thermoZAMP::contractWithConvectedTangent(const ivector& v1, const ivector& v2, itensor& T) const
{
    itensor4 c;
    convectedTangent(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);
                }
}




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




double thermoZAMP::dissipatedEnergy() const
{
    const double C0 = thethermoFiniteStrainMaterial._C0;
    const double C1 = thethermoFiniteStrainMaterial._C1;
    const double C2 = thethermoFiniteStrainMaterial._C2;
    const double C3 = thethermoFiniteStrainMaterial._C3;
    const double C4 = thethermoFiniteStrainMaterial._C4;
    const double C5 = thethermoFiniteStrainMaterial._C5;
    const double n  = thethermoFiniteStrainMaterial._n;
    
    const double taylorQuinney = thethermoFiniteStrainMaterial.taylorQuinney;

    double sigma_y = 0.0;
    double thermal = 0.0;

    if(dgamma > 0.0)
    {
        thermal = (epdot_c>0.0) ? (C1+C2*sqrt(iso_c))*exp(-C3*currT+C4*currT*log(epdot_c)) : 0.0;
        if(C5!=0.0)
        {
            double strain = C5*pow(iso_c,n);
            sigma_y = (1-Dc)*SQ23*(thermal + strain + C0);
        }
        else
        {
            sigma_y = (1-Dc)*SQ23*(thermal + C0);
        }
    }
    
    //Plastic thermal dissipation excluded, must be reviewed
    double mech_diss = taylorQuinney*dgamma*sigma_y;
    double dissEnergy = mech_diss;

    return dissEnergy;
}




double thermoZAMP::effectiveFreeEnergy() const
{
    double logJ = log(Jc);

    const double alpha  = thethermoFiniteStrainMaterial._thermalExpansion;
    const double c0     = thethermoFiniteStrainMaterial._thermalCapacity;
    const double bulk   = thethermoFiniteStrainMaterial.getProperty(muesli::PR_BULK);
    const double theta0 = thethermoFiniteStrainMaterial.referenceTemperature();
    const double deltaTemp = temp_c - theta0;

    double Psi_coup = - 3.0 * alpha * bulk * deltaTemp * logJ;
    double Psi_ther = c0 * (deltaTemp - temp_c * log(temp_c/theta0));

    return effectiveStoredEnergy() + Psi_coup + Psi_ther;
}




double thermoZAMP::effectiveStoredEnergy() const
{
    //dissipated energy not included
    return storedEnergy();
}




double thermoZAMP::entropy() const
{
    double logJ = log(Jc);
    
    const double alpha  = thethermoFiniteStrainMaterial._thermalExpansion;
    const double bulk   = thethermoFiniteStrainMaterial.getProperty(muesli::PR_BULK);
    const double c0     = thethermoFiniteStrainMaterial._thermalCapacity;
    const double theta0 = thethermoFiniteStrainMaterial.referenceTemperature();
    
    return 3.0 * alpha * bulk * logJ + c0 * log(temp_c/theta0);
}




void thermoZAMP::explicitRadialReturn(const ivector &taudev, double ep, double epdot)
{
    const double mu  = thethermoFiniteStrainMaterial.mu;
    const double C0  = thethermoFiniteStrainMaterial._C0;
    const double C1  = thethermoFiniteStrainMaterial._C1;
    const double C2  = thethermoFiniteStrainMaterial._C2;
    const double C3  = thethermoFiniteStrainMaterial._C3;
    const double C4  = thethermoFiniteStrainMaterial._C4;
    const double C5  = thethermoFiniteStrainMaterial._C5;
    const double Tc  = currT;
    const double n   = thethermoFiniteStrainMaterial._n;
    
    double thermal(0.0),dthermal(1.0);
    
    if(epdot>0.0)
    {
        thermal  = (C1+C2*sqrt(ep))*exp(-C3*Tc+C4*Tc*log(epdot));
    }
    
    // Depending on the nature of the material (bcc or fcc structure), the algorithm selects
    // the terms to calculate the sigma yield value
    if (C1==0.0 && epdot > 0.0)
    {
        dthermal = ((SQ23*C2/2.0)*(1.0/sqrt(ep))*exp(-C3*Tc+C4*Tc*log(epdot)))+((C2*sqrt(ep)*exp(-C3*Tc+C4*Tc*log(epdot)))*C4*Tc/epdot);
        double sigma_y = SQ23 * (thermal + C0);
        dgamma = (taudev.norm() - sigma_y)/(2.0*mu);
        iso_c = ep + SQ23*dgamma*dthermal;
    }
    
    else if (C2==0.0 && epdot > 0.0)
    {
        dthermal = C1*exp(-C3*Tc+C4*Tc*log(epdot))*C4*Tc/epdot;
        double strain  = C5*pow(ep,n);
        double dstrain = SQ23*n*C5*pow(ep,n-1.0);
        double sigma_y = SQ23 * (thermal + strain + C0);
        dgamma = (taudev.norm() - sigma_y)/(2.0*mu);
        iso_c = ep + SQ23*dgamma*(dthermal + dstrain);
    }
    
    else if(C1==0.0 && epdot==0.0)
    {
        dthermal=1.0;
        double sigma_y = SQ23 * (thermal + C0);
        dgamma = (taudev.norm() - sigma_y)/(2.0*mu);
        iso_c = ep + SQ23*dgamma*dthermal;
    }
    
    else if(C2==0.0 && epdot==0.0)
    {
        dthermal=0.0;
        double strain  = C5*pow(ep,n);
        double dstrain = SQ23*n*C5*pow(ep,n-1.0);
        double sigma_y = SQ23 * (thermal + strain + C0);
        dgamma = (taudev.norm() - sigma_y)/(2.0*mu);
        iso_c = ep + SQ23*dgamma*(dthermal + dstrain);
    }
}




void thermoZAMP::firstPiolaKirchhoffStress(itensor &P) const
{
    istensor S;
    secondPiolaKirchhoffStress(S);
    P = Fc*S;
}




double thermoZAMP::freeEnergy() const
{
    const double alpha = thethermoFiniteStrainMaterial._thermalExpansion;
    const double c0    = thethermoFiniteStrainMaterial._thermalCapacity;
    const double bulk  = thethermoFiniteStrainMaterial.getProperty(muesli::PR_BULK);
    const double theta0= thethermoFiniteStrainMaterial.referenceTemperature();
    const double deltaTemp   = temp_c - theta0;
    
    double logJ = log(Jc);
    double Psi_coup = - 3.0 * alpha * bulk * deltaTemp * logJ;
    double Psi_ther = c0 * (deltaTemp - temp_c * log(temp_c/theta0));
    
    return storedEnergy() + Psi_coup + Psi_ther;
}




materialState thermoZAMP::getConvergedState() const
{
    materialState state;
    
    state.theDouble.push_back(Jn);
    state.theTensor.push_back(Fn);
    if(thethermoFiniteStrainMaterial.damageModelActivated) state.theDouble.push_back(Dn);
    
    state.theTime = time_n;
    state.theDouble.push_back(temp_n);
    state.theVector.push_back(GradT_n);
    
    return state;
}




materialState thermoZAMP::getCurrentState() const
{
    materialState state;

    state.theDouble.push_back(Jc);
    state.theTensor.push_back(Fc);
    if(thethermoFiniteStrainMaterial.damageModelActivated) state.theDouble.push_back(Dc);
    
    state.theTime = time_c;
    state.theDouble.push_back(temp_c);
    state.theVector.push_back(GradT_c);

    return state;
}




double thermoZAMP::heatCapacity() const
{
    return thethermoFiniteStrainMaterial._thermalCapacity;
}




double thermoZAMP::kineticPotential() const
{
    const double mu = thethermoFiniteStrainMaterial.mu;
    const double C2 = thethermoFiniteStrainMaterial._C2;
    const double C3 = thethermoFiniteStrainMaterial._C3;
    const double C4 = thethermoFiniteStrainMaterial._C4;
    const double C5 = thethermoFiniteStrainMaterial._C5;
    const double n  = thethermoFiniteStrainMaterial._n;

    // logarithmic principal elastic stretches
    ivector  neigvec[3], lambda2n;
    const ivector vone(1.0, 1.0, 1.0);
    be_n.spectralDecomposition(neigvec, lambda2n);
    ivector epse_c, epse_n;
    for (size_t i=0; i<3; i++)
    {
        epse_c[i] = 0.5*log(lambda2[i]);
        epse_n[i] = 0.5*log(lambda2n[i]);
    }

    ivector devepse_c, devepse_n;
    devepse_c = epse_c - 1.0/3.0 * epse_c.dot(vone)*vone;
    devepse_n = epse_n - 1.0/3.0 * epse_n.dot(vone)*vone;
    
    ivector s = 2.0*mu*devepse_c;
    double q = (1-Dc)*(-C5*pow(iso_c,n) - C2*pow(iso_c,0.5) * exp(-C3*currT + C4*currT*log(epdot_c)));
    
    ivector dt_dp = dgamma*nubarTR;
    double dt_KinPot = dt_dp.dot(s) + (iso_c-iso_n)*q;
    
    double dt = tc - tn;
    double K = (dt == 0.0) ? 0.0 : dt_KinPot/dt;
    return K - thermalPotential();
}




void thermoZAMP::KirchhoffStress(istensor &tau) const
{
    istensor sigma;
    CauchyStress(sigma);
    tau = sigma * Jc;
}




void thermoZAMP::KirchhoffStressVector(double tauv[6]) const
{
    istensor tau;
    KirchhoffStress(tau);
    muesli::tensorToVector(tau, tauv);
}




itensor thermoZAMP::materialCouplingTensor() const
{
    itensor Fcinv = Fc.inverse();
    const ivector vone(1.0, 1.0, 1.0);
    
    const double alpha = thethermoFiniteStrainMaterial._thermalExpansion;
    const double bulk  = thethermoFiniteStrainMaterial.getProperty(muesli::PR_BULK);
    
    itensor mechPKDTheta;
    mechPKDTheta.setZero();
    
    if (this->dgamma>0.0)
    {
        mechFirstPKStressDTheta(mechPKDTheta);
    }

    return temp_c * (mechPKDTheta - 3.0 * alpha * bulk * Fcinv.transpose());
}




itensor4 thermoZAMP::materialCouplingTensorDF() const
{
    itensor Fct = Fc.transpose();
    itensor Fctinv = Fct.inverse();
    itensor Fcinv = Fc.inverse();
    
    const double C1 = thethermoFiniteStrainMaterial._C1;
    const double C2 = thethermoFiniteStrainMaterial._C2;
    const double C3 = thethermoFiniteStrainMaterial._C3;
    const double C4 = thethermoFiniteStrainMaterial._C4;
    const double C5 = thethermoFiniteStrainMaterial._C5;
    const double n  = thethermoFiniteStrainMaterial._n;
    const double Tc = currT;
    const double mu = thethermoFiniteStrainMaterial.mu;
    
    double thermal = (epdot_c>1e-20) ? (C1+C2*sqrt(iso_c))*exp(-C3*Tc+C4*Tc*log(epdot_c)) : 0.0;
    double sigmayDTheta, sigmayDgamma;
    sigmayDTheta = sigmayDgamma = 0.0;
    double dt = tc-tn;
    
    if(iso_c>0.0 && epdot_c>1e-20)
    {
        sigmayDgamma = (1.0-Dc)*sqrt(2.0/3.0)*(0.5*sqrt(2.0/3.0)*C2*pow(iso_c,-0.5)*exp(-C3*Tc+C4*Tc*log(epdot_c)) +
                       thermal*C4*Tc*sqrt(2.0/3.0)/(epdot_c*dt) + C5*n*pow(iso_c,n-1)*sqrt(2.0/3.0));
    }

    if (epdot_c>1e-20)
    {
        sigmayDTheta = (1.0-Dc)*sqrt(2.0/3.0)*thermal*(-C3 + C4*log(epdot_c));
    }
    double dgammaDTheta = -1/(2*mu + sigmayDgamma)*sigmayDTheta;
    ivector tauDtheta = -2*mu*dgammaDTheta*nubarTR;
    
    // Reconstruct derivative wrt theta of Cauchy stress
    const double iJ = 1.0/Jc;
    istensor sigmaDTheta;
    sigmaDTheta.setZero();
    for (size_t i=0; i<3; i++)
    {
        sigmaDTheta.addScaledVdyadicV(iJ * tauDtheta[i], nn[i]);
    }
    
    const double alpha = thethermoFiniteStrainMaterial._thermalExpansion;
    const double bulk  = thethermoFiniteStrainMaterial.getProperty(muesli::PR_BULK);
    const double f1    = -3.0*alpha*bulk*temp_c;
    
    itensor4 F_t_F1, F_t_F2, sigmaDthetaDF, F_t_F3, F_t_F4, dFTdF;
    F_t_F1.setZero();
    F_t_F4 = F_t_F3 = F_t_F2 = F_t_F1;
    
    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++)
                {
                    F_t_F4(i,j,k,l) += -f1 * Fcinv(l,i)*Fcinv(j,k);
                    dFTdF(i,j,k,l) += -Fcinv(l,i)*Fcinv(j,k);
                }
    
    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++)
                    for (unsigned m=0; m<3; m++)
                        F_t_F1(i,j,k,l)=temp_c*Jc*Fctinv(k,l)*sigmaDTheta(i,m)*Fctinv(m,j);
    
    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++)
                    for (unsigned m=0; m<3; m++)
                        F_t_F3(i,j,k,l) += temp_c*Jc*sigmaDTheta(i,m)*dFTdF(m,j,k,l);

    //numerical derivation to find d^2sigma/dTheta dF, as a part of F_t_F2,
    //only for testing purposes (to be deleted once function is solved)
    itensor F = Fc;
    const double inc = 1.0e-4;
    thermoZAMP& theMP = const_cast<thermoZAMP&>(*this);
    itensor num_sigmadthetadF;
    
    for (size_t i=0; i<3; i++)
    {
        for (size_t j=0; j<3; j++)
        {
            const double original = F(i,j);

            F(i,j) = original + inc;
            theMP.updateCurrentState(tc, F, GradT_c, temp_c);
            double thermal1 = (epdot_c>0.0) ? (C1+C2*sqrt(iso_c))*exp(-C3*Tc+C4*Tc*log(epdot_c)) : 0.0;
            double sigmayDTheta1, sigmayDgamma1;
            sigmayDTheta1 = sigmayDgamma1 = 0.0;
            double dt = tc-tn;
            
            if(iso_c>0.0)
            {
                sigmayDgamma1 = (1.0-Dc)*sqrt(2.0/3.0)*(0.5*sqrt(2.0/3.0)*C2*pow(iso_c,-0.5)*exp(-C3*Tc+C4*Tc*log(epdot_c)) +
                               thermal1*C4*Tc*sqrt(2.0/3.0)/(epdot_c*dt) + C5*n*pow(iso_c,n-1)*sqrt(2.0/3.0));
            }

            if (epdot_c>0.0)
            {
                sigmayDTheta1 = (1.0-Dc)*sqrt(2.0/3.0)*thermal*(-C3 + C4*log(epdot_c));
            }
            double dgammaDTheta1 = -1/(2*mu + sigmayDgamma1)*sigmayDTheta1;
            ivector tauDtheta1 = -2*mu*dgammaDTheta1*nubarTR;
            
            // Reconstruct derivative wrt theta of Cauchy stress
            double iJ = 1.0/Jc;
            istensor sigmaDTheta1;
            sigmaDTheta1.setZero();
            for (size_t i=0; i<3; i++)
            {
                sigmaDTheta1.addScaledVdyadicV(iJ * tauDtheta1[i], nn[i]);
            }
            
            

            F(i,j) = original + 2.0*inc;
            theMP.updateCurrentState(tc, F, GradT_c, temp_c);
            thermal1 = (epdot_c>0.0) ? (C1+C2*sqrt(iso_c))*exp(-C3*Tc+C4*Tc*log(epdot_c)) : 0.0;
            sigmayDTheta1 = sigmayDgamma1 = 0.0;
            dt = tc-tn;
            
            if(iso_c>0.0)
            {
                sigmayDgamma1 = (1.0-Dc)*sqrt(2.0/3.0)*(0.5*sqrt(2.0/3.0)*C2*pow(iso_c,-0.5)*exp(-C3*Tc+C4*Tc*log(epdot_c)) +
                               thermal1*C4*Tc*sqrt(2.0/3.0)/(epdot_c*dt) + C5*n*pow(iso_c,n-1)*sqrt(2.0/3.0));
            }

            if (epdot_c>0.0)
            {
                sigmayDTheta1 = (1.0-Dc)*sqrt(2.0/3.0)*thermal*(-C3 + C4*log(epdot_c));
            }
            dgammaDTheta1 = -1/(2*mu + sigmayDgamma1)*sigmayDTheta1;
            tauDtheta1 = -2*mu*dgammaDTheta1*nubarTR;
            
            // Reconstruct derivative wrt theta of Cauchy stress
            iJ = 1.0/Jc;
            istensor sigmaDTheta2;
            sigmaDTheta2.setZero();
            for (size_t i=0; i<3; i++)
            {
                sigmaDTheta2.addScaledVdyadicV(iJ * tauDtheta1[i], nn[i]);
            }

            
            F(i,j) = original - inc;
            theMP.updateCurrentState(tc, F, GradT_c, temp_c);
            thermal1 = (epdot_c>0.0) ? (C1+C2*sqrt(iso_c))*exp(-C3*Tc+C4*Tc*log(epdot_c)) : 0.0;
            sigmayDTheta1 = sigmayDgamma1 = 0.0;
            dt = tc-tn;
            
            if(iso_c>0.0)
            {
                sigmayDgamma1 = (1.0-Dc)*sqrt(2.0/3.0)*(0.5*sqrt(2.0/3.0)*C2*pow(iso_c,-0.5)*exp(-C3*Tc+C4*Tc*log(epdot_c)) +
                               thermal1*C4*Tc*sqrt(2.0/3.0)/(epdot_c*dt) + C5*n*pow(iso_c,n-1)*sqrt(2.0/3.0));
            }

            if (epdot_c>0.0)
            {
                sigmayDTheta1 = (1.0-Dc)*sqrt(2.0/3.0)*thermal*(-C3 + C4*log(epdot_c));
            }
            dgammaDTheta1 = -1/(2*mu + sigmayDgamma1)*sigmayDTheta1;
            tauDtheta1 = -2*mu*dgammaDTheta1*nubarTR;
            
            // Reconstruct derivative wrt theta of Cauchy stress
            iJ = 1.0/Jc;
            istensor sigmaDTheta3;
            sigmaDTheta3.setZero();
            for (size_t i=0; i<3; i++)
            {
                sigmaDTheta3.addScaledVdyadicV(iJ * tauDtheta1[i], nn[i]);
            }

            
            F(i,j) = original - 2.0*inc;
            theMP.updateCurrentState(tc, F, GradT_c, temp_c);
            thermal1 = (epdot_c>0.0) ? (C1+C2*sqrt(iso_c))*exp(-C3*Tc+C4*Tc*log(epdot_c)) : 0.0;
            sigmayDTheta1 = sigmayDgamma1 = 0.0;
            dt = tc-tn;
            
            if(iso_c>0.0)
            {
                sigmayDgamma1 = (1.0-Dc)*sqrt(2.0/3.0)*(0.5*sqrt(2.0/3.0)*C2*pow(iso_c,-0.5)*exp(-C3*Tc+C4*Tc*log(epdot_c)) +
                               thermal1*C4*Tc*sqrt(2.0/3.0)/(epdot_c*dt) + C5*n*pow(iso_c,n-1)*sqrt(2.0/3.0));
            }

            if (epdot_c>0.0)
            {
                sigmayDTheta1 = (1.0-Dc)*sqrt(2.0/3.0)*thermal*(-C3 + C4*log(epdot_c));
            }
            dgammaDTheta1 = -1/(2*mu + sigmayDgamma1)*sigmayDTheta1;
            tauDtheta1 = -2*mu*dgammaDTheta1*nubarTR;
            
            // Reconstruct derivative wrt theta of Cauchy stress
            iJ = 1.0/Jc;
            istensor sigmaDTheta4;
            sigmaDTheta4.setZero();
            for (size_t i=0; i<3; i++)
            {
                sigmaDTheta4.addScaledVdyadicV(iJ * tauDtheta1[i], nn[i]);
            }

            
            // fourth order approximation of the derivative
            num_sigmadthetadF = (-sigmaDTheta2 + 8.0*sigmaDTheta1 - 8.0*sigmaDTheta3 + sigmaDTheta4)/(12.0*inc);

            for (unsigned k=0; k<3; k++)
            {
                for (unsigned l=0; l<3; l++)
                {
                    sigmaDthetaDF(k,l,i,j) = num_sigmadthetadF(k,l);
                }
            }
            F(i,j) = original;
            theMP.updateCurrentState(tc, F, GradT_c, temp_c);
        }
    }
    
    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++)
                    for (unsigned m=0; m<3; m++)
                        F_t_F2(i,j,k,l) += temp_c*Jc*sigmaDthetaDF(i,m,k,l)*Fctinv(m,j);
    
    //return F_t_F1 + F_t_F2 + F_t_F3 + F_t_F4;
    return F_t_F1 + F_t_F3 + F_t_F4;
}




itensor thermoZAMP::materialCouplingTensorDTheta() const
{
    itensor Fcinv = Fc.inverse();
    
    const double alpha = thethermoFiniteStrainMaterial._thermalExpansion;
    const double bulk  = thethermoFiniteStrainMaterial.getProperty(muesli::PR_BULK);
    
    itensor mechPKDTheta, mechPKDDTheta;
    mechPKDTheta.setZero();
    mechPKDDTheta.setZero();
    
    if (this->dgamma>0.0)
    {
        mechFirstPKStressDTheta(mechPKDTheta);
        mechFirstPKStressDDTheta(mechPKDDTheta);
    }
    
    return mechPKDTheta - 3.0*alpha*bulk*Fcinv.transpose() + temp_c*mechPKDDTheta;
}




void thermoZAMP::materialTangent(itensor4& cm) const
{
    istensor S;
    secondPiolaKirchhoffStress(S);

    itensor4 cc;
    convectedTangent(cc);

    cm.setZero();
    for (unsigned a=0; a<3; a++)
        for (unsigned b=0; b<3; b++)
            for (unsigned A=0; A<3; A++)
                for (unsigned B=0; B<3; B++)
                {
                    if (a == b) cm(a,A,b,B) += S(A,B);

                    for (unsigned C=0; C<3; C++)
                        for (unsigned D=0; D<3; D++)
                            cm(a,A,b,B) += Fc(a,C) * Fc(b,D) * cc(C,A,D,B);
                }
}




void thermoZAMP::mechCauchyStress(istensor& sigma) const
{
    const double iJ = 1.0/Jc;
    
    // Reconstruct Cauchy stress
    sigma.setZero();
    for (size_t i=0; i<3; i++)
    {
        sigma.addScaledVdyadicV(iJ * tau[i], nn[i]);
    }
}




void thermoZAMP::mechFirstPiolaKirchhoffStress(itensor &P) const
{
    istensor sigma;
    mechCauchyStress(sigma);
    P = Jc * sigma * Fc.inverse().transpose();
}




void thermoZAMP::mechFirstPKStressDDTheta(itensor &DDP) const
{
    const double C1 = thethermoFiniteStrainMaterial._C1;
    const double C2 = thethermoFiniteStrainMaterial._C2;
    const double C3 = thethermoFiniteStrainMaterial._C3;
    const double C4 = thethermoFiniteStrainMaterial._C4;
    const double C5 = thethermoFiniteStrainMaterial._C5;
    const double Tc = currT;
    const double n  = thethermoFiniteStrainMaterial._n;
    const double mu = thethermoFiniteStrainMaterial.mu;

    double thermal = (epdot_c>1e-20) ? (C1+C2*sqrt(iso_c))*exp(-C3*Tc+C4*Tc*log(epdot_c)) : 0.0;
    double sigmayDTheta, sigmayDgamma, sigmayDgammaDTheta, sigmayDDTheta;
    sigmayDTheta = sigmayDgamma = sigmayDgammaDTheta = sigmayDDTheta = 0.0;
    double dt = tc-tn;
    
    if(iso_c>0.0 && epdot_c>1e-20)
    {
        sigmayDgamma = (1.0-Dc)*sqrt(2.0/3.0)*(0.5*sqrt(2.0/3.0)*C2*pow(iso_c,-0.5)*exp(-C3*Tc+C4*Tc*log(epdot_c)) +
                       thermal*C4*Tc*sqrt(2.0/3.0)/(epdot_c*dt) + C5*n*pow(iso_c,n-1)*sqrt(2.0/3.0));
    }

    if (epdot_c>0.0)
    {
        sigmayDTheta = (1.0-Dc)*sqrt(2.0/3.0)*thermal*(-C3 + C4*log(epdot_c));
    }
    double dgammaDTheta = -1/(2*mu + sigmayDgamma)*sigmayDTheta;
    
    if(iso_c>0.0 && epdot_c>1e-20)
    {
        sigmayDgammaDTheta = (1.0-Dc)*sqrt(2.0/3.0)*(-0.25*2.0/3.0*C2*pow(iso_c,-1.5)*dgammaDTheta*exp(-C3*Tc+C4*Tc*log(epdot_c)) +
                             0.5*sqrt(2.0/3.0)*C2*pow(iso_c,-0.5)*exp(-C3*Tc+C4*Tc*log(epdot_c))*(-C3 + C4*log(epdot_c) + C4*Tc*sqrt(2.0/3.0)*dgammaDTheta/(epdot_c*dt)) +
                             0.5*sqrt(2.0/3.0)*C2*pow(iso_c,-0.5)*dgammaDTheta*exp(-C3*Tc + C4*Tc*log(epdot_c))*C4*Tc*sqrt(2.0/3.0)/(epdot_c*dt) +
                             thermal*(-C3 + C4*log(epdot_c) + C4*Tc*dgammaDTheta*sqrt(2.0/3.0)/(epdot_c*dt))*C4*Tc*sqrt(2.0/3.0)/(epdot_c*dt) +
                             thermal*(C4*sqrt(2.0/3.0)/(epdot_c*dt)-C4*Tc*2.0/3.0*dgammaDTheta/(epdot_c*dt*epdot_c*dt)) +
                             C5*n*(n-1)*pow(iso_c,n-2)*2.0/3.0*dgammaDTheta);
    }

    if(iso_c>0.0 && epdot_c>1e-20)
    {
        sigmayDDTheta = (1.0-Dc)*sqrt(2.0/3.0)*(sqrt(2.0/3.0)*0.5*C2*pow(iso_c,-0.5)*dgammaDTheta*exp(-C3*Tc+C4*Tc*log(epdot_c))*(-C3 + C4*log(epdot_c)) +
                        thermal*(-C3 + C4*log(epdot_c) + C4*Tc*dgammaDTheta*sqrt(2.0/3.0)/(epdot_c*dt))*(-C3 + C4*log(epdot_c)) +
                        thermal*C4*dgammaDTheta*sqrt(2.0/3.0)/(epdot_c*dt));
    }
    
    double dgammaDDTheta = pow(2*mu+sigmayDgamma , -2)*sigmayDgammaDTheta*sigmayDTheta - 1/(2*mu+sigmayDgamma)*sigmayDDTheta;
    ivector tauDDtheta = -2*mu*dgammaDDTheta*nubarTR;
    
    // Reconstruct Cauchy stress
    const double iJ = 1.0/Jc;
    istensor sigmaDDtheta;
    sigmaDDtheta.setZero();
    for (size_t i=0; i<3; i++)
    {
        sigmaDDtheta.addScaledVdyadicV(iJ * tauDDtheta[i], nn[i]);
    }
    
    DDP = Jc * sigmaDDtheta * Fc.inverse().transpose();
}




void thermoZAMP::mechFirstPKStressDTheta(itensor &DP) const
{
    const double C1 = thethermoFiniteStrainMaterial._C1;
    const double C2 = thethermoFiniteStrainMaterial._C2;
    const double C3 = thethermoFiniteStrainMaterial._C3;
    const double C4 = thethermoFiniteStrainMaterial._C4;
    const double C5 = thethermoFiniteStrainMaterial._C5;
    const double Tc = currT;
    const double n  = thethermoFiniteStrainMaterial._n;
    const double mu = thethermoFiniteStrainMaterial.mu;
    
    double thermal = (epdot_c>1e-20) ? (C1+C2*sqrt(iso_c))*exp(-C3*Tc+C4*Tc*log(epdot_c)) : 0.0;
    double sigmayDTheta, sigmayDgamma;
    sigmayDTheta = sigmayDgamma = 0.0;
    double dt = tc-tn;
    
    if(iso_c>0.0 && epdot_c>1e-20)
    {
        sigmayDgamma = (1.0-Dc)*sqrt(2.0/3.0)*(0.5*sqrt(2.0/3.0)*C2*pow(iso_c,-0.5)*exp(-C3*Tc+C4*Tc*log(epdot_c)) +
                       thermal*C4*Tc*sqrt(2.0/3.0)/(epdot_c*dt) + C5*n*pow(iso_c,n-1)*sqrt(2.0/3.0));
    }

    if (epdot_c>1e-20)
    {
        sigmayDTheta = (1.0-Dc)*sqrt(2.0/3.0)*thermal*(-C3 + C4*log(epdot_c));
    }
    double dgammaDTheta = -1/(2*mu + sigmayDgamma)*sigmayDTheta;
    ivector tauDtheta = -2*mu*dgammaDTheta*nubarTR;
    
    // Reconstruct derivative wrt theta of Cauchy stress
    const double iJ = 1.0/Jc;
    istensor sigmaDTheta;
    sigmaDTheta.setZero();
    for (size_t i=0; i<3; i++)
    {
        sigmaDTheta.addScaledVdyadicV(iJ * tauDtheta[i], nn[i]);
    }
    
    DP = Jc * sigmaDTheta * Fc.inverse().transpose();
}




void thermoZAMP::mechConvectedTangent(itensor4& nTg) const
{
    itensor4 A;
    itensor dP, Pp1, Pp2, Pm1, Pm2;

    // numerical differentiation stress
    itensor numP;
    numP.setZero();
    
    itensor Fnp;
    Fnp=Fc;
    
    double theta = this->currT;
    ivector gradTp = this->GradT_c;
    
    const double inc = 1.0e-5*Fc.norm();
    double tn1 = this->tc;
    
    thermoZAMP& theMPc = const_cast<thermoZAMP&>(*this);
    
    for (unsigned i=0; i<3; i++)
    {
        for (unsigned j=0; j<3; j++)
        {
            double original = this->Fc(i,j);
            
            Fnp(i,j) = original + inc;
            theMPc.updateCurrentState(tn1, Fnp, gradTp, theta);
            theMPc.mechFirstPiolaKirchhoffStress(Pp1);

            Fnp(i,j) = original + 2.0*inc;
            theMPc.updateCurrentState(tn1, Fnp, gradTp, theta);
            theMPc.mechFirstPiolaKirchhoffStress(Pp2);

            Fnp(i,j) = original - inc;
            theMPc.updateCurrentState(tn1, Fnp, gradTp, theta);
            theMPc.mechFirstPiolaKirchhoffStress(Pm1);

            Fnp(i,j) = original - 2.0*inc;
            theMPc.updateCurrentState(tn1, Fnp, gradTp, theta);
            theMPc.mechFirstPiolaKirchhoffStress(Pm2);

            // derivative of PK stress
            dP = (-Pp2 + 8.0*Pp1 - 8.0*Pm1 + Pm2)/(12.0*inc);
            for (unsigned k=0; k<3; k++)
                for (unsigned l=0; l<3; l++)
                    A(k,l,i,j) = dP(k,l);

            Fnp(i,j) = original;
            theMPc.updateCurrentState(tn1, Fnp, gradTp, theta);
        }
    }

    // transform A to get the convected tangent
    itensor  J  = Fnp.inverse();
    istensor Cn = istensor::tensorTransposedTimesTensor(Fnp);
    istensor Ci = Cn.inverse();
    istensor S;
    theMPc.mechSecondPiolaKirchhoffStress(S);
    for (unsigned a=0; a<3; a++)
        for (unsigned b=0; b<3; b++)
            for (unsigned c=0; c<3; c++)
                for (unsigned d=0; d<3; d++)
                {
                    nTg(c,a,d,b) = - S(a,b)*Ci(c,d);

                    for (unsigned i=0; i<3; i++)
                        for (unsigned j=0; j<3; j++)
                                nTg(c,a,d,b) += J(c,i)*A(i,a,j,b)*J(d,j);
                }

}




void thermoZAMP::mechSecondPiolaKirchhoffStress(itensor &S) const
{
    istensor sigma;
    mechCauchyStress(sigma);
    S =  Jc * istensor::FSFt(Fc.inverse(), sigma);
}




void thermoZAMP::mechSpatialTangent(itensor4& st) const
{
    st.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++)
                {
                    if (i==j && k==l) st(i,j,k,l) += thethermoFiniteStrainMaterial.lambda;
                    if (i==k && j==l) st(i,j,k,l) += thethermoFiniteStrainMaterial.mu;
                    if (i==l && j==k) st(i,j,k,l) += thethermoFiniteStrainMaterial.mu;
                }
}




void  thermoZAMP::numericalFirstPiolaKirchhoffStressDTheta(itensor& P) const
{
    itensor Pp1, Pp2, Pm1, Pm2;
    thermoZAMP& theMPc = const_cast<thermoZAMP&>(*this);
    
    // numerical differentiation stress
    const double inc = 1.0e-4;
    double temp = 0.0;
    
    double tn1 = this->tc;
    const double original = this->currT;
    itensor Fnp = this->Fc;
    ivector gradTp = this->GradT_c;
    
    temp = original + inc;
    theMPc.updateCurrentState(tn1, Fnp, gradTp, temp);
    mechFirstPiolaKirchhoffStress(Pp1);
    
    temp = original + 2.0*inc;
    theMPc.updateCurrentState(tn1, Fnp, gradTp, temp);
    mechFirstPiolaKirchhoffStress(Pp2);
    
    temp = original - inc;
    theMPc.updateCurrentState(tn1, Fnp, gradTp, temp);
    mechFirstPiolaKirchhoffStress(Pm1);
    
    temp = original - 2.0*inc;
    theMPc.updateCurrentState(tn1, Fnp, gradTp, temp);
    mechFirstPiolaKirchhoffStress(Pm2);
    
    temp = original;
    theMPc.updateCurrentState(tn1, Fnp, gradTp, temp);
    
    P = temp * (-Pp2 + 8.0*Pp1 - 8.0*Pm1 + Pm2)/(12.0*inc);
}




double thermoZAMP::plasticSlip() const
{
    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 thermoZAMP::plasticReturn(const ivector& taubarTR)
{
    const double mu      = thethermoFiniteStrainMaterial.mu;
    const double C0      = thethermoFiniteStrainMaterial._C0;
    const double C1      = thethermoFiniteStrainMaterial._C1;
    const double C2      = thethermoFiniteStrainMaterial._C2;
    const double C3      = thethermoFiniteStrainMaterial._C3;
    const double C4      = thethermoFiniteStrainMaterial._C4;
    const double C5      = thethermoFiniteStrainMaterial._C5;
    const double Tc      = currT;
    const double n       = thethermoFiniteStrainMaterial._n;
    const double dt      = tc - tn;
    
    
    double G, DG;
    double eqpn     = iso_n;
    double ntbar    = taubarTR.norm();
    int flagbrent   = 0;
    int flagbrent2  = 1;
    double x        = dt;
    
    plasticReturnResidual(mu, C0, C1, C2, C3, C4, C5, Tc, n, eqpn, ntbar, dt, x, G);
    
    
    // NR scheme. "x" represents the unknown delta-gamma
    size_t count = 0;
    while ( fabs(G) > J2TOL1 && ++count < GAMMAITER1 && x>0.0)
    {
        plasticReturnTangent(mu, C0, C1, C2, C3, C4, C5, Tc, n, eqpn, ntbar, dt, x, DG);
        x -= G/DG;
        plasticReturnResidual(mu, C0, C1, C2, C3, C4, C5, Tc, n, eqpn, ntbar, dt, x, G);
        
    }
    
    // The "if" condition tests whether the convergence with NR scheme has succeded or not, in which case
    // will set flagbrent to 1.
    if (count >= GAMMAITER1 || std::isnan(G) || x<=0.0)
    {
        flagbrent=1;
    }
    else
    {
        dgamma = x;
        iso_c  = iso_n + SQ23*dgamma;
    }
    
    // Initially, two points with different function values wihtin G must be given to the Brent method
    // to assure its convergence. Ga and Gb are both multiplied/divided until their sign is opposite.
    if(flagbrent==1)
    {
        eqpn = iso_n;
        ntbar = taubarTR.norm();
        double a=1e-12, b=1e-5;
        double Ga, Gb;
        
        plasticReturnResidual(mu, C0, C1, C2, C3, C4, C5, Tc, n, eqpn, ntbar, dt, a, Ga);
        plasticReturnResidual(mu, C0, C1, C2, C3, C4, C5, Tc, n, eqpn, ntbar, dt, b, Gb);
        
        while(Gb<0.0)
        {
            b=b*100.0;
            plasticReturnResidual(mu, C0, C1, C2, C3, C4, C5, Tc, n, eqpn, ntbar, dt, b, Gb);
            if (b>2.0) break;
        }
        
        // The value of "a" can be very close to zero due to the Zerilli-Armstrong equation features
        // and because of this the division loop is stopped, making not a big difference compared
        // to the real value, given that it is really small.
        while(Ga>0.0)
        {
            a=a/150.0;
            plasticReturnResidual(mu, C0, C1, C2, C3, C4, C5, Tc, n, eqpn, ntbar, dt, a, Ga);
            if (a<1e-50)
            {
                dgamma = 1e-50;
                flagbrent2 = 0;
                break;
            }
        }
        
        if(flagbrent == 1 && flagbrent2 == 1)
        {
            // Call to Brent method function in case NR fails and "a" value is high enough
            dgamma = brentroot(a, b, Ga, Gb, eqpn, ntbar, mu, C0, C1, C2, C3, C4, C5, Tc, n, dt);
            iso_c  = iso_n + SQ23*dgamma;
        }
    }
}




// Calculation of G function value given each iteration value of dgamma
void thermoZAMP::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)
{
    double deqp = SQ23*dgamma;
    double eqp  = eqpn + deqp;

    double thermal (0.0);
    if(deqp > 0.0 && dt > 0.0)
    {
        thermal = (C1+C2*sqrt(eqp))*exp(-C3*Tc+C4*Tc*log(deqp/dt));
    }
    
    if(C5!=0.0)
    {
        double strain  = C5*pow(eqp,n);
        G = 2.0*mu*dgamma - tau + SQ23*(thermal+strain+C0);
    }
    
    if(C5==0.0)
    {
        G = 2.0*mu*dgamma - tau + SQ23*(thermal+C0);
    }
}



// Calculation of G function derivative value given each iteration value of dgamma
void  thermoZAMP::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)
{
    double deqp = SQ23*dgamma;
    double eqp  = eqpn + deqp;
    
    if(C1 == 0.0 && dgamma > 0.0 && dt > 0.0)
    {
        double dthermal = ((SQ23*C2/2.0)*(1.0/sqrt(eqp))*exp((-C3*Tc)+(C4*Tc*log(deqp/dt))))+((C2*sqrt(eqp)*exp((-C3*Tc)+(C4*Tc*log(deqp/dt))))*(C4*Tc/dgamma));
        DG = 2.0*mu + SQ23*dthermal;
    }
    
    else if(C2 == 0.0 && dgamma > 0.0 && dt > 0.0)
    {
        double dthermal = C1*exp(-C3*Tc+C4*Tc*log(deqp/dt))*C4*Tc/dgamma;
        double dstrain = SQ23*n*C5*pow(eqp,n-1.0);
        DG = 2.0*mu + SQ23*(dthermal + dstrain);
    }
    
    else if(C1 == 0.0 && dgamma == 0.0)
    {
        double dthermal = 1.0;
        DG = 2.0*mu + SQ23*dthermal;
    }
    
    else if(C2 == 0.0 && dgamma == 0.0)
    {
        double dthermal = 1.0;
        double dstrain = SQ23*n*C5*pow(eqp,n-1.0);
        DG = 2.0*mu + SQ23*(dthermal + dstrain);
    }
}




void thermoZAMP::resetCurrentState()
{
    time_n  = time_c;
    GradT_c = GradT_n;
    temp_c  = temp_n;
    
    iso_c   = iso_n;
    be_c    = be_n;
    epdot_c = epdot_n;
    
    tc = tn;
    Fc = Fn;
    Jc = Jn;
    if(thethermoFiniteStrainMaterial.damageModelActivated) Dc = Dn;
}




void thermoZAMP::secondPiolaKirchhoffStress(istensor& S) const
{
    mechSecondPiolaKirchhoffStress(S);
    
    itensor  Finv = Fc.inverse();
    istensor Cinv = istensor::tensorTimesTensorTransposed(Finv);

    const double alpha  = thethermoFiniteStrainMaterial._thermalExpansion;
    const double bulk   = thethermoFiniteStrainMaterial.getProperty(muesli::PR_BULK);
    const double theta0 = thethermoFiniteStrainMaterial.referenceTemperature();
    const double dTemp  = temp_c - theta0;

    S -= 3.0 * alpha * bulk * dTemp * Cinv;
}




void thermoZAMP::setConvergedState(const double theTime, const itensor& F,
                                   const double iso, const ivector& kine,
                                   const istensor& be)
{
    tn     = theTime;
    Fn     = F;
    Jn     = Fn.determinant();
    iso_n  = iso;
    be_n   = be;
}




void thermoZAMP::setRandom()
{
    iso_c = iso_n = muesli::randomUniform(1.0, 2.0);
    ivector tmp; tmp.setRandom();
    ivector vone(1.0, 1.0, 1.0);
    itensor Fe; Fe.setRandom();
    
    if (Fe.determinant() < 0.0) Fe *= -1.0;
    be_c = be_n = istensor::tensorTimesTensorTransposed(Fe);
    
    itensor Fp; Fp.setRandom();
    
    Fp *= 1.0/cbrt(Fp.determinant());
    Fc = Fn = Fe*Fp;
    Jc = Fc.determinant();
    if (thethermoFiniteStrainMaterial.damageModelActivated)
    {
        Dc = muesli::randomUniform(0.0, 1.0);
    }
}




istensor thermoZAMP::spatialConductivity() const
{
    double Jinv = 1/Jc;
    istensor k = thethermoFiniteStrainMaterial._thermalConductivity*istensor::identity();
    
    return Jinv * istensor::FSFt(Fc, k);
}




ivector thermoZAMP::spatialHeatflux() const
{
    const ivector H = materialHeatflux();

    return (1.0/Jc) * Fc * 0.5 * H;
}




void thermoZAMP::spatialTangent(itensor4& Cs) const
{
    /*COMMENTED JUST TO TEST THAT IMPLICIT SIMS WORK WITH SPATIAL TANGENT!
    itensor4 Cc;
    convectedTangent(Cc);
    
    Cs.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++)
                    for (unsigned m=0; m<3; m++)
                        for (unsigned n=0; n<3; n++)
                            for (unsigned p=0; p<3; p++)
                                for (unsigned q=0; q<3; q++)
                                    Cs(i,j,k,l) += Fc(i,m)*Fc(j,n)*Fc(k,p)*Fc(l,q)*Cc(m,n,p,q);
    
    Cs *= 1.0/Jc;*/
    
    mechSpatialTangent(Cs);
    
    itensor Fcinv = Fc.inverse();

    const double alpha = thethermoFiniteStrainMaterial._thermalExpansion;
    const double bulk  = thethermoFiniteStrainMaterial.getProperty(muesli::PR_BULK);
    const double dtheta= temp_c - thethermoFiniteStrainMaterial.referenceTemperature();
    const double f     = -3.0 * alpha * bulk * dtheta;

    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++)
                {
                    Cs(i,j,k,l) += -f * ( Fcinv(l,i)*Fcinv(j,k) );
                }
}




double thermoZAMP::storedEnergy() const
{
    //tentative stored energy function
    const double lambda = thethermoFiniteStrainMaterial.lambda;
    const double mu     = thethermoFiniteStrainMaterial.mu;
    const double C2     = thethermoFiniteStrainMaterial._C2;
    const double C3     = thethermoFiniteStrainMaterial._C3;
    const double C4     = thethermoFiniteStrainMaterial._C4;
    const double C5     = thethermoFiniteStrainMaterial._C5;
    const double n      = thethermoFiniteStrainMaterial._n;
    
    // Logarithmic principal elastic stretches
    ivector lambda2, xeigvec[3];
    be_c.spectralDecomposition(xeigvec, lambda2);
    ivector epse;
    for (size_t i=0; i<3; i++)
    {
        epse[i] = 0.5*log(lambda2[i]);
    }
    
    const double We = 0.5*lambda*( epse(0) + epse(1) + epse(2) ) * ( epse(0) + epse(1) + epse(2) )
    + mu * epse.squaredNorm();
    
    double Wp = 0.0;
    if (dgamma>0.0)
    {
        Wp += (1-Dc)*(C5/(n+1)*pow(iso_c,(n+1)) + C2/1.5 * pow(iso_c,1.5) * exp(-C3*currT + C4*currT*log(epdot_c)));
    }

    return Wp+We;
}




// coupling tensor M = 2.0 * theta * d^2(psi)/( dC dTheta ) = theta * d[S]/d[Theta]
istensor thermoZAMP::symmetricCouplingTensor() const
{
    istensor C  = istensor::tensorTransposedTimesTensor(Fc);
    istensor Cinv = C.inverse();
    itensor Fcinv = Fc.inverse();

    const double alpha = thethermoFiniteStrainMaterial._thermalExpansion;
    const double bulk  = thethermoFiniteStrainMaterial.getProperty(muesli::PR_BULK);

    return -3.0/2.0 * alpha * bulk * temp_c * Cinv;
}



// D = 4 d^2[psi]/d C^2
itensor4 thermoZAMP::symmetricCouplingTensorDC() const
{
    istensor C = istensor::tensorTransposedTimesTensor(Fc);
    istensor Cinv = C.inverse();
    itensor4 D;

    const double alpha = thethermoFiniteStrainMaterial._thermalExpansion;
    const double bulk  = thethermoFiniteStrainMaterial.getProperty(muesli::PR_BULK);
    const double f     = 6.0 * alpha * bulk * temp_c;

    D.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++)
                {
                    D(i,j,k,l) += f * 0.5 * ( Cinv(i,k)* Cinv(j,l) + Cinv(i,l)* Cinv(j,k) );
                }
    return D;
}




istensor thermoZAMP::symmetricCouplingTensorDTheta() const
{
    istensor C  = istensor::tensorTransposedTimesTensor(Fc);
    istensor Cinv = C.inverse();
    
    const double alpha = thethermoFiniteStrainMaterial._thermalExpansion;
    const double bulk  = thethermoFiniteStrainMaterial.getProperty(muesli::PR_BULK);
    return -3.0/2.0 * alpha * bulk * Cinv;
}




void thermoZAMP::updateCurrentState(const double theTime, const istensor& C,
                                    const ivector& GradT, const double& temp)
{
    itensor U = istensor::squareRoot(C);
    this->updateCurrentState(theTime, U, GradT, temp);
}




void thermoZAMP::updateCurrentState(const double theTime, const itensor& F,
                                    const ivector& GradT, const double& temp)
{
    tc = theTime;
    Fc = F;
    Jc = F.determinant();
    
    // Recover material parameters
    const double mu     = thethermoFiniteStrainMaterial.mu;
    const double kappa  = thethermoFiniteStrainMaterial.bulk;
    
    // Re-assignment of current temperature for thermo-mech formulations
    currT=temp;
    if(currT<thethermoFiniteStrainMaterial._curT)
    {
        currT = thethermoFiniteStrainMaterial._curT;
    }
    
    const ivector vone(1.0, 1.0, 1.0);
    int flagD = 1;
    int flagDc = 0;
    
    // Incremental deformation gradient f = F * Fn^-1 */
    itensor f = F * Fn.inverse();
    
    
    //---------------------------------------------------------------------------------------
    //                                     trial state
    //---------------------------------------------------------------------------------------
    // Trial elastic finger tensor b_e
    istensor beTR = istensor::FSFt(f, be_n);
    
    // Logarithmic principal elastic stretches

    beTR.spectralDecomposition(nn, lambda2TR);
    ivector epseTR;
    ivector devEpseTR;
    for (size_t i=0; i<3; i++)
    {
        epseTR[i] = 0.5*log(lambda2TR[i]);
    }
    double theta = epseTR.dot(vone);
    
    // Trial state with frozen plastic flow
    if (!thethermoFiniteStrainMaterial.damageModelActivated || (thethermoFiniteStrainMaterial.damageModelActivated && Dn!=1.0))
    {
        devEpseTR = epseTR - 1.0/3.0*theta*vone;
    }
    double  isoTR     = iso_n;
    double  DTR       = Dn;
    ivector tauVol    = kappa*theta*vone;
    ivector tauTR;
    
    // Trial deviatoric principal Kirchhoff stress -- Hencky model
    if (thethermoFiniteStrainMaterial.damageModelActivated && Dn==1.0)
    {
        tauTR = tauVol;
        if(theta>0.0)
        {
            tauVol[0]=tauVol[1]=tauVol[2]=0.0;
        }
    }
    else
    {
        tauTR = tauVol + 2.0*mu*devEpseTR;
    }
    
    // Yield function value at trial state
    double phiTR;
    if(thethermoFiniteStrainMaterial.damageModelActivated)
    {
        phiTR = yieldFunctionDMG(tauTR, isoTR, 0.0, DTR);
    }
    else
    {
        phiTR = yieldFunction(tauTR, isoTR, 0.0);
    }

    // Explicit computations, deactivated in this version
    //double edot0 = theElastoplasticMaterial._edot0;
    //double edot = (epdot_n < 1e-4*edot0)  ? 1e-4*edot0 : epdot_n;
    //double phiTR = yieldFunction(tauTR, iso_n, edot);
    
    //---------------------------------------------------------------------------------------
    //                   check trial state and radial return
    //---------------------------------------------------------------------------------------
    ivector devepse_c;
        
    if (phiTR <= 0.0)
    {
        // Elastic step: trial -> n+1
        iso_c     = isoTR;
        devepse_c = devEpseTR;
        dgamma    = 0.0;
        flagD = 0;
        if(thethermoFiniteStrainMaterial.damageModelActivated)
        {
            Dc = DTR;
        }
    }
    
    else
    {
        if (!thethermoFiniteStrainMaterial.damageModelActivated || Dn!=1.0)
        {
            // Deviatoric stress
            devTauTR = 2.0*mu*devEpseTR;
        
            // Plastic step : return mapping in principal stretches space,
            // solution goes into MP variables dgamma and iso_c
            plasticReturn(devTauTR);
        
            //explicitRadialReturn(devTauTR, iso_n, edot);
        
            // Correct the trial quantities
            nubarTR = devTauTR * (1.0/devTauTR.norm());
            devepse_c = devEpseTR - dgamma * nubarTR;
        }
    }
    
    // Update elastic finger tensor
    ivector epse_c;
    if (thethermoFiniteStrainMaterial.damageModelActivated && Dn==1.0)
    {
        epse_c = 1.0/3.0*theta*vone;
    }
    else
    {
        epse_c = devepse_c + 1.0/3.0*theta*vone;
    }
    be_c.setZero();
    for (unsigned i=0; i<3; i++)
    {
        lambda2(i) = exp(2.0*epse_c[i]);
        be_c.addScaledVdyadicV( lambda2(i), nn[i] );
    }
    
    double theta_c = epse_c.dot(vone);
    
    // Update of Kirchhof stress and Damage variable, in case Damage model is in use
    // and the temporal step lies within the plastic regime deformation, considering
    // cases of Dmax = 1.0 or Dmax < 1.0, selected by user
    if(thethermoFiniteStrainMaterial.damageModelActivated && Dn!=1.0)
    {
        if (theta_c>0.0)
        {
            tau = (1.0-Dn) * (tauVol + 2.0 * mu * devepse_c);
        }
        else
        {
            tau =  tauVol + (1.0-Dn) * 2.0 * mu * devepse_c;
        }
    }
    else if(thethermoFiniteStrainMaterial.damageModelActivated && Dn==1.0)
    {
        if (theta_c>0.0)
        {
            tau[0]=tau[1]=tau[2]=0.0;
        }
        else
        {
            tau =  tauVol;
        }
    }
    else if (!thethermoFiniteStrainMaterial.damageModelActivated)
    {
        tau = tauVol + 2.0 * mu * devepse_c;
    }
    
    double dt = tc - tn;
    if(dgamma>0.0)
    {
        epdot_c = dt > 0.0  ? SQ23*dgamma/dt : epdot_n;
    }
    else
    {
        epdot_c = dt > 0.0 ? SQ23*1e-50/dt : epdot_n;
    }
    
    
    // Recalculation of current damage variable in case it is possible that it increases.
    if(thethermoFiniteStrainMaterial.damageModelActivated && flagD==1 && !fullyDamaged)
    {
        tau = tau/(1-Dn);
        
        istensor sigma;
        CauchyStress(sigma);
        
        tau = (1-Dn)*tau;
        
        double strainToFracture = thethermoFiniteStrainMaterial.theDamageModel->calculateStrainToFracture(iso_c, epdot_c, sigma, temp);
        Dc = Dn + (iso_c-iso_n)/strainToFracture;
        
        if(Dc>thethermoFiniteStrainMaterial._maxDamage)
        {
            Dc = thethermoFiniteStrainMaterial._maxDamage;
            if (isnan(Dc))
            {
                Dc=Dn;
            }
            
            fullyDamaged = true;
            
            if (theta_c>0.0)
            {
                tau = (1.0-Dc) * (tauVol + 2.0 * mu * devepse_c);
            }
            else
            {
                tau =  tauVol + (1.0-Dc) * 2.0 * mu * devepse_c;
            }
            
            flagDc = 1;
        }
    }
        
    // Testing of yield function to test dgamma root solution with and without Damage model, DEACTIVATED by default
    /*if(thethermoFiniteStrainMaterial.damageModelActivated)
    {
        if (flagDc == 0)
        {
            double kk = yieldFunctionDMG(tau, iso_c, epdot_c, Dn);
            if (fabs(kk)>0.1 && dgamma>0.0)
            {
                std::cout<<"\n"<<kk;
                std::cout<<"\n"<<dgamma;
                std::cout<<"\n"<<theTime<<"\n";
            }
        }
        if (flagDc == 1)
        {
            double kk = yieldFunctionDMG(tau, iso_c, epdot_c, Dc);
            if (fabs(kk)>0.1 && dgamma>0.0)
            {
                std::cout<<"\n"<<kk;
                std::cout<<"\n"<<dgamma;
                std::cout<<"\n"<<theTime<<"\n";
            }
        }
    }
    else
    {
        double kk = yieldFunction(tau, iso_c, epdot_c);
        if (fabs(kk)>0.1 && dgamma>0.0)
        {
            std::cout<<"\n"<<kk;
            std::cout<<"\n"<<dgamma;
            std::cout<<"\n"<<theTime<<"\n";
        }
    }*/
    
    time_c  = theTime;
    
    //condition to assure that the update of the variables applies, when it is a thermo-mechanical material simulation
    if (thethermoFiniteStrainMaterial._curT!=temp)
    {
        GradT_c = GradT;
        temp_c  = temp;
    }
}




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




// Yield function in principal Kirchhoff space
double thermoZAMP::yieldFunction(const ivector& tau,
                                 const double&  eps,
                                 const double&  epsdot) const
{
    const double C0 = thethermoFiniteStrainMaterial._C0;
    const double C1 = thethermoFiniteStrainMaterial._C1;
    const double C2 = thethermoFiniteStrainMaterial._C2;
    const double C3 = thethermoFiniteStrainMaterial._C3;
    const double C4 = thethermoFiniteStrainMaterial._C4;
    const double C5 = thethermoFiniteStrainMaterial._C5;
    const double Tc = currT;
    const double n  = thethermoFiniteStrainMaterial._n;
    
    const ivector vone(1.0, 1.0, 1.0);
    double thermal = (epsdot>0.0) ? (C1+C2*sqrt(eps))*exp(-C3*Tc+C4*Tc*log(epsdot)) : 0.0;
    
    if(C5!=0.0)
    {
        double strain  = C5*pow(eps,n);
        double sigma_y = SQ23 * (thermal + strain + C0);
        ivector taubar = tau - (1.0/3.0) * tau.dot(vone)*vone;
        
        return  taubar.norm() - sigma_y;
    }
    else
    {
        double sigma_y = SQ23 * (thermal + C0);
        ivector taubar = tau - (1.0/3.0) * tau.dot(vone)*vone;
        
        return  taubar.norm() - sigma_y;
    }
}




// Yield function in principal Kirchhoff space, when Damage model is in use
double thermoZAMP::yieldFunctionDMG(const ivector& tau,
                                    const double&  eps,
                                    const double&  epsdot,
                                    const double& D) const
{
    const double C0 = thethermoFiniteStrainMaterial._C0;
    const double C1 = thethermoFiniteStrainMaterial._C1;
    const double C2 = thethermoFiniteStrainMaterial._C2;
    const double C3 = thethermoFiniteStrainMaterial._C3;
    const double C4 = thethermoFiniteStrainMaterial._C4;
    const double C5 = thethermoFiniteStrainMaterial._C5;
    const double Tc = currT;
    const double n  = thethermoFiniteStrainMaterial._n;
    
    const ivector vone(1.0, 1.0, 1.0);
    double thermal = (epsdot>0.0) ? (C1+C2*sqrt(eps))*exp(-C3*Tc+C4*Tc*log(epsdot)) : 0.0;

    if(C5!=0.0)
    {
        double strain  = C5*pow(eps,n);
        double sigma_y = (1-D)* SQ23 * (thermal + strain + C0);
        ivector taubar = tau - (1.0/3.0) * tau.dot(vone)*vone;
        
        return  taubar.norm() - sigma_y;
    }
    else
    {
        double sigma_y = (1-D)* SQ23 * (thermal + C0);
        ivector taubar = tau - (1.0/3.0) * tau.dot(vone)*vone;
        
        return  taubar.norm() - sigma_y;
    }
}
