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

using namespace std;
using namespace muesli;


mooneyMaterial::mooneyMaterial(const std::string& name,
                 const muesli::materialProperties& mp)
    :
    f_invariants(name, mp),
    a0(0.0), a1(0.0), a2(0.0),
    incompressible(false)
{
    muesli::assignValue(mp, "alpha0",  a0);
    muesli::assignValue(mp, "alpha1",  a1);
    muesli::assignValue(mp, "alpha2",  a2);

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




double mooneyMaterial::characteristicStiffness() const
{

    return 2*a1; //G
}




muesli::finiteStrainMP* mooneyMaterial::createMaterialPoint() const
{
    muesli::finiteStrainMP* mp = new mooneyMP(*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 mooneyMaterial::print(std::ostream &of) const
{
    if (incompressible)
    {
        of  << "\n   Elastic, incompressible Mooney-Rivlin material for finite deformation analysis"
            << "\n   Stored energy function:"
            << "\n            W(I1,I2) = a1 (I1t - 3) + a2 (I2t - 3)"
            << "\n            I1t = I1_dev, I2t = I2_dev";
    }
    
    else
    {
        of  << "\n   Elastic, compressible Mooney-Rivlin material for finite deformation analysis"
            << "\n   Stored energy function:"
            << "\n            W(I1,I2,I3) = a1 (I1t - 3) + a2 (I2t - 3) + a/2(J-1)^2"
            << "\n            I1t = I1_dev, I2t = I2_dev";

        of  << "\n   Bulk coeff:     alpha0 = " << a0;
    }
    
    of  << "\n   Shear coeff:  2*alpha1 = " << 2.0*a1;
    of  << "\n   Third coeff:    alpha2 = " << a2;
    of  << "\n   Density                = " << density();
    of  << "\n   Wave velocities C_p    = " << cp;
    of  << "\n                   C_s    = " << cs;
    of  << "\n";
}




void mooneyMaterial::setRandom()
{
    incompressible = true;
    a1 = muesli::randomUniform(1.0, 0.45);
    a2 = muesli::randomUniform(1.0, 100.0);
    
    bulk = a0;
    mu   = 2.0*a1;
    lambda = bulk - 2.0/3.0*mu;
    
    cp     = sqrt((lambda+2.0*mu)/density());
    cs     = sqrt(2.0*mu/density());
}




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

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

    return isok;
}





mooneyMP::mooneyMP(const mooneyMaterial &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 mooneyMP::setConvergedState(const double time, const itensor& F)
{
    tn          = time;
    Fn          = F;
    Jn          = F.determinant();
    istensor Cn = istensor::tensorTransposedTimesTensor(Fn);

    invar_n[0]   = Cn.invariant1();
    invar_n[1]   = Cn.invariant2();
    invar_n[2]   = Cn.invariant3();

    double I1   = invar_n[0];
    double I2   = invar_n[1];
    double I3   = invar_n[2];
    double iJ   = 1.0/Jn;

    const double a0 = mat->a0;
    const double a1 = mat->a1;
    const double a2 = mat->a2;

    ddW_n.setZero();
    if (mat->incompressible)
    {
        dW_n(0) = a1;
        dW_n(1) = a2;
        dW_n(2) = 0.0;
    }

    else
    {
        const double I3m13 = 1.0/pow(I3,1.0/3.0);
        const double I3m23 = I3m13*I3m13;
        const double I3m43 = I3m23*I3m23;
        const double I3m53 = I3m43*I3m13;
        const double I3m73 = I3m53*I3m23;
        const double I3m83 = I3m73*I3m13;


        dW_n(0) = a1 * I3m13;
        dW_n(1) = a2 * I3m23;
        dW_n(2) = -a1/3.0 * I1 * I3m43 - 2.0/3.0 * a2*I2*I3m53 + a0/2.0*(1.0- iJ);

        ddW_n(0,2) = ddW_n(2,0) = -1.0/3.0 * a1 * I3m43;
        ddW_n(1,2) = ddW_n(2,1) = -2.0/3.0 * a2 * I3m53;
        ddW_n(2,2) = 4.0/9.0* a1* I1 * I3m73 + 10.0/9.0 * a2 * I2 * I3m83 + a0/4.0*iJ*iJ*iJ;
    }

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




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




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




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




double mooneyMP::storedEnergy() const
{
    const double I1 = invar_c[0];
    const double I2 = invar_c[1];
    const double I3 = invar_c[2];
    double w;
    
    if ( mat->incompressible)
    {
        w = mat->a1 * (I1-3.0) + mat->a2 * (I2 - 3.0);
    }
    
    else
    {
        const double J  = sqrt(I3);
        const double I3m13 = 1.0/pow(I3,1.0/3.0);
        const double I3m23 = I3m13*I3m13;
        
        const double I1t = I1 * I3m13;
        const double I2t = I2 * I3m23;
        
        w = mat->a1 * (I1t-3.0) + mat->a2 * (I2t - 3.0) + 0.5* mat->a0*(J-1.0)*(J-1.0);
    }
    
    return w;
}




void mooneyMP::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 I2   = invar_c[1];
    double I3   = invar_c[2];
    double J    = sqrt(I3);
    double iJ   = 1.0/J;

    const double a0    = mat->a0;
    const double a1    = mat->a1;
    const double a2    = mat->a2;

    
    ddW_c.setZero();
    
    if (mat->incompressible)
    {
        dW_c(0) = a1;
        dW_c(1) = a2;
        dW_c(2) = 0.0;        
    }
    
    else
    {
        const double I3m13 = 1.0/pow(I3,1.0/3.0);
        const double I3m23 = I3m13*I3m13;
        const double I3m43 = I3m23*I3m23;
        const double I3m53 = I3m43*I3m13;
        const double I3m73 = I3m53*I3m23;
        const double I3m83 = I3m73*I3m13;
        
        
        dW_c(0) = a1 * I3m13;
        dW_c(1) = a2 * I3m23;
        dW_c(2) = -a1/3.0 * I1 * I3m43 - 2.0/3.0 * a2*I2*I3m53 + a0/2.0*(1.0- iJ);
        
        ddW_c(0,2) = ddW_c(2,0) = -1.0/3.0 * a1 * I3m43;
        ddW_c(1,2) = ddW_c(2,1) = -2.0/3.0 * a2 * I3m53;
        ddW_c(2,2) = 4.0/9.0* a1* I1 * I3m73 + 10.0/9.0 * a2 * I2 * I3m83 + a0/4.0*iJ*iJ*iJ;
    }   
    
    computeGammaCoefficients(invar_c, dW_c, ddW_c, G_c);
}


