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

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

using namespace std;
using namespace muesli;




arrheniusTypeMaterial::arrheniusTypeMaterial(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),
                            _a0 (0.0),_a1(0.0),_a2(0.0),_a3(0.0),_q0(0.0),
                            _q1(0.0), _q2(0.0),_q3(0.0),_al0(0.0),_al1(0.0),
                            _al2(0.0),_al3(0.0),_n(0.0),_curT(1.0)
{
    muesli::assignValue(cl, "young",       E);
    muesli::assignValue(cl, "poisson",     nu);
    muesli::assignValue(cl, "lambda",      lambda);
    muesli::assignValue(cl, "mu",          mu);
    muesli::assignValue(cl, "a0",           _a0);
    muesli::assignValue(cl, "a1",           _a1);
    muesli::assignValue(cl, "a2",           _a2);
    muesli::assignValue(cl, "a3",           _a3);
    muesli::assignValue(cl, "q0",           _q0);
    muesli::assignValue(cl, "q1",           _q1);
    muesli::assignValue(cl, "q2",           _q2);
    muesli::assignValue(cl, "q3",           _q3);
    muesli::assignValue(cl, "alfa0",        _al0);
    muesli::assignValue(cl, "alfa1",        _al1);
    muesli::assignValue(cl, "alfa2",        _al2);
    muesli::assignValue(cl, "alfa3",        _al3);
    muesli::assignValue(cl, "n",             _n);
    muesli::assignValue(cl, "temp",       _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;

    double rho = density();
    if (rho > 0.0)
    {
        cp = sqrt((lambda+2.0*mu)/rho);
        cs = sqrt(mu/rho);
    }
}




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




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




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




// This function is much faster than the one with string property names, because
// avoids string comparisons. It should be used.
double arrheniusTypeMaterial::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 arrheniusTypeMaterial::print(std::ostream &of) const
{
    of  << "\n Arrhenius rate- and temperature-dependent plasticity"
    << "\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   The yield Kirchhoff stress is of the form:"
    << "\n   |tau| - sqrt(2/3) 1/alfa(asinh(Z/A)^1/n)"
    << "\n   alfa = alfa0+alfa1 eps + alfa2 eps^2 + alfa3 eps^3"
    << "\n   Q = Q0 + Q1 eps + Q2 eps^2 + Q3 eps^3"
    << "\n  lnA = A0 + A1 eps + A2 eps^2 + A3 eps^3";
    
    of  << "\n   A0                  : " << _a0
    << "\n   A1                      : " << _a1
    << "\n   A2                      : " << _a2
    << "\n   A3                      : " << _a3
    << "\n   Q0                      : " << _q1
    << "\n   Q1                      : " << _q1
    << "\n   Q2                      : " << _q2
    << "\n   Q3                      : " << _q3
    << "\n   alfa0                   : " << _al0
    << "\n   alfa1                   : " << _al1
    << "\n   alfa2                   : " << _al2
    << "\n   alfa3                   : " << _al3
    << "\n   temp                    : " << _curT
    << "\n   n                       : " << _n;
}




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

    E      = muesli::randomUniform(1.0, 10.0);
    nu     = muesli::randomUniform(0.05, 0.45);
    _a0    = muesli::randomUniform(1.0, 10.0);
    _a1    = muesli::randomUniform(1.0, 10.0);
    _a2    = muesli::randomUniform(1.0, 10.0);
    _a3    = muesli::randomUniform(1.0, 10.0);
    _q0    = muesli::randomUniform(1.0, 10.0);
    _q1    = muesli::randomUniform(1.0, 10.0);
    _q2    = muesli::randomUniform(1.0, 10.0);
    _q3    = muesli::randomUniform(1.0, 10.0);
    _al0   = muesli::randomUniform(1.0, 10.0);
    _al1   = muesli::randomUniform(1.0, 10.0);
    _al2   = muesli::randomUniform(1.0, 10.0);
    _al3   = muesli::randomUniform(1.0, 10.0);
    _curT  = muesli::randomUniform(400.0, 500.0);
    _n     = muesli::randomUniform(60.0, 80.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 arrheniusTypeMaterial::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 arrheniusTypeMaterial::waveVelocity() const
{
    return cp;
}




arrheniusTypeMP::arrheniusTypeMP(const arrheniusTypeMaterial &m)
:
finiteStrainMP(m),
theElastoplasticMaterial(m),
iso_n(0.0), iso_c(0.0),
epdot_n(0.0), epdot_c(0.0),
dgamma(0.0)
{
    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);
}




// Brent method function
double arrheniusTypeMP::brentroot(double a, double b, double Ga,
                                    double Gb, double eqpn, double ntbar,
                                    double mu, double a0, double a1,
                                    double a2,double a3, double q0,
                                    double q1, double q2, double q3,
                                    double al0, double al1, double al2,
                                    double al3, 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, a0, a1, a2, a3, q0, q1, q2, q3, al0, al1, al2, al3, 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 arrheniusTypeMP::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 arrheniusTypeMP::spatialTangent(itensor4& C) const
{
    C.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) C(i,j,k,l) += theElastoplasticMaterial.lambda;
                    if (i==k && j==l) C(i,j,k,l) += theElastoplasticMaterial.mu;
                    if (i==l && j==k) C(i,j,k,l) += theElastoplasticMaterial.mu;
                }
}




