/****************************************************************************
*
*                                 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 <stdio.h>
#include "thermofinitestrain.h"
#include "muesli/Finitestrain/finitestrainlib.h"

#define THETA0 100.0

using namespace muesli;

thermofiniteStrainMaterial::thermofiniteStrainMaterial(const std::string& name,
                                                         const materialProperties& cl)
:
material(name, cl),
theFSMaterial(0),
_heatSupply(0.0),
_thermalExpansion(0.0),
_thermalCapacity(0.0),
_thermalConductivity(0.0)
{
    if      (cl.find("mechanical neohookean") != cl.end())   theFSMaterial = new neohookeanMaterial(name, cl);
    else if (cl.find("mechanical arrudaboyce") != cl.end())  theFSMaterial = new arrudaboyceMaterial(name, cl);
    else if (cl.find("mechanical mooney") != cl.end())       theFSMaterial = new mooneyMaterial(name, cl);
    else if (cl.find("mechanical fplastic") != cl.end())     theFSMaterial = new fplasticMaterial(name, cl);
    else if (cl.find("mechanical svk") != cl.end())          theFSMaterial = new svkMaterial(name, cl);
    else if (cl.find("mechanical yeoh") != cl.end())         theFSMaterial = new yeohMaterial(name, cl);
    else if (cl.find("mechanical johnsoncook") != cl.end())  theFSMaterial = new johnsonCookMaterial(name, cl);
    else if (cl.find("mechanical zerilli") != cl.end())      theFSMaterial = new zerilliArmstrongMaterial(name, cl);
    else if (cl.find("mechanical arrhenius") != cl.end())    theFSMaterial = new arrheniusTypeMaterial(name, cl);
    //else if (cl.find("mechanical extended") != cl.end())     theFSMaterial = new extendedMaterial(name, cl);
    
    muesli::assignValue(cl, "heat_supply", _heatSupply);
    muesli::assignValue(cl, "thermal_expansion", _thermalExpansion);
    muesli::assignValue(cl, "heat_capacity", _thermalCapacity);
    muesli::assignValue(cl, "conductivity", _thermalConductivity);
}




thermofiniteStrainMaterial::thermofiniteStrainMaterial(const std::string& name, const materialProperties& cl,
                                                       const std::string& pureThermo)
:
material(name, cl),
theFSMaterial(0),
_heatSupply(0.0),
_thermalExpansion(0.0),
_thermalCapacity(0.0),
_thermalConductivity(0.0)
{
    muesli::assignValue(cl, "heat_supply", _heatSupply);
    muesli::assignValue(cl, "thermal_expansion", _thermalExpansion);
    muesli::assignValue(cl, "heat_capacity", _thermalCapacity);
    muesli::assignValue(cl, "conductivity", _thermalConductivity);
}




thermofiniteStrainMP* thermofiniteStrainMaterial::createMaterialPoint() const
{
    muesli::thermofiniteStrainMP* mp = new thermofiniteStrainMP(*this);
    return mp;
}




bool thermofiniteStrainMaterial::check() const
{
    bool ret = true;
    
    if (_thermalConductivity <= 0.0)
    {
        ret = false;
    }

    return ret;
}




double thermofiniteStrainMaterial::getProperty(const propertyName p) const
{
    double ret=0.0;
    
    // scan all the possible data
    switch (p)
    {
        case PR_CONDUCTIVITY: ret = _thermalConductivity; break;
        case PR_THERMAL_EXP:  ret = _thermalExpansion; break;
        case PR_THERMAL_CAP:  ret = _thermalCapacity; break;
        case PR_HEAT_SUPPLY:  ret = _heatSupply; break;

        default:
            ret = theFSMaterial->getProperty(p);
    }
    return ret;
}




bool thermofiniteStrainMaterial::hasDamageModel() const
{
    return theFSMaterial->hasDamageModel();
}




bool thermofiniteStrainMaterial::hasFiniteStrainMaterial() const
{
    return true;
}




void thermofiniteStrainMaterial::print(std::ostream &of) const
{
    of  << "\n Thermomechanical material"
        << "\n Mechanical part:";
    theFSMaterial->print(of);
    
    of  << "\n Thermal and coupling properties:"
        << "\n   Conductivity           : " << _thermalConductivity
        << "\n   Heat supply            : " << _heatSupply
        << "\n   Thermal expansion coef : " << _thermalExpansion
        << "\n   Heat capacity          : " << _thermalCapacity
        << "\n";
}




void thermofiniteStrainMaterial::setRandom()
{
    int mattype = discreteUniform(0, 5);
    std::string name = "surrogate finite strain material";
    materialProperties mp;
    
    if (mattype == 0)
        theFSMaterial = new neohookeanMaterial(name, mp);

    else if (mattype == 1)
        theFSMaterial = new fplasticMaterial(name, mp);
    
    else if (mattype == 2)
        theFSMaterial = new svkMaterial(name, mp);

    else if (mattype == 3)
        theFSMaterial = new arrudaboyceMaterial(name, mp);
    
    else if (mattype == 4)
        theFSMaterial = new mooneyMaterial(name, mp);
    
    else if (mattype == 5)
        theFSMaterial = new yeohMaterial(name, mp);
    
    else if (mattype == 6)
        theFSMaterial = new johnsonCookMaterial(name, mp);
    
    else if (mattype == 7)
        theFSMaterial = new zerilliArmstrongMaterial(name, mp);
    
    theFSMaterial->setRandom();
    _thermalConductivity = muesli::randomUniform(1.0, 10.0);
    _thermalExpansion    = 1e-6 * muesli::randomUniform(1.0, 10.0);
    _thermalCapacity     = muesli::randomUniform(10.0, 100.0);
}




bool thermofiniteStrainMaterial::test(std::ostream &of)
{
    bool isok = true;
    setRandom();
    muesli::thermofiniteStrainMP* p = this->createMaterialPoint();
    
    isok = p->testImplementation(of);
    return isok;
}




bool thermofiniteStrainMaterial::toBeDeleted() const
{
    return theFSMaterial->toBeDeleted();
}




double thermofiniteStrainMaterial::waveVelocity() const
{
    return theFSMaterial->waveVelocity();
}




thermofiniteStrainMP::thermofiniteStrainMP(const thermofiniteStrainMaterial& m) :
thethermoFiniteStrainMaterial(m),
theFSMP(0)
{
    if (m.hasFiniteStrainMaterial())
    {
        theFSMP=m.theFSMaterial->createMaterialPoint();
    }

    temp_c = m.referenceTemperature();
    temp_n = temp_c;
    GradT_n.setZero();
    GradT_c.setZero();
}




void thermofiniteStrainMP::CauchyStress(istensor& sigma) const
{
    itensor& Fc = theFSMP->deformationGradient();
    double J    = Fc.determinant();
    double Jinv = 1/J;
    istensor S;
    secondPiolaKirchhoffStress(S);
    sigma = Jinv * istensor::FSFt(Fc, S);
}




void thermofiniteStrainMP::CauchyStressVector(double S[6]) const
{
    istensor sigma;
    CauchyStress(sigma);
    muesli::ContraContraSymTensorToVector(sigma, S);
}




void thermofiniteStrainMP::commitCurrentState()
{
    theFSMP->commitCurrentState();

    time_n  = time_c;
    GradT_n = GradT_c;
    temp_n  = temp_c;
}




void thermofiniteStrainMP::contractWithAllTangents(const ivector &v1,
                                                     const ivector& v2,
                                                     itensor&  Tdev,
                                                     istensor& Tmixed,
                                                     double&   Tvol) const
{
    contractWithDeviatoricTangent(v1, v2, Tdev);
    contractWithMixedTangent(Tmixed);
    Tvol = volumetricStiffness();
}




void thermofiniteStrainMP::contractWithConvectedTangent(const ivector& v1, const ivector& v2, itensor& T) const
{
    itensor4 c;
    convectedTangent(c);
    
    T.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++)
                {
                    T(i,k) += c(i,j,k,l)*v1(j)*v2(l);
                }
}




void thermofiniteStrainMP::contractWithDeviatoricTangent(const ivector& v1, const ivector& v2, itensor& T) const
{
    const istensor id = istensor::identity();
    itensor4 Pdev;
    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++)
                    Pdev(i,j,k,l) = 0.5 * ( id(i,k)*id(j,l) + id(i,l)*id(j,k) ) - 1.0/3.0 * id(i,j)*id(k,l);
    
    itensor4 st;
    spatialTangent(st);
    
    itensor4 cdev;
    cdev.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++)
                    for (unsigned m=0; m<3; m++)
                        for (unsigned n=0; n<3; n++)
                            for (unsigned p=0; p<3; p++)
                                for (unsigned q=0; q<3; q++)
                                    cdev(i,j,p,q) += Pdev(i,j,k,l)*st(k,l,m,n)*Pdev(m,n,p,q);
    
    T.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++)
                {
                    T(i,k) += cdev(i,j,k,l)*v1(j)*v2(l);
                }
}




// CM_ij = P_ijkl*c_klmm
void thermofiniteStrainMP::contractWithMixedTangent(istensor& CM) const
{
    const istensor id = istensor::identity();
    itensor4 Pdev;
    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++)
                    Pdev(i,j,k,l) = 0.5 * ( id(i,k)*id(j,l) + id(i,l)*id(j,k) ) - 1.0/3.0 * id(i,j)*id(k,l);
    
    itensor4 st;
    spatialTangent(st);
    
    CM.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++)
                    for (unsigned m=0; m<3; m++)
                        CM(i,j) += Pdev(i,j,k,l)*st(k,l,m,m);
}




void thermofiniteStrainMP::contractWithSpatialTangent(const ivector& v1, const ivector& v2, itensor& T) const
{
    itensor4 c;
    spatialTangent(c);

    T.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++)
                {
                    T(i,k) += c(i,j,k,l)*v1(j)*v2(l);
                }
}




void thermofiniteStrainMP::contractTangent(const ivector& na, const ivector& nb, double& tg) const
{
    tg = na.dot(thethermoFiniteStrainMaterial._thermalConductivity * nb);
}




void thermofiniteStrainMP::convectedTangent(itensor4& ctg) const
{
    theFSMP->convectedTangent(ctg);
    
    itensor& Fc = theFSMP->deformationGradient();
    istensor C = istensor::tensorTransposedTimesTensor(Fc);
    istensor Cinv = C.inverse();

    const double alpha = thethermoFiniteStrainMaterial._thermalExpansion;
    const double bulk  = thethermoFiniteStrainMaterial.theFSMaterial->getProperty(muesli::PR_BULK);
    const double dtheta= temp_c - thethermoFiniteStrainMaterial.referenceTemperature();
    const double f     = 6.0 * alpha * bulk * dtheta;

    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++)
                {
                    ctg(i,j,k,l) += f * 0.5 * ( Cinv(i,k)* Cinv(l,j) + Cinv(i,l) * Cinv(k,j) );
                }
}




void thermofiniteStrainMP::convectedTangentMatrix(double C[6][6]) const
{
    itensor4 ct;
    convectedTangent(ct);
    muesli::tensorToMatrix(ct, C);
}




void thermofiniteStrainMP::convectedTangentTimesSymmetricTensor(const istensor &M, istensor &CM) const
{
    itensor4 C;
    convectedTangent(C);

	CM.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++)
                {
                    CM(i,j) += C(i,j,k,l)*M(k,l);
                }
}




double thermofiniteStrainMP::dissipatedEnergy() const
{
    return theFSMP->dissipatedEnergy();
}




double thermofiniteStrainMP::dissipation() const
{
    double id =   dissipatedEnergy();
    double td = - materialHeatflux().dot(GradT_c) / temp_c;
    return id + td;
}




istensor thermofiniteStrainMP::dissipationTangent() const
{
    istensor zero;
    return zero; //XX theFSMP->dissipationTangent();
}




double thermofiniteStrainMP::effectiveFreeEnergy() const
{
    itensor& Fc = theFSMP->deformationGradient();
    double J    = Fc.determinant();
    double logJ = log(J);

    const double alpha  = thethermoFiniteStrainMaterial._thermalExpansion;
    const double c0     = thethermoFiniteStrainMaterial._thermalCapacity;
    const double bulk   = thethermoFiniteStrainMaterial.theFSMaterial->getProperty(muesli::PR_BULK);
    const double theta0 = thethermoFiniteStrainMaterial.referenceTemperature();
    const double deltaTemp = temp_c - theta0;

    double Psi_coup = - 3.0 * alpha * bulk * deltaTemp * logJ;
    double Psi_ther = c0 * (deltaTemp - temp_c * log(temp_c/theta0));

    return theFSMP->effectiveStoredEnergy() + Psi_coup + Psi_ther;
}




/*void thermofiniteStrainMP::elasticTangent(itensor4& ce) const
{
    const double E = thethermoFiniteStrainMaterial.theFSMaterial->getProperty(muesli::PR_YOUNG);
    const double nu = thethermoFiniteStrainMaterial.theFSMaterial->getProperty(muesli::PR_NU);
    const double lambda = E*nu/(1.0-2.0*nu)/(1.0+nu);
    const double mu     = E/2.0/(1.0+nu);

    itensor& Fc = theFSMP->deformationGradient();
    itensor4 cc;
    cc.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) cc(i,j,k,l) += lambda;
                    if (i==k && j==l) cc(i,j,k,l) += mu;
                    if (i==l && j==k) cc(i,j,k,l) += mu;
                }
    
    ce.setZero();
    for (unsigned a=0; a<3; a++)
        for (unsigned b=0; b<3; b++)
            for (unsigned A=0; A<3; A++)
                for (unsigned B=0; B<3; B++)
                {
                    for (unsigned C=0; C<3; C++)
                        for (unsigned D=0; D<3; D++)
                            ce(a,A,b,B) += Fc(a,C) * Fc(b,D) * cc(C,A,D,B);
                }
}*/




