/****************************************************************************
*
*                                 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 "elastic.h"
#include "muesli/Math/matrix.h"
#include "muesli/Math/realvector.h"
#include <cmath>

using namespace std;
using namespace muesli;



elasticIsotropicMaterial::elasticIsotropicMaterial(const std::string& name,
                                   const double E, const double nu, const double xrho)
:
smallStrainMaterial(name),
lambda(E*nu/(1.0-2.0*nu)/(1.0+nu)),
mu(E/2.0/(1.0+nu)),
rho(xrho)
{
}




elasticIsotropicMaterial::elasticIsotropicMaterial(const std::string& name,
                                   const materialProperties& cl)
:
smallStrainMaterial(name, cl),
lambda(0.0), mu(0.0)
{
    double E, nu;

    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);
    }
}




bool elasticIsotropicMaterial::check() const
{
    if (mu > 0 && lambda+2.0*mu > 0) return true;
    else return false;
}




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




double elasticIsotropicMaterial::density() const
{
    return rho;
}




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

    const double nu   = lambda / 2.0 / (lambda+mu);
    const double E    = mu*2.0*(1.0+nu);
    const double bulk = lambda + 2.0/3.0 * mu;
    const double cp   = rho > 0.0 ? sqrt((lambda+2.0*mu)/rho) : 0.0;
    const double cs   = rho > 0.0 ? sqrt(mu/rho) : 0.0;

    switch (p)
    {
        case PR_LAMBDA:     ret = lambda;   break;
        case PR_MU:         ret = mu;       break;
        case PR_YOUNG:      ret = E;        break;
        case PR_POISSON:    ret = nu;       break;
        case PR_BULK:       ret = bulk; break;
        case PR_CP:         ret = cp;       break;
        case PR_CS:         ret = cs;       break;

        default:
            std::cout << "\n Error in elasticIsotropicMaterial. Property not defined";
    }
    return ret;
}




/* this function is always called once the material is defined, so apart from
 printing its information, we take the opportunity to clean up some mlog its
 data, in particular, setting all the possible constants
 */
void elasticIsotropicMaterial::print(std::ostream &of) const
{
    std::ostream& mlog = material::getLogger();

    const double nu   = lambda / 2.0 / (lambda+mu);
    const double E    = mu*2.0*(1.0+nu);
    const double bulk = lambda + 2.0/3.0 * mu;

    mlog  << "\n   Small strain, elastic, isotropic material ";
    mlog  << "\n   Young modulus   E      : " << E;
    mlog  << "\n   Poisson ratio   nu     : " << nu;
    mlog  << "\n   Lame constants  lambda : " << lambda;
    mlog  << "\n                   mu     : " << mu;
    mlog  << "\n   Bulk modulus    k      : " << bulk;
    mlog  << "\n   Density                : " << rho;

    if (rho > 0.0)
    {
        double cp = sqrt((lambda+2.0*mu)/rho);
        double cs = sqrt(mu/rho);

        mlog  << "\n   Wave velocities C_p    : " << cp;
        mlog  << "\n                   C_s    : " << cs;
    }

    smallStrainMaterial::print(of);
    mlog << std::flush;
}




void elasticIsotropicMaterial::setRandom()
{
    //double E  = muesli::randomUniform(1.0, 10000.0);
    double E  = muesli::randomUniform(1.0e7, 1.0e9);
    double nu = muesli::randomUniform(0.05, 0.45);
    rho       = muesli::randomUniform(1.0, 100.0);
    lambda    = E*nu/(1.0-2.0*nu)/(1.0+nu);
    mu        = E/2.0/(1.0+nu);
}




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

    smallStrainMP* p = this->createMaterialPoint();

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

    return isok;
}




double elasticIsotropicMaterial::waveVelocity() const
{
    double cp = rho > 0.0 ? sqrt((lambda+2.0*mu)/rho) : 0.0;
    return cp;
}




elasticIsotropicMP::elasticIsotropicMP(const elasticIsotropicMaterial &m) :
muesli::smallStrainMP(m),
theElasticIsotropicMaterial(m)
{
}