void arrheniusTypeMP::commitCurrentState()
{
    finiteStrainMP::commitCurrentState();
    
    be_n    = be_c;
    iso_n   = iso_c;
    epdot_n = epdot_c;
}




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




// Calls generic FiniteStrain numerical convected tangent
void arrheniusTypeMP::convectedTangent(itensor4& C) const
{
    this->numericalConvectedTangent(C);
}




void arrheniusTypeMP::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 arrheniusTypeMP::deviatoricEnergy() const
{
    return 0.0;
}




double arrheniusTypeMP::dissipatedEnergy() const
{
    const ivector vone(1.0, 1.0, 1.0);
    ivector taubar = tau - (1.0/3.0) * tau.dot(vone)*vone;
    
    double dissEnergy;
    if(tc-tn>0.0)
    {
        dissEnergy = dgamma*taubar.norm()/(tc-tn);
    }
    else
    {
        dissEnergy = 0.0;
    }
    
    return dissEnergy;
}




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




// Explicit radial return algorithm not tested within this code
// and not used in the "update state" function
void arrheniusTypeMP::explicitRadialReturn(const ivector &taudev, double ep, double epdot)
{
    const double mu  = theElastoplasticMaterial.mu;
    const double a0  = theElastoplasticMaterial._a0;
    const double a1  = theElastoplasticMaterial._a1;
    const double a2  = theElastoplasticMaterial._a2;
    const double a3  = theElastoplasticMaterial._a3;
    const double q0  = theElastoplasticMaterial._q0;
    const double q1  = theElastoplasticMaterial._q1;
    const double q2  = theElastoplasticMaterial._q2;
    const double q3  = theElastoplasticMaterial._q3;
    const double al0 = theElastoplasticMaterial._al0;
    const double al1 = theElastoplasticMaterial._al1;
    const double al2 = theElastoplasticMaterial._al2;
    const double al3 = theElastoplasticMaterial._al3;
    const double Tc  = theElastoplasticMaterial._curT;
    const double n   = theElastoplasticMaterial._n;
    
    
    double Q = q0+q1*ep+q2*ep*ep+q3*ep*ep*ep;
    double Z (0.0);
    double expterm = exp(Q/(R*Tc));
    if(epdot>0.0)
    {
        Z = epdot*expterm;
    }
    double alfa = al0+al1*ep+al2*ep*ep+al3*ep*ep*ep;
    double A = exp(a0+a1*ep+a2*ep*ep+a3*ep*ep*ep);
    double sigma = 1.0/alfa*asinh(pow(Z/A,1.0/n));
    
    double dalfa = SQ23*al1+SQ23*2.0*al2*ep+SQ23*3.0*al3*ep*ep;
    double dZ = SQ23*expterm+epdot*expterm*((SQ23*q1+SQ23*2.0*q2*ep+SQ23*3.0*q3*ep*ep)/(R*Tc));
    double dA = exp(a0+a1*ep+a2*ep*ep+a3*ep*ep*ep)*(SQ23*a1+SQ23*2.0*a2*ep+SQ23*a3*3.0*ep*ep);
    double dsigma = (-dalfa/(alfa*alfa))*(asinh(pow((Z/A),(1.0/n))))+((1.0/alfa)*(1.0/(sqrt(pow(Z/A,2)+1.0)))*((1.0/n)*pow((Z/A),((1.0-n)/n))*((dZ*A-dA*Z)/(A*A))));

    
    double sigma_y = SQ23 * sigma;
    dgamma = (taudev.norm() - sigma_y)/(2.0*mu);
    iso_c = ep + SQ23*dgamma*dsigma;
}




double arrheniusTypeMP::kineticPotential() const
{
    return 0.0;
}




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