//void thermofiniteStrainMP::energyMomentumTensor(itensor &EM) const
//{
//    theFSMP->energyMomentumTensor(EM);
//}




double thermofiniteStrainMP::entropy() const
{
    itensor& Fc = theFSMP->deformationGradient();
    double J = Fc.determinant();
    double logJ = log(J);
    
    const double alpha  = thethermoFiniteStrainMaterial._thermalExpansion;
    const double bulk   = thethermoFiniteStrainMaterial.theFSMaterial->getProperty(muesli::PR_BULK);
    const double c0     = thethermoFiniteStrainMaterial._thermalCapacity;
    const double theta0 = thethermoFiniteStrainMaterial.referenceTemperature();

    return 3.0 * alpha * bulk * logJ + c0 * log(temp_c/theta0);
}




void thermofiniteStrainMP::firstPiolaKirchhoffStress(itensor &P) const
{
    istensor S;
    secondPiolaKirchhoffStress(S);

    itensor& Fc = theFSMP->deformationGradient();
    P = Fc*S;
}




double thermofiniteStrainMP::freeEnergy() const
{
    itensor& Fc = theFSMP->deformationGradient();
    double J    = Fc.determinant();
    double logJ = log(J);
    
    const double alpha = thethermoFiniteStrainMaterial._thermalExpansion;
    const double c0    = thethermoFiniteStrainMaterial._thermalCapacity;
    const double bulk  = thethermoFiniteStrainMaterial.theFSMaterial->getProperty(muesli::PR_BULK);
    const double theta0= thethermoFiniteStrainMaterial.referenceTemperature();
    const double deltaTemp = temp_c - theta0;

    double Psi_coup = -3.0 * alpha * bulk * deltaTemp * logJ;
    double Psi_ther = c0 * (deltaTemp - temp_c * log(temp_c/theta0));

    return theFSMP->storedEnergy() + Psi_coup + Psi_ther;
}




