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

/*
 * viscoplastic.cpp,
 * infinitesimal kinematics visco-elastoplastic element
 * d. del pozo, i. romero, feb 2015
 *
 *
 *  1) von Mises Model with linear isotropic and kinematic hardenings
 *
 *     See "Numerical analysis and simulation of plasticity", by J.C. Simo,
 *     Handbook of numerical analysis vol VI
 *     table 18.2 and pages close to it
 *
 *     We = kappa/2 (theta)^2 + mu * ||e_elastic||^2
 *     Wp = Hiso/2 xi^2 + 2/3 Hkine/2 ||Xi||^2
 *     f(sigma, q, Q) = || dev(sigma) - Q || - sqrt(2/3) (Y0 - q)
 *  2) Tresca with linear isotropic hardening
 *
 *     We = kappa/2 (theta)^2 + mu * ||e_elastic||^2
 *     Wp = Hiso/2 xi^2
 *     f(sigma, q, Q) = sigma_I - sigma_III - (Y0 - q)
 */


#define PROJTOL     1e-8     // residual tolerance for return-to-yield-surface solutions
#define PROJMAXITER 25      // max num of allowed interations in return-to-yield solutions
#define ETOL        1e-13     // residual tolerance for return-to-yield-surface solutions
#define ETOL1       1e-10
#define EPSILON     1e-5


using namespace std;
using namespace muesli;





viscoplasticMaterial::viscoplasticMaterial(const std::string& name,
                                             const materialProperties& cl)
:
smallStrainMaterial(name, cl),
plasticityType("mises"),
E(0.0), nu(0.0), bulk(0.0), cp(0.0), cs(0.0), lambda(0.0), mu(0.0),
Hiso(0.0), Hkine(0.0), Y0(0.0),eta(0.0), vexponent(1.0), epsrateo(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, "isotropich", Hiso);
    muesli::assignValue(cl, "kinematich", Hkine);
    muesli::assignValue(cl, "yieldstress", Y0);
    muesli::assignValue(cl, "eta", eta);
    muesli::assignValue(cl, "viscousexponent", vexponent);
    muesli::assignValue(cl, "srateref", epsrateo);
    
    std::string mstr;
    muesli::assignValue(cl, "model", mstr);
    
    if      (mstr == "von_mises") plasticityType = "mises";
    else if (mstr == "tresca")    plasticityType = "tresca";
    
    // 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 = lambda / 2.0 / (lambda+mu);
        E  = mu*2.0*(1.0+nu);
    }
    
    // we set all the constants, so that later on all of them can be recovered fast
    bulk = lambda + 2.0/3.0 * mu;
    
    if (rho > 0.0)
    {
        cp = sqrt((lambda+2.0*mu)/rho);
        cs = sqrt(mu/rho);
    }
}




viscoplasticMaterial::viscoplasticMaterial(const std::string& name,
                                             const double xE, const double xnu,
                                             const double xrho, const double xHiso,
                                             const double xHkine, const double xY0,
                                             const std::string xplst, const double xeta)
:
smallStrainMaterial(name),
plasticityType(xplst),
E(xE), nu(xnu), bulk(0.0), cp(0.0), cs(0.0),
lambda(0.0), mu(0.0),
rho(xrho), Hiso(xHiso), Hkine (xHkine), Y0(xY0), eta(xeta),
vexponent(0.0), epsrateo(0.0)
{
    lambda = E*nu/(1.0-2.0*nu)/(1.0+nu);
    mu     = E/2.0/(1.0+nu);
    bulk   = lambda + 2.0/3.0 * mu;
    
    if (rho > 0.0)
    {
        cp   = sqrt((lambda+2.0*mu)/rho);
        cs   = sqrt(2.0*mu/rho);
    }
}




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




smallStrainMP* viscoplasticMaterial::createMaterialPoint() const
{
    smallStrainMP* mp = new viscoplasticMP(*this);
    return mp;
}




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




// this function is much faster than the one with string property names, because
// avoids string comparisons. It should be used.
double viscoplasticMaterial::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_YIELD:      ret = Y0;        break;
        case PR_ISOHARD:    ret = Hiso;      break;
        case PR_KINHARD:    ret = Hkine;     break;
        case PR_ETA    :    ret = eta;       break;
        case PR_SRREF:      ret = epsrateo;  break;
        case PR_VEXPONENT:  ret = vexponent; break;
            
            
        default:
            std::cout << "Error in elastoplasticMaterial. Property not defined";
    }
    return ret;
}




// this is a potential law used for the viscoplastic response
// g(x) = eps0 * x^a
//
// derivative =  0, the function g
// derivative =  1, the function g'
// derivative = -1, the function G = integral g
//
double viscoplasticMaterial::potentialLaw(const double &x, const int &derivative) const
{
    double ret = 0.0;
    const  double a    = this->vexponent;
    
    if      (derivative == -1)  ret = 1.0/(a+1.0)*pow(x, a+1.0);//integral
    else if (derivative ==  0)  ret = pow(x, a);
    else if (derivative ==  1)  ret = a*pow(x, a-1.0);//derivative
    else cout << "\n Error in potentialLaw";
    
    
    return ret;
}



void viscoplasticMaterial::print(std::ostream &of) const
{
    
    if (plasticityType == "mises" )
    {
        of  << "\n Viscoelastoplastic material for small strain kinematics."
        << "\n   von Mises yield criterion with linear isotropic and kinematic hardenings."
        << "\n   Yield stress           : " << Y0
        << "\n   Isotropic hardening    : " << Hiso
        << "\n   Kinematic hardening    : " << Hkine
        << "\n   Viscosity              : " << eta;
    }
    
    if (plasticityType == "tresca" )
    {
        of  << "\n Viscoelastoplastic material for small strain kinematics."
        << "\n   Tresca yield criterion with linear isotropic hardening."
        << "\n   Yield stress           : " << Y0
        << "\n   Isotropic hardening    : " << Hiso
        << "\n   reference value for stress rate     : " << epsrateo
        << "\n   power exponent for viscoplastic law : " << vexponent;
    }
    
    
    
    of  << "\n   Young modulus:  E      : " << E;
    of  << "\n   Poisson ratio:  nu     : " << nu;
    of  << "\n   Lame constants: lambda : " << lambda;
    of  << "\n                   mu     : " << mu;
    of  << "\n   Bulk modulus:   k      : " << bulk;
    of  << "\n   Density                : " << rho;
    
    if (rho > 0.0)
    {
        of  << "\n   Wave velocities c_p    : " << cp;
        of  << "\n                   c_s    : " << cs;
    }
}




void viscoplasticMaterial::setRandom()
{

    E   = muesli::randomUniform(1.0, 10000.0);
    nu  = muesli::randomUniform(0.05, 0.45);
    rho = muesli::randomUniform(1.0, 100.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)/rho);
    cs     = sqrt(2.0*mu/rho);
    bulk   = lambda + 2.0/3.0*mu;
    
    int pt = muesli::discreteUniform(0, 1);
    
    if (pt == 0)
    {
        plasticityType = "mises";
        Hiso  = E*1e-4 * muesli::randomUniform(1.0, 2.0);
        Hkine = E*1e-4 * muesli::randomUniform(1.0, 2.0);
        eta   = E*1e-5 * muesli::randomUniform(0.1, 1.0);
        Y0    = E*1e-3 * muesli::randomUniform(0.5, 1.5);
    }
    else if (pt == 1)
    {
        plasticityType = "tresca";
        Y0        = E * 1e-3 *muesli::randomUniform(0.5, 1.5);
        Hiso      = E * 1e-2 * muesli::randomUniform(1.0, 2.0);
        Hkine     = 0.0;
        vexponent = muesli::randomUniform(1.0,2.0);
        epsrateo  = muesli::randomUniform(0.1,100.0);
    }
    
    
}




