/****************************************************************************
*
*                                 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 "splastic.h"
#include <iomanip>
#include <cassert>

#define PROJTOL     1e-10   // 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


/*
 * splastic.cpp,
 * infinitesimal kinematics 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, p.257 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) Drucker-Praguer Model
 *     linear isotropic hardening
 *
 *     See "Computational methods for plasticity", by S. Neto,
 *
 *     The yield function depends on :
 *        - I: trace of sigma
 *        - J: square root of J2
 *        - q: conjugate stress of the internal isotropic variable
 *
 *        f(I,J,q) = J + alphac*I - (alphac+sqrt(3)/3) * (Y0 - q);
 *        q = -K'(xi),    K(xi) = Hiso/2 xi^2
 *
 *
 *  3) Tresca with linear isotropic hardening
 *
 *     See "Computational methods for plasticity", by S. Neto,
 *     p.266 and close to it
 *     p.286 and close to it
 *     Appendix A
 *
 *
 *     We = kappa/2 (theta)^2 + mu * ||e_elastic||^2
 *     Wp = Hiso/2 xi^2
 *     f(sigma, q, Q) = sigma_I - sigma_III - (Y0 - q)
 *
 *
 *  4) 1D Model
 *
 *        See "Computational Inelasticity", Box 1.5, p.45
 *
 *        Q = Hkine * Xi(0,0);
 *        q = Hiso  * xi(0,0);
 *        f = sigma - Q - (Y0 - q);
 *
 */


using namespace std;
using namespace muesli;


splasticMaterial::splasticMaterial(const std::string& name,
                                     const double xE,     const double xnu,
                                     const double xrho,   const double xHiso,
                                     const double xHkine, const double xY0,
                                     const double xalphac,
                                     const char*  plasType)
:
smallStrainMaterial(name),
E(xE), nu(xnu), bulk(0.0),
cp(0.0), cs(0.0), lambda(0.0), mu(0.0),
Hiso(xHiso), Hkine(xHkine), Y0(xY0),
alphac(xalphac)
{
    setDensity(xrho);
    plasticityType = std::string(plasType);

    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 (xrho > 0.0)
    {
        cp   = sqrt((lambda+2.0*mu)/xrho);
        cs   = sqrt(2.0*mu/xrho);
    }
}




splasticMaterial::splasticMaterial(const std::string& name,
                                     const double xE, const double xnu,
                                     const double xrho, const double xHiso,
                                     const double xHkine, const double xY0,
                                     const double xalphac,
                                     const std::string& xplst)
:
smallStrainMaterial(name),
plasticityType(xplst),
E(xE), nu(xnu), bulk(0.0), cp(0.0), cs(0.0),
lambda(0.0), mu(0.0),
Hiso(xHiso), Hkine (xHkine), Y0(xY0), alphac(xalphac)
{
    setDensity(xrho);
    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 (xrho > 0.0)
    {
        cp   = sqrt((lambda+2.0*mu)/xrho);
        cs   = sqrt(2.0*mu/xrho);
    }
}





splasticMaterial::splasticMaterial(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), alphac(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, "alphac", alphac);

    std::string mstr;
    muesli::assignValue(cl, "model", mstr);

    if      (mstr == "von_mises") plasticityType = "mises";
    else if (mstr == "tresca")    plasticityType = "tresca";
    else if (mstr == "drucker")   plasticityType = "drucker";

    // 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 = lambda / 2.0 / (lambda+mu);
        E  = mu*2.0*(1.0+nu);
    }

    // we set all the constants, so that later on all mlog 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(2.0*mu/density());
    }
}





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




smallStrainMP* splasticMaterial::createMaterialPoint() const
{
    smallStrainMP* mp = new splasticMP(*this);
    return mp;
}




// this function is much faster than the one with string property names, because
// avoids string comparisons. It should be used.
double splasticMaterial::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;

        default:
            std::cout << "Error in elastoplasticMaterial. Property not defined";
    }
    return ret;
}