elasticIsotropicMP::~elasticIsotropicMP()
{

}





void  elasticIsotropicMP::commitCurrentState()
{
    smallStrainMP::commitCurrentState();
}




/* Given the fourth order tensor mlog elasticities C, and two vectors v, w
 * compute the second order tensor T with components
 *   T_ij = C_ipjq v_p w_q
 *
 *  which, for linear isotropic materials is just
 *  T = lambda gradV otimes gradW + mu gradW otimes gradV + mu (gradV * gradW) One
 *
 *  Note that the result is not symmetric.
 */

void elasticIsotropicMP::contractWithTangent(const ivector &v1, const ivector &v2, itensor &T) const
{
    itensor A, B;

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

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

    A *= theElasticIsotropicMaterial.lambda;
    B *= theElasticIsotropicMaterial.mu;

    T += A;
    T += B;
}




void elasticIsotropicMP::contractWithDeviatoricTangent(const ivector &v1, const ivector &v2, itensor &T) const
{
    const double&  mu = theElasticIsotropicMaterial.mu;

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




double elasticIsotropicMP::deviatoricEnergy() const
{
    istensor e = istensor::deviatoricPart(eps_c);
    return theElasticIsotropicMaterial.mu * e.contract(e);
}




void elasticIsotropicMP::deviatoricStress(istensor& s) const
{
    s = (2.0*theElasticIsotropicMaterial.mu)* istensor::deviatoricPart(eps_c);
}




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




//II - Cinverse:Cep for small thermomechanical element
void elasticIsotropicMP::dissipationTangent(itensor4& D) const
{
    D.setZero();
}




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




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




materialState elasticIsotropicMP::getConvergedState() const
{
    materialState st;
    st.theTime = time_n;
    st.theStensor.push_back(eps_n);

    return st;
}




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




materialState elasticIsotropicMP::getCurrentState() const
{
    materialState st;
    st.theTime = time_c;
    st.theStensor.push_back(eps_c);

    return st;
}




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




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




double elasticIsotropicMP::pressure() const
{
    const double bulk = theElasticIsotropicMaterial.lambda + 2.0/3.0 * theElasticIsotropicMaterial.mu;
    double p = -bulk * eps_c.trace();
    return p;
}




void elasticIsotropicMP::resetCurrentState()
{
    smallStrainMP::resetCurrentState();
}




void elasticIsotropicMP::setConvergedState(const double xtn, const istensor& strainn)
{
    time_n = xtn;
    eps_n  = strainn;
}




void elasticIsotropicMP::setRandom()
{
    smallStrainMP::setRandom();
}




double elasticIsotropicMP::storedEnergy() const
{
    double tr = eps_c.trace();
    double W  = theElasticIsotropicMaterial.mu*eps_c.squaredNorm() + 0.5*theElasticIsotropicMaterial.lambda*tr*tr;
    return W;
}




void elasticIsotropicMP::stress(istensor& sigma) const
{
    const double mu     = theElasticIsotropicMaterial.mu;
    const double lambda = theElasticIsotropicMaterial.lambda;
    sigma = 2.0*mu*eps_c + lambda * eps_c.trace() * istensor::identity();
}




void elasticIsotropicMP::tangentMatrix(double C[6][6]) const
{
    const double nu = theElasticIsotropicMaterial.lambda / 2.0 / (theElasticIsotropicMaterial.lambda+theElasticIsotropicMaterial.mu);
    const double E  = theElasticIsotropicMaterial.mu*2.0*(1.0+nu);
    const double k  = 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 elasticIsotropicMP::tangentTensor(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) += theElasticIsotropicMaterial.lambda;
                    if (i==k && j==l) C(i,j,k,l) += theElasticIsotropicMaterial.mu;
                    if (i==l && j==k) C(i,j,k,l) += theElasticIsotropicMaterial.mu;
                }
}