bool viscoplasticMaterial::test(std::ostream  &of)
{
    bool isok=true;
    setRandom();
    
    smallStrainMP* p = this->createMaterialPoint();
    p->setRandom();

    if      (this->plasticityType == "mises")   of << "\n   von Mises type";
    else if (this->plasticityType == "tresca") of << "\n   Tresca type";
    
    isok = p->testImplementation(of);
    delete p;
    
    return isok;
}




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




viscoplasticMP::viscoplasticMP(const viscoplasticMaterial &m) :
muesli::smallStrainMP(m),
theViscoplasticMaterial(m),
dg_n(0.0), xi_n(0.0),
dg_c(0.0), xi_c(0.0)
{
    ep_n.setZero();
    ep_c.setZero();
    Xi_n.setZero();
    Xi_c.setZero();
}




void viscoplasticMP::commitCurrentState()
{

    smallStrainMP::commitCurrentState();
    dg_n     = dg_c;
    ep_n     = ep_c;
    xi_n     = xi_c;
    Xi_n     = Xi_c;
}




void viscoplasticMP::contractWithDeviatoricTangent(const ivector &v1, const ivector &v2, itensor &T) const
{
    if (theViscoplasticMaterial.plasticityType == "mises")
    {
        const double mu  = theViscoplasticMaterial.mu;
        const double eta = theViscoplasticMaterial.eta;
        const double tau = eta/(2.0*mu);
        const double dt  = time_c - time_n;
        
        if (dg_c > 0.0)
        {
            const double   Hiso     = theViscoplasticMaterial.Hiso;
            const double   Hkine    = theViscoplasticMaterial.Hkine;
            const istensor strial   = 2.0*mu*(istensor::deviatoricPart(eps_c) - ep_n);
            const istensor Qtrial   = -2.0/3.0*Hkine*Xi_n;
            
            const istensor backs    = strial - Qtrial;
            const double   betanorm = backs.norm();
            const istensor n        = backs * (1.0/betanorm);
            const double   zeta     = 1.0 - 2.0*mu*dg_c/betanorm;
            const double   zbar     = 2.0*mu/(1.0+tau/dt) / (2.0*mu+2.0/3.0*(Hiso+Hkine)) - (1.0-zeta);
            
            T = - 2.0/3.0*mu*zeta   * itensor::dyadic(v1, v2)
            +   mu*zeta*v1.dot(v2)  * itensor::identity()
            +   mu*zeta             * itensor::dyadic(v2,v1)
            -   2.0*mu*zbar         * itensor::dyadic(n*v1,n*v2);
        }
        else
        {
            T = - 2.0/3.0*mu        * itensor::dyadic(v1, v2)
            +   mu*v1.dot(v2)       * itensor::identity()
            +   mu                  * itensor::dyadic(v2,v1);
        }
    }
    else if (theViscoplasticMaterial.plasticityType == "tresca")
    {
        smallStrainMP::contractWithDeviatoricTangent(v1,v2,T);
    }
}




void viscoplasticMP::contractWithTangent(const ivector &v1, const ivector &v2, itensor &T) const
{
    const double mu  = theViscoplasticMaterial.mu;
    const double k   = theViscoplasticMaterial.bulk;
    const double eta = theViscoplasticMaterial.eta;
    const double Hkine = theViscoplasticMaterial.Hkine;
    const double Hiso = theViscoplasticMaterial.Hiso;
    const double tau = eta/(2.0*mu);
    const double dt  = time_c - time_n;
    
    if (theViscoplasticMaterial.plasticityType == "mises")
    {
        if (dg_c == 0.0)
        {
            T = (k - 2.0/3.0*mu) * itensor::dyadic(v1, v2)
            + mu*v1.dot(v2)      * itensor::identity()
            + mu                 * itensor::dyadic(v2,v1);
        }
        else
        {
            const istensor strial = 2.0*mu*(istensor::deviatoricPart(eps_c) - ep_n);
            const istensor Qtrial = -2.0/3.0*Hkine*Xi_n;
            const istensor backs  = strial - Qtrial;
            const istensor n      = backs * (1.0/backs.norm());
            const double zeta     = 1.0 - 2.0*mu*dg_c/backs.norm();
            const double zbar     = 2.0*mu/(1.0+tau/dt) / (2.0*mu+2.0/3.0*(Hiso+Hkine)) - (1.0-zeta);
            
            T = (k - 2.0/3.0*mu*zeta)  * itensor::dyadic(v1, v2)
            + (mu*zeta*v1.dot(v2)) * itensor::identity()
            + (mu*zeta)            * itensor::dyadic(v2,v1)
            - (2.0*mu*zbar)        * itensor::dyadic(n*v1,n*v2);
        }
    }
    else if (theViscoplasticMaterial.plasticityType == "tresca")
    {
        smallStrainMP::contractWithTangent(v1,v2,T);
    }
    
}




double viscoplasticMP::deviatoricEnergy() const
{
    istensor eps_e = istensor::deviatoricPart(eps_c - ep_c);
    return theViscoplasticMaterial.mu * eps_e.contract(eps_e);
}




void viscoplasticMP::deviatoricStress(istensor& s) const
{
    s = 2.0*theViscoplasticMaterial.mu*(istensor::deviatoricPart(eps_c)-ep_c);
}




double viscoplasticMP::dissipatedEnergy() const
{
    const double k     = theViscoplasticMaterial.bulk;
    const double mu    = theViscoplasticMaterial.mu;
    const double Hiso  = theViscoplasticMaterial.Hiso;
    const double Hkine = theViscoplasticMaterial.Hkine;
    const double Y0    = theViscoplasticMaterial.Y0;
    double ret(0.0);
    
    if (theViscoplasticMaterial.plasticityType=="mises")
    {
        const istensor ee_c = istensor::deviatoricPart(eps_c)-ep_c;
        const istensor s    = 2.0*mu*ee_c;
        const istensor Q    = -2.0/3.0 * Hkine*Xi_c;
        const double   q    = -Hiso*xi_c;
        
        const double dt_Psistar = (ep_c-ep_n).dot(s) + (Xi_c - Xi_n).dot(Q) + (xi_c-xi_n)*q;
        ret = dt_Psistar;
    }
    
    else if (theViscoplasticMaterial.plasticityType=="tresca")
    {
        const istensor e_cedev   = istensor::deviatoricPart(eps_c)-ep_c;  // trial elastic deviatoric strain tensor
        const istensor stressdev = 2.0 * mu * e_cedev;                     // deviatoric trial stress tensor
        const double   q         = -Hiso*xi_c;
        const istensor sigma     = stressdev+k*eps_c.trace()*istensor::identity();
        double dt_Psistar        = (ep_c-ep_n).dot(stressdev) + (xi_c-xi_n)*q;
        
        const double dt    = time_c-time_n;
        const double epso  = theViscoplasticMaterial.epsrateo;
        
        double d = trescaDistance(theViscoplasticMaterial, sigma, q);
        double G = theViscoplasticMaterial.potentialLaw(d/Y0, -1);
        double Psi = epso*G;
        dt_Psistar -= dt*Psi;
        
        ret = dt_Psistar;
    }
    
    return ret;
}




// necessary to compute in smallthermo material the dissipated energy
istensor viscoplasticMP::getConvergedPlasticStrain() const
{
    return ep_n;
}




materialState viscoplasticMP::getConvergedState() const
{
    materialState state_n = smallStrainMP::getConvergedState();
    
    state_n.theDouble.push_back(dg_n);
    state_n.theDouble.push_back(xi_n);
    state_n.theStensor.push_back(ep_n);
    state_n.theStensor.push_back(Xi_n);
    
    return state_n;
}




