/****************************************************************************
*
*                                 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 "arrudaboyce.h"
#include "muesli/material.h"
#include "muesli/Finitestrain/fisotropic.h"
#include "muesli/material.h"
#include "muesli/tensor.h"

using namespace std;
using namespace muesli;


arrudaboyceMaterial::arrudaboyceMaterial(const std::string& name,
                 const muesli::materialProperties& mp)
    :
    f_invariants(name, mp),
    _compressible(false),
    _C1(0.0), _lambdam(0.0), _bulk(0.0)
{
    muesli::assignValue(mp, "c1",  _C1);
    muesli::assignValue(mp, "lambdam",  _lambdam);
    muesli::assignValue(mp, "bulk",  _bulk);

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




arrudaboyceMaterial::arrudaboyceMaterial(const std::string& name,
                                           const double C1,
                                           const double lambdam,
                                           const double bulk,
                                           const bool   compressible)
:
f_invariants(name),
_compressible(compressible),
_C1(C1), _lambdam(lambdam), _bulk(bulk)
{

}




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




muesli::finiteStrainMP* arrudaboyceMaterial::createMaterialPoint() const
{
    muesli::finiteStrainMP* mp = new arrudaboyceMP(*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 arrudaboyceMaterial::print(std::ostream &of) const
{
    if (!_compressible)
    {
        of  << "\n   Elastic, incompressible Arruda-Boyce material for finite deformation analysis"
            << "\n   Stored energy function:"
            << "\n          W(I1) = C1 Sum_i=1^5 alpha_i beta^(i-1) (I_1^i - 3^i)"
            << "\n          beta = lambdam^(-2)"
            << "\n          alpha_1 = 1/2, alpha_2 = 1/20, alpha_3 = 11/1050"
            << "\n          alpha_4 = 19/7000, alpha_5 = 519/673750"
            << "\n          C1      = " << _C1
            << "\n          lambdam = " << _lambdam;
    }

    else
    {
        of  << "\n   Elastic, compressible Arruda-Boyce material for finite deformation analysis"
        << "\n   Stored energy function:"
        << "\n          W(I1,I3) = k/2 ( (J^2-1)/2 - log J) + C1 Sum_i=1^5 alpha_i beta^(i-1) (Itilde_1^i - 3^i)"
        << "\n          beta = lambdam^(-2)"
        << "\n          alpha_1 = 1/2, alpha_2 = 1/20, alpha_3 = 11/1050"
        << "\n          alpha_4 = 19/7000, alpha_5 = 519/673750"
        << "\n          C1      = " << _C1
        << "\n          lambdam = " << _lambdam
        << "\n          k       = " << _bulk;
    }

    of  << "\n";
}




void arrudaboyceMaterial::setRandom()
{
    _C1      = muesli::randomUniform(1.0, 10.0);
    _lambdam = muesli::randomUniform(0.1, 1.0);
    _bulk    = muesli::randomUniform(1.0, 10.0);
    _compressible = muesli::discreteUniform(0, 1) == 0;
}




bool arrudaboyceMaterial::test(std::ostream &of)
{
    bool isok = true;
    setRandom();
    muesli::finiteStrainMP* p = this->createMaterialPoint();
    
    isok = p->testImplementation(of);
    delete p;

    return isok;
}





arrudaboyceMP::arrudaboyceMP(const arrudaboyceMaterial &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 arrudaboyceMP::setConvergedState(const double theTime, const itensor& F)
{
    tn          = theTime;
    Fn          = F;
    istensor Cn = istensor::tensorTransposedTimesTensor(Fn);
    tn          = theTime;

    invar_n[0]   = Cn.invariant1();
    invar_n[1]   = Cn.invariant2();
    invar_n[2]   = Cn.invariant3();
    double I1   = invar_n[0];
    double I3   = invar_n[2];
    double J    = sqrt(I3);

    const double alpha[] = {0.5, 1.0/20.0, 11.0/1050.0, 19.0/7000.0, 519.0/673750.0};
    const double beta    = 1.0/(mat->_lambdam * mat->_lambdam);
    const double C1      = mat->_C1;
    const double k       = mat->_bulk;


    if (mat->_compressible)
    {
        const double Jm13 = std::pow(J,-1.0/3.0);
        const double Jm23 = Jm13*Jm13;

        dW_n.setZero();
        for (unsigned i=1; i<=5; i++)
        {
            dW_n(0) += C1 * alpha[i-1] * i * std::pow(beta*I1*Jm23, i-1) * Jm23;
            dW_n(2) -= C1 * alpha[i-1] * i * std::pow(beta*I1*Jm23, i-1) * Jm23 / (3.0*I3) * I1;
        }
        dW_n(2) += 0.25* mat->bulk * (1.0-1.0/I3);


        ddW_n.setZero();
        for (unsigned i=2; i<=5; i++)
        {
            ddW_n(0,0) += C1 * alpha[i-1] * i * (i-1.0) * std::pow(beta*I1*Jm23, i-1) * Jm23 / I1;
        }
        for (unsigned i=1; i<=5; i++)
        {
            ddW_n(0,2) -= C1 * alpha[i-1] * i * i * std::pow(beta*I1*Jm23, i-1) * Jm23 / (3.0*I3);
        }
        ddW_n(2,0) = ddW_n(0,2);
        for (unsigned i=1; i<=5; i++)
        {
            ddW_n(2,2) += C1 * alpha[i-1] * (i+3.0) * i * std::pow(beta*I1*Jm23, i-1) * I1 * Jm23 / (9.0*I3*I3);
        }
        ddW_n(2,2) += 0.25*k/(I3*I3);
    }

    else
    {
        dW_n(0)    = 0.0;
        ddW_n(0,0) = 0.0;
        for (unsigned i=1; i<=5; i++)
        {
            dW_n(0) += C1 * alpha[i-1] * std::pow(beta, i-1) * i * std::pow(I1, i-1);
        }

        for (unsigned i=2; i<=5; i++)
        {
            ddW_n(0,0) += C1 * alpha[i-1] * std::pow(beta, i-1) * i * (i-1.0) * std::pow(I1, i-2);
        }
    }

    computeGammaCoefficients(invar_n, dW_n, ddW_n, G_n);
}




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




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




// U = lambda/2 (log (stretch) )^2 + mu/2 (stretch^2-3) -mu log(stretch)
double arrudaboyceMP::storedEnergy(const double stretch) const
{
    return 0.0;
}




double arrudaboyceMP::storedEnergy() const
{
    const double alpha[] = {0.5, 1.0/20.0, 11.0/1050.0, 19.0/7000.0, 519.0/673750.0};
    const double beta    = 1.0/(mat->_lambdam * mat->_lambdam);
    const double C1      = mat->_C1;
    const double k       = mat->_bulk;

    const double I1 = invar_c[0];
    const double I3 = invar_c[2];
    const double J  = sqrt(I3);
    const double Jm13 = std::pow(J,-1.0/3.0);
    const double Jm23 = Jm13*Jm13;

    double w  = 0.0;


    if (mat->_compressible)
    {
        w  = 0.5 * k * ( 0.5*(I3-1.0) - log(J) );
        for (unsigned i=1; i<=5; i++)
            w += C1 * alpha[i-1] * std::pow(beta, i-1) * ( std::pow(I1*Jm23, i) - std::pow(3.0, i) );
    }
    else
    {
        for (unsigned i=1; i<=5; i++)
            w += C1 * alpha[i-1] * std::pow(beta, i-1) * ( std::pow(I1, i) - std::pow(3.0, i) );
    }

    return w;
}




void arrudaboyceMP::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();
    double I1   = invar_c[0];
    double I3   = invar_c[2];
    double J    = sqrt(I3);

    const double alpha[] = {0.5, 1.0/20.0, 11.0/1050.0, 19.0/7000.0, 519.0/673750.0};
    const double beta    = 1.0/(mat->_lambdam * mat->_lambdam);
    const double C1      = mat->_C1;
    const double k       = mat->_bulk;


    if (mat->_compressible)
    {
        const double Jm13 = std::pow(J,-1.0/3.0);
        const double Jm23 = Jm13*Jm13;
        
        dW_c.setZero();
        for (unsigned i=1; i<=5; i++)
        {
            dW_c(0) += C1 * alpha[i-1] * i * std::pow(beta*I1*Jm23, i-1) * Jm23;
            dW_c(2) -= C1 * alpha[i-1] * i * std::pow(beta*I1*Jm23, i-1) * Jm23 / (3.0*I3) * I1;
        }
        dW_c(2) += 0.25* mat->bulk * (1.0-1.0/I3);


        ddW_c.setZero();
        for (unsigned i=2; i<=5; i++)
        {
            ddW_c(0,0) += C1 * alpha[i-1] * i * (i-1.0) * std::pow(beta*I1*Jm23, i-1) * Jm23 / I1;
        }
        for (unsigned i=1; i<=5; i++)
        {
            ddW_c(0,2) -= C1 * alpha[i-1] * i * i * std::pow(beta*I1*Jm23, i-1) * Jm23 / (3.0*I3);
        }
        ddW_c(2,0) = ddW_c(0,2);
        for (unsigned i=1; i<=5; i++)
        {
            ddW_c(2,2) += C1 * alpha[i-1] * (i+3.0) * i * std::pow(beta*I1*Jm23, i-1) * I1 * Jm23 / (9.0*I3*I3);
        }
        ddW_c(2,2) += 0.25*k/(I3*I3);
    }

    else
    {
        dW_c(0)    = 0.0;
        ddW_c(0,0) = 0.0;
        for (unsigned i=1; i<=5; i++)
        {
            dW_c(0) += C1 * alpha[i-1] * std::pow(beta, i-1) * i * std::pow(I1, i-1);
        }

        for (unsigned i=2; i<=5; i++)
        {
            ddW_c(0,0) += C1 * alpha[i-1] * std::pow(beta, i-1) * i * (i-1.0) * std::pow(I1, i-2);
        }
    }

    computeGammaCoefficients(invar_c, dW_c, ddW_c, G_c);
}