void splasticMaterial::print(std::ostream &of) const
{

    if (plasticityType == "mises")
    {
        of  << "\n Elastoplastic 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;
    }


    else if (plasticityType == "tresca")
    {
        of  << "\n Elastoplastic material for small strain kinematics."
        << "\n   Tresca yield criterion with linear isotropic and kinematic hardenings."
        << "\n   Yield stress           : " << Y0
        << "\n   Isotropic hardening    : " << Hiso
        << "\n   Kinematic hardening    : " << Hkine;
    }


    else if (plasticityType == "drucker")
    {
        of  << "\n Drucker-Praguer yield criterion with linear isotropic hardenings."
        << "\n Yield function --> F = J + alpha*I - (3.0*alpha+sqrt(3.0))/3.0*(Y0-q)"
        << "\n J = sqrt(J2), I = trace(stresstensor)"
        << "\n Conjugate force --> q = -K'(xi), K(xi) = 1/2 Hiso xi^2 "
        << "\n Yield stress, tensile stress for pressure dependent models: Y0    : " << Y0
        << "\n Drucker Prager's cone semiangle:                            alpha : " << alphac
        << "\n Isotropic hardening:                                        Hiso  : " << Hiso;
    }


    of  << "\n   Young modulus:  E      : " << E;
    of  << "\n   Poisson ratio:  nu     : " << nu;
    of  << "\n   Lame constants: lambda : " << lambda;
    of  << "\n   Shear modulus:  mu     : " << mu;
    of  << "\n   Bulk modulus:   k      : " << bulk;
    of  << "\n   Density                : " << density();
    of  << "\n   Wave velocities c_p    : " << cp;
    of  << "\n                   c_s    : " << cs;
}




void splasticMaterial::setRandom()
{
    E   = muesli::randomUniform(1000.0, 10000.0);
    nu  = muesli::randomUniform(0.05, 0.45);
    setDensity(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)/density());
    cs     = sqrt(2.0*mu/density());
    bulk   = lambda + 2.0/3.0*mu;

    int pt = muesli::discreteUniform(0, 2);

    if (pt == 0)
    {
        plasticityType = "mises";
        Y0             = E* muesli::randomUniform(0.5, 1.5);
        Hiso           = E*1e-4 * muesli::randomUniform(1.0, 2.0);
        Hkine          = E*1e-4 * muesli::randomUniform(1.0, 2.0);
    }

    else if (pt == 1)
    {
        plasticityType = "drucker";

        Y0            = E*1e-3* muesli::randomUniform(0.5, 1.5);
        alphac        = muesli::randomUniform(0.1,0.3 );
        Hiso          = E*1e-4 * muesli::randomUniform(1.0, 2.0);
        Hkine         = 0.0;
    }
    else if (pt == 2)
    {
        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;
    }
    
}