// necessary to compute in smallthermo material the dissipated energy
istensor viscoplasticMP::getCurrentPlasticStrain() const
{
    return ep_c;
}




materialState viscoplasticMP::getCurrentState() const
{
    materialState state_c = smallStrainMP::getCurrentState();
    
    state_c.theDouble.push_back(dg_c);
    state_c.theDouble.push_back(xi_c);
    state_c.theStensor.push_back(ep_c);
    state_c.theStensor.push_back(Xi_c);
    
    return state_c;
}




double viscoplasticMP::effectiveStoredEnergy() const
{
    const double mu    = theViscoplasticMaterial.mu;
    const double k     = theViscoplasticMaterial.bulk;
    const double Hiso  = theViscoplasticMaterial.Hiso;
    const double Hkine = theViscoplasticMaterial.Hkine;
    const double Y0    = theViscoplasticMaterial.Y0;
    double ret   = 0.0;

    
    if (theViscoplasticMaterial.plasticityType == "mises")
    {
        istensor ee = istensor::deviatoricPart(eps_c)-ep_c;
        istensor s  = 2.0*mu*ee;
        istensor Q  = -2.0/3.0 * Hkine*Xi_c;
        double   q  = -Hiso*xi_c;
        double   t  = eps_c.trace();
        
        double We = mu*ee.squaredNorm() + 0.5*k*t*t;
        double Wp = 0.5*Hiso*xi_c*xi_c + 0.5*2.0/3.0*Hkine*Xi_c.squaredNorm();
        
        double dt_Psistar = (ep_c-ep_n).dot(s) + (Xi_c-Xi_n).dot(Q) + (xi_c-xi_n)*q;
        
        ret = We+Wp+dt_Psistar;
    }
    
    if (theViscoplasticMaterial.plasticityType == "tresca")
    {
        istensor ee             = istensor::deviatoricPart(eps_c)-ep_c;//elastic deviatoric part
        istensor s              = 2.0*mu*ee;
        const istensor sigma    = s+k*eps_c.trace()*istensor::identity();
        
        const double   q        = -Hiso*xi_c;
        double         t        = eps_c.trace();
        double       dt_Psistar = (ep_c-ep_n).dot(s) + (xi_c-xi_n)*q;
        
        double   We       = mu*ee.squaredNorm() + 0.5*k*t*t;
        double   Wp       = 0.5*Hiso*xi_c*xi_c;
        const double dt    = time_c-time_n;
        const double epso  = theViscoplasticMaterial.epsrateo;
        
        double d = trescaDistance(theViscoplasticMaterial, sigma, q);
        double G = theViscoplasticMaterial.potentialLaw(d/Y0, -1);
        double Psi = epso*G;
        dt_Psistar -= dt*Psi;
        
        
        ret = We+Wp+dt_Psistar;
    }
    
    return ret;
}




double viscoplasticMP::kineticPotential() const
{
    const double k     = theViscoplasticMaterial.bulk;
    const double mu    = theViscoplasticMaterial.mu;
    const double Hiso  = theViscoplasticMaterial.Hiso;
    const double Hkine = theViscoplasticMaterial.Hkine;
    const double Y0    = theViscoplasticMaterial.Y0;
    double dt_Psistar  = 0.0;
    
    if (theViscoplasticMaterial.plasticityType=="mises")
    {
        const istensor ee_c = istensor::deviatoricPart(eps_c)-ep_c;
        const istensor s    = 2.0*mu*ee_c;
        const istensor Q    = -2.0/3.0 * Hkine*Xi_c;
        const double   q    = -Hiso*xi_c;
        
        dt_Psistar = (ep_c-ep_n).dot(s) + (Xi_c - Xi_n).dot(Q) + (xi_c-xi_n)*q;
    }
    
    else if (theViscoplasticMaterial.plasticityType=="tresca")
    {
        const istensor e_cedev   = istensor::deviatoricPart(eps_c)-ep_c;  // trial elastic deviatoric strain tensor
        const istensor stressdev = 2.0 * mu * e_cedev;                     // deviatoric trial stress tensor
        const double   q         = -Hiso*xi_c;
        const istensor sigma     = stressdev+k*eps_c.trace()*istensor::identity();
        dt_Psistar               = (ep_c-ep_n).dot(stressdev) + (xi_c-xi_n)*q;
        
        const double dt    = time_c-time_n;
        const double epso  = theViscoplasticMaterial.epsrateo;
        
        double d = trescaDistance(theViscoplasticMaterial, sigma, q);
        double G = theViscoplasticMaterial.potentialLaw(d/Y0, -1);
        double Psi = epso*G;
        dt_Psistar -= dt*Psi;
    }

    double dt = time_c - time_n;
    return dt == 0.0 ? 0.0 : dt_Psistar/dt;
}



double viscoplasticMP::plasticSlip() const
{
    return ep_c.norm();
}




// return to yield surface for von Mises model
void viscoplasticMP::radialReturn(const istensor& strain, double& dg)
{
    const double mu    = theViscoplasticMaterial.mu;
    const double Hiso  = theViscoplasticMaterial.Hiso;
    const double Hkine = theViscoplasticMaterial.Hkine;
    const double dt    = time_c - time_n;
    const double eta   = theViscoplasticMaterial.eta;
    const double tau   = eta/(2.0*mu);
    
    // trial state
    //const istensor e_c  = istensor::deviatoricPart(strain);       // current deviatoric strain
    const istensor ee_c   = istensor::deviatoricPart(strain-ep_n);  // current elastic deviatoric strain
    const istensor strial = 2.0 * mu * ee_c;                        // deviatoric trial stress
    
    const istensor Qtrial = -2.0/3.0*Hkine*Xi_n;
    const double   qtrial = -Hiso*xi_n;
    
    //check yield condition
    const double   ftrial = yieldfunction(theViscoplasticMaterial, strial, Qtrial, qtrial);
    
    // compute stress and update internal variables
    if (ftrial <= 1.0e-8)
    {
        // elastic step
        ep_c = ep_n;
        Xi_c = Xi_n;
        xi_c = xi_n;
        dg   = 0.0;
    }
    
    else
    {
        // plastic step: do return mapping
        
        const istensor backs = strial - Qtrial;
        const istensor n     = backs / backs.norm();
        
        // plastic slip
        dg = ftrial / ((1.0+tau/dt)*(2.0*mu + 2.0/3.0*(Hiso+Hkine)));
        
        // correct internal variables
        ep_c  = ep_n + dg*n;
        Xi_c  = Xi_n - dg*n;
        xi_c  = xi_n + sqrt(2.0/3.0)*dg;
    }
}




void viscoplasticMP::resetCurrentState()
{

    smallStrainMP::resetCurrentState();
    dg_c   = dg_n;
    ep_c   = ep_n;
    Xi_c   = Xi_n;
    xi_c   = xi_n;
}




void viscoplasticMP::setConvergedState(const double theTime, const double dgn, const istensor & epn,
                                         const double xin, const istensor & Xin, const istensor & strainn)
{

    time_n    = theTime;
    dg_n      = dgn;
    ep_n      = epn;
    xi_n      = xin;
    eps_n     = strainn;
    Xi_n      = Xin;
    
}




