/****************************************************************************
*
*                                 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/>.
*
****************************************************************************/



#define _USE_MATH_DEFINES


#include <cassert>
#include <cmath>
#include <iomanip>
#include <string.h>
#include "sdamage.h"
//#include "muesli/Math/matrix.h"
//#include "muesli/Math/realvector.h"


/*
 * splastic.cpp,
 * infinitesimal kinematics plastic element
 * e.m. andres, j.segurado, i. romero, oct 2017
 *
 *
 *  1) Gurson-Tvergaard-Needleman model with isotropic damage with no hardening
 *
 *     The yield function depends on:
 *        - t:        trace of sigma
 *        - sigma_eq: (3/2)^0.5*sigma_dev
 *
 *        - fv:       volume fraction of cavities or porosity
 *
 *        where sigma_dev is the square root of the deviatoric stress tensor : deviatoric stress tensor
 *
 *     We = (1-d)*e_elastic:C_e_elastickappa/2 (theta)^2 + mu * ||e_elastic||^2
 *     Wp = NONE (at this moment)
 *     f(sigma, Y(fv)) = (sigma_eq/yield)^2 +
 *                        + 2*q1*fv*cosh((1/2)*q2*t^2/yield) -
 *                        - 1 - (q1*fv)^2;
 *
 */



#define FVMAXITER   20        // max num of allowed iterations to compute fv from d value
#define FVTOL       1e-10     // residual tolerance to compute fv from d value
#define TOL         1e-7      // residual tolerance for return-to-yield-surface solutions
#define MAXITER     1e2       // max num of allowed iterations in return-to-yield solutions


using namespace std;
using namespace muesli;




sdamageMaterial::sdamageMaterial(const std::string& name,
                                     const double xE,     const double xnu,
                                     const double xrho,   const double xq1,
                                     const double xq2, const double xY0,
                                     const char*  damgType)
:
smallStrainMaterial(name),
E(xE), nu(xnu), bulk(0.0),
cp(0.0), cs(0.0), lambda(0.0), mu(0.0), rho(xrho),
q1(xq1), q2(xq2), Y0(xY0)
{
    damageType = std::string(damgType);
    
    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);
    }
}




sdamageMaterial::sdamageMaterial(const std::string& name,
                                     const double xE, const double xnu,
                                     const double xrho, const double xq1,
                                     const double xq2, const double xY0,
                                     const std::string& xdamg)
:
smallStrainMaterial(name),
damageType(xdamg),
E(xE), nu(xnu), bulk(0.0), cp(0.0), cs(0.0),
lambda(0.0), mu(0.0), rho(xrho),
q1(xq1), q2(xq2), Y0(xY0)
{
    
    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);
    }
}




sdamageMaterial::sdamageMaterial(const std::string& name,
                                     const materialProperties& cl)
:
smallStrainMaterial(name, cl),
damageType("gurson"),
E(0.0), nu(0.0), bulk(0.0), cp(0.0), cs(0.0), lambda(0.0), mu(0.0),
q1(0.0), q2(0.0), Y0(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, "q1_gurson",   q1);
    muesli::assignValue(cl, "q2_gurson",   q2);
    muesli::assignValue(cl, "yieldstress", Y0);
    
    std::string mstr;
    muesli::assignValue(cl, "model", mstr);
    
    if      (mstr == "gurson_tn") damageType = "gurson";
    
    
    
    // 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 (rho > 0.0)
    {
        cp   = sqrt((lambda+2.0*mu)/rho);
        cs   = sqrt(2.0*mu/rho);
    }
}




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




smallStrainMP* sdamageMaterial::createMaterialPoint() const
{
    smallStrainMP* mp = new sdamageMP(*this);
    return mp;
}




double sdamageMaterial::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 sdamageMaterial::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_Q1_GURSON:  ret = q1;       break;
        case PR_Q2_GURSON:  ret = q2;       break;
            
        default:
            cout << "Error in damageMaterial. Property not defined";
    }
    return ret;
}




