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

using namespace muesli;



fThermoMechMassMaterial::fThermoMechMassMaterial(const std::string& name,
                                                 const materialProperties& cl)
:
material(name, cl),
_variational(false),
_thermalExpansion(0.0),
_heatCapacity(0.0),
_massExpansion(0.0),
_referenceChemicalPotential(0.0),
_R(8.34)
{
    muesli::assignValue(cl, "thermal_expansion", _thermalExpansion);
    muesli::assignValue(cl, "heat_capacity", _heatCapacity);
    muesli::assignValue(cl, "mass_expansion", _massExpansion);
    muesli::assignValue(cl, "muref", _referenceChemicalPotential);
    muesli::assignValue(cl, "r", _R);
    _variational = muesli::hasKeyword(cl, "variational");
}




bool fThermoMechMassMaterial::check() const
{
    return true;
}




double fThermoMechMassMaterial::getProperty(const propertyName p) const
{
    return 0.0;
}




void fThermoMechMassMaterial::print(std::ostream &of) const
{
    of  << "\n Thermo - chemo - mechanical coupled material."
    << "\n Variational formulation: " << _variational
    << "\n Model constants:"
    << "\n    Thermal expansion coeff. : " << _thermalExpansion
    << "\n    Mass expansion coeff.    : " << _massExpansion
    << "\n    Heat capacity            : " << _heatCapacity
    << "\n    Reference chemical pot.  : " << _referenceChemicalPotential
    << "\n    Gas constant             : " << _R;
}




double fThermoMechMassMaterial::referenceChemicalPotential() const
{
    return _referenceChemicalPotential;
}




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

    _thermalExpansion = muesli::randomUniform(1.0, 5.0)*1e-3;
    _massExpansion    = muesli::randomUniform(1.0, 5.0)*1e-3;
    _heatCapacity     = muesli::randomUniform(1.0, 5.0)*1e-5;
    _R                = muesli::randomUniform(1.0, 10.0);
    _referenceChemicalPotential = muesli::randomUniform(1.0, 10.0);
    _variational      = true;
}




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




AnandIJSS2011Material::AnandIJSS2011Material(const std::string& name,
                                             const materialProperties& cl)
:
fThermoMechMassMaterial(name, cl),
theFSMaterial(nullptr),
_thermalConductivity(0.0),
_massDiffusivity(0.0),
_referenceConcentration(0.0),
_NM(0.0), _bulkModulus(0.0)
{
    theFSMaterial = new neohookeanMaterial(name, cl);
    
    muesli::assignValue(cl, "thermal_conductivity", _thermalConductivity);
    muesli::assignValue(cl, "themal_expansion", _thermalExpansion);
    muesli::assignValue(cl, "heat_capacity", _heatCapacity);

    muesli::assignValue(cl, "mass_diffusivity", _massDiffusivity);
    muesli::assignValue(cl, "mass_expansion", _massExpansion);

    muesli::assignValue(cl, "muref", _referenceChemicalPotential);

    muesli::assignValue(cl, "nm", _NM);
    if (theFSMaterial != nullptr) _bulkModulus = theFSMaterial->getProperty(PR_BULK);

    _referenceConcentration = fconcentration(1.0, referenceTemperature(), 0.0);
}




bool AnandIJSS2011Material::check() const
{
    return true;
}




fThermoMechMassMP* AnandIJSS2011Material::createMaterialPoint() const
{
    return new AnandIJSS2011MP(*this);
}




// computes the value of the concentration for a given state, not necessarily
// the current one of the point
double AnandIJSS2011Material::fconcentration(double J, double temp, double mu) const
{
    const double beta  = _massExpansion;
    const double R     = _R;
    const double NM    = _NM;
    const double kappa = _bulkModulus;
    const double logJ  = log(J);
    const double muref = referenceChemicalPotential();
    return NM*exp((mu - muref + 3.0*kappa*beta*logJ)/(R*temp));
}




double AnandIJSS2011Material::getProperty(const propertyName p) const
{
    return 0.0;
}




void AnandIJSS2011Material::print(std::ostream &of) const
{
    fThermoMechMassMaterial::print(of);
    of  << "\n    Thermal (isotropic) conductivity : " << _thermalConductivity
        << "\n    Mass (isotropic) diffusivity     : " << _massDiffusivity
        << "\n    Moles of host material/unit vol  : " << _NM
        << "\n    Reference concentration          : " << _referenceConcentration
        << "\n    Stored bulk modulus              : " << _bulkModulus
        << "\n";
    theFSMaterial->print(of);
}




void AnandIJSS2011Material::setRandom()
{
    fThermoMechMassMaterial::setRandom();

    int mattype = discreteUniform(0, 0);
    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.0e3, 1e5);
    _massDiffusivity     = muesli::randomUniform(1e-3, 1e-2);
    _referenceConcentration = muesli::randomUniform(1e-3, 1e-2);
    _NM                  = muesli::randomUniform(2.0, 8.0)*1e5;
    _bulkModulus         = theFSMaterial->getProperty(PR_BULK);
    _referenceChemicalPotential = -_R*referenceTemperature()*log(_referenceConcentration/_NM);
}




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

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




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




fThermoMechMassMP::fThermoMechMassMP(const fThermoMechMassMaterial& m) :
theThermoMechMassMaterial(m),
time_n(0.0), temp_n(m.referenceTemperature()), J_n(1.0), mu_n(0.0), c_n(0.0),
time_c(0.0), temp_c(m.referenceTemperature()), J_c(1.0), mu_c(0.0), c_c(0.0)
{
    gradT_n.setZero();
    gradT_c.setZero();
    gradMu_n.setZero();
    gradMu_c.setZero();
    F_n = itensor::identity();
    F_c = itensor::identity();
}