void viscoplasticMP::setRandom()
{

    smallStrainMP::setRandom();
    istensor tmp;
    tmp.setRandom();
    
    if (theViscoplasticMaterial.plasticityType == "mises")
    {
        ep_n = istensor::deviatoricPart(tmp);
        ep_c = ep_n;
        
        Xi_n = -ep_n;
        Xi_c = Xi_n;
        
        xi_n = muesli::randomUniform(1.0, 2.0);
        xi_c = xi_n;
    }
    
    else if (theViscoplasticMaterial.plasticityType == "tresca")
    {
        ep_n = istensor::deviatoricPart(tmp);
        ep_c = ep_n;
        
        xi_n = muesli::randomUniform(1.0, 2.0);
        xi_c = xi_n;
    }
    
    
    else
    {
        cout << "\n Error initializating random viscoplastic material" << flush;
    }
}




double viscoplasticMP::storedEnergy() const
{
    const double mu    = theViscoplasticMaterial.mu;
    const double Hiso  = theViscoplasticMaterial.Hiso;
    const double Hkine = theViscoplasticMaterial.Hkine;
    const double k     = theViscoplasticMaterial.bulk;
    
    double th = eps_c.trace() - ep_c.trace();
    double We = mu*(istensor::deviatoricPart(eps_c-ep_c)).squaredNorm() + 0.5 * k * th*th;
    double Wp = 0.5*Hiso*xi_c*xi_c + 0.5*2.0/3.0*Hkine*Xi_c.squaredNorm();
    
    return  Wp+We;
}




void viscoplasticMP::stress(istensor& sigma) const
{
    if (theViscoplasticMaterial.plasticityType=="mises")
    {
        deviatoricStress(sigma);
        sigma += theViscoplasticMaterial.bulk*eps_c.trace()*istensor::identity();
    }
    else if (theViscoplasticMaterial.plasticityType=="tresca")
    {
        deviatoricStress(sigma);
        sigma += theViscoplasticMaterial.bulk*eps_c.trace()*istensor::identity();
    }
}