void sdamageMaterial::print(std::ostream &of) const
{
    
    if (damageType == "gurson")
    {
        of  << "\n Elastoplastic damaged material for small strain kinematics."
        << "\n   Damage model: Gurson-Tvergaard-Needleman yield model with no hardening."
        << "\n   Yield stress           : " << Y0
        << "\n   q1 GTN parameter       : " << q1
        << "\n   q2 GTN parameter       : " << q2;
    }
    
    else
    {
        cout << "\n Error in damageType. Type of damage not defined" << flush;
    }
    
    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                : " << rho;
    
    if (rho > 0.0)
    {
        of  << "\n   Wave velocities c_p    : " << cp;
        of  << "\n                   c_s    : " << cs;
    }
}




void sdamageMaterial::setRandom()
{

    
    E      = muesli::randomUniform(100.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, 0);
    
    if (pt == 0)
    {
        
        damageType = "gurson";
        Y0         = E * 1e-3 * muesli::randomUniform(0.5, 1.5);
        q1         = muesli::randomUniform(1.0, 2.0);
        q2         = 1e-4 * muesli::randomUniform(0.5, 1.5);
    }
    
    else
    {
        cout << "\n Error in damageType. Type of damage not defined" << flush;
    }
    
}




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

    
    if      (this->damageType == "gurson") of << "\n   Gurson TN type";
    
    else
    {
        cout << "\n Error in damageType. Type of damage not defined" << flush;
    }
    
    bool ok = p->testImplementation(of);
    delete p;
    
    
    return ok;
}




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




sdamageMP::sdamageMP(const sdamageMaterial &m) :
muesli::smallStrainMP(m),
theDamageMaterial(m),
d_n(0.0), dg_n(0.0),
d_c(0.0), dg_c(0.0)
{
    ep_n.setZero();
    ep_c.setZero();
}




void sdamageMP::commitCurrentState()
{

    
    smallStrainMP::commitCurrentState();
    d_n  = d_c;
    dg_n = dg_c;
    ep_n = ep_c;
    
}




void sdamageMP::contractWithDeviatoricTangent(const ivector &v1, const ivector &v2, itensor &T) const
{
    smallStrainMP::contractWithDeviatoricTangent(v1,v2,T);
}




void sdamageMP::contractWithTangent(const ivector &v1, const ivector &v2, itensor &T) const
{
    smallStrainMP::contractWithTangent(v1,v2,T);
}




double sdamageMP::damage() const
{
    return d_c;
}




double sdamageMP::deviatoricEnergy() const
{
    istensor eps_e = istensor::deviatoricPart(eps_c - ep_c);
    return (1.0 - d_c) * theDamageMaterial.mu * eps_e.contract(eps_e);
}




void sdamageMP::deviatoricStress(istensor& s) const
{
    s = (1.0 - d_c) * 2.0 * theDamageMaterial.mu * (istensor::deviatoricPart(eps_c-ep_c));
}




double sdamageMP::dissipatedEnergy() const 
{
    const double k  = theDamageMaterial.bulk;
    const double mu = theDamageMaterial.mu;
    
    if (theDamageMaterial.damageType=="gurson")
    {
        const istensor ee_c       = istensor::deviatoricPart(eps_c-ep_c);
        istensor       ee_vol_c   = (eps_c-ep_c) - ee_c;
        istensor       stress_c   = (1.0 - d_c) * (2.0*mu*ee_c + 3.0*k* ee_vol_c);
        const double   dt_Psistar = istensor(ep_c - ep_n).dot(stress_c);
        
        return dt_Psistar;
    }
    
    else
    {
        return 0.0;
    }
}