double thermofiniteStrainMP::freeEntropy() const
{
    return -freeEnergy()/temp_c;
}




materialState thermofiniteStrainMP::getConvergedState() const
{
    materialState state;

    state = theFSMP->getConvergedState();
    state.theTime = time_n;
    state.theDouble.push_back(temp_n);
    state.theVector.push_back(GradT_n);

    return state;
}




materialState thermofiniteStrainMP::getCurrentState() const
{
    materialState state;
    
    state = theFSMP->getCurrentState();
    state.theTime = time_c;
    state.theDouble.push_back(temp_c);
    state.theVector.push_back(GradT_c);
    
    return state;
}




double thermofiniteStrainMP::getDamage() const
{
    return theFSMP->getDamage();
}




ivector& thermofiniteStrainMP::gradT()
{
    return GradT_c;
}




const ivector& thermofiniteStrainMP::gradT() const
{
    return GradT_c;
}




// c = - theta d^2[Psi]/d[theta]^2 Note: specific, to be used rho * c
double thermofiniteStrainMP::heatCapacity() const
{
    return thethermoFiniteStrainMaterial._thermalCapacity;
}




double thermofiniteStrainMP::internalEnergy() const
{
    return freeEnergy() + temp_c * entropy();
}




bool thermofiniteStrainMP::isFullyDamaged() const
{
    return theFSMP->isFullyDamaged();
}




double thermofiniteStrainMP::kineticPotential() const
{
    return theFSMP->kineticPotential() - thermalPotential();
}