void viscoplasticMP::tangentTensor(itensor4& C) const
{
    C.setZero();
    
    //Elastic tangent tensor
    itensor4 Ce;
    Ce.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) Ce(i,j,k,l) += theViscoplasticMaterial.lambda;
                    if (i==k && j==l) Ce(i,j,k,l) += theViscoplasticMaterial.mu;
                    if (i==l && j==k) Ce(i,j,k,l) += theViscoplasticMaterial.mu;
                }
    
    if (theViscoplasticMaterial.plasticityType=="mises")
    {
        const double   mu     = theViscoplasticMaterial.mu;
        const double   kb     = theViscoplasticMaterial.bulk;
        const double   Hkine  = theViscoplasticMaterial.Hkine;
        const double   Hiso   = theViscoplasticMaterial.Hiso;
        const istensor strial = 2.0*mu*(istensor::deviatoricPart(eps_c) - ep_n);
        const istensor Qtrial = -2.0/3.0*Hkine*Xi_n;
        const istensor backs  = strial - Qtrial;
        const istensor n      = backs * (1.0/backs.norm());
        const double   zeta   = 1.0 - 2.0*mu*dg_c/backs.norm();
        const double   zbar   = 2.0*mu / (2.0*mu+2.0/3.0*(Hiso+Hkine)) - (1.0-zeta);
        
        
        if (dg_c == 0.0)
        {
            C=Ce;
        }
        else
        {
            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) += kb-2.0*mu*zeta/3.0;
                            if (i==k && j==l) C(i,j,k,l) += mu*zeta;
                            if (i==l && j==k) C(i,j,k,l) += mu*zeta;
                            C(i,j,k,l) += -2.0*mu*zbar*n(i,j)*n(k,l);
                        }
        }
    }
    else if (theViscoplasticMaterial.plasticityType=="tresca")
    {
        const istensor Q_trial;
        const double   mu                = theViscoplasticMaterial.mu;
        const double   Hiso              = theViscoplasticMaterial.Hiso;
        const double   kb                = theViscoplasticMaterial.bulk;
        const double   qn                = -Hiso*xi_n;
        istensor       edtrial           = istensor::deviatoricPart(eps_c)-ep_n;//deviatoric elastic trial strain
        const istensor s_trial           = 2.0 * mu * edtrial;
        const double   te_trial          = eps_c.trace();
        const istensor sigma_trial       = s_trial + kb * te_trial * istensor::identity();
        
        double f_ = yieldfunction(theViscoplasticMaterial, sigma_trial, Q_trial, qn);
        if (f_ <= PROJTOL*theViscoplasticMaterial.Y0)
        {
            C = Ce;
        }
        else
        {
            //Dd/De; (de) definition
            int      rettypea;
            const double   delta             = theViscoplasticMaterial.Y0;
            const double   Y0                = theViscoplasticMaterial.Y0;
            const double   s0                = theViscoplasticMaterial.epsrateo;
            const double   dt                = time_c - time_n;
            istensor       e_trial           = eps_c - ep_n;
            
            int flag1=0;//flag to know the order of duplicated eigenvalues
            
            const double   distSigmatrial    = trescaDistance(theViscoplasticMaterial, sigma_trial, qn);
            
            istensor ee              = istensor::deviatoricPart(eps_c)-ep_c;//deviatoric elastic strain
            istensor epsev           = eps_c-ep_c-ee;//volumetric elastic strain
            istensor s               = 2.0 * mu * ee;
            istensor sigmavol        = 3.0*kb*epsev;
            istensor sigmax          = s + sigmavol;//current stress
            const double   qn1       = -Hiso*xi_c;//current conjugate force for isotropic strain-like variable
            const double   distSigma = trescaDistance(theViscoplasticMaterial, sigmax, qn1);//exists outside of the elastic set
            const double   g         = theViscoplasticMaterial.potentialLaw(distSigma/delta, 0);
            const double   gder      = theViscoplasticMaterial.potentialLaw(distSigma/delta, 1);
            
            const double   a         = s0*dt*g/distSigma+1.0;
            const double   b         = s0*dt*gder/delta;
            const double   c         = s0*dt*g/distSigma;
            const double   H         = a + b - c;
            const double   Phi       = s0*(gder*(1.0/delta)*distSigma-g)/(distSigma*distSigma);
            
            //Projection instances
            istensor sigmabar;
            double qbar;
            
            trescaProjector(theViscoplasticMaterial, sigmax, qn1, sigmabar, qbar);
            istensor sigmadif        = sigma_trial - sigmabar;
            istensor devsigmadif     = istensor::deviatoricPart(sigmadif);
            istensor p               = 1.0/(9.0*kb)*sigmadif.trace()* istensor::identity()+1.0/(2.0*mu)*devsigmadif;//Cinv times (sigma_trial-sigmabar)
            
            //Elastoplastic tangent tensor definition for Tresca (SIGMATRIAL AND SIGMABAR)
            itensor4 Cep;
            Cep.setZero();
            
            // eigen-decomposition
            ivector str, sdummy, sig, etr;
            ivector evecstr[3], evecsig[3], evecetr[3];
            ///the order in the eigenvalues vector is increasing with the index [0,1,2]->[3,2,1]
            s_trial.spectralDecomposition(evecstr, str);//eigenvectors and eigenvalues of the deviatoric trial stress tensor
            sigmabar.spectralDecomposition(evecsig, sig);//eigenvectors and eigenvalues of the Cauchy stress
            e_trial.spectralDecomposition(evecetr, etr);//eigenvectors and eigenvalues of the trial strain tensor
            
            double dsdetr[3][3];//derivative of the principal deviatoric stresses with respect to the principal elastic trial strains
            
            //Dummy main plane update
            
            double f_trial = str(2) - str(0) - (Y0+Hiso*xi_n);
            double dg = f_trial/(4.0*mu+Hiso);
            
            sdummy(2) = str(2) - 2.0*mu*dg;
            sdummy(1) = str(1);
            sdummy(0) = str(0) + 2.0*mu*dg;
            
            //Now this main plane projection to be correct must ensure that it is contained in the first sextant, which implies:
            // sdummy(2)>=sdummy(1)>=sdummy(0)
            // If not, then do either left or right corner return
            
            if ((sdummy(2) > sdummy(1) || (fabs(sdummy(2)-sdummy(1))<ETOL))  && (sdummy(1) > sdummy(0) || (fabs(sdummy(1)-sdummy(0))<ETOL)))
            {
                //Define tangent tensor with main plane return
                
                double f = 2.0*mu/(4.0*mu+Hiso);
                dsdetr[0][0] = dsdetr[2][2] = 2.0*mu*(1-f);
                dsdetr[1][0] = dsdetr[0][1] = dsdetr[2][1] = dsdetr[1][2] = 0.0;
                dsdetr[1][1] = 2.0*mu;
                dsdetr[2][0] = dsdetr[0][2] = 2.0*mu*f;
                
                rettypea=3;
            }
            else//bad projection, right and left corner respectively
            {
                if (str(0) + str(2) - 2.0*str(1) > ETOL)
                {
                    // return to right corner
                    
                    double d[2][2];
                    d[0][0]=-4.0*mu-Hiso;
                    d[0][1]=d[1][0]=-2.0*mu-Hiso;
                    d[1][1]=-4.0*mu-Hiso;
                    
                    double detd =  d[0][0]* d[1][1] -  d[0][1]* d[1][0];
                    dsdetr[0][0] = 2.0*mu*(1.0-8.0*mu*mu/detd);
                    dsdetr[1][1] = 2.0*mu*(1.0+2.0*mu*d[0][0]/detd);
                    dsdetr[2][2] = 2.0*mu*(1.0+2.0*mu*d[1][1]/detd);
                    dsdetr[1][0] = dsdetr[2][0] = 8.0*mu*mu*mu/detd;
                    dsdetr[0][1] = 4.0*mu*mu/detd*(d[0][1]-d[0][0]);
                    dsdetr[2][1] = -4.0*mu*mu/detd*(d[0][1]);
                    dsdetr[0][2] = 4.0*mu*mu/detd*(d[1][0]-d[1][1]);
                    dsdetr[1][2] = -4.0*mu*mu/detd*(d[1][0]);
                    
                    rettypea=1;
                }
                else
                {
                    // return to left corner
                    
                    double d[2][2];
                    d[0][0]=-4.0*mu-Hiso;
                    d[0][1]=d[1][0]=-2.0*mu-Hiso;
                    d[1][1]=-4.0*mu-Hiso;
                    
                    double detd =  d[0][0]* d[1][1] -  d[0][1]* d[1][0];
                    dsdetr[0][0] = 2.0*mu*(1.0+2.0*mu*d[1][1]/detd);
                    dsdetr[1][1] = 2.0*mu*(1.0+2.0*mu*d[0][0]/detd);
                    dsdetr[2][2] = 2.0*mu*(1.0-8.0*mu*mu/detd);
                    dsdetr[1][2] = dsdetr[0][2] = 8.0*mu*mu*mu/detd;
                    dsdetr[0][1] = -4.0*mu*mu/detd*(d[0][1]);
                    dsdetr[2][1] = 4.0*mu*mu/detd*(d[0][1]-d[0][0]);
                    dsdetr[2][0] = 4.0*mu*mu/detd*(d[1][0]-d[1][1]);
                    dsdetr[1][0] = -4.0*mu*mu/detd*(d[1][0]);
                    
                    rettypea=2;
                }
            }
            double dsigdetr[3][3];
            //initialization of dsigdetr (derivatives of principal stresses with respect to the principal elastic trial strains) at zero
            for (unsigned i=0; i<3; i++)
                for (unsigned j=0; j<3; j++)
                {
                    dsigdetr[i][j] = 0.0;
                }
            //derivatives of principal stresses with respect to the principal elastic trial strains
            
            for (unsigned i=0; i<3; i++){
                for (unsigned j=0; j<3; j++){
                    for (unsigned k=0; k<3; k++)
                    {
                        if (k==j) {
                            dsigdetr[i][j] += dsdetr[i][k] * (1.0 - 1.0/3.0);
                        }
                        else {
                            dsigdetr[i][j] -= dsdetr[i][k] * (1.0/3.0);
                        }
                    }
                    dsigdetr[i][j] += kb;
                }
            }
            //Rate independent tangent definition starts here
            
            if((fabs(etr(0) - etr(1)) < ETOL1) && (fabs(etr(1) - etr(2)) < ETOL1))
            {
                //Define tangent here, principal stresses equals to each other
                size_t owentgtowen[3];
                owentgtowen[0]  = 2;
                owentgtowen[1]  = 0;
                owentgtowen[2]  = 1;
                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) Cep(i,j,k,l) += dsigdetr[owentgtowen[0]][owentgtowen[1]];
                                if (i==k && j==l) Cep(i,j,k,l) += 0.5*(dsigdetr[owentgtowen[0]][owentgtowen[0]]-dsigdetr[owentgtowen[0]][owentgtowen[1]]);
                                if (i==l && j==k) Cep(i,j,k,l) += 0.5*(dsigdetr[owentgtowen[0]][owentgtowen[0]]-dsigdetr[owentgtowen[0]][owentgtowen[1]]);
                            }
            }
            else if (!(fabs(etr(0) - etr(1))< ETOL1) && !(fabs(etr(1) - etr(2))< ETOL1) && !(fabs(etr(0) - etr(2))< ETOL1))
            {
                //All principals stresses are different
                //cyclic permutations mapping [3,2,1]->[0,1,2]
                
                size_t cyclic[3][3];
                
                //Two mappings are required
                cyclic[0][0] = 0;
                cyclic[1][1] = 1;
                cyclic[2][2] = 2;
                cyclic[1][0] = cyclic[0][1] = 2;
                cyclic[2][0] = cyclic[0][2] = 1;
                cyclic[2][1] = cyclic[1][2] = 0;
                
                
                //mapping between the tangent definition owen indeces and C++ indeces
                size_t owentgtoIRIS[3];
                size_t owentgtowen[3];
                owentgtoIRIS[0] = 0;
                owentgtoIRIS[1] = 2;
                owentgtoIRIS[2] = 1;
                owentgtowen[0]  = 2;
                owentgtowen[1]  = 0;
                owentgtowen[2]  = 1;
                
                double ya,xa,xb,xc;
                istensor Ea, Eb, Ec, Ei1, Ej1;
                
                
                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++){
                                    //Terms to assemble the tangent
                                    ya=sig(cyclic[a][0]);
                                    xa=etr(cyclic[a][0]);
                                    xb=etr(cyclic[a][1]);
                                    xc=etr(cyclic[a][2]);
                                    
                                    double commont = ya/((xa-xb)*(xa-xc));
                                    double term1 = (xa-xb)+(xa-xc);
                                    Ea.setZero();
                                    Eb.setZero();
                                    Ec.setZero();
                                    Ea.addDyadic(evecetr[cyclic[a][0]],evecetr[cyclic[a][0]]);
                                    Eb.addDyadic(evecetr[cyclic[a][1]],evecetr[cyclic[a][1]]);
                                    Ec.addDyadic(evecetr[cyclic[a][2]],evecetr[cyclic[a][2]]);
                                    //Assembling the tangent
                                    if (i==k && j==l) Cep(i,j,k,l) += (commont)*(-0.5*(xb+xc));
                                    if (i==l && j==k) Cep(i,j,k,l) += (commont)*(-0.5*(xb+xc));
                                    if (i==k)         Cep(i,j,k,l) += (commont)*(0.5*e_trial(l,j));
                                    if (i==l)         Cep(i,j,k,l) += (commont)*(0.5*e_trial(k,j));
                                    if (j==l)         Cep(i,j,k,l) += (commont)*(0.5*e_trial(i,k));
                                    if (k==j)         Cep(i,j,k,l) += (commont)*(0.5*e_trial(i,l));
                                    
                                    Cep(i,j,k,l) += -commont*term1*Ea(i,j)*Ea(k,l);
                                    Cep(i,j,k,l) += -commont*(xb-xc)*(Eb(i,j)*Eb(k,l)- Ec(i,j)*Ec(k,l));
                                    
                                }
                                
                                for (unsigned i1=0; i1<3; i1++)
                                    for (unsigned j1=0; j1<3; j1++){//the index i1 and j1 belong to owentg
                                        Ei1.setZero();
                                        Ej1.setZero();
                                        Ei1.addDyadic(evecetr[owentgtoIRIS[i1]],evecetr[owentgtoIRIS[i1]]);
                                        Ej1.addDyadic(evecetr[owentgtoIRIS[j1]],evecetr[owentgtoIRIS[j1]]);
                                        Cep(i,j,k,l) += dsigdetr[owentgtowen[i1]][owentgtowen[j1]]*Ei1(i,j)*Ej1(k,l);
                                    }
                            }
            }
            else
            {
                //Two equal principal stresses
                double s1,s2,s3,s4,s5,s6,ya,xa,xc,yc;
                size_t a,b,c, a_o, b_o, c_o;
                
                
                if ((fabs(etr(0) - etr(1)) < ETOL1))//xa=x2=etr(2)_o=0, xb=x3=etr(1)_o=1, xc=x1=etr(0)_o=2; x2 neq (x3 eq x1)
                {
                    a=2;
                    b=1;
                    c=0;
                    a_o=0;
                    b_o=1;
                    c_o=2;
                    flag1=1;
                }

                //if (!(fabs(etr(0) - etr(1)) < ETOL1))//xa=x1=etr(0)_o=2, xb=x2=etr(2)_o=0, xc=x3=etr(1)_o=1; x1 neq (x2 eq x3)
                else
                {
                    a=0;
                    b=2;
                    c=1;
                    a_o=2;
                    b_o=0;
                    c_o=1;
                    flag1=2;
                }
                
                ya=sig(a);
                yc=sig(c);
                xa=etr(a);
                xc=etr(c);
                
                double xaminxc=(xa-xc);
                double xaplxc=(xa+xc);
                double xaminxcp2=xaminxc*xaminxc;
                double xaminxcp3=(xa-xc)*xaminxcp2;
                double yaminyc=(ya-yc);
                
                s1=yaminyc/xaminxcp2+1.0/xaminxc*(dsigdetr[c_o][b_o]-dsigdetr[c_o][c_o]);
                
                s2=2.0*xc*yaminyc/xaminxcp2+xaplxc/xaminxc*(dsigdetr[c_o][b_o]-dsigdetr[c_o][c_o]);
                
                s3=2.0*yaminyc/xaminxcp3+1.0/xaminxcp2*(dsigdetr[a_o][c_o]+dsigdetr[c_o][a_o]-dsigdetr[a_o][a_o]-dsigdetr[c_o][c_o]);
                
                s4=2.0*xc*yaminyc/xaminxcp3+1.0/xaminxc*(dsigdetr[a_o][c_o]-dsigdetr[c_o][b_o])+xc/xaminxcp2*(dsigdetr[a_o][c_o]+dsigdetr[c_o][a_o]-dsigdetr[a_o][a_o]-dsigdetr[c_o][c_o]);
                
                s5=2.0*xc*yaminyc/xaminxcp3+1.0/xaminxc*(dsigdetr[c_o][a_o]-dsigdetr[c_o][b_o])+xc/xaminxcp2*(dsigdetr[a_o][c_o]+dsigdetr[c_o][a_o]-dsigdetr[a_o][a_o]-dsigdetr[c_o][c_o]);
                
                s6=2.0*xc*xc*yaminyc/xaminxcp3+xa*xc/xaminxcp2*(dsigdetr[a_o][c_o]+dsigdetr[c_o][a_o])-xc*xc/xaminxcp2*(dsigdetr[a_o][a_o]+dsigdetr[c_o][c_o])-xaplxc/xaminxc*dsigdetr[c_o][b_o];
                
                
                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==k)         Cep(i,j,k,l) += s1*(0.5*e_trial(l,j));
                                if (i==l)         Cep(i,j,k,l) += s1*(0.5*e_trial(k,j));
                                if (j==l)         Cep(i,j,k,l) += s1*(0.5*e_trial(i,k));
                                if (k==j)         Cep(i,j,k,l) += s1*(0.5*e_trial(i,l));
                                if (i==k && j==l) Cep(i,j,k,l) += -s2*0.5;
                                if (i==l && j==k) Cep(i,j,k,l) += -s2*0.5;
                                if (k==l)         Cep(i,j,k,l) += s4*e_trial(i,j);
                                if (i==j)         Cep(i,j,k,l) += s5*e_trial(k,l);
                                if (i==j && k==l) Cep(i,j,k,l) += -s6;
                                Cep(i,j,k,l) += -s3*e_trial(i,j)*e_trial(k,l);
                            }
            }
            //End of rate independent tangent definition
            double   I1,I2,I3;
            //EIGENVALUES DERIVATIVES COMPUTATION (EPSILON comes into play due to the known result of the not-differentiability of repeated eigenvalues with respect to its tensor)
            ivector edval=edtrial.eigenvalues();
            if(flag1==1){
                I1 = 0.0;
                I2 = -0.5*((edval(0)-EPSILON)*(edval(0)-EPSILON)+(edval(1)+EPSILON)*(edval(1)+EPSILON)+edval(2)*edval(2));
                I3 = (edval(0)-EPSILON)*(edval(1)+EPSILON)*edval(2);
            }
            else if (flag1==2){
                I1 = 0.0;
                I2 = -0.5*((edval(0))*(edval(0))+(edval(1)-EPSILON)*(edval(1)-EPSILON)+(edval(2)+EPSILON)*(edval(2)+EPSILON));
                I3 = edval(0)*(edval(1)-EPSILON)*(edval(2)+EPSILON);
            }
            else{
                
                I1 = edtrial.trace();
                I2 = edtrial.invariant2();
                I3 = edtrial.invariant3();
            }
            
            
            const double   Q      = ( I1*I1 - 3.0*I2)/9.0;
            const double   R      = ( -2.0*I1*I1*I1+9.0*I1*I2-27.0*I3)/54.0;
            const double   sqrtQ3 = sqrt(Q*Q*Q);
            const double   theta  = acos(R/sqrtQ3);
            
            const itensor  dI1dX     = istensor::identity();
            const itensor  dI2dX     = I1*istensor::identity() - edtrial;
            const itensor  dI3dX     = 1.0/3.0*(I1*(I1*istensor::identity() - edtrial) + I2*istensor::identity()-(edtrial*edtrial).trace()*istensor::identity()-2.0*I1*edtrial+3.0*(edtrial*edtrial));
            const itensor  dsqrtQdX  = 1.0/(18.0*sqrt(Q))*(2.0*I1*dI1dX-3.0*dI2dX);
            const itensor  dRdX      = 1.0/(54.0)*(-6.0*I1*I1*dI1dX + 9.0*I1*dI2dX + 9.0*I2*dI1dX - 27.0*dI3dX);
            const itensor  dsqrtQ3dX = 1.0/18.0*3.0*Q*Q/sqrtQ3*(2.0*I1*dI1dX - 3.0*dI2dX);
            const itensor  dthetadX  = -1.0/sqrt(1.0 - (R/sqrtQ3)*(R/sqrtQ3))*(dRdX*sqrtQ3 - R*dsqrtQ3dX)*1.0/(Q*Q*Q);
            
            const itensor  dx3dX     = -2.0*dsqrtQdX*cos(theta/3.0)+2.0*sqrt(Q)*sin(theta/3.0)*dthetadX*1.0/3.0+1.0/3.0*dI1dX;
            const itensor  dx2dX     = -2.0*dsqrtQdX*cos((theta-2.0*M_PI)/3.0)+2.0*sqrt(Q)*sin((theta-2.0*M_PI)/3.0)*dthetadX*1.0/3.0+1.0/3.0*dI1dX;
            const itensor  dx1dX     = -2.0*dsqrtQdX*cos((theta+2.0*M_PI)/3.0)+2.0*sqrt(Q)*sin((theta+2.0*M_PI)/3.0)*dthetadX*1.0/3.0+1.0/3.0*dI1dX;
            const istensor sdx3dX(dx3dX(0,0), dx3dX(1,1), dx3dX(2,2), dx3dX(1,2), dx3dX(0,2), dx3dX(0,1));
            const istensor sdx2dX(dx2dX(0,0), dx2dX(1,1), dx2dX(2,2), dx2dX(1,2), dx2dX(0,2), dx2dX(0,1));
            const istensor sdx1dX(dx1dX(0,0), dx1dX(1,1), dx1dX(2,2), dx1dX(1,2), dx1dX(0,2), dx1dX(0,1));
            
            
            istensor dqbardEtr;
            if (rettypea==3)
            {//MID PLANE
                
                dqbardEtr=-Hiso*2.0*mu/(4.0*mu+Hiso)*(istensor::deviatoricPart(sdx1dX)-istensor::deviatoricPart(sdx3dX));
            }
            else if (rettypea==2)
            {//LEFT CORNER
                
                double d[2][2];
                d[0][0]=-4.0*mu-Hiso;
                d[0][1]=d[1][0]=-2.0*mu-Hiso;
                d[1][1]=-4.0*mu-Hiso;
                //Define tangent tensor with left corner return
                double detd =  d[0][0]* d[1][1] -  d[0][1]* d[1][0];
                dqbardEtr = Hiso*2.0*mu/detd*(-2.0*mu*istensor::deviatoricPart(sdx1dX)+4.0*mu*istensor::deviatoricPart(sdx3dX)-2.0*mu*istensor::deviatoricPart(sdx2dX));
            }
            else if (rettypea==1)
            {//RIGHT CORNER
                
                double d[2][2];
                d[0][0]=-4.0*mu-Hiso;
                d[0][1]=d[1][0]=-2.0*mu-Hiso;
                d[1][1]=-4.0*mu-Hiso;
                //Define tangent tensor with right corner return
                double detd =  d[0][0]* d[1][1] -  d[0][1]* d[1][0];
                dqbardEtr = Hiso*2.0*mu/detd*(-4.0*mu*istensor::deviatoricPart(sdx1dX)+2.0*mu*istensor::deviatoricPart(sdx2dX)+2.0*mu*istensor::deviatoricPart(sdx3dX));
                
            }
            const istensor de = (1.0/H)*delta/distSigmatrial*(sigmadif - Cep(p) - (qn-qbar)/Hiso*dqbardEtr);
            
            
            C=(1.0/a)*Ce + (c/a)*Cep;
            C.addDyadic(sigma_trial,-1.0/(a*a)*dt*Phi*de);
            C.addDyadic(sigmabar,1.0/(a*a)*dt*Phi*de);
        }
    }
}




