/****************************************************************************
*
*                                 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 <stdio.h>
#include "sthermomech.h"
#include "muesli/Smallstrain/smallstrainlib.h"

using namespace muesli;


sThermoMechMaterial::sThermoMechMaterial(const std::string& name,
                                     const materialProperties& cl)
:
material(name, cl),
theSSMaterial(0),
alpha(0.0),
conductivity(0.0), heat_capacity(0.0),
taylor_quinney(1.0), tref(273.0)
{
    if       (cl.find("subtype elastic") != cl.end())      theSSMaterial = new elasticIsotropicMaterial(name, cl);
    else if  (cl.find("subtype plastic") != cl.end())      theSSMaterial = new splasticMaterial(name, cl);
    else if  (cl.find("subtype viscoelastic") != cl.end()) theSSMaterial = new viscoelasticMaterial(name, cl);
    else if  (cl.find("subtype viscoplastic") != cl.end()) theSSMaterial = new viscoplasticMaterial(name, cl);

    muesli::assignValue(cl, "alpha",   alpha);
    muesli::assignValue(cl, "conductivity", conductivity);
    muesli::assignValue(cl, "heatcapacity", heat_capacity);
    muesli::assignValue(cl, "taylor_quinney", taylor_quinney);
    muesli::assignValue(cl, "tref", tref);
}




sThermoMechMaterial::~sThermoMechMaterial()
{
    if (theSSMaterial != 0) delete(theSSMaterial);
}




bool sThermoMechMaterial::check() const
{
    return theSSMaterial->check();
}




sThermoMechMP* sThermoMechMaterial::createMaterialPoint() const
{
    return new sThermoMechMP(*this);
}




double sThermoMechMaterial::density() const
{
    return theSSMaterial->density();
}




double sThermoMechMaterial::getProperty(const propertyName p) const
{
    double ret=0.0;

    switch (p)
    {
        case PR_THERMAL_EXP:    ret = alpha;     break;
        case PR_CONDUCTIVITY:   ret = conductivity; break;
        case PR_THERMAL_CAP:    ret = heat_capacity; break;
        case PR_TAYLOR_QUINNEY: ret = taylor_quinney; break;

        default:
            ret = theSSMaterial->getProperty(p);
    }
    return ret;
}




void sThermoMechMaterial::print(std::ostream &of) const
{
    of  << "\n Small strain sThermoMech material."
        << "\n Surrogate mechanical model: ";
    theSSMaterial->print(of);

    of  << "\n Reference temperature:                 " << tref
        << "\n Thermal expansion coefficient:         " << alpha
        << "\n Thermal conductivity:                  " << conductivity
        << "\n Heat capacity @ const. vol/unit mass:  " << heat_capacity
        << "\n Taylor-Quinney parameter:              " << taylor_quinney;
}




void sThermoMechMaterial::setRandom()
{
    int mattype = discreteUniform(1, 1);
    std::string name = "surrogate small strain material";
    materialProperties mp;
    
    if (mattype == 0)
        theSSMaterial = new elasticIsotropicMaterial(name, mp);
    
    else if (mattype == 1)
        theSSMaterial = new splasticMaterial(name, mp);
    
    else if (mattype == 2)
        theSSMaterial = new viscoelasticMaterial(name, mp);
    
    else if (mattype == 3)
        theSSMaterial = new viscoplasticMaterial(name, mp);
    
    else if (mattype == 4)
        theSSMaterial = new sdamageMaterial(name, mp);
    
    theSSMaterial->setRandom();
    alpha          = randomUniform(1.0, 2.0);
    tref           = randomUniform(100.0, 400.0);
    heat_capacity  = randomUniform(1.0, 2.0);
    taylor_quinney = randomUniform(0.80, 1.0);
}




bool sThermoMechMaterial::test(std::ostream& of)
{
    bool isok = true;
    setRandom();

    sThermoMechMP* p = this->createMaterialPoint();

    isok = p->testImplementation(of);
    delete p;

    return isok;
}




double sThermoMechMaterial::waveVelocity() const
{
    return theSSMaterial->waveVelocity();
}



sThermoMechMP::sThermoMechMP(const sThermoMechMaterial &m)
:
dtemp_n(0.0), dtemp_c(0.0),
theSSMP(0),
theThermoMechMaterial(m)
{
    theSSMP = m.theSSMaterial->createMaterialPoint();
    gradt_n.setZero();
    gradt_c.setZero();
}




sThermoMechMP::~sThermoMechMP()
{
    if (theSSMP != 0) delete (theSSMP);
}




void sThermoMechMP::commitCurrentState()
{
    dtemp_n = dtemp_c;
    gradt_n = gradt_c;
    theSSMP->commitCurrentState();
}




double sThermoMechMP::contractWithConductivity(const ivector &v1, const ivector &v2) const
{
    return theThermoMechMaterial.conductivity * v1.dot(v2);
}




void sThermoMechMP::contractWithTangent(const ivector &v1, const ivector &v2, itensor &T) const
{
    theSSMP->contractWithTangent(v1, v2, T);
}




void sThermoMechMP::contractWithDeviatoricTangent(const ivector &v1, const ivector &v2, itensor &T) const
{
    theSSMP->contractWithDeviatoricTangent(v1, v2, T);
}




double sThermoMechMP::deviatoricEnergy() const
{
    return theSSMP->deviatoricEnergy();
}




void sThermoMechMP::deviatoricStress(istensor& s) const
{
    theSSMP->deviatoricStress(s);
}




double sThermoMechMP::dissipatedEnergy() const
{
    // mechanical part
    double   dissMech   = theSSMP->dissipatedEnergy();
    
    // heat dissipation -> from coupled part of stress
    const istensor ep_c  = theSSMP->getCurrentPlasticStrain();
    const istensor ep_n  = theSSMP->getConvergedPlasticStrain();
    const istensor eps_c = theSSMP->getCurrentState().theStensor[0];
    const istensor eps_n = theSSMP->getConvergedState().theStensor[0];
    istensor M;
    stressTemperatureTensor(M);

    istensor S = dtemp_c * M; 
    double dissHeat = S.dot(ep_c - ep_n);
    
    return dissMech + dissHeat;
}




// II - Cinverse:Cep for small sThermoMech element
void sThermoMechMP::dissipationTangent(itensor4& D) const
{
    theSSMP -> dissipationTangent(D);
}




double sThermoMechMP::effectiveStoredEnergy() const
{
    const double alpha      = theThermoMechMaterial.alpha;
    const double cv         = theThermoMechMaterial.heat_capacity;
    const double kappa      = theSSMP->volumetricStiffness();
    const double tref       = theThermoMechMaterial.tref;
    const double vol        = theSSMP->getCurrentStrain().trace();
    
    double       dt_Psistar = dissipatedEnergy();
    double       fmec       = theSSMP->storedEnergy();
    double       fcou       = -3.0 * kappa * alpha * dtemp_c * vol;
    double       fthe       = -0.5 * cv/tref * dtemp_c * dtemp_c;
    
    return fmec + fcou + fthe + dt_Psistar;
}




double sThermoMechMP::entropy() const
{
    const double    cv      = theThermoMechMaterial.heat_capacity;
    const double    kappa   = theSSMP->volumetricStiffness();
    const double    beta    = 3.0 * kappa * theThermoMechMaterial.alpha;
    const double    tref    = theThermoMechMaterial.tref;

    materialState   state_c = theSSMP->getCurrentState();
    const istensor& eps_c   = state_c.theStensor[0];

    return  beta * eps_c.trace() + cv/tref * dtemp_c;
}





double sThermoMechMP::freeEnergy() const
{
    const double alpha  = theThermoMechMaterial.alpha;
    const double cv     = theThermoMechMaterial.heat_capacity;
    const double kappa  = theSSMP->volumetricStiffness();
    const double tref   = theThermoMechMaterial.tref;
    const double vol    = theSSMP->getCurrentStrain().trace();

    double       fmec   = theSSMP->storedEnergy();
    double       fcou   = -3.0 * kappa * alpha * dtemp_c * vol;
    double       fthe   = -0.5 * cv/tref * dtemp_c * dtemp_c;

    return fmec + fcou + fthe;
}




istensor sThermoMechMP::getConvergedPlasticStrain() const
{
    return theSSMP->getConvergedPlasticStrain();
}




materialState sThermoMechMP::getConvergedState() const
{
    // first store the state of the ssmp
    materialState mat = theSSMP->getConvergedState();

    // then append the thermal data
    mat.theDouble.push_back(dtemp_n);
    mat.theVector.push_back(gradt_n);

    return mat;
}




istensor sThermoMechMP::getCurrentPlasticStrain() const
{
    return theSSMP->getCurrentPlasticStrain();
}




materialState sThermoMechMP::getCurrentState() const
{
    // first store the state of the ssmp
    materialState mat = theSSMP->getCurrentState();

    // then append the thermal data
    mat.theDouble.push_back(dtemp_c);
    mat.theVector.push_back(gradt_c);

    return mat;
}




istensor& sThermoMechMP::getCurrentStrain()
{
    return theSSMP->getCurrentState().theStensor[0];
}




ivector sThermoMechMP::heatflux() const
{
    return -theThermoMechMaterial.conductivity * gradt_c;
}




double sThermoMechMP::plasticSlip() const
{
    return theSSMP->plasticSlip();
}




double sThermoMechMP::pressure() const
{
    return theSSMP->pressure();
}




void  sThermoMechMP::resetCurrentState()
{
    dtemp_c = dtemp_n;
    gradt_c = gradt_n;
    theSSMP->resetCurrentState();
}




void sThermoMechMP::setRandom()
{
    dtemp_c = randomUniform(0.06, 0.13) * theThermoMechMaterial.tref;
    gradt_c.setRandom();
    theSSMP->setRandom();
}




void sThermoMechMP::stress(istensor& sigma) const
{
    const double alpha = theThermoMechMaterial.alpha;
    const double kappa = theSSMP->volumetricStiffness();

    theSSMP->stress(sigma);
    sigma -= 3.0 * kappa * alpha * dtemp_c * istensor::identity();
}




void  sThermoMechMP::stressTemperatureTensor(istensor& M) const
{
    const double alpha = theThermoMechMaterial.alpha;
    const double kappa = theSSMP->volumetricStiffness();
    M = -3.0 * kappa * alpha * istensor::identity();
}




void sThermoMechMP::tangentElasticities(itensor4& c) const
{
    theSSMP->tangentTensor(c);
}




double& sThermoMechMP::temperature()
{
    return dtemp_c;
}




const double& sThermoMechMP::temperature() const
{
    return dtemp_c;
}




bool sThermoMechMP::testImplementation(std::ostream& os) const
{
    bool isok = true;
    
    // set a random update in the material
    istensor eps;
    eps.setZero();

    double temp;
    temp = randomUniform(0.7, 1.3) * theThermoMechMaterial.tref;

    ivector gradt;
    gradt.setRandom();

    sThermoMechMP* ssmp = const_cast<sThermoMechMP*>(this);
    ssmp->updateCurrentState(0.0, eps, temp, gradt);
    ssmp->commitCurrentState();

    double tn1 = muesli::randomUniform(0.1,1.0);
    eps.setRandom();
    sThermoMechMP& theMP = const_cast<sThermoMechMP&>(*this);
    theMP.updateCurrentState(tn1, eps, temp, gradt);

    // programmed tangent of elasticities
    itensor4 tg;
    tangentElasticities(tg);

    // (1)  compare DEnergy with the derivative of Energy
    if (true)
    {
        // programmed stress
        istensor sigma;
        this->stress(sigma);

        // numerical differentiation stress
        istensor numSigma;
        numSigma.setZero();
        const double   inc = 1.0e-3;

        for (size_t i=0; i<3; i++)
        {
            for (size_t j=i; j<3; j++)
            {
                double original = eps(i,j);

                eps(i,j) = eps(j,i) = original + inc;
                theMP.updateCurrentState(tn1, eps, dtemp_c, gradt_c);
                double Wp1 = effectiveStoredEnergy();

                eps(i,j) = eps(j,i) = original + 2.0*inc;
                theMP.updateCurrentState(tn1, eps, dtemp_c, gradt_c);
                double Wp2 = effectiveStoredEnergy();

                eps(i,j) = eps(j,i) = original - inc;
                theMP.updateCurrentState(tn1, eps, dtemp_c, gradt_c);
                double Wm1 = effectiveStoredEnergy();

                eps(i,j) = eps(j,i) = original - 2.0*inc;
                theMP.updateCurrentState(tn1, eps, dtemp_c, gradt_c);
                double Wm2 = effectiveStoredEnergy();

                // fourth order approximation of the derivative
                double der = (-Wp2 + 8.0*Wp1 - 8.0*Wm1 + Wm2)/(12.0*inc);
                numSigma(i,j) = der;
                if (i != j) numSigma(i,j) *= 0.5;

                numSigma(j,i) = numSigma(i,j);

                eps(i,j) = eps(j,i) = original;
                theMP.updateCurrentState(tn1, eps, dtemp_c, gradt_c);
            }
        }

        // relative error less than 0.01%
        istensor error = numSigma - sigma;
        isok = (error.norm()/sigma.norm() < 1e-4);

        os << "\n   1. Comparing stress with DWeff.";
        if (isok)
        {
            os << " Test passed.";
        }
        else
        {
            os << "\n      Test failed.";
            os << "\n      Relative error in DWeff computation: " <<  error.norm()/sigma.norm();
            os << "\n      Stress:\n" << sigma;
            os << "\n      Numeric stress:\n" << numSigma;
        }
    }


    // (2) compare tensor c with derivative of stress
    if (true)
    {
        // numeric C
        itensor4 nC;
        nC.setZero();

        // numerical differentiation sigma
        istensor dsigma, sigmap1, sigmap2, sigmam1, sigmam2;
        double   inc = 1.0e-3;

        for (unsigned i=0; i<3; i++)
        {
            for (unsigned j=i; j<3; j++)
            {
                double original = eps(i,j);

                eps(i,j) = eps(j,i) = original + inc;
                theMP.updateCurrentState(tn1, eps, dtemp_c, gradt_c);
                stress(sigmap1);

                eps(i,j) = eps(j,i) = original + 2.0*inc;
                theMP.updateCurrentState(tn1, eps, dtemp_c, gradt_c);
                stress(sigmap2);

                eps(i,j) = eps(j,i) = original - inc;
                theMP.updateCurrentState(tn1, eps, dtemp_c, gradt_c);
                stress(sigmam1);

                eps(i,j) = eps(j,i) = original - 2.0*inc;
                theMP.updateCurrentState(tn1, eps, dtemp_c, gradt_c);
                stress(sigmam2);

                // fourth order approximation of the derivative
                dsigma = (-sigmap2 + 8.0*sigmap1 - 8.0*sigmam1 + sigmam2)/(12.0*inc);

                if (i != j) dsigma *= 0.5;


                for (unsigned k=0; k<3; k++)
                    for (unsigned l=0; l<3; l++)
                    {
                        nC(k,l,i,j) = dsigma(k,l);
                        nC(k,l,j,i) = dsigma(k,l);
                    }

                eps(i,j) = original;
                eps(j,i) = original;
                theMP.updateCurrentState(tn1, eps, dtemp_c, gradt_c);
            }
        }

        // relative error less than 0.01%
        double error = 0.0;
        double norm = 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++)
                    {
                        error += pow(nC(i,j,k,l)-tg(i,j,k,l),2);
                        norm  += pow(tg(i,j,k,l),2);
                    }
        error = sqrt(error);
        norm = sqrt(norm);
        isok = (error/norm < 1e-4);

        os << "\n   2. Comparing tensor C with DStress.";
        if (isok)
        {
            os << " Test passed.";
        }
        else
        {
            os << "\n      Test failed.";
            os << "\n      Relative error in DWeff computation: " <<  error/norm;
        }
    }


    // (3) compare stress and voigt stress


    // (4) compare contract tangent with cijkl
    if (true)
    {
        itensor Tvw;
        ivector v, w;
        v.setRandom();
        w.setRandom();
        contractWithTangent(v, w, Tvw);

        istensor S; S.setRandom();

        itensor nTvw; nTvw.setZero();
        for (unsigned i=0; i<3; i++)
            for (unsigned j=0; j<3; j++)
                for (unsigned k=0; k<3; k++)
                    for (unsigned l=0; l<3; l++)
                    {
                        nTvw(i,k) += tg(i,j,k,l)*v(j)*w(l);
                    }

        // relative error less than 0.01%
        itensor error = Tvw - nTvw;
        isok = (error.norm()/Tvw.norm() < 1e-4);

        os << "\n   4. Comparing contract tangent with C_ijkl v_j w_l.";
        if (isok)
        {
            os << " Test passed.";
        }
        else
        {
            os << "\n      Test failed.";
            os << "\n      Relative error: " << error.norm()/Tvw.norm();
            os << "\n      C{a,b} \n" << Tvw;
            os << "\n      C_ijkl a_j b_l:\n" << nTvw;
        }
    }


    // (5) compare volumetric and deviatoric tangent contraction


    // (6) compare tangent with volumetric+deviatoric


    return isok;
}




void sThermoMechMP::updateCurrentState(const double t, const istensor& strain, const double temp, const ivector& gradt)
{
    dtemp_c  = temp;
    gradt_c  = gradt;
    theSSMP->updateCurrentState(t, strain);
}




double sThermoMechMP::volumetricEnergy() const
{
    return theSSMP->volumetricEnergy();
}




double sThermoMechMP::volumetricStiffness() const
{
    return theSSMP->volumetricStiffness();
}