bool splasticMaterial::test(std::ostream  &of)
{
    setRandom();
    smallStrainMP* p = this->createMaterialPoint();
    p->setRandom();
    if      (this->plasticityType == "mises")   of << "\n   von Mises type";
    else if (this->plasticityType == "drucker") of << "\n   Drucker-Prager type";
    else if (this->plasticityType == "tresca") of << "\n    Tresca type";

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




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




splasticMP::splasticMP(const splasticMaterial &m) :
muesli::smallStrainMP(m),
theElastoplasticMaterial(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 splasticMP::commitCurrentState()
{
    smallStrainMP::commitCurrentState();
    dg_n   = dg_c;
    ep_n   = ep_c;
    xi_n   = xi_c;
    Xi_n   = Xi_c;
}




void splasticMP::contractWithDeviatoricTangent(const ivector &v1, const ivector &v2, itensor &T) const
{
    const double mu  = theElastoplasticMaterial.mu;
    const double k   = theElastoplasticMaterial.bulk;
    const double Hiso = theElastoplasticMaterial.Hiso;


    if (theElastoplasticMaterial.plasticityType == "mises")
    {
        if (dg_c > 0.0)
        {
            const double   Hkine    = theElastoplasticMaterial.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/(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 (theElastoplasticMaterial.plasticityType == "drucker")
    {
        if (dg_c == 0.0)
        {
            T = - 2.0/3.0*mu        * itensor::dyadic(v1, v2)
            +   mu*v1.dot(v2)       * itensor::identity()
            +   mu                  * itensor::dyadic(v2,v1);
        }
        else
        {
            const istensor edtrial     = istensor::deviatoricPart(eps_c-ep_n);//deviatoric elastic trial strain
            //const istensor evoltrial   = eps_c-ep_n-edtrial;//trial elastic volumetric strain
            const istensor D           = edtrial * (1.0/edtrial.norm());
            const double   alphac      = theElastoplasticMaterial.alphac;
            const double   normed      = sqrt(2.0)*edtrial.norm();
            const double   beta        = (3.0*alphac+sqrt(3.0))/3.0;
            const double   A           = 1.0/(mu+9.0*k*alphac*alphac+beta*beta*Hiso);
            const istensor s_trial     = 2.0 * mu * edtrial;
            const double   J_trial     = sqrt(s_trial.J2());

            bool isInCone = J_trial-mu*dg_c >= 0.0;
            if (isInCone)//case within cone's definition
            {
                T = mu*(1.0-dg_c/normed)* v1.dot(v2)* itensor::identity()
                + mu*(1.0-dg_c/normed)*itensor::dyadic(v2,v1)
                - 2.0*mu*(1.0-dg_c/normed)/3.0*itensor::dyadic(v1, v2)
                + 2.0*mu*(dg_c/normed)*itensor::dyadic(D*v1, D*v2)
                - 2.0*mu*mu*A*itensor::dyadic(D*v1, D*v2);
            }
            else //case beyond the apex
            {
                istensor Q; Q.setZero();

                T = Q;
            }
        }
    }

    else if (theElastoplasticMaterial.plasticityType == "tresca")
    {
        smallStrainMP::contractWithDeviatoricTangent(v1,v2,T);
    }
}




void splasticMP::contractWithTangent(const ivector &v1, const ivector &v2, itensor &T) const
{
    const double mu  = theElastoplasticMaterial.mu;
    const double k   = theElastoplasticMaterial.bulk;
    const double Hkine = theElastoplasticMaterial.Hkine;
    const double Hiso = theElastoplasticMaterial.Hiso;



    if (theElastoplasticMaterial.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 / (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 (theElastoplasticMaterial.plasticityType == "drucker")
    {
        // Elastic part
        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);
        }
        // Elastoplastic part
        else
        {
            const istensor edtrial     = istensor::deviatoricPart(eps_c-ep_n);//deviatoric elastic trial strain
            const istensor D           = edtrial * (1.0/edtrial.norm());
            const double   alphac      = theElastoplasticMaterial.alphac;
            const double   normed      = sqrt(2.0)*edtrial.norm();
            const double   root2       = sqrt(2.0);
            const double   beta        = (3.0*alphac+sqrt(3.0))/3.0;
            const double   A           = 1.0/(mu+9.0*k*alphac*alphac+beta*beta*Hiso);
            const istensor s_trial     = 2.0 * mu * edtrial;
            const double   J_trial     = sqrt(s_trial.J2());

            bool isInCone = J_trial-mu*dg_c >= 0.0;
            if (isInCone)//case within cone's definition
            {
                T =

                //Deviatoric Part
                mu*(1.0-dg_c/normed)* v1.dot(v2)* itensor::identity()
                + mu*(1.0-dg_c/normed)*itensor::dyadic(v2,v1)
                - 2.0/3.0*mu*(1.0-dg_c/normed)*itensor::dyadic(v1, v2)
                + 2.0*mu*(dg_c/normed)*itensor::dyadic(D*v1, D*v2)
                - 2.0*mu*mu*A*itensor::dyadic(D*v1, D*v2)
                - A*root2*mu*3.0*k*alphac*itensor::dyadic(D*v1, v2)

                //Volumetric Part
                + k*itensor::dyadic(v1, v2)
                - A*root2*mu*3.0*k*alphac*itensor::dyadic(v1, D*v2)
                - alphac*alphac*9.0*k*k*A*itensor::dyadic(v1, v2);
            }
            else //case beyond the apex
            {
                const double   Q = 1.0/(3.0*k*alphac+beta*beta/(3.0*alphac)*Hiso);
                T = k*(1.0 - 3.0*k*Q*alphac) * itensor::dyadic(v1, v2);
            }
        }
    }

    else if (theElastoplasticMaterial.plasticityType == "tresca")
    {smallStrainMP::contractWithTangent(v1,v2,T);}
}




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




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




double splasticMP::dissipatedEnergy() const
{
    double ret = 0.0;


    if (theElastoplasticMaterial.plasticityType == "mises")
        ret = sqrt(2.0/3.0) * theElastoplasticMaterial.Y0 * dg_c;

    else if (theElastoplasticMaterial.plasticityType == "drucker")
        ret = (time_c - time_n)*kineticPotential();

    else if (theElastoplasticMaterial.plasticityType == "tresca")
        ret = theElastoplasticMaterial.Y0 * dg_c;

    return ret;
}




// when trial point lies beyond the apex, this return mapping must be used
// return value is the increment in volumetric plastic strain Delta-eps-p-vol
// dpvol is 3*trace(eps^p)

double splasticMP::druckerApexReturn(const splasticMaterial& m, const istensor& sigma, const double& q,
                                       istensor& sigmabar, double& qbar)
{
    const double   alpha = m.alphac;
    const double   Hiso  = m.Hiso;
    const double   k     = m.bulk;
    const double   Y0    = m.Y0;
    const double   beta  = (3.0*alpha+sqrt(3.0))/3.0;

    // starting point, trial state
    const double   Itrial = sigma.invariant1();
    double dpvol  = 0.0;


    // Only valid for linear hardening
    ivector b, c, d;
    d(0) = 9.0*alpha*beta*k/(Hiso*beta*beta+9.0*alpha*alpha*k);
    d(1) = Hiso*beta*beta/(Hiso*beta*beta+9.0*alpha*alpha*k);
    d(2) = -d(0);

    c(0) = d(1);
    c(1) = -Hiso*alpha*beta/(Hiso*beta*beta+9.0*alpha*alpha*k);
    c(2) = 9.0*alpha*alpha*k/(Hiso*beta*beta+9.0*alpha*alpha*k);

    double Ibar;
    Ibar = d(0)*Y0 + d(1)*Itrial + d(2)*q;
    qbar = c(0)*Y0 + c(1)*Itrial + c(2)*q;

    sigmabar = Ibar/3.0*istensor::identity();
    dpvol    = (Itrial-Ibar)/(9.0*k);


    return dpvol;
}




// mapping return when trial point lies within the normal to the smooth part
// of the cone
// Drucker-Prager close point projection
// return delta gamma
// Note: only valid for linear isotropic hardening

double splasticMP::druckerConeProjection(const splasticMaterial& m, const istensor& sigma, const double& q)
{
    const double   mu     = m.mu;
    const double   alphac = m.alphac;
    const double   k      = m.bulk;
    const double   Hiso   = m.Hiso;
    const double   Y0     = m.Y0;
    const double   beta   = (3.0*alphac+sqrt(3.0))/3.0;

    const double   I      = sigma.invariant1();
    const istensor s      = istensor::deviatoricPart(sigma);
    const double   J      = sqrt(s.J2());

    double dg = (J + alphac*I - beta*Y0 + beta*q)/(mu+9.0*alphac*alphac*k+beta*beta*Hiso);

    return dg;
}




// this function updates the current plastic variables for a given strain
void splasticMP::druckerReturn(const istensor& strain, double& dg)
{
    const double mu     = theElastoplasticMaterial.mu;
    const double k      = theElastoplasticMaterial.bulk;
    const double alphac = theElastoplasticMaterial.alphac;
    const double beta   = (3.0*alphac+sqrt(3.0))/3.0;
    const double Hiso   = theElastoplasticMaterial.Hiso;
    const double Y0     = theElastoplasticMaterial.Y0;

    // trial state
    istensor ee_trial    = istensor::deviatoricPart(strain-ep_n);  // trial elastic deviatoric strain tensor
    double   te_trial    = strain.trace() - ep_n.trace();
    istensor s_trial     = 2.0 * mu * ee_trial;                    // deviatoric trial stress tensor
    istensor sigma_trial = s_trial + k* te_trial * istensor::identity();
    double   q_trial     = -Hiso*xi_n;
    istensor Q_trial;
    istensor sigma;
    double q;


    //check yield condition

    const double f_trial = yieldfunction(theElastoplasticMaterial, sigma_trial, Q_trial, q_trial);

    // compute stress and update internal variables
    if (f_trial <= PROJTOL*Y0)
    {
        // elastic step

        ep_c = ep_n;
        xi_c = xi_n;
        dg   = 0.0;
    }

    else  // plastic step: do return mapping
    {
        //closest point projection algorithm
        dg = druckerConeProjection(theElastoplasticMaterial, sigma_trial, q_trial);
        double  J_trial = sqrt(s_trial.J2());

        if (J_trial-mu*dg >= 0.0)
        {
            //correct internal variables
            ep_c   = ep_n + dg*s_trial/(2.0*J_trial) + dg*alphac * istensor::identity();
            xi_c   = (Hiso > 0.0) ? xi_n + beta*dg : 0.0;
            toApex = false;
        }
        else
        {
            //apex dpvol is alphac*Delta gamma

            toApex = true;
            double dpvol = druckerApexReturn(theElastoplasticMaterial, sigma_trial, q_trial, sigma, q);
            ep_c  = ep_n + dpvol* (istensor::identity()) + ee_trial;
            xi_c  = xi_n + dpvol/alphac*beta;
        }
    }

    //Compute the deviatoric and volumetric elastic strain tensor in order to obtain
    //sigma for "assert" function
    istensor ee_c, epsve_c;
    ee_c    = istensor::deviatoricPart(strain-ep_c); //current elastic deviatoric strain tensor
    epsve_c = strain-ep_c-ee_c;                   //current elastic volumetric strain tensor

    assert(yieldfunction(theElastoplasticMaterial,
                         2.0*mu*ee_c+3.0*k*epsve_c, Q_trial, -Hiso*xi_c)/Y0 < 1e-3);
}




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




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




materialState splasticMP::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;
}




materialState splasticMP::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;
}




double splasticMP::kineticPotential() const
{
    const double mu    = theElastoplasticMaterial.mu;
    const double Hiso  = theElastoplasticMaterial.Hiso;
    const double Hkine = theElastoplasticMaterial.Hkine;
    const double k     = theElastoplasticMaterial.bulk;
    double dt_Psistar  = 0.0;

    if (theElastoplasticMaterial.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 (theElastoplasticMaterial.plasticityType == "drucker")
    {
        const istensor e_cedev        = istensor::deviatoricPart(eps_c-ep_c);
        const istensor e_cevol        = eps_c-ep_c-e_cedev;
        const istensor stressdev      = 2.0 * mu * e_cedev;
        const istensor stressvol      = 3.0 * k * e_cevol;
        istensor       sigma          = stressdev+stressvol;
        const double   q              = -Hiso*xi_c;

        dt_Psistar = (ep_c-ep_n).dot(sigma) + (xi_c-xi_n)*q;
    }

    else if (theElastoplasticMaterial.plasticityType == "tresca")
    {
        const istensor e_cedev        = istensor::deviatoricPart(eps_c)-ep_c;
        const istensor stressdev      = 2.0 * mu * e_cedev;
        const double   q              = -Hiso*xi_c;

        dt_Psistar = (ep_c-ep_n).dot(stressdev) + (xi_c-xi_n)*q;
    }

    else
    {
        dt_Psistar = 0.0;
    }

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




double splasticMP::plasticSlip() const
{
    return xi_c;
}




// return to yield surface for von Mises model
void splasticMP::radialReturn(const istensor& strain, double& dg)
{
    const double mu    = theElastoplasticMaterial.mu;
    const double Hiso  = theElastoplasticMaterial.Hiso;
    const double Hkine = theElastoplasticMaterial.Hkine;

    // trial state
    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(theElastoplasticMaterial, 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  / (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 splasticMP::resetCurrentState()
{
    smallStrainMP::resetCurrentState();

    dg_c   = dg_n;
    ep_c   = ep_n;
    Xi_c   = Xi_n;
    xi_c   = xi_n;
}




void splasticMP::setConvergedState(const double theTime, const istensor & strainn,
                                     const double gamman,
                                     const istensor &epn, const double xin, const istensor &Xin)
{
    time_n    = theTime;
    dg_n      = gamman;
    ep_n      = epn;
    xi_n      = xin;
    eps_n     = strainn;
    Xi_n      = Xin;
}




void splasticMP::setRandom()
{
    smallStrainMP::setRandom();

    istensor tmp;
    tmp.setRandom();

    if (theElastoplasticMaterial.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 (theElastoplasticMaterial.plasticityType == "drucker")
    {
        ep_n = istensor::symmetricPartOf(tmp);
        ep_c = ep_n;

        xi_n = muesli::randomUniform(1.0, 2.0);
        xi_c = xi_n;
    }

    else if (theElastoplasticMaterial.plasticityType == "tresca")
    {

        ep_n = istensor::deviatoricPart(tmp);
        ep_c = ep_n;
        xi_n = muesli::randomUniform(1.0, 2.0);
        xi_c = xi_n;

    }
}




double splasticMP::storedEnergy() const
{
    const double mu    = theElastoplasticMaterial.mu;
    const double Hiso  = theElastoplasticMaterial.Hiso;
    const double Hkine = theElastoplasticMaterial.Hkine;
    const double k     = theElastoplasticMaterial.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 splasticMP::stress(istensor& sigma) const
{
    if (theElastoplasticMaterial.plasticityType=="mises")
    {
        deviatoricStress(sigma);
        sigma += theElastoplasticMaterial.bulk*eps_c.trace()*istensor::identity();
    }

    else if (theElastoplasticMaterial.plasticityType=="drucker")
    {
        deviatoricStress(sigma);
        sigma += theElastoplasticMaterial.bulk*(eps_c-ep_c).trace()*istensor::identity();
    }

    else if (theElastoplasticMaterial.plasticityType=="tresca")
    {
        deviatoricStress(sigma);
        sigma += theElastoplasticMaterial.bulk*eps_c.trace()*istensor::identity();
    }
}




void splasticMP::tangentTensor(itensor4& C) const
{
    C.setZero();

    if (theElastoplasticMaterial.plasticityType=="drucker")
    {
        if (dg_c == 0.0)
        {
            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;
                        }
        }
        else
        {
            const double   mu          = theElastoplasticMaterial.mu;
            const double   kb          = theElastoplasticMaterial.bulk;
            const double   Hiso        = theElastoplasticMaterial.Hiso;
            const double   alphac      = theElastoplasticMaterial.alphac;
            const istensor edtrial     = istensor::deviatoricPart(eps_c-ep_n);//deviatoric elastic trial strain
            const istensor D           = edtrial * (1.0/edtrial.norm());
            const double   normed      = sqrt(2.0)*edtrial.norm();
            const double   root2       = sqrt(2.0);
            const double   beta        = (3.0*alphac+sqrt(3.0))/3.0;
            const double   A           = 1.0/(mu+9.0*kb*alphac*alphac+beta*beta*Hiso);
            const istensor s_trial     = 2.0 * mu * edtrial;
            const double   J_trial     = sqrt(s_trial.J2());
            const double   Q           = 1.0/(3.0*kb*alphac+beta*beta/(3.0*alphac)*Hiso);


            bool isInCone = J_trial-mu*dg_c >= 0.0;
            if (isInCone)//case within cone's definition
            {
                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) += -2.0*mu*(1.0-dg_c/normed)/3.0+kb-9.0*kb*kb*alphac*alphac*A;
                                if (i==k && j==l) C(i,j,k,l) += mu*(1.0-dg_c/normed);
                                if (i==l && j==k) C(i,j,k,l) += mu*(1.0-dg_c/normed);
                                if (l==k)         C(i,j,k,l) += - A*2.0/root2*mu*3.0*kb*alphac*D(i,j);
                                if (i==j)         C(i,j,k,l) += - A*2.0/root2*mu*3.0*kb*alphac*D(k,l);
                                C(i,j,k,l) += (2.0*mu*(dg_c/normed)-2.0*mu*mu*A)*D(i,j)*D(k,l);
                            }
            }
            else//case returned to the apex
            {
                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-3.0*kb*kb*Q*alphac;
                            }
            }
        }
    }
    else if (theElastoplasticMaterial.plasticityType=="mises")
    {
        if (dg_c == 0.0)
        {
            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;
                        }
        }
        else
        {
            const double   mu     = theElastoplasticMaterial.mu;
            const double   kb     = theElastoplasticMaterial.bulk;
            const double   Hkine  = theElastoplasticMaterial.Hkine;
            const double   Hiso   = theElastoplasticMaterial.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);

            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 (theElastoplasticMaterial.plasticityType=="tresca")
    {
        const double mu     = theElastoplasticMaterial.mu;
        const double kb     = theElastoplasticMaterial.bulk;
        const double Hiso   = theElastoplasticMaterial.Hiso;
        const double Y0     = theElastoplasticMaterial.Y0;

        // compute stress and update internal variables
        if (dg_c == 0.0)
        {
            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;
                        }
        }
        else  // plastic step: do return mapping
        {
            // trial state
            istensor ee_trial    = istensor::deviatoricPart(eps_c-ep_n);  // trial elastic deviatoric strain tensor
            istensor s_trial     = 2.0 * mu * ee_trial;                   // deviatoric trial stress tensor
            istensor e_trial     = eps_c-ep_n;                            // elastic trial strain
            //double   tetr        = eps_c.trace();
            //istensor sigma_trial = s_trial + kb*tetr*istensor::identity();

            // current state
            istensor ee          = istensor::deviatoricPart(eps_c-ep_c);  // elastic deviatoric strain tensor
            double   te          = eps_c.trace();
            istensor stensor     = 2.0 * mu * ee;                         // deviatoric stress tensor
            istensor sigma       = stensor + kb * te * istensor::identity();

            // 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
            sigma.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;
            }
            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]);
                }
                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]);
                }

            }
            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;
                }
            }
            //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) C(i,j,k,l) += dsigdetr[owentgtowen[0]][owentgtowen[1]];
                                if (i==k && j==l) C(i,j,k,l) += 0.5*(dsigdetr[owentgtowen[0]][owentgtowen[0]]-dsigdetr[owentgtowen[0]][owentgtowen[1]]);
                                if (i==l && j==k) C(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++){
                        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]]);

                        if (i==k && j==l) C(i,j,k,l) += (commont)*(-0.5*(xb+xc));
                        if (i==l && j==k) C(i,j,k,l) += (commont)*(-0.5*(xb+xc));
                        if (i==k)         C(i,j,k,l) += (commont)*(0.5*e_trial(l,j));
                        if (i==l)         C(i,j,k,l) += (commont)*(0.5*e_trial(k,j));
                        if (j==l)         C(i,j,k,l) += (commont)*(0.5*e_trial(i,k));
                        if (k==j)         C(i,j,k,l) += (commont)*(0.5*e_trial(i,l));

                        C(i,j,k,l) += -commont*term1*Ea(i,j)*Ea(k,l);
                        C(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]]);
                        C(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;
                }
                //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;
                }

                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)         C(i,j,k,l) += s1*(0.5*e_trial(l,j));
                    if (i==l)         C(i,j,k,l) += s1*(0.5*e_trial(k,j));
                    if (j==l)         C(i,j,k,l) += s1*(0.5*e_trial(i,k));
                    if (k==j)         C(i,j,k,l) += s1*(0.5*e_trial(i,l));
                    if (i==k && j==l) C(i,j,k,l) += -s2*0.5;
                    if (i==l && j==k) C(i,j,k,l) += -s2*0.5;
                    if (k==l)         C(i,j,k,l) += s4*e_trial(i,j);
                    if (i==j)         C(i,j,k,l) += s5*e_trial(k,l);
                    if (i==j && k==l) C(i,j,k,l) += -s6;
                    C(i,j,k,l) += -s3*e_trial(i,j)*e_trial(k,l);
                }
            }
        }
    }
}




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