thPotentials viscoplasticMP::thermodynamicPotentials() const
{
    thPotentials tp;
    return tp;
}




double viscoplasticMP::trescaDistance(const viscoplasticMaterial& m, const istensor& sigma, const double& q)
{
    istensor Qjunk;
    
    const double Y0    = m.Y0;
    if (yieldfunction(m, sigma, Qjunk, q) <= PROJTOL*Y0) return 0.0;
    
    const double mu    = m.mu;
    const double k     = m.bulk;
    const double Hiso  = m.Hiso;
    const double delta = m.Y0;
    
    istensor sigmabar;
    double   qbar;
    
    trescaProjector(m, sigma, q, sigmabar, qbar);
    
    istensor dsigma = sigma-sigmabar;
    double   dq = q - qbar;
    istensor ds = istensor::deviatoricPart(dsigma);
    double   dp = dsigma.trace()/3.0;
    
    // generalized distance
    double  d1;
    if (Hiso > 0.0)
        d1 = sqrt(delta * ( ds.squaredNorm()/(2.0*mu) + dp*dp/k + dq*dq/Hiso ) );
    else
        d1 = sqrt(delta * ( ds.squaredNorm()/(2.0*mu) + dp*dp/k) );
    
    return d1;
}





// this function computes the projecton of stresses (sigma,q)->(sigmabar,qbar) onto the elastic domain of Tresca
// returns the delta-gamma of the plastic step