double fThermoMechMassMP::chemicalPotential() const
{
    return mu_c;
}




void fThermoMechMassMP::CauchyStress(istensor& sigma) const
{
    istensor S;
    secondPiolaKirchhoffStress(S);
    sigma = istensor::FSFt(F_c, S);
    sigma *= 1.0/J_c;
}




void fThermoMechMassMP::commitCurrentState()
{
    time_n   = time_c;
    temp_n   = temp_c;
    J_n      = J_c;
    mu_n     = mu_c;
    c_n      = c_c;
    gradMu_n = gradMu_c;
    F_n      = F_c;
}




double fThermoMechMassMP::concentration() const
{
    return c_c;
}




void fThermoMechMassMP::contractWithAllTangents(const ivector &v1,
                                                const ivector& v2,
                                                itensor&  Tdev,
                                                istensor& Tmixed,
                                                double&   Tvol) const
{

}




void fThermoMechMassMP::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 fThermoMechMassMP::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 fThermoMechMassMP::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 fThermoMechMassMP::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 fThermoMechMassMP::contractTangent(const ivector& na, const ivector& nb, double& tg) const
{

}




void fThermoMechMassMP::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);
                }
}




void fThermoMechMassMP::energyMomentumTensor(itensor &EM) const
{

}




void fThermoMechMassMP::firstPiolaKirchhoffStress(itensor &P) const
{
    istensor S;
    secondPiolaKirchhoffStress(S);
    P = F_c*S;
}




double fThermoMechMassMP::freeEntropy() const
{
    double tref = theThermoMechMassMaterial.referenceTemperature();
    return -freeEnergy()/tref;
}



// Legendre transform of the free energy wrt the concentration
// gives potential that is function of F, theta, mu
double fThermoMechMassMP::grandCanonicalPotential() const
{
    return freeEnergy() - c_c * mu_c;
}




double fThermoMechMassMP::heatCapacityPerUnitVolume() const
{
    return theThermoMechMassMaterial._heatCapacity;
}




double fThermoMechMassMP::internalEnergy() const
{
    return grandCanonicalPotential() + temp_c*entropy() + c_c*mu_c;
}




void fThermoMechMassMP::KirchhoffStress(istensor& tau) const
{
    istensor sigma;
    CauchyStress(sigma);
    tau = sigma * J_c;
}




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

    itensor4 cc;
    convectedTangent(cc);

    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) += F_c(a,C) * F_c(b,D) * cc(C,A,D,B);
                }
}




const fThermoMechMassMaterial& fThermoMechMassMP::parentMaterial() const
{
    return theThermoMechMassMaterial;
}




void fThermoMechMassMP::resetCurrentState()
{
    time_c   = time_n;
    temp_c   = temp_n;
    J_c      = J_n;
    mu_c     = mu_n;
    c_c      = c_n;
    gradMu_c = gradMu_n;
    F_c      = F_n;
}




void fThermoMechMassMP::spatialTangent(itensor4& Cs) const
{
    const itensor& Fc = F_c;
    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_c;
}




// coupling tensor M = 2.0 * theta * d^2(psi)/( dC dTheta ) = theta * d[S]/d[Theta]
istensor fThermoMechMassMP::symmetricMaterialStressTemperatureTensor() const
{
    itensor  Finv = F_c.inverse();
    istensor Cinv = istensor::tensorTimesTensorTransposed(Finv);

    const double alpha = theThermoMechMassMaterial._thermalExpansion;
    const double bulk  = theThermoMechMassMaterial.getProperty(muesli::PR_BULK);

    std::cout << "\n Needs review";

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




double fThermoMechMassMP::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;
}




double fThermoMechMassMP::waveVelocity() const
{
    return theThermoMechMassMaterial.waveVelocity();
}




AnandIJSS2011MP::AnandIJSS2011MP(const AnandIJSS2011Material& m):
fThermoMechMassMP(m),
theAnandMaterial(m)
{
    theFSMP = theAnandMaterial.theFSMaterial->createMaterialPoint();
    c_n = c_c = theAnandMaterial._referenceConcentration;
}




// d^2 [GCP] / d mu^2 = - d [ concentration] / d [ chemical potential ]
double AnandIJSS2011MP::chemicalTangent() const
{
    const double R = theAnandMaterial._R;
    
    return -concentration()/(R*temp_c);
}




void AnandIJSS2011MP::commitCurrentState()
{
    fThermoMechMassMP::commitCurrentState();
    theFSMP->commitCurrentState();
}




void AnandIJSS2011MP::convectedTangent(itensor4& ctg) const
{
    const double alpha = theAnandMaterial._thermalExpansion;
    const double beta  = theAnandMaterial._massExpansion;
    const double cref  = theAnandMaterial._referenceConcentration;
    const double kappa = theAnandMaterial._bulkModulus;
    const double tref  = theAnandMaterial.referenceTemperature();
    const double R     = theAnandMaterial._R;

    theFSMP->convectedTangent(ctg);

    istensor C = istensor::tensorTransposedTimesTensor(F_c);
    istensor Cinv = C.inverse();

    const double dtheta= temp_c - tref;
    const double f     = 6.0 * alpha * kappa * dtheta;

    const double dcon  = c_c - cref;
    const double g     = 6.0 * beta * kappa * dcon;
    const double gg    = 9.0 * beta * beta * kappa * kappa * c_c / R /temp_c;

    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+g) * 0.5 * ( Cinv(i,k)* Cinv(l,j) + Cinv(i,l) * Cinv(k,j) )
                                 - gg * Cinv(i,j) * Cinv(k,l);
                }
}