double sdamageMP::effectiveStoredEnergy() const
{
    const double k         = theDamageMaterial.bulk;
    const double mu        = theDamageMaterial.mu;
    
    if (theDamageMaterial.damageType == "gurson")
    {
        istensor ee_c       = istensor::deviatoricPart(eps_c-ep_c);
        istensor ee_vol_c   = (eps_c-ep_c) - ee_c;
        istensor stress_c   = (1.0 - d_c) * (2.0*mu*ee_c + 3.0*k* ee_vol_c);
        double   t          = (eps_c-ep_c).trace();
        
        double   We         = (1.0 - d_c) * (mu*ee_c.squaredNorm() + 0.5*k*t*t);
        double   Wp         = 0.0;
        
        double   dt_Psistar = istensor(ep_c-ep_n).dot(stress_c);
        
        return We+Wp+dt_Psistar;
    }
    
    else
    {
        return 0.0;
    }
}




double sdamageMP::eq_stress() const
{
    const double mu        = theDamageMaterial.mu;
    
    if (theDamageMaterial.damageType == "gurson")
    {
        istensor ee_c      = istensor::deviatoricPart(eps_c-ep_c);
        istensor stress_c  = (1.0 - d_c) * 2.0*mu*ee_c;
        return sqrt(1.5)*stress_c.norm();
    }
    
    else
    {
        return 0.0;
    }
}




double sdamageMP::fvFunction(const sdamageMaterial& m, const double& d, const double& t)
{
    if (m.damageType == "gurson")
    {
        const double   q1        = m.q1;
        const double   q2        = m.q2;
        
        double         aux_cosh  = 0.5*q2*t;
        double         fv        = 0.0;
        double         fv_check  = (1.0-d)*(1.0-d) - 1.0 - (q1*fv)*(q1*fv) + 2.0*q1*fv*cosh(aux_cosh);
        double         dfv_check = 0.0;
        size_t         count     = 0;
        
        
        while (abs(fv_check) >= FVTOL && ++count < FVMAXITER)
        {
            fv_check  = (1.0-d)*(1.0-d) - 1.0 - (q1*fv)*(q1*fv) + 2.0*q1*fv*cosh(aux_cosh);
            dfv_check = 2.0*q1*cosh(aux_cosh) - 2.0*q1*q1*fv;

            fv = fv - fv_check/dfv_check;
        }
        
        if (count == FVMAXITER)
        {
            cout << "\n   computation for fv as a function of d not converged:\n";
        }
        
        return fv;
    }
    
    else
        
    {
        return false;
    }
}




istensor sdamageMP::getConvergedPlasticStrain() const
{
    return ep_n;
}




materialState sdamageMP::getConvergedState() const
{
    materialState state_n = smallStrainMP::getConvergedState();
    
    state_n.theDouble.push_back(d_n);
    state_n.theDouble.push_back(dg_n);
    state_n.theStensor.push_back(ep_n);
    
    return state_n;
}




istensor sdamageMP::getCurrentPlasticStrain() const
{
    return ep_c;
}




materialState sdamageMP::getCurrentState() const
{
    materialState state_c = smallStrainMP::getCurrentState();
    
    state_c.theDouble.push_back(d_c);
    state_c.theDouble.push_back(dg_c);
    state_c.theStensor.push_back(ep_c);
    //state_c.theStensor.push_back(fv_c);
    
    return state_c;
}