double viscoplasticMP::trescaProjector(const viscoplasticMaterial& m, const istensor& sigma, const double& q, istensor& sigmabar, double& qbar)
{
    const double   mu     = m.mu;
    const double   Hiso   = m.Hiso;
    const double   Y0     = m.Y0;
    
    
    
    istensor Qdummy;
    
    //check yield condition
    double f_trial = yieldfunction(m, sigma, Qdummy, q);
    
    // compute stress and update internal variables
    if (f_trial <= PROJTOL*Y0)
    {
        // elastic step
        sigmabar = sigma;
        qbar = q;
        return  0.0;
    }
    
    else  // plastic step: do return mapping
    {
        double dg;
        ivector str, s;
        ivector evec[3];
        (istensor::deviatoricPart(sigma)).spectralDecomposition(evec, str);//eigenvalues and eigenvectors
        //
        while (true)//mid plane return by default
        {
            dg = f_trial/(4.0*mu+Hiso);
            
            //update 8.16
            s(2) = str(2) - 2.0*mu*dg;
            s(1) = str(1);
            s(0) = str(0) + 2.0*mu*dg;
            
            break;
        }
        
        if (!((s(2) > s(1) || (fabs(s(2)-s(1))<ETOL))  && (s(1) > s(0) || (fabs(s(1)-s(0))<ETOL))))//bad projection, right and left corner respectively
        {
            if (str(0) + str(2) - 2.0*str(1) > 0.0)
            {
                // return to right corner
                double dga(0.0), dgb(0.0);
                double sa = str(2) - str(0);
                double sb = str(2) - str(1);
                double res[2];
                
                res[0] = sa - 2.0*mu*(2.0*dga + dgb)-(Y0+Hiso*( dga + dgb)-q);
                res[1] = sb - 2.0*mu*(dga + 2.0*dgb)-(Y0+Hiso*( dga + dgb)-q);
                
                double detH = (4.0*mu+Hiso)*(4.0*mu+Hiso)-(2.0*mu+Hiso)*(2.0*mu+Hiso);
                double Hinv[2][2];
                Hinv[0][0] = Hinv[1][1] = (-4.0*mu-Hiso)/detH;
                Hinv[1][0] = Hinv[0][1] = (2.0*mu+Hiso)/detH;
                
                dga -= Hinv[0][0]*res[0] + Hinv[0][1]*res[1];
                dgb -= Hinv[1][0]*res[0] + Hinv[1][1]*res[1];
                dg = dga + dgb;
                s(2) = str(2) - 2.0*mu*dg;
                s(1) = str(1) + 2.0*mu*dgb;
                s(0) = str(0) + 2.0*mu*dga;
            }
            else
            {
                // return to left corner
                double dga(0.0), dgb(0.0);
                double sa = str(2) - str(0);
                double sb = str(1) - str(0);
                double res[2];
                
                res[0] = sa - 2.0*mu*(2.0*dga + dgb)-(Y0+Hiso*(dga + dgb)-q);
                res[1] = sb - 2.0*mu*(dga + 2.0*dgb)-(Y0+Hiso*(dga + dgb)-q);
                
                double detH = (4.0*mu+Hiso)*(4.0*mu+Hiso)-(2.0*mu+Hiso)*(2.0*mu+Hiso);
                double Hinv[2][2];
                Hinv[0][0] = Hinv[1][1] = (-4.0*mu-Hiso)/detH;
                Hinv[1][0] = Hinv[0][1] = (2.0*mu+Hiso)/detH;
                
                dga -= Hinv[0][0]*res[0] + Hinv[0][1]*res[1];
                dgb -= Hinv[1][0]*res[0] + Hinv[1][1]*res[1];
                dg = dga + dgb;
                
                s(2) = str(2) - 2.0*mu*dga;
                s(1) = str(1) - 2.0*mu*dgb;
                s(0) = str(0) + 2.0*mu*dg;
            }
        }
        ///Projections of sigma onto the elastic set
        qbar     = q - Hiso*dg;
        sigmabar = sigma-istensor::deviatoricPart(sigma);
        for (unsigned a=0; a<3; a++) sigmabar.addScaledVdyadicV(s(a), evec[a]);
        
        double test = yieldfunction(m, sigmabar, Qdummy, qbar);
        if (test/Y0 > 1e-6)
        {
            std::cout << "\n Tresca warning: projector not projecting!!"<<endl;
        }
        return dg;
    }
}

