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


using namespace std;
using namespace muesli;


JCDamageModel::JCDamageModel(const std::string& name,
                             const materialProperties& cl)
:
edot0(1.0), modelRefTemp(1.0), meltT(1.0), D1(0.0),
D2(0.0), D3(0.0), D4(0.0), D5(0.0)
{
    muesli::assignValue(cl, "edot0",       edot0);
    muesli::assignValue(cl, "modelreftemp",modelRefTemp);
    muesli::assignValue(cl, "melttemp",    meltT);
    muesli::assignValue(cl, "d1",          D1);
    muesli::assignValue(cl, "d2",          D2);
    muesli::assignValue(cl, "d3",          D3);
    muesli::assignValue(cl, "d4",          D4);
    muesli::assignValue(cl, "d5",          D5);
}




// Alternative material creation method, for MP individual tests
JCDamageModel::JCDamageModel(const std::string& name,
                             const double x_edot0, const double x_refT,
                             const double x_meltT, const double x_D1,
                             const double x_D2, const double x_D3,
                             const double x_D4, const double x_D5)
:
edot0(x_edot0), modelRefTemp(x_refT), meltT(x_meltT),
D1(x_D1), D2(x_D2), D3(x_D3),
D4(x_D4), D5(x_D5)
{
}




double JCDamageModel::calculateStrainToFracture(double eps, double epsdot,
                                      istensor sigma, double temp) const
{
    double epsFracture = 0.0;
    const double tempterm  = (temp/meltT <= 1.0) ? (temp-modelRefTemp)/(meltT-modelRefTemp) : 0.0;
    double hydro = (sigma(0,0) + sigma(1,1) + sigma(2,2))/3;
    double vonMisesStress=sqrt((((sigma(0,0) - sigma(1,1))*(sigma(0,0) - sigma(1,1))
                         + (sigma(1,1) - sigma(2,2))*(sigma(1,1) - sigma(2,2))
                         + (sigma(2,2) - sigma(0,0))*(sigma(2,2) - sigma(0,0)))/2)
                         + 3*((sigma(0,1)*sigma(0,1)) + (sigma(1,2)*sigma(1,2)) + (sigma(2,0)*sigma(2,0))));
    double triax = hydro/vonMisesStress;
    if (hydro==0.0 || std::isnan(triax))
    {
        triax = 0.0;
    }
    
    if (triax<1.5)
    {
        epsFracture = (D1 + D2*exp(D3*triax))*(1.0 + D4*log(epsdot/edot0))*(1.0 + D5*tempterm);
    }
    else
    {
        epsFracture = (D1 + D2*exp(D3*1.5))*(1.0 + D4*log(epsdot/edot0))*(1.0 + D5*tempterm);
    }
    
    return epsFracture;
}




bool JCDamageModel::check() const
{
    bool ok = theJCDamageModel->check();
    
    return ok;
}




void JCDamageModel::print(std::ostream &of) const
{
    of  << "\n Johnson - Cook original damage model";

    of  << "\n   The damage variable is obtained as:"
    << "\n    Dc = Dn + (iso_c - iso_n)/epsfract"
    << "\n    (for an instantaneous triaxiality < 1.5)"
    << "\n    epsfract = (D1 + D2*exp(D3*triax))*(1 + D4*log(edot/edot0))*(1 + D5*theta)"
    << "\n    (for an instantaneous triaxiality > 1.5)"
    << "\n    epsfract = (D1 + D2*exp(D3*1.5))*(1 + D4*log(edot/edot0))*(1 + D5*theta)"
    << "\n    theta = (Temp_c - T0)/(Tm - T0)";
    
    of  <<"\n  D1                    : " << D1
    <<"\n  D2                    : " << D2
    <<"\n  D3                    : " << D3
    <<"\n  D4                    : " << D4
    <<"\n  D5                    : " << D5

    << "\n   dot{eps}_0             : " << edot0
    << "\n   T0                     : " << modelRefTemp
    << "\n   Tm                     : " << meltT;
}