thPotentials elasticIsotropicMP::thermodynamicPotentials() const
{
    thPotentials tp;

    const double mu   = theElasticIsotropicMaterial.mu;
    const double bulk = theElasticIsotropicMaterial.lambda + 2.0/3.0 * mu;

    istensor e = istensor::deviatoricPart(eps_c);
    double theta = eps_c.trace();

    tp.volumetric     = 0.5 * bulk * theta * theta;
    tp.deviatoric     = mu * e.contract(e);
    tp.internalEnergy = tp.volumetric+tp.deviatoric;
    tp.freeEntropy    = tp.internalEnergy;

    return tp;
}



void elasticIsotropicMP::updateCurrentState(const double t, const istensor& strain)
{
    smallStrainMP::updateCurrentState(t, strain);
}




double elasticIsotropicMP::volumetricStiffness() const
{
    const double bulk = theElasticIsotropicMaterial.lambda + 2.0/3.0 * theElasticIsotropicMaterial.mu;
    return bulk;
}




double elasticIsotropicMP::volumetricEnergy() const
{
    const double bulk = theElasticIsotropicMaterial.lambda + 2.0/3.0 * theElasticIsotropicMaterial.mu;
    double theta = eps_c.trace();
    double V     = 0.5*bulk * theta * theta;
    return V;
}




///////////////////////////////////////////////////////////////////////////////////////////

elasticAnisotropicMaterial::elasticAnisotropicMaterial(const std::string& name)
:
 muesli::smallStrainMaterial(name),
 rho(1.0)
{
    cc.resize(6,6);
    cc.setZero();
}




elasticAnisotropicMaterial::elasticAnisotropicMaterial(const std::string& name,
                                                         const double cv[21],
                                                         const double xrho)
:
smallStrainMaterial(name),
rho(xrho)
{
    cc.resize(6,6);
    unsigned k = 0;
    for (unsigned i=0; i<6; i++)
        for (unsigned j=i; j<6; j++)
            cc(i,j) = cc(j,i) = cv[k++];
}




elasticAnisotropicMaterial::elasticAnisotropicMaterial(const std::string& name,
                                                     const materialProperties& cl)
:
smallStrainMaterial(name, cl),
rho(1.0)
{
    double cv[21];
    muesli::assignValue(cl, "c1111",   cv[0]);
    muesli::assignValue(cl, "c1122",   cv[1]);
    muesli::assignValue(cl, "c1133",   cv[2]);
    muesli::assignValue(cl, "c1123",   cv[3]);
    muesli::assignValue(cl, "c1113",   cv[4]);
    muesli::assignValue(cl, "c1112",   cv[5]);

    muesli::assignValue(cl, "c2222",   cv[6]);
    muesli::assignValue(cl, "c2233",   cv[7]);
    muesli::assignValue(cl, "c2223",   cv[8]);
    muesli::assignValue(cl, "c2213",   cv[9]);
    muesli::assignValue(cl, "c2212",   cv[10]);

    muesli::assignValue(cl, "c3333",   cv[11]);
    muesli::assignValue(cl, "c3323",   cv[12]);
    muesli::assignValue(cl, "c3313",   cv[13]);
    muesli::assignValue(cl, "c3312",   cv[14]);

    muesli::assignValue(cl, "c2323",   cv[15]);
    muesli::assignValue(cl, "c2313",   cv[16]);
    muesli::assignValue(cl, "c2312",   cv[17]);

    muesli::assignValue(cl, "c1313",   cv[18]);
    muesli::assignValue(cl, "c1312",   cv[19]);

    muesli::assignValue(cl, "c1212",   cv[20]);


    cc.resize(6,6);
    unsigned k = 0;
    for (unsigned i=0; i<6; i++)
        for (unsigned j=i; j<6; j++)
            cc(i,j) = cc(j,i) = cv[k++];
}




bool elasticAnisotropicMaterial::check() const
{
    // we should check at least that the stiffness matrix is positive definite
    matrix evec(6,6);
    complexvector ev = evec.eigenvalues();

    double maxev = 0.0;
    for (unsigned a=0; a<6; a++)
        if ( ev.real()[a] > maxev ) maxev = ev.real()[a];

    bool isok = true;
    if (maxev <= 0.0) isok = false;

    return isok;
}




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