int viscoplasticMP::rettype(const viscoplasticMaterial& m, const istensor& sigma, const double& q)
{
    const double   mu     = m.mu;
    const double   Hiso   = m.Hiso;
    const double   Y0     = m.Y0;
    int ret;
    istensor Qdummy;
    
    //check yield condition
    double f_trial = yieldfunction(m, sigma, Qdummy, q);
    
    // compute stress and update internal variables
    if (f_trial <= PROJTOL*Y0)
    {
        ret=0;
    }
    
    else  // plastic step: do return mapping
    {
        double dg;
        ivector str, s;
        ivector evec[3];
        (istensor::deviatoricPart(sigma)).spectralDecomposition(evec, str);//eigenvalues and eigenvectors
        //
        while (true)//mid plane return by default
        {
            dg = f_trial/(4.0*mu+Hiso);
            
            //update 8.16
            s(2) = str(2) - 2.0*mu*dg;
            s(1) = str(1);
            s(0) = str(0) + 2.0*mu*dg;
            ret=3;
            
            //
            break;
        }
        
        if (!((s(2) > s(1) || (fabs(s(2)-s(1))<ETOL))  && (s(1) > s(0) || (fabs(s(1)-s(0))<ETOL))))//bad projection, right and left corner respectively
        {
            if (str(0) + str(2) - 2.0*str(1) > 0.0)
            {
                // return to right corner
                ret=1;
            }
            else
            {
                // return to left corner
                ret=2;
            }
        }
    }
    return ret;
}




void viscoplasticMP::updateCurrentState(const double theTime, const istensor& strain)
{

    smallStrainMP::updateCurrentState(theTime, strain);
    
    
    if (theViscoplasticMaterial.plasticityType == "mises")
    {
        radialReturn(strain, dg_c);
        
    }
    
    else if (theViscoplasticMaterial.plasticityType == "tresca")
    {
        viscousTrescaReturn(strain);
        
    }
    
    eps_c = strain;
    
}




// vicoplastic update.
// Given the current state of the point and the argument 'strain'
// this function updates the state ep_c, xi_c
void viscoplasticMP::viscousTrescaReturn(const istensor& strain)
{
    const double mu    = theViscoplasticMaterial.mu;
    const double k     = theViscoplasticMaterial.bulk;
    const double Hiso  = theViscoplasticMaterial.Hiso;
    const double Y0    = theViscoplasticMaterial.Y0;
    const double epso  = theViscoplasticMaterial.epsrateo;
    const double dt    = time_c-time_n;
    
    
    // Compute trial state
    const istensor ee_trial    = istensor::deviatoricPart(strain)-ep_n;
    const double   te_trial    = strain.trace();
    const istensor s_trial     = 2.0*mu * ee_trial;
    const istensor sigma_trial = s_trial + k * te_trial * istensor::identity();
    const istensor Q_trial;
    const double   q_trial = -Hiso*xi_n;
    
    
    //Check yield condition
    const double   f_trial = yieldfunction(theViscoplasticMaterial, sigma_trial, Q_trial, q_trial);
    
    double omega;
    istensor sigma;
    double q, dn1;
    double gn1 = 1.0;
    
    
    if (f_trial <= PROJTOL*Y0) // elastic step
    {
        ep_c = ep_n;
        xi_c = xi_n;
    }
    
    else // viscoplastic step
    {
        istensor sigma_bar;
        double q_bar;
        
        // Newton-Rapson solution for d = distance(sigma)
        
        double d_trial = trescaDistance(theViscoplasticMaterial, sigma_trial, q_trial);
        dn1 = d_trial;
        
        size_t count  = 0;
        double r      = dn1 + dt*epso * theViscoplasticMaterial.potentialLaw(dn1/Y0, 0) - d_trial;
        while (fabs(r) > PROJTOL*Y0 && ++count < PROJMAXITER)
        {
            double tg = 1.0 + dt*epso/Y0*theViscoplasticMaterial.potentialLaw(dn1/Y0, 1);
            dn1 -= r/tg;
            r  = dn1 + dt*epso * theViscoplasticMaterial.potentialLaw(dn1/Y0, 0) - d_trial;
        }
        gn1   = theViscoplasticMaterial.potentialLaw(dn1/Y0, 0);
        omega = 1.0 + dt*epso * gn1 / dn1;
        if (count == PROJMAXITER) cout << "\n   NR for distance not converged:\n" ;
        
        // stress update
        dg_c  = trescaProjector(theViscoplasticMaterial, sigma_trial, q_trial, sigma_bar, q_bar);//the current plastic slip "dg_c" refers to the rate independent solution
        sigma = sigma_bar + 1.0/omega * (sigma_trial-sigma_bar);
        q     = q_bar     + 1.0/omega * (q_trial - q_bar);
        
        // internal variable update
        istensor dsigma = sigma - sigma_bar;
        istensor ds     = istensor::deviatoricPart(dsigma);
        double   dp     = dsigma.trace()/3.0;
        double   dq     = q - q_bar;
        ep_c = ep_n + dt*epso*gn1/dn1 * ( ds/(2.0*mu) + dp/(3.0*k)*istensor::identity() );
        xi_c = (Hiso > 0.0) ? xi_n + dt*epso * gn1 / dn1 *   dq / Hiso : xi_n;
    }
}




double viscoplasticMP::volumetricEnergy() const
{
    double th = (eps_c-ep_c).trace();
    double V  = 0.5*theViscoplasticMaterial.bulk * th * th;
    return V;
}




double viscoplasticMP::yieldfunction(const viscoplasticMaterial& m, const istensor& sigma, const istensor& Q, const double& q)
{
    const istensor s = istensor::deviatoricPart(sigma);
    const double   Y0 = m.Y0;
    
    if (m.plasticityType == "mises")
    {
        istensor backs = s - Q;
        return backs.norm() - sqrt(2.0/3.0)*(Y0 - q);
    }
    
    else if (m.plasticityType == "tresca")
    {
        
        ivector eval = s.eigenvalues();
        
        return eval(2) - eval(0) - (Y0 - q);
    }
    
    else
    {
        return 0.0;
    }
}