// Gurson mapping return: return to yield surface for Gurson-Tvergaard-Needleman model
// Non linear system of 5 equations with 5 unknowns: tr(sigma_c),||s_c||, d_c, fv_c, dg_c
void sdamageMP::gursonReturn(double &dg, const istensor& strain)
{
    const double k  = theDamageMaterial.bulk;
    const double mu = theDamageMaterial.mu;
    
    if (theDamageMaterial.damageType == "gurson")
        
    {
        const double   q1   = theDamageMaterial.q1;
        const double   q2   = theDamageMaterial.q2;
        const double   Y0   = theDamageMaterial.Y0;
        
        // trial state
        
        double   dTrial     = d_n;                                            // trial value of damage variable
        double   dgTrial    = dg;
        istensor eTrial     = strain-ep_n;                                    // trial elastic strain tensor
        istensor eeTrial    = istensor::deviatoricPart(strain-ep_n);        // trial elastic deviatoric strain tensor
        istensor eVolTrial  = eTrial - eeTrial;                               // trial elastic volumetric strain tensor
        istensor sigmaTrial = (1.0-d_n) * (3.0*k*eVolTrial + 2.0*mu*eeTrial); // trial stress tensor
        istensor sTrial     = istensor::deviatoricPart(sigmaTrial);         // trial deviatoric stress tensor
        double   ssTrial    = sTrial.norm()/Y0;                               // norm of the trial deviatoric stress tensor over Y0
        double   tsTrial    = sigmaTrial.trace()/Y0;                          // trace of the trial stress tensor over Y0
        
        double   fvTrial    =  fvFunction(theDamageMaterial,dTrial,tsTrial);  // trial volume fraction of voids

        
        //check yield condition
        const double fTrial = yieldfunction(theDamageMaterial, fvTrial, ssTrial, tsTrial);

        
        // compute stress and update internal variables
        
        // elastic step
        if (fTrial <= TOL)
        {
            
            d_c  = d_n;
            ep_c = ep_n;
            dg   = 0.0;
        }

        
        // plastic step: do return mapping
        else
        {
            // initial shot for Newton-Raphson solver
            
            
            double gursonSolution[5];  // array containing the trial value of the variables in the system of equations for Gurson model
            
            gursonSolution[0] = tsTrial;
            gursonSolution[1] = ssTrial;
            gursonSolution[2] = dTrial;
            gursonSolution[3] = fvTrial;
            gursonSolution[4] = dgTrial;
            
            realvector gursonIndependentTerm = systemEqFunction(theDamageMaterial,gursonSolution,ssTrial,tsTrial,dTrial); // value of the system of equations for Gurson model for the trial variables
            
            matrix gursonJacobian;
            gursonJacobian.setZero();
            
            
            // apply Newton-Raphson method
            
            double error = gursonIndependentTerm.norm();
            
            size_t count = 0;

            while (error > TOL && ++count < MAXITER)
            {
                gursonIndependentTerm = systemEqFunction(theDamageMaterial,gursonSolution,ssTrial,tsTrial,dTrial); // array containing the value of the system of equations of the Gurson model for the new computed value of the variables
                
                gursonJacobian = jacobianFunction(theDamageMaterial,gursonSolution,dTrial,ssTrial,tsTrial);  // Jacobian matrix of the system of equations for Gurson model
                
                // solve Ax=b (gursonJacobian * gursonSolution = gursonIndependentTerm), overwriting b with x
                gursonJacobian.solveFull(gursonIndependentTerm);
                
                size_t j = gursonIndependentTerm.size();
                
                for (size_t k=0; k<j; k++)
                {
                    gursonSolution[k] += - gursonIndependentTerm[k];
                }
                
                error = systemEqFunction(theDamageMaterial,gursonSolution,ssTrial,tsTrial,dTrial).norm();
                realvector error_vector = systemEqFunction(theDamageMaterial,gursonSolution,ssTrial,tsTrial,dTrial);
            }
            
            if (count == MAXITER) cout << "\n   NR for gursonReturn not converged \n" ;
            
            istensor s_c = (gursonSolution[1]/(ssTrial))*sTrial;
            
            // correct internal variables
            d_c  = gursonSolution[2];
            dg   = gursonSolution[4];
            ep_c = ep_n + dg*(3.0*s_c/(Y0*Y0)+(q1*q2*gursonSolution[3]/Y0)*sinh(0.5*q2*gursonSolution[0])*istensor::identity());
        }
    }
}




