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

using namespace muesli;


smallStrainMaterial::smallStrainMaterial(const std::string& name) :
material(name),
theFailureCriterium(0)
{
}




smallStrainMaterial::smallStrainMaterial(const std::string& name,
                                           const materialProperties& cl)
:
material(name, cl),
theFailureCriterium(0)
{}




void smallStrainMaterial::print(std::ostream &of) const
{
    if (theFailureCriterium != 0)
    {
        theFailureCriterium->print(of);
    }
}




void smallStrainMaterial::setFailureCriterium(const ssFailureCriterium& fc)
{
    theFailureCriterium = &fc;
}




smallStrainMP::smallStrainMP(const smallStrainMaterial &m)
:
theSmallStrainMaterial(m),
theFailurePoint(0),
time_n(0.0), time_c(0.0)
{
    eps_n.setZero();
    eps_c.setZero();

    if (m.theFailureCriterium != 0) theFailurePoint = m.theFailureCriterium->createFailurePoint();
}




smallStrainMP::~smallStrainMP()
{
    if (theFailurePoint != 0) delete theFailurePoint;
}




void smallStrainMP::commitCurrentState()
{
    time_n = time_c;
    eps_n  = eps_c;

    if (theFailurePoint != 0) theFailurePoint->commitCurrentState();
}




void smallStrainMP::contractWithTangent(const ivector &v1, const ivector &v2, itensor &T) const
{
    itensor4 C;
    tangentTensor(C);

    T.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++)
                    T(i,k) += C(i,j,k,l)*v1[j]*v2[l];
}




/* let Pdev the 4th order tensor that projects sym tensors to deviatoric sym tensors
 Compute Cdev = Pdev^T C Pdev, then
 compute the contraction (Cdev)_{ijkl} v_j w_l
 */
void smallStrainMP::contractWithDeviatoricTangent(const ivector &v, const ivector &w, itensor &T) const
{
    itensor4 C;
    tangentTensor(C);

    double Cpp[3][3];
    for (unsigned i=0; i<3; i++)
        for (unsigned j=0; j<3; j++)
        {
            Cpp[i][j] = 0.0;
            for (unsigned p=0; p<3; p++)
                Cpp[i][j] += C(p,p,i,j);
        }

    double Cppqq=0.0;
    for (unsigned q=0; q<3; q++)
    {
        Cppqq += Cpp[q][q];
    }


    T.setZero();
    const double th = 1.0/3.0;
    for (unsigned i=0; i<3; i++)
    {
        for (unsigned k=0; k<3; k++)
        {
            for (unsigned j=0; j<3; j++)
            {
                for (unsigned l=0; l<3; l++)
                {
                    T(i,k) += C(i,j,k,l)*v[j]*w[l];
                }

                T(i,k) -= th*Cpp[k][j]*v[i]*w[j]
                        + th*Cpp[i][j]*v[j]*w[k];
            }
            T(i,k) += th*th*Cppqq *v[i]*w[k];
        }
    }
}




// CM_ij = 1.0/3.0 * P_ijkl*c_klmm
void smallStrainMP::contractWithMixedTangent(istensor& CM) const
{
    const istensor id = istensor::identity();
    itensor4 Pdev;
    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++)
                    Pdev(i,j,k,l) = 0.5 * ( id(i,k)*id(j,l) + id(i,l)*id(j,k) ) - 1.0/3.0 * id(i,j)*id(k,l);
    
    itensor4 st;
    tangentTensor(st);

    CM.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++)
                    for (unsigned m=0; m<3; m++)
                            CM(i,j) += 1.0/3.0*Pdev(i,j,k,l)*st(k,l,m,m);

}