void thermofiniteStrainMP::KirchhoffStress(istensor &tau) const
{
    theFSMP->KirchhoffStress(tau);
}




void thermofiniteStrainMP::KirchhoffStressVector(double tauv[6]) const
{
    theFSMP->KirchhoffStressVector(tauv);
}




istensor thermofiniteStrainMP::materialConductivity() const
{
    return istensor::scaledIdentity(thethermoFiniteStrainMaterial._thermalConductivity);
}




// coupling tensor M = Theta * d^2(psi)/(d_F d_Theta) = Theta d_P/d_Theta
itensor thermofiniteStrainMP::materialCouplingTensor() const
{
    itensor& Fc = theFSMP->deformationGradient();
    itensor Fcinv = Fc.inverse();

    const double alpha = thethermoFiniteStrainMaterial._thermalExpansion;
    const double bulk  = thethermoFiniteStrainMaterial.theFSMaterial->getProperty(muesli::PR_BULK);

    return -3.0 * alpha * bulk * temp_c * Fcinv.transpose();
}




itensor thermofiniteStrainMP::materialCouplingTensorDTheta() const
{
    itensor& Fc = theFSMP->deformationGradient();
    itensor Fcinv = Fc.inverse();

    const double alpha = thethermoFiniteStrainMaterial._thermalExpansion;
    const double bulk  = thethermoFiniteStrainMaterial.theFSMaterial->getProperty(muesli::PR_BULK);

    return -3.0 * alpha * bulk * Fcinv.transpose();
}




itensor4 thermofiniteStrainMP::materialCouplingTensorDF() const
{
    itensor& Fc = theFSMP->deformationGradient();
    itensor Fcinv = Fc.inverse();

    const double alpha = thethermoFiniteStrainMaterial._thermalExpansion;
    const double bulk  = thethermoFiniteStrainMaterial.theFSMaterial->getProperty(muesli::PR_BULK);
    const double f     = 3.0 * alpha * bulk * temp_c;

    itensor4 F_t_F;
    F_t_F.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++)
                {
                    F_t_F(i,j,k,l) += f*Fcinv(l,i)*Fcinv(j,k);
                }

    return F_t_F;
}




ivector thermofiniteStrainMP::materialHeatflux() const
{
    return -thethermoFiniteStrainMaterial._thermalConductivity * GradT_c;
}




void thermofiniteStrainMP::materialTangent(itensor4& cm) const
{
    istensor S;
    secondPiolaKirchhoffStress(S);

    itensor4 cc;
    convectedTangent(cc);

    const itensor& Fc = theFSMP->deformationGradient();
    cm.setZero();
    for (unsigned a=0; a<3; a++)
        for (unsigned b=0; b<3; b++)
            for (unsigned A=0; A<3; A++)
                for (unsigned B=0; B<3; B++)
                {
                    if (a == b) cm(a,A,b,B) += S(A,B);

                    for (unsigned C=0; C<3; C++)
                        for (unsigned D=0; D<3; D++)
                            cm(a,A,b,B) += Fc(a,C) * Fc(b,D) * cc(C,A,D,B);
                }
}




double thermofiniteStrainMP::plasticSlip() const
{
    return theFSMP->plasticSlip();
}




void thermofiniteStrainMP::resetCurrentState()
{
    theFSMP->resetCurrentState();

    GradT_c = GradT_n;
    temp_c  = temp_n;
}




void thermofiniteStrainMP::secondPiolaKirchhoffStress(istensor& S) const
{
    theFSMP->secondPiolaKirchhoffStress(S);

    itensor  Finv = theFSMP->deformationGradient().inverse();
    istensor Cinv = istensor::tensorTimesTensorTransposed(Finv);

    const double alpha  = thethermoFiniteStrainMaterial._thermalExpansion;
    const double bulk   = thethermoFiniteStrainMaterial.theFSMaterial->getProperty(muesli::PR_BULK);
    const double theta0 = thethermoFiniteStrainMaterial.referenceTemperature();
    const double dTemp  = temp_c - theta0;

    S -= 3.0 * alpha * bulk * dTemp * Cinv;
}




void  thermofiniteStrainMP::secondPiolaKirchhoffStressVector(double Sv[6]) const
{
    istensor S;
    secondPiolaKirchhoffStress(S);
    muesli::tensorToVector(S, Sv);
}




istensor thermofiniteStrainMP::spatialConductivity() const
{
    itensor& Fc = theFSMP->deformationGradient();
    double J = Fc.determinant();
    double Jinv = 1.0/J;
    istensor k = thethermoFiniteStrainMaterial._thermalConductivity*istensor::identity();
    return Jinv * istensor::FSFt(Fc, k);
}




ivector thermofiniteStrainMP::spatialHeatflux() const
{
    itensor& Fc     = theFSMP->deformationGradient();
    const ivector H = materialHeatflux();
        
    double J = Fc.determinant();
    return (1.0/J) * Fc * 0.5 * H;
}




void thermofiniteStrainMP::spatialTangent(itensor4& Cs) const
{
    /*COMMENTED JUST TO TEST THAT IMPLICIT SIMS WORK WITH SPATIAL TANGENT!
    itensor& Fc = theFSMP->deformationGradient();
    double J = Fc.determinant();
    itensor4 Cc;
    convectedTangent(Cc);
    
    Cs.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++)
                    for (unsigned m=0; m<3; m++)
                        for (unsigned n=0; n<3; n++)
                            for (unsigned p=0; p<3; p++)
                                for (unsigned q=0; q<3; q++)
                                    Cs(i,j,k,l) += Fc(i,m)*Fc(j,n)*Fc(k,p)*Fc(l,q)*Cc(m,n,p,q);
    
    Cs *= 1.0/J;*/
    
    theFSMP->spatialTangent(Cs);
    
    itensor& Fc = theFSMP->deformationGradient();
    itensor Fcinv = Fc.inverse();

    const double alpha = thethermoFiniteStrainMaterial._thermalExpansion;
    const double bulk  = thethermoFiniteStrainMaterial.getProperty(muesli::PR_BULK);
    const double dtheta= temp_c - thethermoFiniteStrainMaterial.referenceTemperature();
    const double f     = -3.0 * alpha * bulk * dtheta;

    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++)
                {
                    Cs(i,j,k,l) += -f * ( Fcinv(l,i)*Fcinv(j,k) );
                }
}


 