double sdamageMP::kineticPotential() const
{
    const double k  = theDamageMaterial.bulk;
    const double mu = theDamageMaterial.mu;
    double dt_Psistar = 0.0;

    if (theDamageMaterial.damageType=="gurson")
    {
        const istensor ee_c       = istensor::deviatoricPart(eps_c-ep_c);
        istensor       ee_vol_c   = (eps_c-ep_c) - ee_c;
        istensor       stress_c   = (1.0 - d_c) * (2.0*mu*ee_c + 3.0*k* ee_vol_c);
        dt_Psistar = istensor(ep_c - ep_n).dot(stress_c);
    }

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




// Function which provides the Jacobian matrix of the system of equations for the damage model
matrix sdamageMP::jacobianFunction (const sdamageMaterial& m, double *variables, const double d, const double s, const double t)
{
    const double mu       = m.mu;
    const double k        = m.bulk;
    
    matrix matrix_jacobian(5,5);
    matrix_jacobian.setZero();
    
    // Function which provides the Jacobian matrix of the system of equations for Gurson model
    if (m.damageType == "gurson")
    {
        const double  q1           = m.q1;
        const double  q2           = m.q2;
        const double  Y0           = m.Y0;
        
        const double  aux_A        = 1.0 + (1.0 - variables[2])*6.0*mu*variables[4]/(Y0*Y0);
        const double  aux_B        = (1.0 - variables[2])/(1.0-d);
        const double  aux_cosh     = 0.5*q2*variables[0];
        
        matrix_jacobian(0,0) = 1.0 + 9.0*k*q1*q2*q2*(1.0 - variables[2])*variables[3]*variables[4]*cosh(aux_cosh)/(2.0*Y0*Y0);
        matrix_jacobian(0,1) = 0.0;
        matrix_jacobian(0,2) = t/(1.0 - d) - 9.0*k*q1*q2*variables[3]*variables[4]*sinh(aux_cosh)/(Y0*Y0);
        matrix_jacobian(0,3) = (1.0 - variables[2])*9.0*k*q1*q2*variables[4]*sinh(aux_cosh)/(Y0*Y0);
        matrix_jacobian(0,4) = (1.0 - variables[2])*9.0*k*q1*q2*variables[3]*sinh(aux_cosh)/(Y0*Y0);
        
        matrix_jacobian(1,0) = 0.0;
        matrix_jacobian(1,1) = 2.0*aux_A*aux_A*variables[1];
        matrix_jacobian(1,2) = -12.0*mu*aux_A*variables[1]*variables[1]*variables[4]/(Y0*Y0) + 2.0*s*s*aux_B/(1.0-d);
        matrix_jacobian(1,3) = 0.0;
        matrix_jacobian(1,4) = 12.0*mu*variables[1]*variables[1]*aux_A*(1.0 - variables[2])/(Y0*Y0);

        matrix_jacobian(2,0) = -variables[4]*9.0*k*q1*q2*q2*variables[3]*(1.0 - variables[2])*(1.0 - variables[2])*cosh(aux_cosh)/(2.0*Y0*Y0*variables[0]) + variables[4]*9.0*k*q1*q2*variables[3]*(1.0 - variables[2])*(1.0 - variables[2])*sinh(aux_cosh)/(Y0*Y0*variables[0]*variables[0]);
        matrix_jacobian(2,1) = 0.0;
        matrix_jacobian(2,2) = 1.0 + variables[4]*18.0*k*q1*q2*variables[3]*(1.0 - variables[2])*sinh(aux_cosh)/(Y0*Y0*variables[0]);
        matrix_jacobian(2,3) = -variables[4]*9.0*k*q1*q2*(1.0 - variables[2])*(1.0 - variables[2])*sinh(aux_cosh)/(Y0*Y0*variables[0]);
        matrix_jacobian(2,4) = -variables[3]*9.0*k*q1*q2*(1.0 - variables[2])*(1.0 - variables[2])*sinh(aux_cosh)/(Y0*Y0*variables[0]);
        
        matrix_jacobian(3,0) = q1*q2*variables[3]*sinh(aux_cosh);
        matrix_jacobian(3,1) = 3.0*variables[1];
        matrix_jacobian(3,2) = 0.0;
        matrix_jacobian(3,3) = 2.0*q1*cosh(aux_cosh) - 2.0*q1*q1*variables[3];
        matrix_jacobian(3,4) = 0.0;

        matrix_jacobian(4,0) = q1*q2*variables[3]*sinh(aux_cosh);
        matrix_jacobian(4,1) = 0.0;
        matrix_jacobian(4,2) = -2.0*(1.0 - variables[2]);
        matrix_jacobian(4,3) = -2.0*q1*q1*variables[3] + 2.0*q1*cosh(aux_cosh);
        matrix_jacobian(4,4) = 0.0;
    }
    return matrix_jacobian;
}




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




void sdamageMP::resetCurrentState()
{

    
    smallStrainMP::resetCurrentState();
    
    d_c    = d_n;
    dg_c   = dg_n;
    ep_c   = ep_n;
    
}




void sdamageMP::setConvergedState(const double theTime, const double &dn, const double dgn,
                                     const istensor &epn, const istensor &strainn)
{

    
    time_n    = theTime;
    d_n       = dn;
    dg_n      = dgn;
    ep_n      = epn;
    eps_n     = strainn;
    
}




void sdamageMP::setRandom()
{

    
    smallStrainMP::setRandom();
    
    istensor tmp;
    tmp.setRandom();
    
    if (theDamageMaterial.damageType == "gurson")
    {
        ep_n = tmp;
        ep_c = ep_n;
        
        d_n  = muesli::randomUniform(0.0, 0.5);
        d_c  = d_n;
    }

    else
    {
        cout << "\n Error initializating random damage material" << flush;
    }

}




double sdamageMP::storedEnergy() const
{
    const double k     = theDamageMaterial.bulk;
    const double mu    = theDamageMaterial.mu;
    
    double       th    = eps_c.trace() - ep_c.trace();
    double       We    = (1.0 - d_c) * (mu*(istensor::deviatoricPart(eps_c-ep_c)).squaredNorm() + 0.5 * k * th*th);
    double       Wp    = 0.0;
    
    return We+Wp;
}




void sdamageMP::stress(istensor& sigma) const
{
    if (theDamageMaterial.damageType=="gurson")
    {
        deviatoricStress(sigma);
        sigma += (1.0 - d_c)*theDamageMaterial.bulk*(eps_c-ep_c).trace()*istensor::identity();
    }
}

    
    
// Function which provides the system of equations for the damage model
realvector sdamageMP::systemEqFunction(const sdamageMaterial& m, double *variables, const double s, const double t, const double d)
{
    const double mu       = m.mu;
    const double k        = m.bulk;
        
    // Function which provides the system of equations for Gurson model
    if (m.damageType == "gurson")
    {
        const double   q1       = m.q1;
        const double   q2       = m.q2;
        const double   Y0       = m.Y0;
        
        const double   aux_A    = 1.0 + (1.0 - variables[2])*6.0*mu*variables[4]/(Y0*Y0);
        const double   aux_B    = (1.0 - variables[2])/(1.0-d);
        const double   aux_cosh = 0.5*q2*variables[0];
        
        realvector     system_eq(5);
        
        system_eq[0] = variables[0] - (1.0-variables[2])*t/(1.0-d) + (1.0-variables[2])*variables[4]*9.0*k*q1*q2*variables[3]*sinh(aux_cosh)/(Y0*Y0);

        system_eq[1] = variables[1]*variables[1]*aux_A*aux_A - aux_B*aux_B*s*s;
        
        system_eq[2] = variables[2] - d - variables[4]*variables[3]*9.0*k*q1*q2*(1.0 - variables[2])*(1.0 - variables[2])*sinh(aux_cosh)/(Y0*Y0*variables[0]);
        
        system_eq[3] = 1.5*variables[1]*variables[1] + 2.0*q1*variables[3]*cosh(aux_cosh) - 1.0 - q1*q1*variables[3]*variables[3];
        
        system_eq[4] = (1.0 - variables[2])*(1.0 - variables[2]) - 1.0 - q1*q1*variables[3]*variables[3] + 2.0*q1*variables[3]*cosh(aux_cosh);
        
        return system_eq;
    }
    else
    {
        cout << "\n Error create system of equations for random damage material" << flush;
        return realvector(5);
    }
}




thPotentials sdamageMP::thermodynamicPotentials() const
{
    thPotentials tp;

    return tp;
}




void sdamageMP::tangentTensor(itensor4& C) const
{
    C.setZero();
    
    const double   lambda      = theDamageMaterial.lambda;
    const double   mu          = theDamageMaterial.mu;
    
    if (theDamageMaterial.damageType=="gurson")
    {
        // STEP 1: build the fourth order elastic stiffness tensor
        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) += lambda;
                        if (i==k && j==l) C(i,j,k,l) += mu;
                        if (i==l && j==k) C(i,j,k,l) += mu;
                    }

       if (dg_c == 0.0)
        {
            // STEP 2: build the fourth order damaged elastic stiffness tensor
            C = (1.0-d_n)*C;
        }
        
                
       else
        {
            // STEP 2: build the fourth order damaged elastic stiffness tensor
            // PENDING TO BE CORRECTED (e.andres 20171024)
            C = (1.0-d_n)*C;
        }
    }
}