// if variational, return omega-bar, otherwise, omega
double AnandIJSS2011MP::diffusionPotential() const
{
    const double diff = theAnandMaterial._massDiffusivity;

    ivector D;
    if (theAnandMaterial._variational)
    {
        D = -temp_n/temp_c * gradMu_c;
    }
    else
    {
        D = -gradMu_c;
    }

    return -0.5*diff*D.squaredNorm();
}




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




double AnandIJSS2011MP::effectiveFreeEnergy() const
{
    return freeEnergy();
}




double AnandIJSS2011MP::entropy() const
{
    const double alpha = theAnandMaterial._thermalExpansion;
    const double cap   = theAnandMaterial._heatCapacity;
    const double R     = theAnandMaterial._R;
    const double NM    = theAnandMaterial._NM;
    const double kappa = theAnandMaterial._bulkModulus;
    const double tref  = theAnandMaterial.referenceTemperature();

    double s = cap*log(temp_c/tref) + 3.0*kappa*alpha*log(J_c) - R*c_c*(log(c_c/NM)-1.0);

    return s;
}



// free energy as a function of F, theta, c
double AnandIJSS2011MP::freeEnergy() const
{
    const double alpha = theAnandMaterial._thermalExpansion;
    const double beta  = theAnandMaterial._massExpansion;
    const double cap   = theAnandMaterial._heatCapacity;
    const double cref  = theAnandMaterial._referenceConcentration;
    const double R     = theAnandMaterial._R;
    const double muref = theAnandMaterial._referenceChemicalPotential;
    const double NM    = theAnandMaterial._NM;
    const double kappa = theAnandMaterial._bulkModulus;
    const double tref  = theAnandMaterial.referenceTemperature();
    const double dTemp = temp_c - tref;
    const double dCon  = c_c-cref;

    double logJ = log(J_c);

    double psi = theFSMP->effectiveStoredEnergy()
        - 3.0*kappa*alpha*dTemp*logJ
        - 3.0*kappa*beta*dCon*logJ
        + cap*dTemp - cap*temp_c*log(temp_c/tref)
        + muref*c_c + R*temp_c*c_c*(log(c_c/NM)-1.0);

    return psi;
}




materialState AnandIJSS2011MP::getConvergedState() const
{
    materialState state;
    
    state.theTime = time_n;
    state.theDouble.push_back(temp_n);
    state.theDouble.push_back(J_n);
    state.theDouble.push_back(mu_n);
    state.theDouble.push_back(c_n);
    state.theVector.push_back(gradMu_n);
    state.theTensor.push_back(F_n);
    
    return state;
}




materialState AnandIJSS2011MP::getCurrentState() const
{
    materialState state;
    
    state.theTime = time_c;
    state.theDouble.push_back(temp_c);
    state.theDouble.push_back(J_c);
    state.theDouble.push_back(mu_c);
    state.theDouble.push_back(c_c);
    state.theVector.push_back(gradMu_c);
    state.theTensor.push_back(F_c);

    return state;
}




double AnandIJSS2011MP::kineticPotential() const
{
    return diffusionPotential() + thermalPotential();
}




istensor AnandIJSS2011MP::materialConductivity() const
{
    return istensor::scaledIdentity(theAnandMaterial._thermalConductivity*temp_n);
}




// -d[J]/d[D]
istensor AnandIJSS2011MP::materialDiffusivity() const
{
    const double diff = theAnandMaterial._massDiffusivity;

    istensor K;
    if (theAnandMaterial._variational)
    {
        K = istensor::scaledIdentity(diff*temp_n/temp_c);
    }
    else
    {
        K = istensor::scaledIdentity(diff);
    }
    return K;
}




ivector AnandIJSS2011MP::materialHeatFlux() const
{
    ivector G = gradT_c/(-temp_c);
    return theAnandMaterial._thermalConductivity*temp_n * G;
}




ivector AnandIJSS2011MP::materialMassFlux() const
{
    const double diff = theAnandMaterial._massDiffusivity;

    ivector J;
    if (theAnandMaterial._variational)
    {
        J = -diff * temp_n /temp_c * gradMu_c;
    }
    else
    {
        J = -diff*gradMu_c;
    }
    return J;
}




ivector AnandIJSS2011MP::materialMassFluxDTheta() const
{
    const double diff = theAnandMaterial._massDiffusivity;

    ivector Jprime;
    if (theAnandMaterial._variational)
    {
        Jprime = 2.0 * diff * temp_n * temp_n /(temp_c * temp_c * temp_c) * gradMu_c;
    }
    else
    {
        Jprime.setZero();
    }
    return Jprime;
}




// coupling tensor M = d_P/d_mu
itensor AnandIJSS2011MP::materialStressChemicalTensor() const
{
    itensor Fcinv = F_c.inverse();

    const double beta  = theAnandMaterial._massExpansion;
    const double kappa = theAnandMaterial._bulkModulus;
    const double R     = theAnandMaterial._R;

    return -c_c * 3.0 * beta * kappa /(R*temp_c) * Fcinv.transpose();
}




// coupling tensor M = d^2(GCP)/(d_F d_Theta) = d_P/d_Theta
itensor AnandIJSS2011MP::materialStressTemperatureTensor() const
{
    itensor Fcinv = F_c.inverse();

    const double alpha = theAnandMaterial._thermalExpansion;
    const double kappa = theAnandMaterial._bulkModulus;
    const double beta  = theAnandMaterial._massExpansion;
    const double NM    = theAnandMaterial._NM;

    return (log(c_c/NM)*3.0*beta*kappa*c_c/temp_c-3.0*alpha*kappa)*Fcinv.transpose();
}




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




