/****************************************************************************
*
*                                 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 "neohook.h"
#define TESTNEOHOOK 0

using namespace std;
using namespace muesli;


neohookeanMaterial::neohookeanMaterial(const std::string& name,
                                         const materialProperties& cl)
:
  f_invariants(name, cl),
  E(0.0), nu(0.0), lambda(0.0), mu(0.0),
  regularized(false)
{
    muesli::assignValue(cl, "young",   E);
    muesli::assignValue(cl, "poisson", nu);
    muesli::assignValue(cl, "lambda",  lambda);
    muesli::assignValue(cl, "mu",      mu);

    if ( cl.find("subtype regularized") != cl.end()) regularized = true;

    // E and nu have priority. If they are defined, compute lambda and mu
    if (E*E  > 0)
    {
        lambda = E*nu/(1.0-2.0*nu)/(1.0+nu);
        mu     = E/2.0/(1.0+nu);
    }
    else
    {
        nu = lambda / 2.0 / (lambda+mu);
        E  = mu*2.0*(1.0+nu);
    }

    // we set all the constants, so that later on all of them can be recovered fast
    double rho = density();
    if (rho > 0.0)
    {
        cp = sqrt((lambda+2.0*mu)/rho);
        cs = sqrt(mu/rho);
    }
    bulk = lambda + 2.0/3.0 * mu;
}




neohookeanMaterial::neohookeanMaterial(const std::string& name,
                                         const double xE, const double xnu, const double rho)
:
f_invariants(name),
  E(xE), nu(xnu), lambda(0.0), mu(0.0), regularized(false)
{
    setDensity(rho);
    lambda = E*nu/(1.0-2.0*nu)/(1.0+nu);
    mu     = E/2.0/(1.0+nu);
}




double neohookeanMaterial::characteristicStiffness() const
{
    return E;
}




finiteStrainMP* neohookeanMaterial::createMaterialPoint() const
{
    muesli::finiteStrainMP* mp = new neohookeanMP(*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 neohookeanMaterial::print(std::ostream &of) const
{
    if (!regularized)
    {
        of  << "\n   Elastic Neo-Hookean material for finite deformation analysis"
        << "\n   Stored energy function:"
        << "\n            W(I1,J) = U(J) + mu/2 ( I1_C - 3) - mu log J"
        << "\n            U(J)    = lambda/2 * (log J)^2\n";
    }
    else
    {
        of  << "\n   Regularized elastic Neo-Hookean material for finite deformation analysis"
        << "\n   Stored energy function:"
        << "\n            W(I1,J) = U(J) + mu/2 (I1^tilde_C - 3)"
        << "\n            U(J)    = kappa/2 * (log J)^2\n";
    }

    of  << "\n   Young modulus:  E      = " << E;
    of  << "\n   Poisson ratio:  nu     = " << nu;
    of  << "\n   Lame constants: Lambda = " << lambda;
    of  << "\n                   Mu     = " << mu;
    of  << "\n   Bulk modulus:   k      = " << bulk;
    of  << "\n   Density                = " << density();
    of  << "\n   Wave velocities C_p    = " << cp;
    of  << "\n                   C_s    = " << cs;
    of  << "\n";
}




void neohookeanMaterial::setRandom()
{
    material::setRandom();
    E   = muesli::randomUniform(1.0, 10.0);
    nu  = muesli::randomUniform(0.05, 0.45);
    //setDensity(muesli::randomUniform(1.0, 100.0));

    lambda = E*nu/(1.0-2.0*nu)/(1.0+nu);
    mu     = E/2.0/(1.0+nu);
    cp     = sqrt((lambda+2.0*mu)/density());
    cs     = sqrt(2.0*mu/density());
    bulk   = lambda + 2.0/3.0 * mu;

    double s = muesli::randomUniform(-1.0, 1.0);
    regularized = s > 0.0 ? true : false;
    regularized = false;
}




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

    muesli::finiteStrainMP* p = this->createMaterialPoint();

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

    return isok;
}





neohookeanMP::neohookeanMP(const neohookeanMaterial&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];
}




// volumetric and deviatoric responses are decoupled
void neohookeanMP::contractWithMixedTangent(istensor& CM) const
{
    CM.setZero();
}




void neohookeanMP::convectedTangent(itensor4& c) const
{
    fisotropicMP::convectedTangent(c);
}




void neohookeanMP::convectedTangentTimesSymmetricTensor(const istensor& T,istensor& result) const
{
    const istensor C  = istensor::tensorTransposedTimesTensor(Fc);
    istensor Cinv = C.inverse();
    double J = sqrt(invar_c[2]);
    double logJ = log(J);

    if (mat->regularized)
    {
        fisotropicMP::convectedTangentTimesSymmetricTensor(T, result);
    }
    else
    {
        // result using compact notation
        result = mat->lambda * Cinv.dot(T)*Cinv
                + 2.0*(mat->mu - mat->lambda*logJ) * istensor::FSFt(Cinv, T);
    }
}




void neohookeanMP::secondPiolaKirchhoffStress(istensor &S) const
{
    fisotropicMP::secondPiolaKirchhoffStress(S);
}




void neohookeanMP::setConvergedState(const double time, const itensor& F)
{
    tn          = time;
    Fn          = F;
    Jn          = Fn.determinant();
    istensor Cn = istensor::tensorTransposedTimesTensor(Fn);

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

    double I3   = invar_n[2];
    double J    = sqrt(I3);
    double logJ = log(J);
    double iJ   = 1.0/J;

    ddW_n.setZero();

    if (mat->regularized)
    {
        double I1  = invar_n[0];
        double iJ4 = iJ*iJ*iJ*iJ;
        double Jm23= pow(J, -2.0/3.0);
        double Jm83= Jm23*Jm23*Jm23*Jm23;
        double Jm143= Jm83*Jm23*Jm23*Jm23;

        dW_n(0) = 0.5 * mat->mu * Jm23;
        dW_n(1) = 0.0;
        dW_n(2) = -1.0/6.0*mat->mu * I1 * Jm83 + 0.5 * mat->bulk * logJ * iJ * iJ;

        ddW_n(0,2) = ddW_n(2,0) = -mat->mu/6.0* Jm83;
        ddW_n(2,2) = 2.0/9.0*mat->mu*I1* Jm143 + 0.25*mat->bulk*iJ4 - 0.50* mat->bulk * logJ *iJ4;
    }

    else
    {
        dW_n(0) = mat->mu * 0.5;
        dW_n(1) = 0.0;
        dW_n(2) = -0.5*mat->mu * iJ*iJ + 0.5 * mat->lambda *logJ*iJ*iJ;
        ddW_n(2,2) = 1.0/( I3 * I3) * (0.5*mat->mu - 0.5*mat->lambda * logJ + 0.25*mat->lambda);
    }
    
    computeGammaCoefficients(invar_n, dW_n, ddW_n, G_n);
}




double neohookeanMP::storedEnergy() const
{
    double w(0.0);
    double J = sqrt( invar_c[2] );
    double I1 = invar_c[0];
    double logJ = log(J);

    // regularized  ->  W = U(J) + mu/2 ( I1tilde_C-3), U(J) = kappa/2 (log J)^2
    if (mat->regularized)
    {
        double U = 0.5*mat->bulk * logJ * logJ;
        double I1tilde = I1 * pow(J, -2.0/3.0);
        w = U + 0.5 * mat->mu * (I1tilde -3.0);
    }

    else
    {
        double U = 0.5*mat->lambda * logJ*logJ;
        w = U + 0.5*mat->mu* (I1 - 3.0) - mat->mu * logJ;
    }

    return w;
}




/* compressible neo hookean material.
 * stored energy function ->  W = U(J) + mu/2 ( I1_C - 3) - mu log J,  U(J) = lambda/2 * (log J)^2
 *           regularized  ->  W = U(J) + mu/2 ( I1tilde_C-3), U(J) = kappa/2 (log J)^2
 */