// 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 ARR equation), the function jumps to a Brent scheme, whose convergence is assured.
void arrheniusTypeMP::plasticReturn(const ivector& taubarTR)
{
    const double mu      = theElastoplasticMaterial.mu;
    const double a0      = theElastoplasticMaterial._a0;
    const double a1      = theElastoplasticMaterial._a1;
    const double a2      = theElastoplasticMaterial._a2;
    const double a3      = theElastoplasticMaterial._a3;
    const double q0      = theElastoplasticMaterial._q0;
    const double q1      = theElastoplasticMaterial._q1;
    const double q2      = theElastoplasticMaterial._q2;
    const double q3      = theElastoplasticMaterial._q3;
    const double al0     = theElastoplasticMaterial._al0;
    const double al1     = theElastoplasticMaterial._al1;
    const double al2     = theElastoplasticMaterial._al2;
    const double al3     = theElastoplasticMaterial._al3;
    const double Tc      = theElastoplasticMaterial._curT;
    const double n       = theElastoplasticMaterial._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, a0, a1, a2, a3, q0, q1, q2, q3, al0, al1, al2, al3, 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, a0, a1, a2, a3, q0, q1, q2, q3, al0, al1, al2, al3, Tc, n, eqpn, ntbar, dt, x, DG);
        x -= G/DG;
        if (std::isnan(x)) break;
        plasticReturnResidual(mu, a0, a1, a2, a3, q0, q1, q2, q3, al0, al1, al2, al3, 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(x) || 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, a0, a1, a2, a3, q0, q1, q2, q3, al0, al1, al2, al3, Tc, n, eqpn, ntbar, dt, a, Ga);
        plasticReturnResidual(mu, a0, a1, a2, a3, q0, q1, q2, q3, al0, al1, al2, al3, Tc, n, eqpn, ntbar, dt, b, Gb);
        
        
        while(Gb<0.0)
        {
            b=b*100.0;
            plasticReturnResidual(mu, a0, a1, a2, a3, q0, q1, q2, q3, al0, al1, al2, al3, 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 Arrhenius equation features
        // and because of this the division loop could be conditionally stopped. The error induced
        // is thought to be assumable
        while(Ga>0.0)
        {
            a=a/1000.0;
            plasticReturnResidual(mu, a0, a1, a2, a3, q0, q1, q2, q3, al0, al1, al2, al3, 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 big enough
        dgamma = brentroot(a, b, Ga, Gb, eqpn, ntbar, mu, a0, a1, a2, a3, q0, q1, q2, q3, al0, al1, al2, al3, Tc, n, dt);
        iso_c  = iso_n + SQ23*dgamma;
        }
    }
}




// Calculation of G function value given each iteration value of dgamma
void arrheniusTypeMP::plasticReturnResidual(double mu, double a0, double a1, double a2,
                                              double a3, double q0, double q1, double q2,
                                              double q3, double al0, double al1, double al2,
                                              double al3, double Tc, double n, double eqpn,
                                              double tau, double dt, double dgamma, double& G)
{
    double deqp = SQ23*dgamma;
    double eqp  = eqpn + deqp;
    
    double alfa = al0+al1*eqp+al2*eqp*eqp+al3*eqp*eqp*eqp;
    double Q = q0+q1*eqp+q2*eqp*eqp+q3*eqp*eqp*eqp;
    double A = exp(a0+a1*eqp+a2*eqp*eqp+a3*eqp*eqp*eqp);
    double expterm = exp(Q/(R*Tc));
    double Z (0.0);
    
    if(deqp>0.0)
    {
        Z = deqp/dt*expterm;
    }
    
    double sigma = 1.0/alfa*asinh(pow(Z/A,1.0/n));
    
    G = 2.0*mu*dgamma - tau + SQ23*sigma;
}




// Calculation of G function derivative value given each iteration value of dgamma
void  arrheniusTypeMP::plasticReturnTangent(double mu, double a0, double a1, double a2,
                                              double a3, double q0, double q1, double q2,
                                              double q3, double al0, double al1, double al2,
                                              double al3, double Tc,   double n, double eqpn,
                                              double tau, double dt, double dgamma, double& DG)
{
    double deqp = SQ23*dgamma;
    double eqp  = eqpn + deqp;
    
    double Q = q0+q1*eqp+q2*eqp*eqp+q3*eqp*eqp*eqp;;
    double expterm = exp(Q/(R*Tc));
    double Z (0.0);
    if(deqp>0.0)
    {
        Z = deqp/dt*expterm;
    }
    double alfa = al0+al1*eqp+al2*eqp*eqp+al3*eqp*eqp*eqp;
    double A = exp(a0+a1*eqp+a2*eqp*eqp+a3*eqp*eqp*eqp);
    
    double dalfa = SQ23*al1+SQ23*2.0*al2*eqp+SQ23*3.0*al3*eqp*eqp;
    double dZ = SQ23/dt*expterm+deqp/dt*expterm*((SQ23*q1+SQ23*2.0*q2*eqp+SQ23*3.0*q3*eqp*eqp)/(R*Tc));
    double dA = exp(a0+a1*eqp+a2*eqp*eqp+a3*eqp*eqp*eqp)*(SQ23*a1+SQ23*2.0*a2*eqp+SQ23*a3*3.0*eqp*eqp);
    double dsigma = (-dalfa/(alfa*alfa))*(asinh(pow(Z/A,1.0/n)))+((1.0/alfa)*(1.0/(sqrt(((Z/A)*(Z/A))+1.0)))*((1.0/n)*pow((Z/A),((1.0-n)/n))*((dZ*A-dA*Z)/(A*A))));
    DG = 2.0*mu + SQ23*dsigma;
}




void arrheniusTypeMP::resetCurrentState()
{
    finiteStrainMP::resetCurrentState();
    
    iso_c   = iso_n;
    be_c    = be_n;
    epdot_c = epdot_n;
}




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




double arrheniusTypeMP::storedEnergy() const
{
    const double lambda = theElastoplasticMaterial.lambda;
    const double mu     = theElastoplasticMaterial.mu;
    
    
    // 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 = dissipatedEnergy()*0.1/0.9;;
    return Wp+We;
}




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




void arrheniusTypeMP::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 ivector vone(1.0, 1.0, 1.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;
    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
    ivector devEpseTR = epseTR - 1.0/3.0*theta*vone;
    double isoTR     = iso_n;
    
    // Trial deviatoric principal Kirchhoff strees -- Hencky model
    const ivector tauTR = kappa*theta * vone + 2.0*mu* devEpseTR;
    
    // Yield function value at trial state
    double 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;
    }
    else
    {
        // Deviatoric stress
       const ivector devTauTR = 2.0*mu*devEpseTR;
        
        // Plastic step : return mapping in principal stretches space.
        // Solution goes into MP variables dgamma and iso_c
        plasticReturn(devTauTR);
     
        // Correct the trial quantities
        if(devTauTR.norm()==0.0)
        {
            nubarTR = devTauTR*0.0;
            devepse_c = devEpseTR;
        }
        else
        {
            nubarTR = devTauTR * (1.0/devTauTR.norm());
            devepse_c = devEpseTR - dgamma * nubarTR;
        }
    }
    
    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 dt = tc - tn;
    epdot_c = dt > 0.0 ? SQ23*dgamma/dt : epdot_n;
    
    // Testing of yield function to test dgamma root solution, deactivated by default
    //double kk = yieldFunction(tau, iso_c, epdot_c);
    //if (fabs(kk)>0.1 && dgamma>0.0)
    //{
    //   std::cout<<"\n"<<kk;
    //   std::cout<<"\n"<<dgamma;
    //}
}




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