void JCDamageModel::setRandom()
{
    D1 = muesli::randomUniform(0.0, 10.0);
    D2 = muesli::randomUniform(0.0, 10.0);
    D3 = muesli::randomUniform(0.0, 10.0);
    D4 = muesli::randomUniform(0.0, 10.0);
    D5 = muesli::randomUniform(0.0, 10.0);
    
    edot0 = muesli::randomUniform(1.0, 10.0);
    modelRefTemp = muesli::randomUniform(1.0, 10.0);
    meltT = muesli::randomUniform(260.0, 400.0);
}




JCCustomDamageModel::JCCustomDamageModel(const std::string& name,
                                         const materialProperties& cl)
:
edot0(1.0), modelRefTemp(1.0), meltT(1.0), D1(0.0),
D2(0.0), D3(0.0), D4(0.0), D5(0.0), D6(0.0)
{
    muesli::assignValue(cl, "edot0",       edot0);
    muesli::assignValue(cl, "modelreftemp",modelRefTemp);
    muesli::assignValue(cl, "melttemp",    meltT);
    muesli::assignValue(cl, "d1",          D1);
    muesli::assignValue(cl, "d2",          D2);
    muesli::assignValue(cl, "d3",          D3);
    muesli::assignValue(cl, "d4",          D4);
    muesli::assignValue(cl, "d5",          D5);
    muesli::assignValue(cl, "d6",          D6);
    muesli::assignValue(cl, "d7",          D7);
    muesli::assignValue(cl, "d8",          D8);
}




// Alternative material creation method, for MP individual tests
JCCustomDamageModel::JCCustomDamageModel(const std::string& name,
                                         const double x_edot0, const double x_refT,
                                         const double x_meltT,
                                         const double x_D1, const double x_D2,
                                         const double x_D3, const double x_D4,
                                         const double x_D5, const double x_D6)
:
edot0(x_edot0), modelRefTemp(x_refT), meltT(x_meltT),
D1(x_D1), D2(x_D2), D3(x_D3),
D4(x_D4), D5(x_D5), D6(x_D6)
{
}




double JCCustomDamageModel::calculateStrainToFracture(double eps, double epsdot,
                                                      istensor sigma, double temp) const
{
    double epsFracture = 0.0;
    const double tempterm  = (temp/meltT <= 1.0) ? (temp-modelRefTemp)/(meltT-modelRefTemp) : 0.0;
    double hydro = (sigma(0,0) + sigma(1,1) + sigma(2,2))/3;
    double vonMisesStress=sqrt((((sigma(0,0) - sigma(1,1))*(sigma(0,0) - sigma(1,1))
                         + (sigma(1,1) - sigma(2,2))*(sigma(1,1) - sigma(2,2))
                         + (sigma(2,2) - sigma(0,0))*(sigma(2,2) - sigma(0,0)))/2)
                         + 3*((sigma(0,1)*sigma(0,1)) + (sigma(1,2)*sigma(1,2)) + (sigma(2,0)*sigma(2,0))));
    double triax = hydro/vonMisesStress;
    if (hydro==0.0 || std::isnan(triax))
    {
        triax = 0.0;
    }
    
    if (triax<0.0)
    {
        epsFracture = (D1 + D2*exp(D3*triax))*(1.0 + D4*log(epsdot/edot0))*(1.0 + D5*tempterm);
    }
    else if (triax>=0.0 && triax<1.5)
    {
        epsFracture = (D1 + D2*exp(D3*D6*sin(D7 + triax*D8)))*(1.0 + D4*log(epsdot/edot0))*(1.0 + D5*tempterm);
    }
    else if (triax>=1.5)
    {
        epsFracture = (D1 + D2*exp(D3*D6*sin(D7 + 1.5*D8)))*(1.0+D4*log(epsdot/edot0))*(1.0 + D5*tempterm);
    }
    
    return epsFracture;
}




bool JCCustomDamageModel::check() const
{
    bool ok = theJCCustomDamageModel->check();
    
    return ok;
}




