/****************************************************************************
*
*                                 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 <math.h>
#include <string.h>

#include "yeoh.h"
#include "muesli/material.h"
#include "muesli/Finitestrain/fisotropic.h"
#include "muesli/material.h"
#include "muesli/tensor.h"

using namespace std;
using namespace muesli;


static void computeEnergyDerivatives(const bool compressible, const ivector& invar,
                                     const double C1, const double C2, const double C3, const double k,
                                     ivector& dW, istensor& ddW);



yeohMaterial::yeohMaterial(const std::string& name,
                 const muesli::materialProperties& mp)
    :
    f_invariants(name, mp),
    _compressible(false),
    _C1(0.0), _C2(0.0), _C3(0.0), _bulk(0.0)
{
    muesli::assignValue(mp, "c1",  _C1);
    muesli::assignValue(mp, "c2",  _C2);
    muesli::assignValue(mp, "c3",  _C3);
    muesli::assignValue(mp, "bulk",_bulk);

    if ( mp.find("compressible") != mp.end() ) _compressible = true;
}




yeohMaterial::yeohMaterial(const std::string& name,
                             const double C1,
                             const double C2,
                             const double C3,
                             const double bulk,
                             const bool   compressible)
:
f_invariants(name),
_compressible(compressible),
_C1(C1), _C2(C2), _C3(C3), _bulk(bulk)
{

}




double yeohMaterial::characteristicStiffness() const
{
    return _C1;
}




muesli::finiteStrainMP* yeohMaterial::createMaterialPoint() const
{
    muesli::finiteStrainMP* mp = new yeohMP(*this);
    return mp;
}




/* this function is always called once the material is defined, so apart from
printing its information, we take the opportunity to clean up some of its
data, in particular, setting all the possible constants
*/
void yeohMaterial::print(std::ostream &of) const
{
    if (!_compressible)
    {
        of  << "\n   Elastic, incompressible Yeoh material for finite deformation analysis"
            << "\n   Stored energy function:"
            << "\n          W(I1) = C1 (I_1 - 3) + C2 (I_1-3)^2 + C3 (I1-3)^3"
            << "\n          C1      = " << _C1
            << "\n          C2      = " << _C2
            << "\n          C3      = " << _C3;
    }

    else
    {
        of  << "\n   Elastic, compressible Yeoh material for finite deformation analysis"
        << "\n   Stored energy function:"
        << "\n          W(I1,I3) = k/2 ( (J^2-1)/2 - log J) + C1 (I_1 - 3) + C2 (I_1-3)^2 + C3 (I1-3)^3"
        << "\n          C1      = " << _C1
        << "\n          C2      = " << _C2
        << "\n          C3      = " << _C3
        << "\n          k       = " << _bulk;
    }

    of  << "\n";
}




void yeohMaterial::setRandom()
{
    _C1      = muesli::randomUniform(1.0, 10.0);
    _C2      = muesli::randomUniform(1.0, 10.0);
    _C3      = muesli::randomUniform(1.0, 10.0);
    _bulk    = muesli::randomUniform(1.0, 10.0);
    _compressible = muesli::discreteUniform(0, 1) == 0;
}




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

    muesli::finiteStrainMP* p = this->createMaterialPoint();
    isok = p->testImplementation(of);
    delete p;

    return isok;
}





yeohMP::yeohMP(const yeohMaterial &m) :
    fisotropicMP(m),
    mat(&m)
{
    itensor F = itensor::identity();
    updateCurrentState(0.0, F);
    
    tn      = tc;
    Fn      = Fc;
    dW_n    = dW_c;
    ddW_n   = ddW_c;
    invar_n = invar_c;
    for (size_t a=0; a<8; a++) G_n[a] = G_c[a];
}




void yeohMP::setConvergedState(const double theTime, const itensor& F)
{
    tn = theTime;
    Fn = F;
    istensor Cn = istensor::tensorTransposedTimesTensor(Fn);
    
    invar_n[0] = Cn.invariant1();
    invar_n[1] = Cn.invariant2();
    invar_n[2] = Cn.invariant3();

    const double C1 = mat->_C1;
    const double C2 = mat->_C2;
    const double C3 = mat->_C3;
    const double k  = mat->_bulk;

    computeEnergyDerivatives(mat->_compressible, invar_n, C1, C2, C3, k,
                             dW_n, ddW_n);
    computeGammaCoefficients(invar_n, dW_n, ddW_n, G_n);
}




