/****************************************************************************
*
*                                 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 "fplastic.h"
#include <string.h>
#include <cmath>

#define J2TOL     1e-10
#define GAMMAITER 15

using namespace std;
using namespace muesli;


fplasticMaterial::fplasticMaterial(const std::string& name,
                     const materialProperties& cl)
:
  finiteStrainMaterial(name, cl),
  E(0.0), nu(0.0), lambda(0.0), mu(0.0), bulk(0.0), cp(0.0), cs(0.0),
  Y0(0.0), Yinf(-1.0), Yexp(0.0), soft(0.0), Hiso(0.0), Hkine(0.0)
{
    muesli::assignValue(cl, "young",       E);
    muesli::assignValue(cl, "poisson",     nu);
    muesli::assignValue(cl, "lambda",      lambda);
    muesli::assignValue(cl, "mu",          mu);
    muesli::assignValue(cl, "isotropich",  Hiso);
    muesli::assignValue(cl, "kinematich",  Hkine);
    muesli::assignValue(cl, "yieldstress", Y0);
    muesli::assignValue(cl, "yieldinf",    Yinf);
    muesli::assignValue(cl, "hardexp",     Yexp);
    muesli::assignValue(cl, "softening",   soft);

    if (Yinf == -1.0) Yinf = Y0;

    // E and nu have priority. If they are defined, define lambda and mu
    if (E*E > 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 (density() > 0.0)
    {
        cp = sqrt((lambda+2.0*mu)/density());
        cs = sqrt(mu/density());
    }
}




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




bool fplasticMaterial::check() const
{
    if (mu > 0 && lambda+2.0*mu > 0) return true;
    else return false;
}




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




// this function is much faster than the one with string property names, because
// avoids string comparisons. It should be used.
double fplasticMaterial::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 fplasticMaterial::print(std::ostream &of) const
{
    of  << "\n Elastoplastic material for finite strain kinematics."
        << "\n Von Mises yield criterion with linear kinematic hardening"
        << "\n and saturation isotropic hardening."
        << "\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
        << "\n   Yield stress           : " << Y0
        << "\n   Yield softening        : " << soft
        << "\n   Isotropic hardening    : " << Hiso
        << "\n   Kinematic hardening    : " << Hkine
        << "\n   Saturation yield stress: " << Yinf
        << "\n   Hardening exponent     : " << Yexp
        << "\n   Reference temperature  : " << referenceTemperature();
}




void fplasticMaterial::setRandom()
{
    E      = muesli::randomUniform(1.0, 10.0);
    nu     = muesli::randomUniform(0.05, 0.45);
    setDensity(muesli::randomUniform(1.0, 100.0));
    Y0     = 0.001*E*muesli::randomUniform(0.5, 1.5);
    soft   = randomUniform(0.0, 1.0);
    Yinf   = Y0*muesli::randomUniform(1.5, 2.5);
    Yexp   = muesli::randomUniform(0.9, 1.1);
    Hiso   = E*1e-4 * muesli::randomUniform(1.0, 2.0);
    Hkine  = E*1e-4 * muesli::randomUniform(1.0, 2.0);

    lambda = E*nu/(1.0-2.0*nu)/(1.0+nu);
    mu     = E/2.0/(1.0+nu);
    cp     = sqrt((lambda+2.0*mu)/density());
    cs     = sqrt(2.0*mu/density());
    bulk   = lambda + 2.0/3.0 * mu;
}




bool fplasticMaterial::test(std::ostream &of)
{
    bool isok=true;
    setRandom();

    muesli::finiteStrainMP* p = this->createMaterialPoint();
    p->setRandom();

    isok = p->testImplementation(of);
    delete p;
    return isok;
}




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




fplasticMP::fplasticMP(const fplasticMaterial &m)
:
    finiteStrainMP(m),
    theElastoplasticMaterial(m),
    iso_n(0.0), iso_c(0.0),
    dgamma(0.0),
    tgs(0.0), tgdelta(0.0)
{
    be_n = be_c = istensor::identity();
    kine_n.setZero();
    kine_c.setZero();
    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);
}




void fplasticMP::CauchyStress(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 fplasticMP::commitCurrentState()
{
    finiteStrainMP::commitCurrentState();

    be_n   = be_c;
    kine_n = kine_c;
    iso_n  = iso_c;
}




void fplasticMP::contractWithConvectedTangent(const ivector &v1, const ivector& v2, itensor &T) const
{
    const itensor G = Fc.inverse();
    itensor TT;
    contractWithSpatialTangent(G.transpose()*v1, G.transpose()*v2, TT);

    T = Jc * G * TT * G.transpose();
}




void fplasticMP::contractWithSpatialTangent(const ivector &v1, const ivector &v2, itensor &T) const
{
    // recover material parameters
    const double mu     = theElastoplasticMaterial.mu;
    const double kappa  = theElastoplasticMaterial.bulk;
    const ivector vone(1.0, 1.0, 1.0);

    istensor cep;
    if (dgamma <= 0.0)
    {
        cep = 2.0*mu*istensor::identity();
        cep.addScaledVdyadicV(kappa-2.0/3.0*mu, vone);
    }

    else
    {
        cep = 2.0*mu*tgs*istensor::identity();
        cep.addScaledVdyadicV(kappa-2.0/3.0*mu*tgs, vone);
        cep.addScaledVdyadicV(-2.0*mu*tgdelta, nubarTR);
    }

    T.setZero();
    // Part i: tangent in ppal strains
    for (unsigned a=0; a<3; a++)
    {
        for (unsigned b=0; b<3; b++)
        {
            const double f = cep(a,b) * nn[a].dot(v1) * nn[b].dot(v2);
            T.addDyadic( f *            nn[a],          nn[b]);
        }
    }


    // Part ii: g
    double gab = 0.0;
    for (unsigned a=0; a<3; a++)
    {
        const double na_v1 = nn[a].dot(v1);
        const double na_v2 = nn[a].dot(v2);

        const double fa = -2.0 * tau(a) * na_v1 * na_v2;
        T.addDyadic(fa*nn[a], nn[a]);

        for (unsigned b=0; b<3; b++)
        {
            if (a != b)
            {
                const double nb_v1 = nn[b].dot(v1);
                const double nb_v2 = nn[b].dot(v2);

                if ( fabs(lambda2TR[a] - lambda2TR[b]) > 1e-14 )
                {
                    gab = (tau[a]*lambda2TR[b] - tau[b]*lambda2TR[a])/(lambda2TR[a] - lambda2TR[b]);
                }
                else
                {
                    gab = 0.5*cep(a,a) - tau(a) - 0.5*cep(b,a);
                }

                const double h = 0.5;

                // symmetrized tangent
                const double f1 = h * gab * nb_v1 * nb_v2;
                const double f2 = h * gab * nb_v1 * na_v2;

                T.addDyadic( nn[a],  f1 * nn[a] + f2 * nn[b]);


                const double f3 = h * gab * na_v1 * nb_v2;
                const double f4 = h * gab * na_v1 * na_v2;

                T.addDyadic( nn[b],  f3 * nn[a] + f4 * nn[b]);
            }
        }
    }

    const double iJ = 1.0/Jc;
    T *= iJ;

    return;
}




void fplasticMP::convectedTangent(itensor4& C) const
{
    itensor4 Cs;
    spatialTangent(Cs);

    C.setZero();
    itensor Finv = Fc.inverse();
    for (size_t i=0; i < 3; i++)
        for (size_t j=0; j < 3; j++)
            for (size_t k=0; k < 3; k++)
                for (size_t l=0; l < 3; l++)
                    for (size_t m=0; m < 3; m++)
                        for (size_t n=0; n < 3; n++)
                            for (size_t p=0; p < 3; p++)
                                for (size_t q=0; q < 3; q++)
                                    C(i,j,k,l) +=  Jc * Finv(i,m)*Finv(j,n)*Finv(k,p)*Finv(l,q)*Cs(m,n,p,q);
}




double fplasticMP::getDamage() const 
{
    const double soft   = theElastoplasticMaterial.soft;
    const double temp0  = theElastoplasticMaterial.referenceTemperature();
    const double Y00    = theElastoplasticMaterial.Y0;
    const double Y0     = Y00*exp(-soft*(tempc-temp0));
    const double Hiso0  = theElastoplasticMaterial.Hiso;
    const double Hiso   = Hiso0*exp(-soft*(tempc-temp0));

    return Y0 + Hiso*iso_c;
}




void fplasticMP::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);
                }
}




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




double fplasticMP::dissipatedEnergy() const
{
    double dt = tc - tn;
    return dt * kineticPotential();
}




double fplasticMP::effectiveStoredEnergy() const
{
    return storedEnergy() + dissipatedEnergy();
}




void fplasticMP::firstPiolaKirchhoffStress(itensor &P) const
{
    // reconstruct Kirchhoff stress
    istensor kir;
    kir.setZero();
    for (size_t i=0; i<3; i++)
    {
        kir.addScaledVdyadicV(tau[i], nn[i]);
    }

    itensor G = Fc.inverse();
    P = kir*G.transpose();
}




double fplasticMP::kineticPotential() const
{
    const ivector vone(1.0, 1.0, 1.0);
    const double temp0 = theElastoplasticMaterial.referenceTemperature();
    const double mu    = theElastoplasticMaterial.mu;
    const double Hkine = theElastoplasticMaterial.Hkine;
    const double Y00   = theElastoplasticMaterial.Y0;
    const double soft  = theElastoplasticMaterial.soft;
    const double Y0    = Y00*exp(-soft*(tempc-temp0));
    const double Hiso0 = theElastoplasticMaterial.Hiso;
    const double Hiso  = Hiso0*exp(-soft*(tempc-temp0));
    const double Yinf  = theElastoplasticMaterial.Yinf;
    const double Yexp  = theElastoplasticMaterial.Yexp;


    // logarithmic principal elastic stretches
    ivector  neigvec[3], lambda2n;;
    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;
    ivector Q = -2.0/3.0 * Hkine * kine_c;
    double  q = -Hiso*iso_c - (Yinf - Y0)*(1.0-exp(-Yexp*iso_c));
    ivector dt_dp = dgamma*nubarTR;
    double  dt_Psistar = dt_dp.dot(s) + (kine_c - kine_n).dot(Q) + (iso_c - iso_n)*q;

    double dt = tc - tn;

    return (dt == 0.0) ? 0.0 : dt_Psistar/dt;
}




double fplasticMP::plasticSlip() const
{
    return iso_c;
}




// beta: deviatoric, shifted principal Kirchhoff stresses
// xi  : isotropic hardening strain-like internal variable
double fplasticMP::radialReturn(const ivector& beta, const double& xi) const
{
    const double sq23  = sqrt(2.0/3.0);
    const double temp0 = theElastoplasticMaterial.referenceTemperature();
    const double mu    = theElastoplasticMaterial.mu;
    const double Hkine = theElastoplasticMaterial.Hkine;
    const double Y00   = theElastoplasticMaterial.Y0;
    const double soft  = theElastoplasticMaterial.soft;
    const double Y0    = Y00*exp(-soft*(tempc-temp0));
    const double Hiso0 = theElastoplasticMaterial.Hiso;
    const double Hiso  = Hiso0*exp(-soft*(tempc-temp0));
    const double Yinf  = theElastoplasticMaterial.Yinf;
    const double Yexp  = theElastoplasticMaterial.Yexp;
    const double bnorm = beta.norm();

    double dg   = 0.0;
    double G    = 1.0;
    double iso  = xi;
    size_t count = 0;

    double Kp = 0.0, Kpp = 0.0, DG = 0.0;
    const double tol = J2TOL;

    while (G > (J2TOL*Y0) && ++count < GAMMAITER)
    {
        Kp  = Hiso*iso + (Yinf - Y0)*(1.0-exp(-Yexp*iso));
        Kpp = Hiso+(Yinf-Y0)*Yexp*exp(-Yexp*iso);
        G   = bnorm - sq23*(Y0 + Kp)*(1.0+tol) - dg*(2.0*mu+2.0/3.0*Hkine);
        DG  = -2.0*mu -2.0/3.0*(Hkine + Kpp);
        dg -= G/DG;
        iso = xi + sq23 * dg;
    }

    if (count == GAMMAITER)
    {
	  std::cout << "\n zero radial return";
	  dg = 0.0;
    }

    return dg;
}




void fplasticMP::resetCurrentState()
{
    finiteStrainMP::resetCurrentState();

    iso_c  = iso_n;
    kine_c = kine_n;
    be_c   = be_n;
}




void fplasticMP::secondPiolaKirchhoffStress(istensor &S) const
{
    // reconstruct Kirchhoff stress
    istensor kirchhoff;
    kirchhoff.setZero();
    for (size_t i=0; i<3; i++)
    {
        kirchhoff.addScaledVdyadicV(tau[i], nn[i]);
    }

    S = istensor::FtCF((Fc.inverse()).transpose(), kirchhoff);
}





void fplasticMP::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;
    kine_n = kine;
    be_n   = be;
}



void fplasticMP::setRandom()
{
    iso_c = iso_n = muesli::randomUniform(1.0, 2.0);

    ivector tmp; tmp.setRandom();
    ivector vone(1.0, 1.0, 1.0);
    kine_c = kine_n = tmp - 1.0/3.0* tmp.dot(vone) * vone;

    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();
}




void fplasticMP::spatialTangent(itensor4& Cs) const
{
    Cs.setZero();

    // recover material parameters
    const double mu     = theElastoplasticMaterial.mu;
    const double kappa  = theElastoplasticMaterial.bulk;
    const ivector vone(1.0, 1.0, 1.0);

    istensor cep;
    if (dgamma <= 0.0)
    {
        cep = 2.0*mu*istensor::identity();
        cep.addScaledVdyadicV(kappa-2.0/3.0*mu, vone);
    }

    else
    {
        cep = 2.0*mu*tgs*istensor::identity();
        cep.addScaledVdyadicV(kappa-2.0/3.0*mu*tgs, vone);
        cep.addScaledVdyadicV(-2.0*mu*tgdelta, nubarTR);
    }

    const double iJ = 1.0/Jc;

    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 a=0; a<3; a++)
                        for(unsigned b=0; b<3; b++)
                        {
                            Cs(i,j,k,l) +=  iJ * cep(a,b) * nn[a][i]*nn[a][j]*nn[b][k]*nn[b][l];
                        }

                    double gab = 0.0;
                    for (unsigned a=0; a<3; a++)
                    {

                        Cs(i,j,k,l) += -2.0 * iJ * tau(a) * nn[a][i]*nn[a][j]*nn[a][k]*nn[a][l];

                        for (unsigned b=0; b<3; b++)
                        {
                            if (a != b)
                            {
                                if ( fabs(lambda2TR[a] - lambda2TR[b]) > 1e-14 )
                                {
                                    gab = (tau[a]*lambda2TR[b] - tau[b]*lambda2TR[a])/(lambda2TR[a] - lambda2TR[b]);
                                }
                                else
                                {
                                    gab = 0.5*cep(a,a) - tau(a) - 0.5*cep(b,a);
                                }

                                Cs(i,j,k,l) += 0.5 * iJ * gab * nn[a][i]*nn[b][j]*nn[a][k]*nn[b][l] + 0.5 * iJ * gab * nn[a][i]*nn[b][j]*nn[b][k]*nn[a][l];
                                Cs(i,j,k,l) += 0.5 * iJ * gab * nn[b][i]*nn[a][j]*nn[a][k]*nn[b][l] + 0.5 * iJ * gab * nn[b][i]*nn[a][j]*nn[b][k]*nn[a][l];
                            }
                        }
                    }
                }
}




double fplasticMP::storedEnergy() const
{
    const double lambda= theElastoplasticMaterial.lambda;
    const double mu    = theElastoplasticMaterial.mu;
    const double Hkine = theElastoplasticMaterial.Hkine;
    const double soft  = theElastoplasticMaterial.soft;
    const double temp0 = theElastoplasticMaterial.referenceTemperature();
    const double Y00   = theElastoplasticMaterial.Y0;
    const double Y0    = Y00*exp(-soft*(tempc-temp0));
    const double Hiso0 = theElastoplasticMaterial.Hiso;
    const double Hiso  = Hiso0*exp(-soft*(tempc-temp0));
    const double Yinf  = theElastoplasticMaterial.Yinf;
    const double Yexp  = theElastoplasticMaterial.Yexp;

    ivector epse;
    for (size_t i=0; i<3; i++)
    {
        epse[i] = 0.5*log(lambda2[i]);
    }

    double trace = epse(0) + epse(1) + epse(2);
    double We = 0.5*lambda* trace*trace + mu * epse.squaredNorm();
    double Wp = 0.5 * Hkine * kine_c.squaredNorm() + 0.5 * Hiso * iso_c*iso_c;

    if (Yexp != 0.0)
    {
        Wp += (Yinf-Y0) * (iso_c + exp(-Yexp*iso_c)/ Yexp - 1.0/Yexp) ;
    }

    return Wp+We;
}




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




void fplasticMP::updateCurrentState(const double theTime, const itensor& F)
{
    finiteStrainMP::updateCurrentState(theTime, F);

    // recover material parameters
    const double mu     = theElastoplasticMaterial.mu;
    const double kappa  = theElastoplasticMaterial.bulk;
    const double soft   = theElastoplasticMaterial.soft;
    const double temp0  = theElastoplasticMaterial.referenceTemperature();
    const double Y00    = theElastoplasticMaterial.Y0;
    const double Y0     = Y00*exp(-soft*(tempc-temp0));
    const double Hkin   = theElastoplasticMaterial.Hkine;
    const double Hiso0  = theElastoplasticMaterial.Hiso;
    const double Hiso   = Hiso0*exp(-soft*(tempc-temp0));
    const double Yinf   = theElastoplasticMaterial.Yinf;
    const double Yexp   = theElastoplasticMaterial.Yexp;
    const ivector vone(1.0, 1.0, 1.0);


    // Incremental deformation gradient f = F * Fn^-1 */
    const itensor f = F * Fn.inverse();


    //---------------------------------------------------------------------------------------
    //                                     trial state
    //---------------------------------------------------------------------------------------
    // trial elastic finger tensor b_e
    const istensor beTR = istensor::FSFt(f, be_n);

    // logarithmic principal elastic stretches
    beTR.spectralDecomposition(nn, lambda2TR);
    ivector epseTR;
    for (size_t i=0; i<3; i++)
    {
        epseTR[i] = 0.5*log(lambda2TR[i]);
    }
    const double theta = epseTR.dot(vone);

    // trial state with frozen plastic flow
    const ivector devepseTR = epseTR - 1.0/3.0*theta*vone;
    const double  isoTR     = iso_n;
    const ivector kineTR    = kine_n;

    // trial stress-like hardening variables
    const double  Kp  =  Hiso*isoTR + (Yinf - Y0)*(1.0-exp(-Yexp*isoTR));
    const double  qTR = -Kp;
    const ivector QTR = -2.0/3.0*Hkin*kineTR;

    // trial deviatoric principal Kirchhoff strees -- Hencky model
    const ivector tauTR = kappa*theta * vone + 2.0*mu* devepseTR;


    //---------------------------------------------------------------------------------------
    //                   check trial state and radial return
    //---------------------------------------------------------------------------------------
    ivector devepse_c;

    double phiTR = yieldFunction(tauTR, QTR, qTR);
    if (phiTR <= 0.0)
    {
        // elastic step: trial -> n+1
        kine_c    = kineTR;
        iso_c     = isoTR;
        devepse_c = devepseTR;
        dgamma    = 0.0;

        // factors for tangent
        tgs       = 1.0;
        tgdelta   = 1.0;
    }

    else
    {
        // deviatoric shifted stress
        const ivector devbetaTR = 2.0*mu*devepseTR - QTR;
        nubarTR = devbetaTR * (1.0/devbetaTR.norm());

        // plastic step : return mapping in principal stretches space
        dgamma = radialReturn(devbetaTR, isoTR);

        // correct the trial quantities
        devepse_c = devepseTR - dgamma * nubarTR;
        kine_c    = kineTR    - dgamma * nubarTR;
        iso_c     = isoTR     + sqrt(2.0/3.0) * dgamma;

        // factors for tangent
        double Kpp = Hiso + (Yinf-Y0)*Yexp*exp(-Yexp*iso_c);
        tgs     = 1.0 - 2.0*mu*dgamma/devbetaTR.norm();
        tgdelta = 1.0/(1.0 + (Kpp+Hkin)/(3.0*mu)) - (1.0-tgs);
    }


    tau = kappa * theta * vone + 2.0 * mu * devepse_c;

    // update elastic finger tensor
    ivector 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 fplasticMP::volumetricEnergy() const
{
    return 0.0;
}




// yield function in principal Kirchhoff space
//
//      phi(tau, Q, q) = || dev(tau) - Qbar || - sqrt(2/3) (Yo - q)
//
//      tau : principal Kirchhoff stresses
//      Q   : kinematic hardening stress-like internal variable
//      q   : isotropic hardening stress-like internal variable

double fplasticMP::yieldFunction(const ivector& tau,
                                 const ivector& Q,
                                 const double&  q) const
{
    const double Y00   = theElastoplasticMaterial.Y0;
    const double soft  = theElastoplasticMaterial.soft;
    const double temp0 = theElastoplasticMaterial.referenceTemperature();
    const double Y0    = Y00*exp(-soft*(tempc-temp0));

    const ivector vone(1.0, 1.0, 1.0);
    const ivector taubar = tau - (1.0/3.0) * tau.dot(vone)*vone;

    const double phi = (taubar - Q).norm() - sqrt(2.0/3.0)*(Y0 - q);
    return phi;
}