void thermofiniteStrainMP::spatialTangentMatrix(double c[6][6]) const
{
    itensor4 st;
    spatialTangent(st);
    muesli::tensorToMatrix(st, c);
}




// coupling tensor M = 2.0 * theta * d^2(psi)/( dC dTheta ) = theta * d[S]/d[Theta]
istensor thermofiniteStrainMP::symmetricCouplingTensor() const
{
    itensor& Fc = theFSMP->deformationGradient();
    istensor C  = istensor::tensorTransposedTimesTensor(Fc);
    istensor Cinv = C.inverse();
    itensor Fcinv = Fc.inverse();

    const double alpha = thethermoFiniteStrainMaterial._thermalExpansion;
    const double bulk  = thethermoFiniteStrainMaterial.theFSMaterial->getProperty(muesli::PR_BULK);

    return -3.0/2.0 * alpha * bulk * temp_c * Cinv;
}




// D = 4 d^2[psi]/d C^2
itensor4 thermofiniteStrainMP::symmetricCouplingTensorDC() const
{
    itensor& Fc = theFSMP->deformationGradient();
    istensor C = istensor::tensorTransposedTimesTensor(Fc);
    istensor Cinv = C.inverse();

    const double alpha = thethermoFiniteStrainMaterial._thermalExpansion;
    const double bulk  = thethermoFiniteStrainMaterial.theFSMaterial->getProperty(muesli::PR_BULK);
    const double f     = 6.0 * alpha * bulk * temp_c;

    itensor4 D;
    D.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++)
                {
                    D(i,j,k,l) += f * 0.5 * ( Cinv(i,k)* Cinv(j,l) + Cinv(i,l)* Cinv(j,k) );
                }
    return D;
}




istensor thermofiniteStrainMP::symmetricCouplingTensorDTheta() const
{
    itensor& Fc = theFSMP->deformationGradient();
    istensor C  = istensor::tensorTransposedTimesTensor(Fc);
    istensor Cinv = C.inverse();

    const double alpha = thethermoFiniteStrainMaterial._thermalExpansion;
    const double bulk  = thethermoFiniteStrainMaterial.theFSMaterial->getProperty(muesli::PR_BULK);
    return -3.0/2.0 * alpha * bulk * Cinv;
}




double& thermofiniteStrainMP::temperature()
{
    return temp_c;
}




const double& thermofiniteStrainMP::temperature() const
{
    return temp_c;
}