void AnandIJSS2011MP::resetCurrentState()
{
    fThermoMechMassMP::resetCurrentState();
    theFSMP->resetCurrentState();
}




// 4 d^2 GCP / dC^2
void AnandIJSS2011MP::secondPiolaKirchhoffStress(istensor& S) const
{
    theFSMP->secondPiolaKirchhoffStress(S);

    itensor  Finv = F_c.inverse();
    istensor Cinv = istensor::tensorTimesTensorTransposed(Finv);

    const double alpha = theAnandMaterial._thermalExpansion;
    const double kappa = theAnandMaterial._bulkModulus;
    const double tref  = theAnandMaterial.referenceTemperature();
    const double dTemp = temp_c - tref;

    const double beta  = theAnandMaterial._massExpansion;
    const double cref  = theAnandMaterial._referenceConcentration;
    const double dCon  = c_c - cref;

    S -= 3.0 * kappa * (alpha * dTemp + beta * dCon) * Cinv;
}



// - d chi / d theta = d^ GCP / (dtheta dmu)
double AnandIJSS2011MP::temperatureChemicalCoupling() const
{
    const double NM = theAnandMaterial._NM;
    return c_c/temp_c * log(c_c/NM);
}




double AnandIJSS2011MP::thermalPotential() const
{
    ivector G = gradT_c/(-temp_c);
    const double k  = theAnandMaterial._thermalConductivity*temp_n;
    return -0.5 * k * G.squaredNorm();
}




// d^2 GCP / d(theta)^2
double AnandIJSS2011MP::temperatureTangent() const
{
    const double R  = theAnandMaterial._R;
    const double NM = theAnandMaterial._NM;
    
    return -heatCapacityPerUnitVolume()/temp_c
        - temperatureChemicalCoupling() * R * log(c_c/NM);
}




void AnandIJSS2011MP::updateCurrentState(const double theTime, const itensor& F,
                                         const double temp, const ivector& gradT,
                                         const double mu, const ivector& gradMu)
{
    time_c   = theTime;
    F_c      = F;
    J_c      = F.determinant();
    temp_c   = temp;
    gradT_c  = gradT;
    mu_c     = mu;
    gradMu_c = gradMu;
    theFSMP->updateCurrentState(theTime, F);

    c_c = theAnandMaterial.fconcentration(J_c, temp_c, mu_c);
    if (isnan(c_c))
        c_c = theAnandMaterial.fconcentration(J_c, temp_c, mu_c);
}




double AnandIJSS2011MP::volumeFraction() const
{
    return c_c/theAnandMaterial._NM;
}