double sdamageMP::uniaxialStiffness() const
{
    //Need to be defined
    return 0.0;
}




double sdamageMP::uniaxialStress() const
{
    //Need to be defined
    return 0.0;
}




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

    
    smallStrainMP::updateCurrentState(theTime, strain);
    
    if (theDamageMaterial.damageType == "gurson")
    {
        gursonReturn(dg_c, strain);
    }
    
    eps_c  = strain;
    
}




double sdamageMP::pressure() const
{
    istensor sigma;
    stress(sigma);
    double p = -sigma.trace()/3.0;
    return p;
}




double sdamageMP::void_fraction() const
{
    const double k        = theDamageMaterial.bulk;
    const double Y0       = theDamageMaterial.Y0;
    
    if (theDamageMaterial.damageType == "gurson")
    {
        double   d        = d_c;
        istensor ee_dev_c = istensor::deviatoricPart(eps_c-ep_c);
        istensor ee_vol_c = (eps_c-ep_c) - ee_dev_c;
        double   t        = ((1.0 - d_c) * (3.0*k*ee_vol_c)).trace()/Y0;
        
        return fvFunction(theDamageMaterial,d,t);
    }
    
    else
    {
        return 0.0;
    }
}




double sdamageMP::vol_pslip() const
{
    
    if (theDamageMaterial.damageType == "gurson")
    {
        istensor ep_dev_c  = istensor::deviatoricPart(ep_c);
        istensor ep_vol_c  = ep_c - ep_dev_c;
        
        return ep_vol_c.trace();
    }
    
    else
    {
        return 0.0;
    }
}




double sdamageMP::vol_stress() const
{
    const double k             = theDamageMaterial.bulk;
    
    if (theDamageMaterial.damageType == "gurson")
    {
        istensor ee_c          = istensor::deviatoricPart(eps_c-ep_c);
        istensor ee_vol_c      = (eps_c-ep_c) - ee_c;
        istensor stress_vol_c  = (1.0 - d_c) * 3.0*k*ee_vol_c;
        
        return   stress_vol_c.trace();
    }
    
    else
    {
        return 0.0;
    }
}




double sdamageMP::volumetricEnergy() const
{
    double th = (eps_c-ep_c).trace();
    double V  = (1.0 - d_c) * 0.5 * theDamageMaterial.bulk * th * th;

    return V;
}

    
    
    
double sdamageMP::yieldfunction(const sdamageMaterial& m, const double& fv, const double& s, const double& t)
{
    if (m.damageType == "gurson")
    {
        const double q1       = m.q1;
        const double q2       = m.q2;
        
        const double aux_cosh = 0.5*q2*t;
        
        return 1.5*s*s + 2.0*q1*fv*cosh(aux_cosh) - 1.0 - (q1*fv)*(q1*fv);
    }
    else return false;
}
