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

using namespace std;
using namespace muesli;


/* the input line of a viscoelastic material should set the elastic constants

 material, label = 1, type = viscoelastic,
 E = 1000.0, nu = 0.3,       # this sets the purely elastic constants

 and then as many pairs (eta,tau) as desired. these can be given in any
 order, but eta_i <-> tau_i are paired

 eta = 1.0, tau = 1.3, eta = 2.0, tau = 0.7, ....
 */
viscoelasticMaterial::viscoelasticMaterial(const std::string& name,
                                             const materialProperties& cl)
:
smallStrainMaterial(name, cl),
E(0.0), nu(0.0), lambda(0.0), mu(0.0), bulk(0.0), cp(0.0), cs(0.0)
{
    muesli::assignValue(cl, "young",   E);
    muesli::assignValue(cl, "poisson", nu);
    muesli::assignValue(cl, "lambda",  lambda);
    muesli::assignValue(cl, "mu",      mu);
    muesli::assignValue(cl, "eta",     eta);
    muesli::assignValue(cl, "tau",     tau);

    // E and nu have priority. If they are defined, define 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 mlog them can be recovered fast
    bulk = lambda + 2.0/3.0 * mu;
    double rho = density();
    if (rho > 0.0)
    {
        cp   = sqrt((lambda+2.0*mu)/rho);
        cs   = sqrt(mu/rho);
    }
}




viscoelasticMaterial::viscoelasticMaterial(const std::string& name,
                                             const double xE, const double xnu, const double rho,
                                             const size_t nvisco, const double* xeta, const double *xtau)
:
    smallStrainMaterial(name),
    E(xE), nu(xnu), lambda(0.0), mu(0.0), bulk(0.0), cp(0.0), cs(0.0)
{
    setDensity(rho);

    // E and nu have priority. If they are defined, define 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 mlog them can be recovered fast
    bulk = lambda + 2.0/3.0 * mu;
    if (rho > 0.0)
    {
        cp   = sqrt((lambda+2.0*mu)/rho);
        cs   = sqrt(mu/rho);
    }

    if (nvisco > 0 && xeta != 0 && xtau != 0 )
    {
        for (size_t k=0; k<nvisco; k++)
        {
            eta.push_back(xeta[k]);
            tau.push_back(xtau[k]);
        }
    }
}




bool viscoelasticMaterial::check() const
{
    bool ret = true;

    if ( eta.size() != tau.size() )  ret = false;

    return ret;
}




/* an object of the type "viscoelastic" creates a material point of type viscoelastic. The
 material point holds information that is not part of the material itself but
 that is particular of the specific (physical) point.
 */
smallStrainMP* viscoelasticMaterial::createMaterialPoint() const
{
    smallStrainMP* mp = new viscoelasticMP(*this);
    return mp;
}




// this function is much faster than the one with string property names, because
// avoids string comparisons. It should be used.
double viscoelasticMaterial::getProperty(const propertyName p) const
{
    return 0.0;
}





/* 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 viscoelasticMaterial::print(std::ostream &of) const
{
    of  << "\n   Small strain, viscoelastic, isotropic material. ";
    of  << "\n   In 3D problems, the viscoelastic response affects only the";
    of  << "\n   deviatoric strain. In 1D problems, the whole respose is";
    of  << "\n   viscoelastic and the limit value of the stiffness is E.\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;

    for (size_t a=0; a<eta.size(); a++)
        of  << "\n   (eta_i, tau_i) pair    : (" << eta[a] << ", " << tau[a] << ")";
}




void viscoelasticMaterial::setRandom()
{
    material::setRandom();

    E      = muesli::randomUniform(1.0, 10000.0);
    nu     = muesli::randomUniform(0.05, 0.45);
    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;

    int n = muesli::discreteUniform(1, 5);
    for (size_t a=0; a<n; a++)
    {
        eta.push_back( muesli::randomUniform(1.0,10.0) );
        tau.push_back( muesli::randomUniform(0.1,2.0) );
    }
}




bool viscoelasticMaterial::test(std::ostream  &of)
{
    bool isok=true;
    this->setRandom();

    smallStrainMP* p = this->createMaterialPoint();
    p->setRandom();

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

    return isok;
}




double viscoelasticMaterial::waveVelocity() const
{
    return cp;
}




viscoelasticMP::viscoelasticMP(const viscoelasticMaterial &m) :
smallStrainMP(m),
theta_n(0.0),
theta_c(0.0),
theViscoelasticIsotropicMaterial(m)
{
    size_t v = theViscoelasticIsotropicMaterial.eta.size();
    epsvisco_n.resize(v);
    epsvisco_c.resize(v);

    for (size_t a=0; a<v; a++)
    {
        epsvisco_n[a].setZero();
        epsvisco_c[a].setZero();
    }
}




viscoelasticMP::~viscoelasticMP()
{
    epsvisco_n.clear();
    epsvisco_c.clear();
}




void viscoelasticMP::commitCurrentState()
{


    smallStrainMP::commitCurrentState();
    epsvisco_n = epsvisco_c;
    edev_n     = edev_c;
    theta_n    = theta_c;

}




/* Given the fourth order tensor of viscoelasticities C, and two vectors v, w
 * compute the second order tensor T with components
 *   T_ij = C_ipjq v_p w_q
 *
 *  Note that the result is not symmetric.
 */