void JCCustomDamageModel::print(std::ostream &of) const
{
    of  << "\n Johnson - Cook custom damage model";

    of  << "\n   The damage variable is obtained as:"
    << "\n    Dc = Dn + (iso_c - iso_n)/epsfract"
    << "\n    (for an instantaneous triaxiality < 1.5)"
    << "\n    epsfract = (D1 + D2*exp(triax))*(1 + D4*log(edot/edot0))*(1 + D5*theta)"
    << "\n    (for an instantaneous triaxiality < 1.5)"
    << "\n    epsfract = (D1 + D2*exp(D3*D6*sin(D7 + triax*D8)))*(1 + D4*log(edot/edot0))*(1 + D5*theta)"
    << "\n    (for an instantaneous triaxiality > 1.5)"
    << "\n    epsfract = (D1 + D2*exp(D3*D6*sin(D7+1.5*D8)))*(1 + D4*log(edot/edot0))*(1 + D5*theta)"
    << "\n    theta = (Temp_c - T0)/(Tm - T0)";
    
    of  <<"\n  D1                    : " << D1
    <<"\n  D2                    : " << D2
    <<"\n  D3                    : " << D3
    <<"\n  D4                    : " << D4
    <<"\n  D5                    : " << D5
    <<"\n  D6                    : " << D6
    <<"\n  D7                    : " << D7
    <<"\n  D8                    : " << D8

    << "\n   dot{eps}_0             : " << edot0
    << "\n   T0                     : " << modelRefTemp
    << "\n   Tm                     : " << meltT;
}




void JCCustomDamageModel::setRandom()
{
    D1 = muesli::randomUniform(0.0, 10.0);
    D2 = muesli::randomUniform(0.0, 10.0);
    D3 = muesli::randomUniform(0.0, 10.0);
    D4 = muesli::randomUniform(0.0, 10.0);
    D5 = muesli::randomUniform(0.0, 10.0);
    D6 = muesli::randomUniform(0.0, 10.0);
    D7 = muesli::randomUniform(0.0, 10.0);
    D8 = muesli::randomUniform(0.0, 10.0);
    
    edot0 = muesli::randomUniform(1.0, 10.0);
    modelRefTemp = muesli::randomUniform(1.0, 10.0);
    meltT = muesli::randomUniform(260.0, 400.0);
}



JCCustom2DamageModel::JCCustom2DamageModel(const std::string& name,
                                         const materialProperties& cl)
:
edot0(1.0), modelRefTemp(1.0), meltT(1.0), D1(0.0),
D2(0.0), D3(0.0), D4(0.0), D5(0.0), D6(0.0),
D7(0.0), D8(0.0)
{
    muesli::assignValue(cl, "edot0",       edot0);
    muesli::assignValue(cl, "modelreftemp",modelRefTemp);
    muesli::assignValue(cl, "melttemp",    meltT);
    muesli::assignValue(cl, "d1",          D1);
    muesli::assignValue(cl, "d2",          D2);
    muesli::assignValue(cl, "d3",          D3);
    muesli::assignValue(cl, "d4",          D4);
    muesli::assignValue(cl, "d5",          D5);
    muesli::assignValue(cl, "d6",          D6);
    muesli::assignValue(cl, "d7",          D7);
    muesli::assignValue(cl, "d8",          D8);
}




// Alternative material creation method, for MP individual tests
JCCustom2DamageModel::JCCustom2DamageModel(const std::string& name,
                                         const double x_edot0, const double x_refT,
                                         const double x_meltT,
                                         const double x_D1, const double x_D2,
                                         const double x_D3, const double x_D4,
                                         const double x_D5, const double x_D6)
:
edot0(x_edot0), modelRefTemp(x_refT), meltT(x_meltT),
D1(x_D1), D2(x_D2), D3(x_D3),
D4(x_D4), D5(x_D5), D6(x_D6)
{
}