double elasticAnisotropicMaterial::density() const
{
    return rho;
}




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




/* this function is always called once the material is defined, so apart from
 printing its information, we take the opportunity to clean up some mlog its
 data, in particular, setting all the possible constants
 */
void elasticAnisotropicMaterial::print(std::ostream &of) const
{
    std::ostream& mlog = material::getLogger();

    mlog << "\n   Small strain, elastic, anisotropic material ";
    mlog << "\n   Density                : " << rho;
    mlog << "\n   Stiffness matrix in Voigt notation: ";
    mlog << "\n   (11, 22, 33, 23, 13, 12)";
    cc.print(mlog);
    mlog << std::flush;
}




void elasticAnisotropicMaterial::setRandom()
{
    cc.setZero();
    for (unsigned a=0; a<6; a++)
    {
        double eval = muesli::randomUniform(1.0, 10.0);

        realvector evec(6);
        evec.setRandom();
        evec.normalize();

        for (unsigned i=0; i<6; i++)
            for (unsigned j=0; j<6; j++)
                cc(i,j) += eval * evec[i] * evec[j];
    }
}




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

    smallStrainMP* p = this->createMaterialPoint();

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

    return isok;
}




double elasticAnisotropicMaterial::waveVelocity() const
{
    double c = 0.0;

    istensor K;
    K(0,0) = cc(0,0);
    K(0,1) = K(1,0) = cc(0,5);
    K(0,2) = K(2,0) = cc(0,4);
    K(1,1) = cc(5,5);
    K(1,2) = K(2,1) = cc(4,5);
    K(2,2) = cc(4,4);

    istensor M;
    M = rho * istensor::identity();
    istensor sqrtM = istensor::squareRoot(M);

    istensor A = istensor::FSFt(sqrtM, K);
    ivector  ev = A.eigenvalues();

    c = ev.max();

    return c;
}




elasticAnisotropicMP::elasticAnisotropicMP(const elasticAnisotropicMaterial &m) :
muesli::smallStrainMP(m),
theElasticAnisotropicMaterial(m)
{
}




elasticAnisotropicMP::~elasticAnisotropicMP()
{

}





void  elasticAnisotropicMP::commitCurrentState()
{
    smallStrainMP::commitCurrentState();
}




double elasticAnisotropicMP::deviatoricEnergy() const
{
    istensor edev = istensor::deviatoricPart(eps_c);
    itensor4 c; tangentTensor(c);
    istensor sdev = c(edev);

    return 0.5 * edev.contract(sdev);
}




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




//II - Cinverse:Cep for small thermomechanical element
void elasticAnisotropicMP::dissipationTangent(itensor4& D) const
{
    D.setZero();
}




materialState elasticAnisotropicMP::getConvergedState() const
{
    materialState st;
    st.theTime = time_n;
    st.theStensor.push_back(eps_n);

    return st;
}




materialState elasticAnisotropicMP::getCurrentState() const
{
    materialState st;
    st.theTime = time_c;
    st.theStensor.push_back(eps_c);

    return st;
}




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




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





void elasticAnisotropicMP::resetCurrentState()
{
    smallStrainMP::resetCurrentState();
}




void elasticAnisotropicMP::setConvergedState(const double xtn, const istensor& strainn)
{
    time_n = xtn;
    eps_n  = strainn;
}




void elasticAnisotropicMP::setRandom()
{
    smallStrainMP::setRandom();
}




double elasticAnisotropicMP::storedEnergy() const
{
    istensor sigma;
    stress(sigma);

    return 0.5 * eps_c.contract(sigma);
}




void elasticAnisotropicMP::stress(istensor& sigma) const
{
    realvector strain(6);
    strain[0] = eps_c(0,0);
    strain[1] = eps_c(1,1);
    strain[2] = eps_c(2,2);
    strain[3] = 2.0*eps_c(1,2);
    strain[4] = 2.0*eps_c(2,0);
    strain[5] = 2.0*eps_c(0,1);

    realvector stress(6);
    stress = theElasticAnisotropicMaterial.cc*strain;
    double S[6];
    for (unsigned i=0; i<6; i++) S[i] = stress[i];
    vectorToContraContraSymTensor(S, sigma);
}