void viscoelasticMP::contractWithTangent(const ivector &v1, const ivector &v2,
                                           itensor &T) const
{
    // shear response
    const double dt    = time_c - time_n;
    const double muinf = theViscoelasticIsotropicMaterial.mu;
    double       mu_eq = muinf;

    for (size_t a=0; a < epsvisco_n.size(); a++)
    {
        double   taua  = theViscoelasticIsotropicMaterial.tau[a];
        double   mua   = theViscoelasticIsotropicMaterial.eta[a]/(2.0*taua);
        mu_eq         +=  mua * taua/(dt+taua);
    }

    // bulk response
    const double& kinf     = theViscoelasticIsotropicMaterial.bulk;
    const double lambda_eq = kinf - 2.0/3.0*mu_eq;

    itensor A, B;

    T = itensor::identity();
    T *= v1.dot(v2) * mu_eq;

    A = itensor::dyadic(v1, v2);
    B = A.transpose();

    A *= lambda_eq;
    B *= mu_eq;

    T += A;
    T += B;
}




void viscoelasticMP::contractWithDeviatoricTangent(const ivector &v1, const ivector &v2, itensor &T) const
{
    const double dt = time_c - time_n;

    // shear response
    double   muinf = theViscoelasticIsotropicMaterial.mu;
    double   mu_eq = muinf;
    for (size_t a=0; a < epsvisco_n.size(); a++)
    {
        double   tau   = theViscoelasticIsotropicMaterial.tau[a];
        double   mua   = theViscoelasticIsotropicMaterial.eta[a]/(2.0*tau);

        mu_eq +=  mua * tau/(dt+tau);
    }

    T = mu_eq * v1.dot(v2) * itensor::identity()
    +   mu_eq *              itensor::dyadic(v2, v1)
    -   2.0/3.0*mu_eq *      itensor::dyadic(v1, v2);
}




double viscoelasticMP::deviatoricEnergy() const
{
    return theViscoelasticIsotropicMaterial.mu * edev_c.contract(edev_c);
}




void viscoelasticMP::deviatoricStress(istensor& s) const
{
    const double dt    = time_c - time_n;
    const double muinf = theViscoelasticIsotropicMaterial.mu;

    // contribution from purely elastic terms
    s = 2.0*muinf*edev_c;

    // viscous contributions from elements
    for (size_t a=0; a < epsvisco_n.size(); a++)
    {
        double   taua  = theViscoelasticIsotropicMaterial.tau[a];
        double   mua   = theViscoelasticIsotropicMaterial.eta[a]/(2.0*taua);
        istensor einc  = edev_c - epsvisco_n[a];

        s += 2.0 * mua * taua/(dt+taua) * einc;
    }
}




double viscoelasticMP::dissipatedEnergy() const
{
    double diss  = 0.0;
    double dt = time_c - time_n;

    // elastic potential
    for (size_t a=0; a < epsvisco_n.size(); a++)
    {
        double   taua  = theViscoelasticIsotropicMaterial.tau[a];
        double   etaa  = theViscoelasticIsotropicMaterial.eta[a];
        double   mua   = etaa/(2.0*taua);
        istensor einc  = edev_c - epsvisco_n[a];

        diss += mua * (taua/dt)/(dt+taua) * einc.squaredNorm();
    }

    return diss;
}




double viscoelasticMP::effectiveStoredEnergy() const
{
    const double dt = time_c - time_n;
    return storedEnergy() + dt*kineticPotential();
}




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




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

    state_n.theDouble.push_back(theta_n);
    state_n.theStensor.push_back(edev_n);
    state_n.theStensor.insert( state_n.theStensor.end(), epsvisco_n.begin(), epsvisco_n.end() );

    return state_n;
}




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




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

    state_c.theDouble.push_back(theta_c);
    state_c.theStensor.push_back(edev_c);
    state_c.theStensor.insert( state_c.theStensor.end(), epsvisco_c.begin(), epsvisco_c.end() );

    return state_c;
}