// II - Cinverse:Cep for small thermomechanical element
void smallStrainMP::dissipationTangent(itensor4 &D) const
{
    const double kappa = volumetricStiffness();
    const double mu    = theSmallStrainMaterial.getProperty(muesli::PR_MU);
    
    const itensor4 idSymm = itensor4::identitySymm();
    const itensor4 Pdev   = itensor4::deviatoricIdentity();
    
    // First step: obtain fourth order consistent tangent tensor, Cep
    itensor4 Cep;
    tangentTensor(Cep);
    
    // Second step: build the inverse of the fourth order elastic tensor, Cinverse
    itensor4 ioi;
    ioi.setZero();
    ioi.addDyadic(istensor::identity(), istensor::identity());

    itensor4 Cinverse = 1.0/(9.0 * kappa) * ioi + 1.0/(2.0 * mu) * Pdev;
    
    // Third step: compute D = II - Cinverse*Cep
    itensor4 aux = Cinverse*Cep;
    D.setZero();
    D = idSymm - aux;
}




void smallStrainMP::deviatoricStress(istensor& s) const
{
    istensor sigma;
    this->stress(sigma);
    s = istensor::deviatoricPart(sigma);
}




double smallStrainMP::effectiveStoredEnergy() const
{
    return storedEnergy() + dissipatedEnergy();
}




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




materialState smallStrainMP::getConvergedState() const
{
    materialState state_n;

    state_n.theTime = time_n;
    state_n.theStensor.push_back(eps_n);

    return state_n;
}




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




materialState smallStrainMP::getCurrentState() const
{
    materialState state_c;

    state_c.theTime = time_c;
    state_c.theStensor.push_back(eps_c);

    return state_c;
}




istensor& smallStrainMP::getCurrentStrain()
{
    return eps_c;
}



const istensor& smallStrainMP::getCurrentStrain() const
{
    return eps_c;
}




double smallStrainMP::logLife() const
{
    return (theFailurePoint == 0) ? 0.0 : theFailurePoint->remainingLife();
}




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




void smallStrainMP::resetCurrentState()
{
    time_c = time_n;
    eps_c = eps_n;

    if (theFailurePoint != 0) theFailurePoint->resetCurrentState();
}




void smallStrainMP::setRandom()
{
    time_c = muesli::randomUniform(0.01, 0.1);
    eps_c.setRandom();
}




void smallStrainMP::tangentMatrix(double C[6][6]) const
{
    itensor4 tg;
    tangentTensor(tg);
    muesli::tensorToMatrix(tg, C);
}