// this function updates the current plastic variables for a given strain
void splasticMP::trescaReturn(const istensor& strain, double& dg)
{
    const double mu     = theElastoplasticMaterial.mu;
    const double k      = theElastoplasticMaterial.bulk;
    const double Hiso   = theElastoplasticMaterial.Hiso;
    const double Y0     = theElastoplasticMaterial.Y0;

    // trial state
    istensor ee_trial    = istensor::deviatoricPart(strain-ep_n);  // trial elastic deviatoric strain tensor
    double   te_trial    = strain.trace();
    istensor s_trial     = 2.0 * mu * ee_trial;                    // deviatoric trial stress tensor
    istensor sigma_trial = s_trial + k * te_trial * istensor::identity();
    double   q_trial     = -Hiso*xi_n;
    istensor sigma;
    istensor Qdummy;

    //check yield condition
    double f_trial = yieldfunction(theElastoplasticMaterial, sigma_trial, Qdummy, q_trial);//check consistency condition

    // compute stress and update internal variables

    if (f_trial <= PROJTOL*Y0)
    {
        // elastic step
        ep_c = ep_n;
        xi_c = xi_n;
        dg   = 0.0;
    }

    else  // plastic step: do return mapping, f_trial > PROJTOL*Y0
    {

        ivector str, s;
        ivector evec[3];
        s_trial.spectralDecomposition(evec, str);//eigenvalues and eigenvectors

        //
        while (true)//main plane return by default
        {
            //f_trial = str(2) - str(0) - (Y0+Hiso*xi_n);
            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*(xi_n + dga + dgb));
                res[1] = sb - 2.0*mu*(dga + 2.0*dgb)-(Y0+Hiso*(xi_n + dga + dgb));


                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*(xi_n + dga + dgb));
                res[1] = sb - 2.0*mu*(dga + 2.0*dgb)-(Y0+Hiso*(xi_n + dga + dgb));
                
                
                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;
            }
            
        }
        
        
        
        //Update of internal variables in plastic step
        xi_c  = xi_n + dg;
        ep_c = istensor::deviatoricPart(strain);
        for (unsigned a=0; a<3; a++) ep_c.addScaledVdyadicV(-s(a)/(2.0*mu), evec[a]);
        
        
    }
    
    //Compute the deviatoric and volumetric elastic strain tensor in order to obtain
    //sigma for "assert" function
    istensor ee_c, epsve_c;
    ee_c    = istensor::deviatoricPart(strain)-ep_c; //current elastic deviatoric strain tensor
    epsve_c = strain-ep_c-ee_c;                   //current elastic volumetric strain tensor
    
    assert(yieldfunction(theElastoplasticMaterial,
                         2.0*mu*ee_c+3.0*k*epsve_c, Qdummy, -Hiso*xi_c)/Y0 < 1e-5);
}




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

    
    smallStrainMP::updateCurrentState(theTime, strain);
    
    if (theElastoplasticMaterial.plasticityType == "mises")
    {
        radialReturn(strain, dg_c);
    }
    
    else if (theElastoplasticMaterial.plasticityType == "drucker")
    {
        druckerReturn(strain, dg_c);
    }
    
    else if (theElastoplasticMaterial.plasticityType == "tresca")
    {
        trescaReturn(strain, dg_c);
    }
    
    eps_c = strain;
    
}