double viscoelasticMP::kineticPotential() const
{
    double dt = time_c - time_n;

    double psi = 0.0;
    for (size_t a=0; a < epsvisco_n.size(); a++)
    {
        double   etaa  = theViscoelasticIsotropicMaterial.eta[a];
        istensor einc  = epsvisco_c[a] - epsvisco_n[a];

        psi += 0.5 * etaa * einc.squaredNorm() / (dt*dt);
    }
    return psi;
}




double viscoelasticMP::plasticSlip() const
{
    return 0.0;
}




void viscoelasticMP::resetCurrentState()
{

    smallStrainMP::resetCurrentState();

    epsvisco_c = epsvisco_n;
    edev_c     = edev_n;
    theta_c    = theta_n;
}




void viscoelasticMP::setConvergedState(const double& tn, const istensor& epsn,
                                         const std::vector<istensor>& epsv, const istensor& epsdev, const double& theta)
{


    time_n     = tn;
    eps_n      = epsn;
    epsvisco_n = epsv;
    edev_n     = epsdev;
    theta_n    = theta;

}




void viscoelasticMP::setRandom()
{
    smallStrainMP::setRandom();
    theta_c = theta_n = muesli::randomUniform(0.01, 0.1);
    istensor tmp;
    tmp.setRandom();
    edev_c  = edev_n  = istensor::deviatoricPart(tmp);

    for (size_t a=0; a < epsvisco_n.size(); a++)
    {
        tmp.setRandom();
        epsvisco_n[a] = istensor::deviatoricPart(tmp);
    }
    epsvisco_c = epsvisco_n;

}




double viscoelasticMP::storedEnergy() const
{
    // elastic potential
    const double  muinf = theViscoelasticIsotropicMaterial.mu;
    const double  kinf  = theViscoelasticIsotropicMaterial.bulk;

    double W = 0.5*kinf*theta_c*theta_c + muinf*edev_c.squaredNorm();

    for (size_t a=0; a < epsvisco_n.size(); a++)
    {
        double   taua  = theViscoelasticIsotropicMaterial.tau[a];
        double   mua   = theViscoelasticIsotropicMaterial.eta[a]/(2.0*taua);
        istensor einc  = edev_c - epsvisco_c[a];
        W += mua * einc.squaredNorm();
    }
    return W;
}




void viscoelasticMP::stress(istensor& sigma) const
{
    deviatoricStress(sigma);
    sigma += theViscoelasticIsotropicMaterial.bulk * theta_c * istensor::identity();
}




void viscoelasticMP::tangentTensor(itensor4& C) const
{
    // shear response
    const double dt    = time_c - time_n;
    const double muinf = theViscoelasticIsotropicMaterial.mu;

    double mu_eq = muinf;
    for (size_t a=0; a < epsvisco_n.size(); a++)
    {
        double   taua  = theViscoelasticIsotropicMaterial.tau[a];
        double   mua   = theViscoelasticIsotropicMaterial.eta[a]/(2.0*taua);

        mu_eq +=  mua * taua/(dt+taua);
    }

    // bulk response
    const double& kinf  = theViscoelasticIsotropicMaterial.bulk;
    double lambda_eq = kinf - 2.0/3.0*mu_eq;


    C.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++)
                {
                    if (i==j && k==l) C(i,j,k,l) += lambda_eq;
                    if (i==k && j==l) C(i,j,k,l) += mu_eq;
                    if (i==l && j==k) C(i,j,k,l) += mu_eq;
                }
}




thPotentials viscoelasticMP::thermodynamicPotentials() const
{
    thPotentials tp;
    return tp;
}




void viscoelasticMP::updateCurrentState(const double theTime, const istensor& strain)
{


    smallStrainMP::updateCurrentState(theTime, strain);
    edev_c  = istensor::deviatoricPart(strain);
    theta_c = strain.trace();

    // viscous contributions from prony elements
    const double dt = time_c - time_n;
    for (size_t a=0; a < epsvisco_n.size(); a++)
    {
        const double& taua = theViscoelasticIsotropicMaterial.tau[a];
        epsvisco_c[a]      = epsvisco_n[a] + dt/(dt+taua)* (edev_c - epsvisco_n[a]);
    }

}




double viscoelasticMP::volumetricEnergy() const
{
    double V = 0.5* theViscoelasticIsotropicMaterial.bulk * theta_c * theta_c;
    return V;
}




double viscoelasticMP::volumetricStiffness() const
{
    return theViscoelasticIsotropicMaterial.bulk;
}

