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

using namespace std;
using namespace muesli;


svkMaterial::svkMaterial(const std::string& name,
                           const materialProperties& cl)
:
finiteStrainMaterial(name, cl),
E(0.0), nu(0.0), lambda(0.0), mu(0.0)
{
    muesli::assignValue(cl, "young",   E);
    muesli::assignValue(cl, "poisson", nu);
    muesli::assignValue(cl, "lambda",  lambda);
    muesli::assignValue(cl, "mu",      mu);

    // 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 of them can be recovered fast
    if (density() > 0.0)
    {
        cp = sqrt((lambda+2.0*mu)/density());
        cs   = sqrt(mu/density());
    }
    bulk = lambda + 2.0/3.0 * mu;
}



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




finiteStrainMP* svkMaterial::createMaterialPoint() const
{
    finiteStrainMP *mp = new svkMP(*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 svkMaterial::print(std::ostream &of) const
{
    of  << "\n   Elastic Saint Venant-Kirchhoff material for finite deformation analysis";
    of  << "\n   Young modulus:  E      = " << E;
    of  << "\n   Poisson ratio:  nu     = " << nu;
    of  << "\n   Lame constants: Lambda = " << lambda;
    of  << "\n                   Mu     = " << mu;
    of  << "\n   Density                = " << density();
    of  << "\n   Wave velocities C_p    = " << cp;
    of  << "\n                   C_s    = " << cs;
}




void svkMaterial::setRandom()
{
    E   = muesli::randomUniform(1.0, 10000.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;
}




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

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




svkMP::svkMP(const svkMaterial &m) :
finiteStrainMP(m),
mat(&m)
{
    tn = tc = 0.0;
    Fn = itensor::identity();
    Fc = itensor::identity();
    Jn = 1.0;
    Jc = 1.0;
}




void svkMP::CauchyStress(istensor &sigma) const
{
    istensor S;
    secondPiolaKirchhoffStress(S);
    sigma = 1.0/Jc * istensor::FSFt(Fc, S);
}




void  svkMP::commitCurrentState()
{
    finiteStrainMP::commitCurrentState();
}




// this operation is the one with index representation C_abcd Vb Wc, which is NOT the
// one appearing in the material tangent
void svkMP::contractWithConvectedTangent(const ivector &V, const ivector &W, itensor &T) const
{
    T = mat->lambda * itensor::dyadic(V,W)
    +   mat->mu * itensor::dyadic(W,V)
    +   mat->mu *  V.dot(W) * itensor::identity();
}




// spatial elasticity tangent is
// cijkl = J^(-1) (lambda Bij Bkl + mu (Bik Bjl + Bil Bjk))
//
/* Given the fourth order tensor of elasticities C, and two vectors v, w
 * compute the second order tensor T with components
 *   T_ij = C_ipqj v_p w_q
 *
 *   T = J^(-1) * (lambda Bv otimes Bw + mu Bw otimes Bv + mu (v * Bw) B )
 *
 *  Note that the result is not symmetric.
 */
void svkMP::contractWithSpatialTangent(const ivector &v, const ivector &w, itensor &T) const
{
    istensor B = istensor::tensorTimesTensorTransposed(Fc);
    double   iJ = 1.0/sqrt(B.determinant());
    ivector  Bv(B*v), Bw(B*w);

    itensor D( itensor::dyadic(Bv, Bw) );

    T = (iJ*mat->lambda)*D + (iJ*mat->mu)*D.transpose() + (iJ*mat->mu*v.dot(Bw))*B;
}




void svkMP::convectedTangent(itensor4& C) const
{
    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) += mat->lambda;
                    if (i==k && j==l) C(i,j,k,l) += mat->mu;
                    if (i==l && j==k) C(i,j,k,l) += mat->mu;
                }
}




void svkMP::convectedTangentMatrix(double C[6][6]) const
{
    const double nu = mat->nu;
    const double k  = mat->E/(1.0+nu)/(1.0-2.0*nu);

    for (unsigned i=0; i<6; i++)
        for (unsigned j=0; j<6; j++)
            C[i][j] = 0.0;

    for (unsigned i=0; i<3; i++)
    {
        C[i][i]     = (1.0-nu)*k;
        C[i+3][i+3] = 0.5*(1.0-2.0*nu)*k;
        for (unsigned j=i+1; j<3; j++)
            C[i][j] = C[j][i] = nu*k;
    }
}




void svkMP::convectedTangentTimesSymmetricTensor(const istensor& M, istensor& CM) const
{
    CM = mat->lambda * M.trace() * istensor::identity() + 2.0*mat->mu*M;
}




double svkMP::dissipatedEnergy() const
{
    return 0.0;
}




double svkMP::effectiveStoredEnergy() const
{
    return storedEnergy();
}




double svkMP::kineticPotential() const
{
    return 0.0;
}




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




void  svkMP::resetCurrentState()
{
    finiteStrainMP::resetCurrentState();
}




void svkMP::secondPiolaKirchhoffStress(istensor &S) const
{
    istensor Cc = istensor::tensorTransposedTimesTensor(Fc);
    istensor E = 0.5*(Cc - istensor::identity() );

    S = (mat->lambda * E.trace() )*istensor::identity() + (2.0 * mat->mu) * E;
}




void svkMP::setConvergedState(const double time, const itensor& xFn)
{
    tn = time;
    Fn = xFn;
    Jn = Fn.determinant();
}




void svkMP::setRandom()
{
    finiteStrainMP::setRandom();
}




double svkMP::stiffness(const double stretch) const
{
    return 2.0*mat->mu*stretch + mat->lambda*stretch;
}




double svkMP::storedEnergy() const
{
    istensor Cc = istensor::tensorTransposedTimesTensor(Fc);
    istensor E = 0.5*(Cc - istensor::identity() );
    double   trE = E.trace();

    double W = 0.5* mat->lambda * trE*trE + mat->mu * E.contract(E);
    return W;
}




double svkMP::storedEnergy(const double stretch) const
{
    double E ( 0.5*(stretch*stretch - 1.0) );
    return mat->mu*E*E + 0.5*mat->lambda*E*E;
}




// for 1-dimensional problems
double svkMP::stress(const double stretch) const
{
    double E ( 0.5*(stretch*stretch - 1.0) );
    return 2.0*mat->mu*E+mat->lambda*E;
}




void  svkMP::updateCurrentState(const double theTime, const itensor& F)
{
    finiteStrainMP::updateCurrentState(theTime, F);
}