bool smallStrainMP::testImplementation(std::ostream& os, const bool testDE, const bool testDDE) const
{
    bool isok = true;
    smallStrainMP& theMP = const_cast<smallStrainMP&>(*this);

    // set a random update in the material
    istensor eps;
    eps.setZero();

    theMP.updateCurrentState(0.0, eps);
    theMP.commitCurrentState();
    
    double tn1 = muesli::randomUniform(0.1,1.0);
    eps.setRandom();
    eps *= 0.01;

    theMP.updateCurrentState(tn1, eps);
    
    // programmed tangent of elasticities
    itensor4 tg;
    tangentTensor(tg);


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

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

        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);
                double Wp1 = effectiveStoredEnergy();

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

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

                eps(i,j) = eps(j,i) = original - 2.0*inc;
                theMP.updateCurrentState(tn1, eps);
                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);
            }
        }

        // 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 (testDDE)
    {
        // numeric C
        itensor4 nC;
        nC.setZero();

        // numerical differentiation sigma
        istensor dsigma, sigmap1, sigmap2, sigmam1, sigmam2;
        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);
                stress(sigmap1);

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

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

                eps(i,j) = eps(j,i) = original - 2.0*inc;
                theMP.updateCurrentState(tn1, eps);
                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) = nC(k,l,j,i) = dsigma(k,l);

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

        // 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;
        }
    }



    // (4) compare contract tangent with cijkl
    if (testDDE)
    {
        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-6) compare volumetric and deviatoric tangent contraction
    if (testDDE)
    {
        double kappa = volumetricStiffness();
        double nkappa = 0.0;

        for (unsigned i=0; i<3; i++)
            for (unsigned j=0; j<3; j++)
                nkappa += 1.0/9.0 * tg(i,i,j,j);

        // relative error less than 0.01%
        double error = kappa - nkappa;
        isok = (error/kappa < 1e-4);

        os << "\n   5. Comparing volumetric tangent with (1/9) C_iijj.";
        if (isok)
        {
            os << " Test passed.";
        }
        else
        {
            os << "\n      Test failed.";
            os << "\n      Relative error: " << error/kappa;
            os << "\n      kappa: " << kappa;
            os << "\n      (1/9)C_iijj: " << nkappa;
        }


        ivector v, w;
        v.setRandom();
        w.setRandom();

        itensor Tvw;
        contractWithDeviatoricTangent(v, w, Tvw);

        itensor Dvw;
        smallStrainMP::contractWithDeviatoricTangent(v, w, Dvw);


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

        os << "\n   6. Comparing specific and generic dev tangent contraction";
        if (isok)
        {
            os << " Test passed.";
        }
        else
        {
            os << "\n      Test failed.";
            os << "\n      Relative error: " << Terror.norm()/Tvw.norm();
            os << "\n      C{a,b} \n" << Tvw;
            os << "\n      C_ijkl a_j b_l:\n" << Dvw;
        }
    }


    // (7) compare tangent with volumetric+deviatoric+mixed
    if (testDDE)
    {
        const double kappa = volumetricStiffness();
        
        itensor Tvw, Tdevvw, Tvolvw;
        ivector v, w;
        v.setRandom();
        w.setRandom();
        contractWithTangent(v, w, Tvw);
        contractWithDeviatoricTangent(v, w, Tdevvw);
        Tvolvw = kappa * itensor::dyadic(v, w);

        istensor CM;
        contractWithMixedTangent(CM);
        itensor Tmixed =  itensor::dyadic(v, CM * w) + itensor::dyadic(CM * v, w);

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

        os << "\n   7. Comparing tangent with deviatoric/volumetric/mixed split.";
        if (isok)
        {
            os << " Test passed.";
        }
        else
        {
            os << "\n      Test failed.";
            os << "\n      Relative error: " << error.norm()/Tvw.norm();
            os << "\n      Tvw: " << Tvw;
            os << "\n      Tvw_split" << Tdevvw+Tvolvw;
        }
    }


    // (8) compare voigt matrix and elasticity tensor
    if (testDDE)
    {
        double cm[6][6];
        tangentMatrix(cm);
        //extern const unsigned voigt[2][6];


        double error = 0.0;
        double norm  = 0.0;
        for (unsigned i=0; i<6; i++)
        {
            error += pow(cm[i][i] - tg( voigt(0,i), voigt(1,i), voigt(0,i), voigt(1,i)) ,2);
            norm  += pow(cm[i][i],2);
            for (unsigned j=i+1; j<6; j++)
            {
                error += pow(cm[i][j] - tg( voigt(0,i), voigt(1,i), voigt(0,j), voigt(1,j) ),2);
                error += pow(cm[j][i] - tg( voigt(0,j), voigt(1,j), voigt(0,i), voigt(1,i) ),2);
                norm  += pow(cm[i][j],2);
                norm  += pow(cm[j][i],2);
            }
        }
        error = sqrt(error);
        norm  = sqrt(norm);


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

        os << "\n   8. Comparing tangent tensor and Voigt matrix.";
        if (isok)
        {
            os << " Test passed.";
        }
        else
        {
            os << "\n      Test failed.";
            os << "\n      Relative error: " << error/norm;
        }
    }


    return isok;
}





void smallStrainMP::updateCurrentState(const double theTime, const istensor& strain)
{
    time_c = theTime;
    eps_c  = strain;

    if (theFailurePoint != 0)
    {
        istensor sigma;
        this->stress(sigma);
        theFailurePoint->updateCurrentState(theTime, strain, sigma);
    }
}




double smallStrainMP::volumetricStiffness() const
{
    double kappa = 0.0;
    const double i19 = 1.0/9.0;
    
    itensor4 C;
    tangentTensor(C);
    
    for (unsigned i=0; i<3; i++)
        for (unsigned j=0; j<3; j++)
            kappa += i19 * C(i,i,j,j);
    
    return kappa;
}