void elasticAnisotropicMP::tangentMatrix(double C[6][6]) const
{
    const matrix& cc = theElasticAnisotropicMaterial.cc;
    for (unsigned i=0; i<6; i++)
        for (unsigned j=0; j<6; j++)
            C[i][j] = cc(i,j);
}




void elasticAnisotropicMP::tangentTensor(itensor4& C) const
{
    double Cmat[6][6];
    tangentMatrix(Cmat);
    matrixToTensor(Cmat, C);
}




thPotentials elasticAnisotropicMP::thermodynamicPotentials() const
{
    thPotentials tp;

    return tp;
}






void elasticAnisotropicMP::updateCurrentState(const double t, const istensor& strain)
{
    smallStrainMP::updateCurrentState(t, strain);
}




double elasticAnisotropicMP::volumetricEnergy() const
{
    return 0.0;
}




elasticOrthotropicMaterial::elasticOrthotropicMaterial(const std::string& name,
                           const double c[9],
                           const double rhox)
:
elasticAnisotropicMaterial(name)
{
    cc(0,0) = c[0];
    cc(0,1) = cc(1,0) = c[1];
    cc(0,2) = cc(2,0) = c[2];
    cc(1,1) = c[3];
    cc(1,2) = cc(2,1) = c[4];
    cc(2,2) = c[5];
    cc(3,3) = c[6];
    cc(4,4) = c[7];
    cc(5,5) = c[8];

    rho = rhox;
}




elasticOrthotropicMaterial::elasticOrthotropicMaterial(const std::string& name,
                                                         const materialProperties& cl)
:
elasticAnisotropicMaterial(name, cl)
{

}




/* an object of the type "elasticAnisotropic" creates a material point of type elastic. The
 material point holds information that is not part of the material itself but
 that is particular mlog the specific (physical) point.
 */
smallStrainMP* elasticOrthotropicMaterial::createMaterialPoint() const
{
    smallStrainMP* mp = new elasticOrthotropicMP(*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 mlog its
 data, in particular, setting all the possible constants
 */
void elasticOrthotropicMaterial::print(std::ostream &of) const
{
    std::ostream& mlog = material::getLogger();

    mlog << "\n   Small strain, elastic, orthotropic material ";
    mlog << "\n   Density                : " << rho;
    mlog << "\n   Stiffness matrix in Voigt notation: ";
    mlog << "\n   (11, 22, 33, 23, 13, 12)";
    cc.print(mlog);
    mlog << std::flush;
}




void elasticOrthotropicMaterial::setRandom()
{
    cc.setZero();
    const double E1 = randomUniform(1.0, 10.0);
    const double E2 = randomUniform(1.0, 10.0);
    const double E3 = randomUniform(1.0, 10.0);

    const double nu23 = randomUniform(0.0, 0.49);
    const double nu13 = randomUniform(0.0, 0.49);
    const double nu12 = randomUniform(0.0, 0.49);

    const double G23 = randomUniform(1.0, 10.0);
    const double G13 = randomUniform(1.0, 10.0);
    const double G12 = randomUniform(1.0, 10.0);

    cc(0,0) = 1.0/E1;   cc(0,1) = -nu12/E1; cc(0,2) = -nu13/E1;
    cc(1,0) = -nu12/E1; cc(1,1) = 1.0/E2;   cc(1,2) = -nu23/E2;
    cc(2,0) = -nu13/E1; cc(2,1) = -nu23/E2; cc(2,2) = 1.0/E3;
    cc(3,3) = 0.5/G23;
    cc(4,4) = 0.5/G13;
    cc(5,5) = 0.5/G12;

    rho = randomUniform(1.0, 2.0);
}




elasticOrthotropicMP::elasticOrthotropicMP(const elasticOrthotropicMaterial &m)
:
muesli::elasticAnisotropicMP(m)
{

}




void elasticOrthotropicMP::setRandom()
{
    smallStrainMP::setRandom();
}