double JCCustom2DamageModel::calculateStrainToFracture(double eps, double epsdot,
                                                      istensor sigma, double temp) const
{
    double epsFracture = 0.0;
    const double tempterm  = (temp/meltT <= 1.0) ? (temp-modelRefTemp)/(meltT-modelRefTemp) : 0.0;
    double hydro = (sigma(0,0) + sigma(1,1) + sigma(2,2))/3;
    double vonMisesStress=sqrt((((sigma(0,0) - sigma(1,1))*(sigma(0,0) - sigma(1,1))
                         + (sigma(1,1) - sigma(2,2))*(sigma(1,1) - sigma(2,2))
                         + (sigma(2,2) - sigma(0,0))*(sigma(2,2) - sigma(0,0)))/2)
                         + 3*((sigma(0,1)*sigma(0,1)) + (sigma(1,2)*sigma(1,2)) + (sigma(2,0)*sigma(2,0))));
    double triax = hydro/vonMisesStress;
    if (hydro==0.0 || std::isnan(triax))
    {
        triax = 0.0;
    }
    
    if (triax<0.0)
    {
        epsFracture = (D1 + D2*exp(D3*triax))*(1.0 + D4*log(epsdot/edot0))*(1.0 + D5*tempterm);
    }
    else if (triax>=0.0 && triax<1.5)
    {
        epsFracture = (D6 + D7*sin(triax*D8))*(1.0 + D4*log(epsdot/edot0))*(1.0 + D5*tempterm);
    }
    else if (triax>=1.5)
    {
        epsFracture = (D6 + D7*sin(1.5*D8))*(1.0 + D4*log(epsdot/edot0))*(1.0 + D5*tempterm);
    }
    
    return epsFracture;
}




bool JCCustom2DamageModel::check() const
{
    bool ok = theJCCustom2DamageModel->check();
    
    return ok;
}




void JCCustom2DamageModel::print(std::ostream &of) const
{
    of  << "\n Johnson - Cook custom damage model";

    of  << "\n   The damage variable is obtained as:"
    << "\n    Dc = Dn + (iso_c - iso_n)/epsfract"
    << "\n    (for an instantaneous triaxiality < 1.5)"
    << "\n    epsfract = (D1 + D2*exp(triax))*(1 + D4*log(edot/edot0))*(1 + D5*theta)"
    << "\n    (for an instantaneous triaxiality < 1.5)"
    << "\n    epsfract = (D1 + D2*exp(D3*D6*sin(D7 + triax*D8)))*(1 + D4*log(edot/edot0))*(1 + D5*theta)"
    << "\n    (for an instantaneous triaxiality > 1.5)"
    << "\n    epsfract = (D1 + D2*exp(D3*D6*sin(D7+1.5*D8)))*(1 + D4*log(edot/edot0))*(1 + D5*theta)"
    << "\n    theta = (Temp_c - T0)/(Tm - T0)";
    
    of  <<"\n  D1                    : " << D1
    <<"\n  D2                    : " << D2
    <<"\n  D3                    : " << D3
    <<"\n  D4                    : " << D4
    <<"\n  D5                    : " << D5
    <<"\n  D6                    : " << D6
    <<"\n  D7                    : " << D7
    <<"\n  D8                    : " << D8

    << "\n   dot{eps}_0             : " << edot0
    << "\n   T0                     : " << modelRefTemp
    << "\n   Tm                     : " << meltT;
}




void JCCustom2DamageModel::setRandom()
{
    D1 = muesli::randomUniform(0.0, 10.0);
    D2 = muesli::randomUniform(0.0, 10.0);
    D3 = muesli::randomUniform(0.0, 10.0);
    D4 = muesli::randomUniform(0.0, 10.0);
    D5 = muesli::randomUniform(0.0, 10.0);
    D6 = muesli::randomUniform(0.0, 10.0);
    D7 = muesli::randomUniform(0.0, 10.0);
    D8 = muesli::randomUniform(0.0, 10.0);
    
    edot0 = muesli::randomUniform(1.0, 10.0);
    modelRefTemp = muesli::randomUniform(1.0, 10.0);
    meltT = muesli::randomUniform(260.0, 400.0);
}