bool thermofiniteStrainMP::testImplementation(std::ostream& of, const bool testDE, const bool testDDE) const
{
    bool isok = true;
    const double inc = 1.0e-4;
    thermofiniteStrainMP& theMP = const_cast<thermofiniteStrainMP&>(*this);

    // set a random update in the material
    itensor F;  F.setRandom();
    if (F.determinant() < 0.0) F *= -1.0;
    ivector gradT; gradT.setRandom();
    double temp = randomUniform(2.0, 2.1) * thethermoFiniteStrainMaterial.referenceTemperature();

    theMP.updateCurrentState(0.0, F, gradT, temp);
    theMP.commitCurrentState();

    double tn1 = muesli::randomUniform(0.1,1.0);
    F.setRandom();
    if (F.determinant() < 0.0) F *= -1.0;
    gradT.setRandom();
    temp = randomUniform(2.1, 2.2) * thethermoFiniteStrainMaterial.referenceTemperature();
    theMP.updateCurrentState(tn1, F, gradT, temp);

    // check programmed quantities by numerically differentiating the free energy
    theMP.resetCurrentState();

    // Derivatives with respect to F
    // Compute numerical value of P, 1PK stress tensor, = d Psi / d F
    // Compute numerical value of A, material tangent A_{iAjB} = d (P_iA) / d F_jB
    // Compute numerical value of d M / d F
    itensor num_P;
    itensor4 num_A, num_dMdF;
    {
        itensor dP, Pp1, Pp2, Pm1, Pm2, dM;

        for (size_t i=0; i<3; i++)
        {
            for (size_t j=0; j<3; j++)
            {
                const double original = F(i,j);

                F(i,j) = original + inc;
                theMP.updateCurrentState(tn1, F, gradT, temp);
                double Wp1 = effectiveFreeEnergy();
                firstPiolaKirchhoffStress(Pp1);
                itensor dMdF_p1 = materialCouplingTensor();

                F(i,j) = original + 2.0*inc;
                theMP.updateCurrentState(tn1, F, gradT, temp);
                double Wp2 = effectiveFreeEnergy();
                firstPiolaKirchhoffStress(Pp2);
                itensor dMdF_p2 = materialCouplingTensor();

                F(i,j) = original - inc;
                theMP.updateCurrentState(tn1, F, gradT, temp);
                double Wm1 = effectiveFreeEnergy();
                firstPiolaKirchhoffStress(Pm1);
                itensor dMdF_m1 = materialCouplingTensor();

                F(i,j) = original - 2.0*inc;
                theMP.updateCurrentState(tn1, F, gradT, temp);
                double Wm2 = effectiveFreeEnergy();
                firstPiolaKirchhoffStress(Pm2);
                itensor dMdF_m2 = materialCouplingTensor();

                // fourth order approximation of the derivative
                num_P(i,j) = (-Wp2 + 8.0*Wp1 - 8.0*Wm1 + Wm2)/(12.0*inc);

                // derivative of PK stress
                dP = (-Pp2 + 8.0*Pp1 - 8.0*Pm1 + Pm2)/(12.0*inc);
                dM = (-dMdF_p2 + 8.0*dMdF_p1 - 8.0*dMdF_m1 + dMdF_m2)/(12.0*inc);
                for (unsigned k=0; k<3; k++)
                    for (unsigned l=0; l<3; l++)
                    {
                        num_A(k,l,i,j) = dP(k,l);
                        num_dMdF(k,l,i,j) = dM(k,l);
                    }

                F(i,j) = original;
                theMP.updateCurrentState(tn1, F, gradT, temp);
            }
        }
    }


    // Derivatives with respect to theta
    // Compute numerical value of entropy = - d Psi / d theta
    // Compute numerical value of coupling tensor M = theta * d P / d theta
    // Compute numerical value of heat capacity =  theta * (d s / d theta)
    // Compute numerical value of d M / d theta
    double num_entropy;
    itensor num_coupling;
    itensor num_dMdt;
    double num_c;
    {
        const double original = temp;

        temp = original + inc;
        theMP.updateCurrentState(tn1, F, gradT, temp);
        double Psip1 = effectiveFreeEnergy();
        itensor PK_p1; firstPiolaKirchhoffStress(PK_p1);
        double s_p1 = entropy();
        itensor M_p1 = materialCouplingTensor();

        temp = original + 2.0*inc;
        theMP.updateCurrentState(tn1, F, gradT, temp);
        double Psip2 = effectiveFreeEnergy();
        itensor PK_p2; firstPiolaKirchhoffStress(PK_p2);
        double s_p2 = entropy();
        itensor M_p2 = materialCouplingTensor();

        temp = original - inc;
        theMP.updateCurrentState(tn1, F, gradT, temp);
        double Psim1 = effectiveFreeEnergy();
        itensor PK_m1; firstPiolaKirchhoffStress(PK_m1);
        double s_m1 = entropy();
        itensor M_m1 = materialCouplingTensor();

        temp = original - 2.0*inc;
        theMP.updateCurrentState(tn1, F, gradT, temp);
        double Psim2 = effectiveFreeEnergy();
        itensor PK_m2; firstPiolaKirchhoffStress(PK_m2);
        double s_m2 = entropy();
        itensor M_m2 = materialCouplingTensor();

        temp = original;
        theMP.updateCurrentState(tn1, F, gradT, temp);

        double theta = temp_c;

        // fourth order approximation of the derivative
        num_entropy  = -(-Psip2 + 8.0*Psip1 - 8.0*Psim1 + Psim2)/(12.0*inc);
        num_coupling = theta * (-PK_p2 + 8.0*PK_p1 - 8.0*PK_m1 + PK_m2)/(12.0*inc);
        num_c        = theta * (-s_p2 + 8.0*s_p1 - 8.0*s_m1 + s_m2)/(12.0*inc);
        num_dMdt     = (-M_p2 + + 8.0*M_p1 - 8.0*M_m1 + M_m2)/(12.0*inc);
    }


    // compare 1st PK stress with derivative of free energy wrt F
    if (testDE)
    {
        itensor pr_P;
        firstPiolaKirchhoffStress(pr_P);

        itensor errorP = num_P - pr_P;
        isok = (errorP.norm()/pr_P.norm() < 1e-4);
        of << "\n   1. Comparing P with derivative [d Psi / d F].";
        if (isok)
        {
            of << " Test passed.";
        }
        else
        {
            of << "\n Relative error in 1st PK computation %e. Test failed." << errorP.norm()/pr_P.norm();
            of << "\n " << pr_P;
            of << "\n " << num_P;
        }
    }
    else
    {
        of << "\n   2. Comparing P with derivative of Psi ::: not run for this material";
    }
    
    // test the consistency of the stress tensors sigma and P
    {
        istensor sigma, S;
        itensor pr_P;
        CauchyStress(sigma);
        firstPiolaKirchhoffStress(pr_P);
        secondPiolaKirchhoffStress(S);
        
        itensor P1 = F.determinant() * sigma * F.inverse().transpose();
        itensor P2 = F*S;
        
        itensor errorP1 = P1 - pr_P;
        itensor errorP2 = P2 - pr_P;
        double  error = errorP1.norm() + errorP2.norm();
        isok = error/pr_P.norm() < 1e-4;
        of << "\n   3. Checking consistency of the 1PK, 2PK, and Cauchy stresses.";
        if (isok)
        {
            of << " Test passed.";
        }
        else
        {
            of << "\n Test failed. Relative error: " << error/pr_P.norm();
            of << "\n P: " << pr_P;
            of << "\n P1: " << P1;
            of << "\n P2: " << P2;
            of << "\n numP: " << num_P;
            of << "\n Cauchy: " << sigma;
        }
        of << std::flush;
    }

    // compare entropy with derivative of Psi wrt theta (-)
    {
        double pr_entropy = entropy();

        double error = num_entropy - pr_entropy;
        isok = (fabs(error)/fabs(pr_entropy) < 1e-4);
        of << "\n   4. Comparing entropy with  [- d Psi / d theta ].";
        if (isok)
        {
            of << " Test passed.";
        }
        else
        {
            of << "\n Relative error entropy computation %e. Test failed." << fabs(error)/fabs(pr_entropy);
            of << "\n " << pr_entropy;
            of << "\n " << num_entropy;
        }
    }


    // compare material tangent with derivative of the 1st PK w.r.t. F
    if (testDDE)
    {
        // programmed material tangent
        itensor4 pr_A;
        materialTangent(pr_A);

        // relative error
        itensor4 errorA = num_A - pr_A;
        double error = errorA.norm();
        double norm = pr_A.norm();
        isok  = (error/norm < 1e-4);

        of << "\n   5. Comparing material tangent with [d 1PK / d F ].";
        if (isok)
        {
            of << " Test passed.";
        }
        else
        {
            of << "\n      Test failed.";
            of << "\n      Relative error in DStress computation: " <<  error/norm;
            of << "\n      Error 4th order tensor:\n";
            of << errorA;
            of << "\n      Programmed tangent:\n";
            of << pr_A;
            of << "\n      Numerical tangent:\n";
            of << num_A;
        }
    }
    
    // test tangent as derivative of the stress
    if (testDDE)
    {
        // programmed convected tangent
        itensor4 tg;
        convectedTangent(tg);
        
        // numeric convected tangent
        itensor4 nTg;
        nTg.setZero();
        
        // transform num_A to get the convected tangent
        itensor  J  = F.inverse();
        istensor C  = istensor::tensorTransposedTimesTensor(F);
        istensor Ci = C.inverse();
        istensor S;
        secondPiolaKirchhoffStress(S);
        for (unsigned a=0; a<3; a++)
            for (unsigned b=0; b<3; b++)
                for (unsigned c=0; c<3; c++)
                    for (unsigned d=0; d<3; d++)
                    {
                        nTg(c,a,d,b) = - S(a,b)*Ci(c,d);
                        
                        for (unsigned i=0; i<3; i++)
                            for (unsigned j=0; j<3; j++)
                                nTg(c,a,d,b) += J(c,i)*num_A(i,a,j,b)*J(d,j);
                    }
        
        // relative
        double error = 0.0;
        double norm = 0.0;
        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++)
                    {
                        error += pow(nTg(i,j,k,l)-tg(i,j,k,l),2);
                        norm  += pow(tg(i,j,k,l),2);
                    }
        
        error = sqrt(error);
        norm = sqrt(norm);
        isok = (error/norm < 1e-4);
        
        of << "\n   6. Comparing convected tangent with derivative of stress.";
        if (isok)
        {
            of << " Test passed.";
        }
        else
        {
            of << "\n      Test failed.";
            of << "\n      Relative error in DStress computation: " <<  error/norm;
        }
        of << std::flush;
    }



    // test coupling tensor with theta * d P / d theta
    {
        itensor pr_coupling = materialCouplingTensor();

        itensor error = num_coupling - pr_coupling;
        isok = (error.norm()/pr_coupling.norm() < 1e-4);
        of << "\n   7. Comparing M with [theta d P / d theta].";
        if (isok)
        {
            of << " Test passed.";
        }
        else
        {
            of << "\n Relative error in coupling tensor computation. Test failed." << error.norm()/pr_coupling.norm();
            of << "\n " << pr_coupling;
            of << "\n " << num_coupling;
        }
    }


    // compare heat capacity with temperature * derivative of entropy wrt theta
    {
        double pr_c = heatCapacity();

        double error = num_c - pr_c;
        isok = (fabs(error)/fabs(pr_c) < 1e-4);
        of << "\n   8. Comparing heat capacity with  [theta * d s / d theta].";
        if (isok)
        {
            of << " Test passed.";
        }
        else
        {
            of << "\n Relative error c computation %e. Test failed." << fabs(error)/fabs(pr_c);
            of << "\n " << pr_c;
            of << "\n " << num_c;
        }
    }


    // compare derivative of coupling tensor
    {
        itensor dMdt = materialCouplingTensorDTheta();
        itensor error = num_dMdt - dMdt;
        isok = (error.norm()/dMdt.norm() < 1e-4);

        of << "\n   9. Comparing DM/dtheta with derivative of M.";
        if (isok)
        {
            of << " Test passed.";
        }
        else
        {
            of << "\n Relative error in DM/Dtheta. Test failed." << error.norm()/dMdt.norm();
            of << "\n " << dMdt;
            of << "\n " << num_dMdt;
        }
    }
    
    
    // compare derivative of coupling tensor w.r.t. F
    {
        itensor4 dMdF = materialCouplingTensorDF();
        itensor4 error = num_dMdF - dMdF;
        isok = (error.norm()/dMdF.norm() < 1e-4);
        of << "\n   10. Comparing DM/dF with derivative of M.";
        if (isok)
        {
            of << " Test passed.";
        }
        else
        {
            of << "\n Relative error in DM/DF. Test failed." << error.norm()/dMdF.norm();
            of << "\n " << dMdF;
            of << "\n " << num_dMdF;
        }
    }
    
    // tangent contracted functions
    if ((true))
    {
        itensor4 tg;
        convectedTangent(tg);
        
        istensor Csym;
        istensor sym;  sym.setZero();
        ivector r1; r1.setRandom();
        ivector r2; r2.setRandom();
        sym.addSymmetrizedDyadic(r1, r2);
        convectedTangentTimesSymmetricTensor(sym, Csym);
        
        istensor refCsym;   refCsym.setZero();
        itensor refCuv;     refCuv.setZero();
        itensor refCuv1;     refCuv1.setZero();
        itensor refCuv2;     refCuv2.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++)
                    {
                        // tangent times symmetric tensor
                        refCsym(i,j) += tg(i,j,k,l)*sym(k,l);
                        
                        // inner contraction of tangent with two vectors
                        refCuv(i,k) += tg(i,j,k,l)*r1(j)*r2(l);
                    }
        
        // relative error less than 0.01%
        istensor error = refCsym - Csym;
        isok = (error.norm()/Csym.norm() < 1e-4);
        of << "\n   11. Checking convected tangent times symmetric tensor.";
        if (isok)
        {
            of << " Test passed.";
        }
        else
        {
            of << "\n   Test failed.";
            of << "\n   Relative error: " << error.norm()/Csym.norm();
            of << "\n   C*sym function \n" << Csym;
            of << "\n   C*sym component-wise:\n" << refCsym;
        }

        
        // relative error less than 0.01%
        itensor Cuv;
        contractWithConvectedTangent(r1, r2, Cuv);
        itensor Cuverror = refCuv - Cuv;
        isok = (Cuverror.norm()/Cuv.norm() < 1e-4);
        
        of << "\n   12. Checking contract convected tangent.";
        if (isok)
        {
            of << " Test passed.";
        }
        else
        {
            of << "\n   Test failed.";
            of << "\n   Relative error: " << Cuverror.norm()/Cuv.norm();
            of << "\n   Cuv \n" << Cuv;
            of << "\n   Reference Cuv: \n" << refCuv;
        }
        of << std::flush;


        // checking other contractions
        itensor Cdev_uv, ref_Cdev_uv;
        contractWithDeviatoricTangent(r1, r2, Cdev_uv);
        thermofiniteStrainMP::contractWithDeviatoricTangent(r1, r2, ref_Cdev_uv);
        itensor Cdev_uverror = ref_Cdev_uv - Cdev_uv;
        double relerror = Cdev_uv.norm() > 1e-4 ? Cdev_uverror.norm()/Cdev_uv.norm() : Cdev_uverror.norm();
        isok = relerror < 1e-4;
        of << "\n   13. Checking contract with deviatoric tangent.";
        if (isok)
        {
            of << " Test passed.";
        }
        else
        {
            of << "\n   Test failed.";
            of << "\n   Relative error: " << relerror;
            of << "\n   Cuv \n" << Cdev_uv;
            of << "\n   Reference Cuv: \n" << ref_Cdev_uv;
        }
        of << std::flush;


        // checking other contractions
        istensor CM, CM_ref;
        contractWithMixedTangent(CM);
        thermofiniteStrainMP::contractWithMixedTangent(CM_ref);
        itensor CM_error = CM_ref - CM;
        relerror = CM.norm() > 1e-4 ? CM_error.norm()/CM.norm() : CM_error.norm();
        isok = (relerror < 1e-4);
        of << "\n   14. Checking contract with mixed tangent.";
        if (isok)
        {
            of << " Test passed.";
        }
        else
        {
            of << "\n   Test failed.";
            of << "\n   Relative error: " << relerror;
            of << "\n   Cuv \n" << CM;
            of << "\n   Reference Cuv: \n" << CM_ref;
        }
        of << std::flush;


        double kv = volumetricStiffness();
        double kvref = thermofiniteStrainMP::volumetricStiffness();
        double kv_error = fabs(kv-kvref);
        isok = (kv_error/fabs(kvref) < 1e-4);
        of << "\n   15. Checking volumetric stiffness.";
        if (isok)
        {
            of << " Test passed.";
        }
        else
        {
            of << "\n   Test failed.";
            of << "\n   Relative error: " << kv_error/fabs(kvref);
            of << "\n   k          : " << kv;
            of << "\n   Reference k: " << kvref;
        }
        of << std::flush;

    }


    {
        // check voigt maps
        istensor sigma;
        CauchyStress(sigma);
        double S[6];
        CauchyStressVector(S);
        
        double error = 0.0;
        for (unsigned i=0; i<6; i++)
        {
            error += pow( S[i] - sigma( voigt(0,i) , voigt(1,i) ), 2);
        }

        // relative error less than 0.01%
        isok = (error/sigma.norm() < 1e-4);

        of << "\n   16. Comparing Cauchy stress tensor and Voigt Cauchy stress.";
        if (isok)
        {
            of << " Test passed.";
        }
        else
        {
            of << "\n      Test failed.";
            of << "\n      Relative error in stress comparison: " << error/sigma.norm();
        }
        of << std::flush;
    }


    {
        itensor4 tg;
        convectedTangent(tg);
        
        double cm[6][6];
        convectedTangentMatrix(cm);
        
        double error = 0.0;
        double norm  = 0.0;
        for (unsigned i=0; i<6; i++)
        {
            error += pow(cm[i][i] - tg( voigt(0,i), voigt(1,i), voigt(0,i), voigt(1,i) ),2);
            norm  += pow(cm[i][i],2);
            for (unsigned j=i+1; j<6; j++)
            {
                error += pow(cm[i][j] - tg( voigt(0,i), voigt(1,i), voigt(0,j), voigt(1,j) ),2);
                error += pow(cm[j][i] - tg( voigt(0,i), voigt(1,i), voigt(0,j), voigt(1,j) ),2);
                norm  += pow(cm[i][j],2);
                norm  += pow(cm[j][i],2);
            }
        }
        error = sqrt(error);
        norm  = sqrt(norm);
        
        
        // relative error less than 0.01%
        isok = (error/norm < 1e-4);
        of << "\n   17. Comparing convected tangent tensor and Voigt matrix.";
        if (isok)
        {
            of << " Test passed.";
        }
        else
        {
            of << "\n      Test failed.";
            of << "\n      Relative error: " << error/norm;
        }
        of << std::flush;
    }

    return isok;
}




double thermofiniteStrainMP::thermalPotential() const
{
    double k = thethermoFiniteStrainMaterial._thermalConductivity;
    return 0.5 * k * GradT_c.squaredNorm();
}




double thermofiniteStrainMP::volumetricStiffness() const
{
    itensor4 tg;
    spatialTangent(tg);
    
    double vs = 0.0;
    for (unsigned i=0; i<3; i++)
    {
        for (unsigned j=0; j<3; j++)
        {
            vs += tg(i,i,j,j);
        }
    }
    return vs/9.0;
}




void thermofiniteStrainMP::updateCurrentState(const double theTime, const itensor& F, const ivector& GradT, const double& temp)
{
    theFSMP->setTemperature(temp_n);
    theFSMP->updateCurrentState(theTime, F);

    time_c  = theTime;
    GradT_c = GradT;
    temp_c  = temp;
}