double yeohMP::stress(const double stretch) const
{
    return 0.0;
}




double yeohMP::stiffness(const double stretch) const
{
    return 0.0;
}




double yeohMP::storedEnergy(const double stretch) const
{
    return 0.0;
}




double yeohMP::storedEnergy() const
{
    const double I1 = invar_c[0];
    const double C1 = mat->_C1;
    const double C2 = mat->_C2;
    const double C3 = mat->_C3;
    double w;

    if ( !mat->_compressible)
    {
        const double I1m3 = I1 - 3.0;
        w = C1*I1m3 + C2*I1m3*I1m3 + C3*I1m3*I1m3*I1m3;
    }

    else
    {
        const double I3 = invar_c[2];
        const double k  = mat->_bulk;
        const double J  = sqrt(I3);
        const double I3rm13 = 1.0/std::cbrt(I3);
        const double I1t = I1 * I3rm13;
        const double I1tm3 = I1t - 3.0;
        w = C1*I1tm3 + C2*I1tm3*I1tm3 + C3*I1tm3*I1tm3*I1tm3 + 0.5*k*log(J)*log(J);
    }

    return w;
}





void yeohMP::updateCurrentState(const double theTime, const itensor& F)
{
    finiteStrainMP::updateCurrentState(theTime, F);
    istensor Cc = istensor::tensorTransposedTimesTensor(Fc);
    invar_c[0]  = Cc.invariant1();
    invar_c[1]  = Cc.invariant2();
    invar_c[2]  = Cc.invariant3();

    const double C1 = mat->_C1;
    const double C2 = mat->_C2;
    const double C3 = mat->_C3;
    const double k  = mat->_bulk;

    computeEnergyDerivatives(mat->_compressible, invar_c, C1, C2, C3, k,
                             dW_c, ddW_c);
    computeGammaCoefficients(invar_c, dW_c, ddW_c, G_c);
}



void computeEnergyDerivatives(const bool compressible, const ivector& invar,
                              const double C1, const double C2, const double C3, const double k,
                              ivector& dW, istensor& ddW)
{


    const double I1 = invar[0];

    if (compressible)
    {
        const double I3     = invar[2];
        const double J      = sqrt(I3);
        const double I3r13  = std::cbrt(I3);
        const double I3r23  = I3r13*I3r13;
        const double I3rm13 = 1.0/I3r13;
        const double I3rm43 = I3rm13*I3rm13*I3rm13*I3rm13;
        const double I1t    = I1 * I3rm13;
        const double I1tm3  = I1t - 3.0;
        const double logJ   = log(J);

        dW(0) = C1*I3rm13 + 2.0*C2*I1tm3*I3rm13 +3.0*C3*I1tm3*I1tm3*I3rm13;
        dW(1) = 0.0;
        dW(2) = -C1*I1*I3rm43/3.0 - 2.0/3.0*C2*I1*I1tm3*I3rm43
        - C3 * I1 * I1tm3 * I1tm3 * I3rm43 + 0.5*k*log(J)/I3;


        ddW.setZero();
        ddW(0,0) = (6.0*C3*I1+2.0*(C2-9.0*C3)*I3r13)/I3;
        ddW(0,2) = -(9.0*C3*I1*I1+4.0*(C2-9.0*C3)*I1*I3r13+(C1-6.0*C2+27.0*C3)*I3r23)/(3.0*I3*I3);
        ddW(2,0) = ddW(0,2);
        ddW(2,2) = (8.0*I1*(9*C3*I1*I1+5.0*(C2-9.0*C3)*I1*I3r13+2.0*(C1-6.0*C2+27*C3)*I3r23)+9.0*I3*k
                      -9.0*I3*k*2.0*logJ)/(36.0*I3*I3*I3);
    }

    else
    {
        const double I1m3 = I1 - 3.0;
        dW(0) = C1 + 2.0*C2*I1m3 + 3.0*C3*I1m3*I1m3;
        dW(1) = 0.0;
        dW(2) = 0.0;

        ddW.setZero();
        ddW(0,0) = 2.0*C2 + 6.0*C3*(I1-3.0);
    }
}