void neohookeanMP::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 I3   = invar_c[2];
    double J    = sqrt(I3);
    double logJ = log(J);
    double iJ   = 1.0/J;

    ddW_c.setZero();
    if (mat->regularized)
    {
        double I1  = invar_c[0];
        double iJ4 = iJ*iJ*iJ*iJ;
        double Jm23= pow(J, -2.0/3.0);
        double Jm83= Jm23*Jm23*Jm23*Jm23;
        double Jm143= Jm83*Jm23*Jm23*Jm23;

        dW_c(0) = 0.5 * mat->mu * Jm23;
        dW_c(1) = 0.0;
        dW_c(2) = -1.0/6.0*mat->mu * I1 * Jm83 + 0.5 * mat->bulk * logJ * iJ * iJ;

        ddW_c(0,2) = ddW_c(2,0) = -mat->mu/6.0* Jm83;
        ddW_c(2,2) = 2.0/9.0*mat->mu*I1* Jm143 + 0.25*mat->bulk*iJ4 - 0.50* mat->bulk * logJ *iJ4;
    }

    else
    {
        dW_c(0) = mat->mu * 0.5;
        dW_c(1) = 0.0;
        dW_c(2) = -0.5*mat->mu * iJ*iJ + 0.5 * mat->lambda *logJ*iJ*iJ;
        ddW_c(2,2) = 1.0/( I3 * I3) * (0.5*mat->mu - 0.5*mat->lambda * logJ + 0.25*mat->lambda);
    }

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