// Yield function in principal Kirchhoff space
double arrheniusTypeMP::yieldFunction(const ivector& tau,
                                        const double&  eps,
                                        const double&  epsdot) const
{
    const double a0   = theElastoplasticMaterial._a0;
    const double a1   = theElastoplasticMaterial._a1;
    const double a2   = theElastoplasticMaterial._a2;
    const double a3   = theElastoplasticMaterial._a3;
    const double q0   = theElastoplasticMaterial._q0;
    const double q1   = theElastoplasticMaterial._q1;
    const double q2   = theElastoplasticMaterial._q2;
    const double q3   = theElastoplasticMaterial._q3;
    const double al0  = theElastoplasticMaterial._al0;
    const double al1  = theElastoplasticMaterial._al1;
    const double al2  = theElastoplasticMaterial._al2;
    const double al3  = theElastoplasticMaterial._al3;
    const double Tc   = theElastoplasticMaterial._curT;
    const double n    = theElastoplasticMaterial._n;
    
    const ivector vone(1.0, 1.0, 1.0);
    double alfa = al0+al1*eps+al2*eps*eps+al3*eps*eps*eps;
    double Q = q0+q1*eps+q2*eps*eps+q3*eps*eps*eps;
    double A = exp(a0+a1*eps+a2*eps*eps+a3*eps*eps*eps);
    double Z = epsdot*exp(Q/(R*Tc));
    
    double sigma = 1.0/alfa*asinh(pow(Z/A,1.0/n));
    double sigma_y = SQ23 * sigma;
    ivector taubar = tau - (1.0/3.0) * tau.dot(vone)*vone;
    return  taubar.norm() - sigma_y;
}