bool fThermoMechMassMP::testImplementation(std::ostream& of, const bool testDE, const bool testDDE) const
{
    bool isok = true;
    fThermoMechMassMP& theMP = const_cast<fThermoMechMassMP&>(*this);

    // set a random state in the material
    double t = randomUniform(0.0, 1.0);
    itensor F; F.setRandom(); F *= 2e-1;
    F += itensor::identity();
    if (F.determinant() < 0.0) F *= -1.0;

    double temp = theMP.parentMaterial().referenceTemperature()*randomUniform(0.5, 0.8);
    ivector gradT; gradT.setRandom(); gradT *= 1000.0;

    double  mu = randomUniform(-5000.0, -2000.0);
    ivector gradMu; gradMu.setRandom(); gradMu *= 1000.0;

    theMP.updateCurrentState(t, F, temp, gradT, mu, gradMu);
    theMP.commitCurrentState();

    // current state
    t = muesli::randomUniform(0.1,1.0);
    F.setRandom(); F *= 2e-1;
    F += itensor::identity();
    if (F.determinant() < 0.0) F *= -1.0;

    temp = theMP.parentMaterial().referenceTemperature()*randomUniform(2.0, 3.0);
    gradT.setRandom();  gradT *= 100.0;

    mu = randomUniform(-5000.0, -4000.0);
    gradMu.setRandom(); gradMu *= 100.0;

    theMP.updateCurrentState(t, F, temp, gradT, mu, gradMu);

    if (testDE)
    {
        {
            itensor pr_P; firstPiolaKirchhoffStress(pr_P);

            class gcpDF : public muesli::NumDiff
            {
            public:
                gcpDF(unsigned ii, unsigned jj, fThermoMechMassMP& xp, double xt, itensor& xF, double xtemp, ivector& xgradT, double xmu, ivector xgradMu)
                : i(ii), j(jj), point(xp), _t(xt), _temp(xtemp), _mu(xmu), _gradT(xgradT), _gradMu(xgradMu), _F(xF) {}

                virtual double eval(){return point.grandCanonicalPotential();}
                virtual void update(double dx)
                {
                    itensor F = _F;
                    F(i,j) += dx;
                    point.updateCurrentState(_t, F, _temp, _gradT, _mu, _gradMu);
                }

            private:
                unsigned i, j;
                fThermoMechMassMP& point;
                double _t, _temp, _mu;
                ivector _gradT, _gradMu;
                itensor _F;
            };

            itensor numP;
            for (unsigned i=0; i<3; i++)
            {
                for (unsigned j=0; j<3; j++)
                {
                    gcpDF dd(i, j, theMP, t, F, temp, gradT, mu, gradMu);
                    numP(i,j) = dd();
                }
            }

            itensor errorP = numP - pr_P;
            isok = (errorP.norm()/pr_P.norm() < 1e-4);
            of << "\n   1. Comparing P with derivative [d GCP / d F].";
            if (isok)
            {
                of << " Test passed.";
            }
            else
            {
                of << "\n Test failed. Relative error in 1st PK computation: " << errorP.norm()/pr_P.norm();
                of << "\n " << pr_P;
                of << "\n " << numP;
            }
        }

        {
            class gcpDTheta : public muesli::NumDiff
            {
            public:
                gcpDTheta(fThermoMechMassMP& xp, double xt, itensor& xF, double xtemp, ivector& xgradT, double xmu, ivector xgradMu)
                : point(xp), _t(xt), _temp(xtemp), _mu(xmu), _gradT(xgradT), _gradMu(xgradMu), _F(xF) {}

                virtual double eval(){return point.grandCanonicalPotential();}
                virtual void update(double dx) {point.updateCurrentState(_t, _F, _temp+dx, _gradT, _mu, _gradMu);}

            private:
                fThermoMechMassMP& point;
                double _t, _temp, _mu;
                ivector _gradT, _gradMu;
                itensor _F;
            };

            gcpDTheta dd(theMP, t, F, temp, gradT, mu, gradMu);
            double numEntropy = -dd();
            double error = numEntropy - entropy();

            isok = fabs(error)/fabs(entropy()) < 1e-4;
            of << "\n   2. Comparing entropy with derivative [- d GCP / d theta].";
            if (isok)
            {
                of << " Test passed.";
            }
            else
            {
                of << "\n Test failed. Relative error in entropy computation: " << error;
                of << "\n " << entropy();
                of << "\n " << numEntropy;
            }
        }

        {
            class gcpDmu : public muesli::NumDiff
            {
            public:
                gcpDmu(fThermoMechMassMP& xp, double xt, itensor& xF, double xtemp, ivector& xgradT, double xmu, ivector xgradMu)
                :point(xp), _t(xt), _temp(xtemp), _mu(xmu), _gradT(xgradT), _gradMu(xgradMu), _F(xF){}
                virtual double eval(){return point.grandCanonicalPotential();}
                virtual void update(double dx){point.updateCurrentState(_t, _F, _temp, _gradT, _mu+dx, _gradMu);}

            private:
                fThermoMechMassMP& point;
                const double _t, _temp, _mu;
                const ivector _gradT, _gradMu;
                const itensor _F;
            };

            gcpDmu dd(theMP, t, F, temp, gradT, mu, gradMu);
            double numConcentration = -dd();
            double error = numConcentration - c_c;

            isok = fabs(error)/c_c < 1e-4;
            of << "\n   3. Comparing concentration with derivative [- d GCP / d mu].";
            if (isok)
            {
                of << " Test passed.";
            }
            else
            {
                of << "\n Test failed. Relative error in concentration computation: " << error;
                of << "\n " << c_c;
                of << "\n " << numConcentration;
            }
        }

        {
            class cDmu : public muesli::NumDiff
            {
            public:
                cDmu(fThermoMechMassMP& xp, double xt, itensor& xF, double xtemp, ivector& xgradT, double xmu, ivector xgradMu)
                : point(xp), _t(xt), _temp(xtemp), _mu(xmu), _gradT(xgradT), _gradMu(xgradMu), _F(xF){}
                virtual double eval(){return point.concentration();}
                virtual void update(double dx){point.updateCurrentState(_t, _F, _temp, _gradT, _mu+dx, _gradMu);}

            private:
                fThermoMechMassMP& point;
                const double _t, _temp, _mu;
                const ivector _gradT, _gradMu;
                const itensor _F;
            };

            cDmu dd(theMP, t, F, temp, gradT, mu, gradMu);
            double numChemicalTangent = -dd();
            double error = numChemicalTangent - chemicalTangent();

            isok = fabs(error)/chemicalTangent() < 1e-4;
            of << "\n   4. Comparing chemical tangent with [- d c / d mu].";
            if (isok)
            {
                of << " Test passed.";
            }
            else
            {
                of << "\n Test failed. Relative error in chemical tangent: " << error;
                of << "\n " << chemicalTangent();
                of << "\n " << numChemicalTangent;
            }
        }

        {
            class conDgradT : public muesli::NumDiff
            {
            public:
                conDgradT(unsigned xi, fThermoMechMassMP& xp, double xt, itensor& xF, double xtemp, ivector& xgradT, double xmu, ivector xgradMu)
                : i(xi), point(xp), _t(xt), _temp(xtemp), _mu(xmu), _gradT(xgradT), _gradMu(xgradMu), _F(xF){}
                virtual double eval(){return point.thermalPotential();}
                virtual void update(double dx){ivector g=_gradT; g(i)+=dx; point.updateCurrentState(_t, _F, _temp, g, _mu, _gradMu);}

            private:
                const unsigned i;
                fThermoMechMassMP& point;
                const double _t, _temp, _mu;
                const ivector _gradT, _gradMu;
                const itensor _F;
            };

            ivector numHeat;
            for (unsigned i=0; i<3; i++)
            {
                conDgradT dd(i, theMP, t, F, temp, gradT, mu, gradMu);
                numHeat(i) = -dd();
            }

            numHeat *= -temp_c;
            ivector heat = materialHeatFlux();
            double error = (numHeat - heat).norm();

            isok = error/heat.norm() < 1e-4;
            of << "\n   5. Comparing heat flux with temp*[d Potential / d gradT].";
            if (isok)
            {
                of << " Test passed.";
            }
            else
            {
                of << "\n Test failed. Relative error in heat flux: " << error;
                of << "\n " << heat;
                of << "\n " << numHeat;
            }
        }

        {
            class heatDgradT : public muesli::NumDiff
            {
            public:
                heatDgradT(unsigned xi, unsigned xj, fThermoMechMassMP& xp, double xt, itensor& xF, double xtemp, ivector& xgradT, double xmu, ivector xgradMu)
                : i(xi), j(xj), point(xp), _t(xt), _temp(xtemp), _mu(xmu), _gradT(xgradT), _gradMu(xgradMu), _F(xF){}
                virtual double eval(){return point.materialHeatFlux()(i);}
                virtual void update(double dx){ivector g=_gradT; g(j)+=dx; point.updateCurrentState(_t, _F, _temp, g, _mu, _gradMu);}

            private:
                const unsigned i, j;
                fThermoMechMassMP& point;
                const double _t, _temp, _mu;
                const ivector _gradT, _gradMu;
                const itensor _F;
            };

            istensor numK;
            for (unsigned i=0; i<3; i++)
            {
                for (unsigned j=0; j<3; j++)
                {
                    heatDgradT dd(i, j, theMP, t, F, temp, gradT, mu, gradMu);
                    numK(i,j) = dd();
                }
            }
            numK *= -temp_c;

            istensor K = materialConductivity();
            double error = (numK - K).norm();

            isok = error/K.norm() < 1e-4;
            of << "\n   6. Comparing conductivity with (-temp) [d heat / d gradT].";
            if (isok)
            {
                of << " Test passed.";
            }
            else
            {
                of << "\n Test failed. Relative error in heat flux: " << error;
                of << "\n " << K;
                of << "\n " << numK;
            }
        }

        {
            class conDgradT : public muesli::NumDiff
            {
            public:
                conDgradT(unsigned xi, fThermoMechMassMP& xp, double xt, itensor& xF, double xtemp, ivector& xgradT, double xmu, ivector xgradMu)
                : i(xi), point(xp), _t(xt), _temp(xtemp), _mu(xmu), _gradT(xgradT), _gradMu(xgradMu), _F(xF){}
                virtual double eval(){return point.diffusionPotential();}
                virtual void update(double dx){ivector g=_gradMu; g(i)-=dx; point.updateCurrentState(_t, _F, _temp, _gradT, _mu, g);}

            private:
                const unsigned i;
                fThermoMechMassMP& point;
                const double _t, _temp, _mu;
                const ivector _gradT, _gradMu;
                const itensor _F;
            };

            ivector numFlux;
            for (unsigned i=0; i<3; i++)
            {
                conDgradT dd(i, theMP, t, F, temp, gradT, mu, gradMu);
                numFlux(i) = -dd();
            }

            ivector flux = materialMassFlux();
            double error = (numFlux - flux).norm();

            isok = error/flux.norm() < 1e-4;
            of << "\n   7. Comparing mass flux with  -[d Potential / d (-gradMu)].";
            if (isok)
            {
                of << " Test passed.";
            }
            else
            {
                of << "\n Test failed. Relative error in mass flux: " << error;
                of << "\n " << flux;
                of << "\n " << numFlux;
            }
        }

        {
            class heatDgradT : public muesli::NumDiff
            {
            public:
                heatDgradT(unsigned xi, unsigned xj, fThermoMechMassMP& xp, double xt, itensor& xF, double xtemp, ivector& xgradT, double xmu, ivector xgradMu)
                : i(xi), j(xj), point(xp), _t(xt), _temp(xtemp), _mu(xmu), _gradT(xgradT), _gradMu(xgradMu), _F(xF){}
                virtual double eval(){return point.materialMassFlux()(i);}
                virtual void update(double dx){ivector g=_gradMu; g(j)+=dx; point.updateCurrentState(_t, _F, _temp, _gradT, _mu, g);}
                virtual void reset(){point.updateCurrentState(_t, _F, _temp, _gradT, _mu, _gradMu);}

            private:
                const unsigned i, j;
                fThermoMechMassMP& point;
                const double _t, _temp, _mu;
                const ivector _gradT, _gradMu;
                const itensor _F;
            };

            istensor numK;
            for (unsigned i=0; i<3; i++)
            {
                for (unsigned j=0; j<3; j++)
                {
                    heatDgradT dd(i, j, theMP, t, F, temp, gradT, mu, gradMu);
                    numK(i,j) = -dd();
                }
            }

            istensor K = materialDiffusivity();
            double error = (numK - K).norm();

            isok = error/K.norm() < 1e-4;
            of << "\n   8. Comparing diffusivity with  -[d flux / d gradMu].";
            if (isok)
            {
                of << " Test passed.";
            }
            else
            {
                of << "\n Test failed. Relative error in diffusivity: " << error;
                of << "\n " << K;
                of << "\n " << numK;
            }
        }

        {
            class stressDmu : public muesli::NumDiff
            {
            public:
                stressDmu(unsigned xi, unsigned xj, fThermoMechMassMP& xp, double xt, itensor& xF, double xtemp, ivector& xgradT, double xmu, ivector xgradMu)
                : i(xi), j(xj), point(xp), _t(xt), _temp(xtemp), _mu(xmu), _gradT(xgradT), _gradMu(xgradMu), _F(xF){}
                virtual double eval(){itensor P; point.firstPiolaKirchhoffStress(P); return P(i,j);}
                virtual void update(double dx){point.updateCurrentState(_t, _F, _temp, _gradT, _mu+dx, _gradMu);}

            private:
                const unsigned i, j;
                fThermoMechMassMP& point;
                const double _t, _temp, _mu;
                const ivector _gradT, _gradMu;
                const itensor _F;
            };

            itensor numK;
            for (unsigned i=0; i<3; i++)
            {
                for (unsigned j=0; j<3; j++)
                {
                    stressDmu dd(i, j, theMP, t, F, temp, gradT, mu, gradMu);
                    numK(i,j) = dd();
                }
            }

            itensor K = materialStressChemicalTensor();
            double error = (numK - K).norm();

            isok = error/K.norm() < 1e-4;
            of << "\n   9. Comparing stressChemicalTensor with  [d P / d mu].";
            if (isok)
            {
                of << " Test passed.";
            }
            else
            {
                of << "\n Test failed. Relative error in stressChemical tensor: " << error;
                of << "\n " << K;
                of << "\n " << numK;
            }
        }

        {
            class stressDtheta : public muesli::NumDiff
            {
            public:
                stressDtheta(unsigned xi, unsigned xj, fThermoMechMassMP& xp, double xt, itensor& xF, double xtemp, ivector& xgradT, double xmu, ivector xgradMu)
                : i(xi), j(xj), point(xp), _t(xt), _temp(xtemp), _mu(xmu), _gradT(xgradT), _gradMu(xgradMu), _F(xF){}
                virtual double eval(){itensor P; point.firstPiolaKirchhoffStress(P); return P(i,j);}
                virtual void update(double dx){point.updateCurrentState(_t, _F, _temp+dx, _gradT, _mu, _gradMu);}

            private:
                const unsigned i, j;
                fThermoMechMassMP& point;
                const double _t, _temp, _mu;
                const ivector _gradT, _gradMu;
                const itensor _F;
            };

            itensor numK;
            for (unsigned i=0; i<3; i++)
            {
                for (unsigned j=0; j<3; j++)
                {
                    stressDtheta dd(i, j, theMP, t, F, temp, gradT, mu, gradMu);
                    numK(i,j) = dd();
                }
            }

            theMP.updateCurrentState(t, F, temp, gradT, mu, gradMu);
            itensor K = materialStressTemperatureTensor();
            double error = (numK - K).norm();

            isok = error/K.norm() < 1e-4;
            of << "\n   10. Comparing stressTemperatureTensor with  [d P / d theta].";
            if (isok)
            {
                of << " Test passed.";
            }
            else
            {
                of << "\n Test failed. Relative error in stress-temp tensor: " << error;
                of << "\n " << K;
                of << "\n " << numK;
            }
        }

        {
            class minuscDtheta : public muesli::NumDiff
            {
            public:
                minuscDtheta(fThermoMechMassMP& xp, double xt, itensor& xF, double xtemp, ivector& xgradT, double xmu, ivector xgradMu)
                : point(xp), _t(xt), _temp(xtemp), _mu(xmu), _gradT(xgradT), _gradMu(xgradMu), _F(xF){}
                virtual double eval(){return -point.concentration();}
                virtual void update(double dx){point.updateCurrentState(_t, _F, _temp+dx, _gradT, _mu, _gradMu);}

            private:
                fThermoMechMassMP& point;
                const double _t, _temp, _mu;
                const ivector _gradT, _gradMu;
                const itensor _F;
            };

            minuscDtheta dd(theMP, t, F, temp, gradT, mu, gradMu);
            double numK = dd();

            theMP.updateCurrentState(t, F, temp, gradT, mu, gradMu);
            double K = temperatureChemicalCoupling();
            double error = fabs(K-numK);

            isok = error/fabs(K) < 1e-4;
            of << "\n   11. Comparing temperatureChemicalTensor with  -[d c / d theta].";
            if (isok)
            {
                of << " Test passed.";
            }
            else
            {
                of << "\n Test failed. Relative error in temp-chemical tensor: " << error;
                of << "\n " << K;
                of << "\n " << numK;
            }
        }

        {
            class minuscDF : public muesli::NumDiff
            {
            public:
                minuscDF(unsigned ii, unsigned jj, fThermoMechMassMP& xp, double xt, itensor& xF, double xtemp, ivector& xgradT, double xmu, ivector xgradMu)
                : i(ii), j(jj), point(xp), _t(xt), _temp(xtemp), _mu(xmu), _gradT(xgradT), _gradMu(xgradMu), _F(xF){}
                virtual double eval(){return -point.concentration();}
                virtual void update(double dx){itensor F(_F); F(i,j)+=dx; point.updateCurrentState(_t, F, _temp, _gradT, _mu, _gradMu);}

            private:
                unsigned i, j;
                fThermoMechMassMP& point;
                const double _t, _temp, _mu;
                const ivector _gradT, _gradMu;
                const itensor _F;
            };

            itensor numK;
            for (unsigned i=0; i<3; i++)
            {
                for (unsigned j=0; j<3; j++)
                {
                    minuscDF dd(i, j, theMP, t, F, temp, gradT, mu, gradMu);
                    numK(i,j) = dd();
                }
            }

            theMP.updateCurrentState(t, F, temp, gradT, mu, gradMu);
            itensor K = materialStressChemicalTensor();
            double error = (K-numK).norm();

            isok = error/K.norm() < 1e-4;
            of << "\n   12. Comparing stressChemicalTensor with  -[d c / d F].";
            if (isok)
            {
                of << " Test passed.";
            }
            else
            {
                of << "\n Test failed. Relative error in stress-chemical tensor: " << error;
                of << "\n " << K;
                of << "\n " << numK;
            }
        }
    }

    if (testDDE)
    {
        itensor4 A, numA;
        materialTangent(A);

        class stressDF : public muesli::NumDiff
        {
        public:
            stressDF(unsigned xi, unsigned xj, unsigned xk, unsigned xl,
                     fThermoMechMassMP& xp, double xt, itensor& xF, double xtemp, ivector& xgradT, double xmu, ivector xgradMu)
            : i(xi), j(xj), k(xk), l(xl), point(xp), _t(xt), _temp(xtemp), _mu(xmu), _gradT(xgradT), _gradMu(xgradMu), _F(xF){}
            virtual double eval(){itensor P; point.firstPiolaKirchhoffStress(P); return P(i,j);}
            virtual void update(double dx){itensor F(_F); F(k,l)+= dx; point.updateCurrentState(_t, F, _temp, _gradT, _mu, _gradMu);}

        private:
            const unsigned i, j, k, l;
            fThermoMechMassMP& point;
            const double _t, _temp, _mu;
            const ivector _gradT, _gradMu;
            const itensor _F;
        };

        itensor numK;
        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++)
                    {
                        stressDF dd(i, j, k, l, theMP, t, F, temp, gradT, mu, gradMu);
                        numA(i,j,k,l) = dd();
                    }
                }
            }
        }

        double error = (A - numA).norm();

        isok = error/A.norm() < 1e-4;
        of << "\n   13. Comparing materialTangent with  [d P / d F].";
        if (isok)
        {
            of << " Test passed.";
        }
        else
        {
            of << "\n Test failed. Relative error in materialTangent tensor: " << error;
            of << "\n " << A;
            of << "\n " << numA;
        }
    }


    /*


     // 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-3);

     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 << "NUMÉRICA numA:  " << num_A;
     of << "ANALÍTICA prA:  " << pr_A;
     of << "DIFERENCIA ENTRE AMBAS:  " << errorA;
     }
     }

     // 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-3);

     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 << "     " ;
     of << "NUMERICA nTg" << nTg;
     of << "ANALITICA tg" << tg;
     itensor4 errorTg = nTg - tg;
     of << "DIFERENCIA ENTRE AMBAS" << errorTg;
     }
     of << std::flush;
     }



     // test coupling tensor with d P / d mu
     {
     itensor pr_coupling = materialCouplingTensor();
     itensor error = num_coupling - pr_coupling;
     isok = (error.norm()/pr_coupling.norm() < 2e-3);
     of << "\n   7. Comparing M with [d P / d mu].";
     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 chemical tangent
     {
     double pr_kmu = chemicalTangent();

     double error = num_kmu - pr_kmu;
     isok = (fabs(error)/fabs(pr_kmu) < 1e-3);
     of << "\n   8. Comparing chemical tangent  -[d c / d mu].";
     if (isok)
     {
     of << " Test passed.";
     }
     else
     {
     of << "\n Relative error kmu computation %e. Test failed." << fabs(error)/fabs(pr_kmu);
     of << "\n " << pr_kmu;
     of << "\n " << num_kmu;
     }
     }


     // Derivatives with respect to G = Grad[mu]
     inc = 1e-5*theThermoMechMassMaterial._mu0;

     // approximation of derivative
     const size_t nnumder = 4.0;
     const double ndtimes[] = {+1.0, +2.0, -1.0, -2.0};
     const double ndfact[]  = {+8.0, -1.0, -8.0, +1.0};
     const double ndden = 12.0;


     // Compute numerical value of flux j = d Omega / d G
     {
     ivector num_j;
     for (unsigned i=0; i<3; i++)
     {
     double der = 0.0;
     for (unsigned q=0; q<nnumder; q++)
     {
     ivector gmu = gradMu_n;
     gmu[i] = gradMu_n[i] + inc*ndtimes[q];
     theMP.updateCurrentState(tn1, F_n, mu_n, gmu);
     der += ndfact[q]*diffusionPotential();
     theMP.resetCurrentState();
     }
     der /= inc*ndden;
     num_j(i) = der;
     }

     // compare mass flus
     {
     ivector pr_j = materialMassFlux();

     double error = (num_j - pr_j).norm();
     isok = (fabs(error)/pr_j.norm() < 1e-3);
     of << "\n   9. Comparing mass flux.";
     if (isok)
     {
     of << " Test passed.";
     }
     else
     {
     of << "\n Test failed.";
     of << "\n Relative error mass flux computation:" << fabs(error)/pr_j.norm();
     of << "\n " << pr_j;
     of << "\n " << num_j;
     }
     }
     }


     // Compute numerical value of hessian h = - d j / d G
     {
     istensor num_h;
     for (unsigned i=0; i<3; i++)
     {
     for (unsigned j=0; j<3; j++)
     {
     double der = 0.0;
     for (unsigned q=0; q<nnumder; q++)
     {
     ivector gmu = gradMu_n;
     gmu[j] = gradMu_n[j] + inc*ndtimes[q];
     theMP.updateCurrentState(tn1, F_n, mu_n, gmu);
     der += ndfact[q]*materialMassFlux()(i);
     theMP.resetCurrentState();
     }
     der /= inc*ndden;
     num_h(i,j) = -der;
     }
     }

     // compare mass hessian
     {
     istensor pr_h = diffusionTangent();

     double error = (num_h - pr_h).norm();
     isok = (fabs(error)/pr_h.norm() < 1e-3);
     of << "\n   10. Comparing mass diffusion hessian.";
     if (isok)
     {
     of << " Test passed.";
     }
     else
     {
     of << "\n Test failed.";
     of << "\n Relative error:" << fabs(error)/pr_h.norm();
     of << "\n " << pr_h;
     of << "\n " << num_h;
     }
     }
     }
     */
    return true;
}