double splasticMP::volumetricStiffness() const
{
    double ret = 0.0;
    double k   = theElastoplasticMaterial.bulk;
    if (theElastoplasticMaterial.plasticityType == "mises" || theElastoplasticMaterial.plasticityType == "tresca")
    {
        ret = k;
    }
    
    if (theElastoplasticMaterial.plasticityType == "drucker")
    {
        const double alphac = theElastoplasticMaterial.alphac;
        const double mu     = theElastoplasticMaterial.mu;
        const double beta   = (3.0*alphac+sqrt(3.0))/3.0;
        const double Hiso   = theElastoplasticMaterial.Hiso;
        const double A      = 1.0/(mu+9.0*k*alphac*alphac+beta*beta*Hiso);
        const double Q      = 1.0/(3.0*k*alphac+beta*beta/(3.0*alphac)*Hiso);
        if (dg_c == 0.0)
        ret = k;
        else
        {
            const           istensor e_cedev        = istensor::deviatoricPart(eps_c-ep_n);
            const           istensor stressdevtrial = 2.0 * mu * e_cedev;
            const double    J_trial                 = sqrt(stressdevtrial.J2());
            
            if (J_trial-mu*dg_c >= 0.0)
            {
                ret = k-9.0*k*k*alphac*alphac*A;
            }
            else
            {
                ret = (1.0-3.0*k*Q*alphac)*k;
            }
        }
    }
    
    return ret;
}




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




double splasticMP::yieldfunction(const splasticMaterial& m, const istensor& sigma, const istensor& Q, const double& q)
{
    const istensor s = istensor::deviatoricPart(sigma);
    const double   Y0 = m.Y0;
    
    if (m.plasticityType == "drucker")
    {
        const double   I      = sigma.invariant1();
        const double   J      = sqrt(s.J2());
        const double   alphac = m.alphac;
       
        //std::cout << "splastic, yieldF, I, J, alphac, res  = " << I << " , " << J << " , " << alphac << " , " << J + alphac*I - (3.0*alphac+sqrt(3.0))/3.0*(Y0-q) << std::endl;
        return J + alphac*I - (3.0*alphac+sqrt(3.0))/3.0*(Y0-q);
    }
    
    else 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;
    }
}
