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

using namespace muesli;

fMechMassMaterial::fMechMassMaterial(const std::string& name,
                                     const materialProperties& cl)
:
material(name, cl),
_variational(false),
_Omega(0.0),
_G(0.0),
_K(0.0),
_R(0.0),
_mu0(0.0),
_theta0(298.0),
_chi(0.0),
_diff(0.0)
{
    muesli::assignValue(cl, "omega", _Omega);
    muesli::assignValue(cl, "g", _G);
    muesli::assignValue(cl, "k", _K);
    muesli::assignValue(cl, "r", _R);
    muesli::assignValue(cl, "mu0", _mu0);
    muesli::assignValue(cl, "tref", _theta0);
    muesli::assignValue(cl, "chi", _chi);
    muesli::assignValue(cl, "diffusivity", _diff);
    _variational = muesli::hasKeyword(cl, "variational");
}




fMechMassMP* fMechMassMaterial::createMaterialPoint() const
{
    muesli::fMechMassMP* mp = new fMechMassMP(*this);
    return mp;
}




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




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




void fMechMassMaterial::print(std::ostream &of) const
{
    of  << "\n Mechanical - mass transport coupled material."
        << "\n Model based on the work by Chester & Anand."
        << "\n Variational formulation: " << _variational
        << "\n Model constants:"
        << "\n   Omega:  " << _Omega
        << "\n   G:      " << _G
        << "\n   K:      " << _K
        << "\n   R:      " << _R
        << "\n   mu0:    " << _mu0
        << "\n   theta0: " << _theta0
        << "\n   chi:    " << _chi
        << "\n   Diffus: " << _diff;
}




void fMechMassMaterial::setRandom()
{
     _Omega = muesli::randomUniform(1.0, 5.0)*1e-4;
     _G     = muesli::randomUniform(1.0, 10.0)*1e3;
     _K     = muesli::randomUniform(1.0, 10.0)*1e3;
     _R     = muesli::randomUniform(1.0, 10.0);
     _theta0= muesli::randomUniform(250.0, 350.0);
     _chi   = muesli::randomUniform(0.1, 0.2);
     _diff  = muesli::randomUniform(1.0, 10.0)*1e-4;
    
    double J   = 1.0;
    double c   = 1.0;
    double Js  = 1.0 + _Omega*c;
    double Je  = J/Js;
    double phi = 1.0/Js;
    double lje = log(Je);

    _mu0 = -(_R* _theta0 *(log(1.0-phi) + phi + _chi*phi*phi)
            - _Omega * _K * lje + 0.5 * _K * _Omega * lje * lje);
}




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




fMechMassMP::fMechMassMP(const fMechMassMaterial& m) :
theMechMassMaterial(m),
time_n(0.0),
mu_n(0.0),
c_n(0.0),
J_n(1.0),
time_c(0.0),
mu_c(0.0),
c_c(0.0),
J_c(1.0)
{
    gradMu_n.setZero();
    gradMu_c.setZero();
    F_n = itensor::identity();
    F_c = itensor::identity();
    c_c = c_n = fConcentration(F_c, mu_c, 10.0);
}




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




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




// (minus) d [ concentration] / d [ chemical potential ]
double fMechMassMP::chemicalTangent() const
{
    const double inc = std::max<double>(1e-5, 1e-5*theMechMassMaterial._mu0);
    double Xip1 = fConcentration(F_c, mu_c+inc, c_c);
    double Xip2 = fConcentration(F_c, mu_c+2.0*inc, c_c);
    double Xim1 = fConcentration(F_c, mu_c-inc, c_c);
    double Xim2 = fConcentration(F_c, mu_c-2.0*inc, c_c);

    // fourth order approximation of the (negative) derivative
    return -(-Xip2 + 8.0*Xip1 - 8.0*Xim1 + Xim2)/(12.0*inc);
}



/*
itensor fMechMassMP::der_cr_F() const
{
    fMechMassMP& theMP = const_cast<fMechMassMP&>(*this);
    double inc = 1.0e-3;
    itensor F = F_c;
    itensor DcF;
    
    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(time_c, F, mu_c, gradMu_c);
            double Xp1 = fConcentration(F, mu_c, c_c);
            
            F(i,j) = original + 2.0*inc;
            theMP.updateCurrentState(time_c, F, mu_c, gradMu_c);
            double Xp2 = fConcentration(F, mu_c, c_c);
            
            F(i,j) = original - inc;
            theMP.updateCurrentState(time_c, F, mu_c, gradMu_c);
            double Xm1 = fConcentration(F, mu_c, c_c);
            
            F(i,j) = original -2.0*inc;
            theMP.updateCurrentState(time_c, F, mu_c, gradMu_c);
            double Xm2 = fConcentration(F, mu_c, c_c);
            
            DcF(i,j) = (-Xp2 + 8.0*Xp1 - 8.0*Xm1 + Xm2)/(12.0*inc);
            
            F(i,j) = original;
            theMP.updateCurrentState(time_c, F, mu_c, gradMu_c);
        }
    }
    
    return DcF;
}
*/



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




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




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




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

}




void fMechMassMP::convectedTangent(itensor4& ctg) const
{
    istensor C  = istensor::tensorTransposedTimesTensor(F_c);
    istensor Ci = C.inverse();

    itensor4 mtang;
    materialTangent(mtang);

	itensor  Fi = F_c.inverse();
    istensor SS;
    secondPiolaKirchhoffStress(SS);

    ctg.setZero();
    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++)
                {
                    ctg(c,a,d,b) = -SS(a,b)*Ci(c,d);
                    
                    for (unsigned i=0; i<3; i++)
                        for (unsigned j=0; j<3; j++)
                            ctg(c,a,d,b) += Fi(c,i)*mtang(i,a,j,b)*Fi(d,j);
                }
}




void fMechMassMP::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 fMechMassMP::convergedConcentration() const
{
    return c_n;
}




double fMechMassMP::diffusionPotential() const
{
    const double diff   = theMechMassMaterial._diff;
    const double R      = theMechMassMaterial._R;
    const double theta0 = theMechMassMaterial._theta0;

    double cJ = theMechMassMaterial._variational ? c_n/J_n : c_c/J_c;
    return -0.5 * cJ*diff/(R*theta0) * gradMu_c.squaredNorm();
}




// -d[J]/d[D]
istensor fMechMassMP::diffusionTangent() const
{
    const double diff   = theMechMassMaterial._diff;
    const double R      = theMechMassMaterial._R;
    const double theta0 = theMechMassMaterial._theta0;

    double cJ = theMechMassMaterial._variational ? c_n/J_n : c_c/J_c;
    return istensor::scaledIdentity(diff*cJ/(R*theta0));
}




double fMechMassMP::dissipation() const
{
    return 0.0;
}




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




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

}



// this functio re-computes the chemical potential as function of F and c,
// not the current state of the point
double fMechMassMP::fChemicalPotential(const itensor& F, double c) const
{
    const double K = theMechMassMaterial._K;
    const double O = theMechMassMaterial._Omega;
    const double mu0 = theMechMassMaterial._mu0;
    const double R = theMechMassMaterial._R;
    const double chi = theMechMassMaterial._chi;
    const double tref = theMechMassMaterial._theta0;
    double J    = F.determinant();
    double Js   = 1.0 + O*c;
    double Je   = J/Js;
    double phi  = 1.0/Js;
    double lje  = log(Je);
    double mu = mu0 + R*tref*(log(1.0-phi) + phi + chi*phi*phi)
        - O*K*lje + 0.5*K*O*lje*lje;

    return mu;
}




double fMechMassMP::fConcentration(const itensor&F, const double mu, const double c0) const
{
    const double K = theMechMassMaterial._K;
    const double O = theMechMassMaterial._Omega;
    const double R = theMechMassMaterial._R;
    const double chi = theMechMassMaterial._chi;
    const double tref = theMechMassMaterial._theta0;
    const double mu0 = theMechMassMaterial._mu0;
    double c = c0;
    double error = mu - fChemicalPotential(F, c);
    unsigned count = 0;
    const double tol = std::max<double>(1e-9*mu0, 1e-9);
    double J = F.determinant();

    while (fabs(error) > tol && count++ < 50)
    {
        double Js  = 1.0 + O*c;
        double Je  = J/Js;
        double phi = 1.0/Js;
        double lje = log(Je);

        double tg = -K*O*O/Js*(1.0 - lje) + R*tref*O/(Js*Js)*( -phi/(1.0-phi) + 2.0*chi*phi);
        c -= error/tg;
        c = std::max<double>(1e-7, c);
        error =  mu - fChemicalPotential(F, c);
    }

    if (count == 50)
        std::cout << "\n fConcentration not converged.";
    
    return c;
}




itensor fMechMassMP::fFirstPiolaKirchhoffStress(const itensor& F, double mu) const
{
    const double G = theMechMassMaterial._G;
    const double K = theMechMassMaterial._K;
    const double O = theMechMassMaterial._Omega;

    double c  = fConcentration(F, mu, c_c);
    double J  = F.determinant();
    itensor iFt = F.inverse().transpose();

    double Js = 1.0 + O*c;
    double Je = J/Js;

    return G*(F - iFt) + K*Js*log(Je)*iFt;
}




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




double fMechMassMP::freeEnergy() const
{
    const double G = theMechMassMaterial._G;
    const double K = theMechMassMaterial._K;
    const double O = theMechMassMaterial._Omega;
    const double mu0 = theMechMassMaterial._mu0;
    const double R = theMechMassMaterial._R;
    const double chi = theMechMassMaterial._chi;
    const double tref = theMechMassMaterial._theta0;

    const istensor C = istensor::tensorTransposedTimesTensor(F_c);
    double lbar  = sqrt(C.trace()/3.0);
    double Js    = 1.0 + O*c_c;
    double Je    = J_c/Js;
    double logJe = log(Je);


    double psi = mu0*c_c + R*tref*c_c*( log(O*c_c/Js) + chi/Js)
              + 0.5*G*(3.0*(lbar*lbar - 1.0) - 2.0*log(J_c))
           + Js*0.5*K*logJe*logJe;
    return psi;
}




double fMechMassMP::freeEntropy() const
{
    double tref = theMechMassMaterial._theta0;
    return -freeEnergy()/tref;
}




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




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

    return state;
}




double fMechMassMP::grandCanonicalPotential() const
{
    return freeEnergy() - c_c * mu_c;
}




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




double fMechMassMP::internalEnergy() const
{
    return freeEnergy();
}




double fMechMassMP::kineticPotential() const
{
    return diffusionPotential();
}




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




// coupling tensor M = d P / d mu
itensor fMechMassMP::materialCouplingTensor() const
{
    const double inc = std::max<double>(1e-5*theMechMassMaterial._mu0, 1e-5);
    itensor Pp1 = fFirstPiolaKirchhoffStress(F_c, mu_c+inc);
    itensor Pp2 = fFirstPiolaKirchhoffStress(F_c, mu_c+2.0*inc);
    itensor Pm1 = fFirstPiolaKirchhoffStress(F_c, mu_c-inc);
    itensor Pm2 = fFirstPiolaKirchhoffStress(F_c, mu_c-2.0*inc);

    // fourth order approximation of the derivative
    itensor M = (-Pp2 + 8.0*Pp1 - 8.0*Pm1 + Pm2)/(12.0*inc);

    return M;
}




ivector fMechMassMP::materialMassFlux() const
{
    const double diff   = theMechMassMaterial._diff;
    const double R      = theMechMassMaterial._R;
    const double theta0 = theMechMassMaterial._theta0;

    double cJ = theMechMassMaterial._variational ? c_n/J_n : c_c/J_c;
    return -cJ*diff/(R*theta0) * gradMu_c;
}




void fMechMassMP::materialTangent(itensor4& cm) const
{
    const double G = theMechMassMaterial._G;
    const double K = theMechMassMaterial._K;
    const double O = theMechMassMaterial._Omega;
    
    double Js    = 1.0 + O*c_c;
    double Je    = J_c/Js;
    double logJe = log(Je);
    
    itensor Fi  = F_c.inverse();
    itensor Fit = Fi.transpose();
    
    itensor d_P_cr = K*O*(logJe-1.0)*Fit;
    itensor d_cr_F = -materialCouplingTensor(); //der_cr_F();

    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,k,l) += G*Fi(l,i)*Fi(j,k) + Js*K* Fi(j,i)*Fi(l,k) - Js*K*logJe*Fi(l,i)*Fi(j,k);
                    cm(i,j,k,l) += d_P_cr(i,j) * d_cr_F(k,l);
                }
    cm += G*itensor4::identity(); //check with Ángel
}




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




double fMechMassMP::referenceChemicalPotential() const
{
    return theMechMassMaterial._mu0;
}




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




void fMechMassMP::secondPiolaKirchhoffStress(istensor& S) const
{
    const double G = theMechMassMaterial._G;
    const double K = theMechMassMaterial._K;
    const double O = theMechMassMaterial._Omega;

    istensor Cinv = istensor::tensorTimesTensorTransposed(F_c.inverse());
    double Js = 1.0 + O*c_c;
    double Je = J_c/Js;

    S = G*(istensor::identity()-Cinv) + K*Js*log(Je)*Cinv;
}




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




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

    // set a random update in the material
    
    //*****RANDOM
    itensor F_n;
    F_n.setRandom();
    F_n *= 1e-2;
    F_n += itensor::identity();
    if (F_n.determinant() < 0.0) F_n *= -1.0;

    ivector gradMu_n;
    gradMu_n.setRandom();
    
    double  mu_n;
    mu_n = randomUniform(-5000.0, -5000.0);

    theMP.updateCurrentState(0.0, F_n, mu_n, gradMu_n);
    theMP.commitCurrentState();

    double tn1 = muesli::randomUniform(0.1,1.0);
    
    itensor F;
    F.setRandom();
    F *= 1e-1;
    F += itensor::identity();
    if (F.determinant() < 0.0) F *= -1.0;
    
    ivector gradMu;
    gradMu.setRandom();
    
    double mu;
    mu = randomUniform(-5000.0, -5000.0);
    theMP.updateCurrentState(tn1, F, mu, gradMu);


    itensor4 ctg;
    theMP.convectedTangent(ctg);
    itensor4 ctg2;

    // 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
    itensor num_P;
    itensor4 num_A;
    {
        itensor dP, Pp1, Pp2, Pm1, Pm2;

        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, mu, gradMu);
                double Wp1 = grandCanonicalPotential();
                firstPiolaKirchhoffStress(Pp1);

                F(i,j) = original + 2.0*inc;
                theMP.updateCurrentState(tn1, F, mu, gradMu);
                double Wp2 = grandCanonicalPotential();
                firstPiolaKirchhoffStress(Pp2);

                F(i,j) = original - inc;
                theMP.updateCurrentState(tn1, F, mu, gradMu);
                double Wm1 = grandCanonicalPotential();
                firstPiolaKirchhoffStress(Pm1);

                F(i,j) = original - 2.0*inc;
                theMP.updateCurrentState(tn1, F, mu, gradMu);
                double Wm2 = grandCanonicalPotential();
                firstPiolaKirchhoffStress(Pm2);

                // 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);
                for (unsigned k=0; k<3; k++)
                    for (unsigned l=0; l<3; l++)
                        num_A(k,l,i,j) = dP(k,l);

                F(i,j) = original;
                theMP.updateCurrentState(tn1, F, mu, gradMu);
            }
        }
    }



    // Derivatives with respect to mu
    // Compute numerical value of concentration c =  -(d GCP / d mu)
    // Compute numerical value of coupling tensor M = d P / d mu
    // Compute numerical value of chemical tangent =  -(d c / d mu)
    inc = 1e-5*theMechMassMaterial._mu0;
    double  num_c;
    itensor num_coupling;
    double  num_kmu;
    {
        theMP.updateCurrentState(tn1, F, mu+inc, gradMu);
        double  GCP_p1 = grandCanonicalPotential();
        itensor PK_p1; firstPiolaKirchhoffStress(PK_p1);
        double  c_p1 = concentration();

        theMP.updateCurrentState(tn1, F, mu+2.0*inc, gradMu);
        double  GCP_p2 = grandCanonicalPotential();
        itensor PK_p2; firstPiolaKirchhoffStress(PK_p2);
        double  c_p2 = concentration();

        theMP.updateCurrentState(tn1, F, mu-inc, gradMu);
        double  GCP_m1 = grandCanonicalPotential();
        itensor PK_m1; firstPiolaKirchhoffStress(PK_m1);
        double  c_m1 = concentration();

        theMP.updateCurrentState(tn1, F, mu-2.0*inc, gradMu);
        double  GCP_m2 = grandCanonicalPotential();
        itensor PK_m2; firstPiolaKirchhoffStress(PK_m2);
        double  c_m2 = concentration();

        theMP.updateCurrentState(tn1, F, mu, gradMu);

        // fourth order approximation of the derivative
        num_c        = -(-GCP_p2 + 8.0*GCP_p1 - 8.0*GCP_m1 + GCP_m2)/(12.0*inc);
        num_coupling =  (-PK_p2 + 8.0*PK_p1 - 8.0*PK_m1 + PK_m2)/(12.0*inc);
        num_kmu      = -(-c_p2 + 8.0*c_p1  - 8.0*c_m1  + c_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 Test failed. Relative error in 1st PK computation: " << errorP.norm()/pr_P.norm();
            of << "\n " << pr_P;
            of << "\n " << num_P;
        }
    }
    else
    {
        of << "\n   2. Comparing 1PK with derivative of Psi";
    }
    
    // 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 concentration with derivative of GCP wrt mu
    {
        double pr_c = concentration();

        double error = num_c - pr_c;
        isok = (fabs(error)/fabs(pr_c) < 1e-4);
        of << "\n   4. Comparing concentration potential with  -[d GCP / d mu ].";
        if (isok)
        {
            of << " Test passed.";
        }
        else
        {
            of << "\n Test failed. Relative error " << fabs(error)/fabs(pr_c);
            of << "\n Programmed c: " << pr_c;
            of << "\n Numeric c:    " << num_c;
        }
    }


    // 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*theMechMassMaterial._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;
}




double fMechMassMP::volumeFraction() const
{
    const double O = theMechMassMaterial._Omega;

    double Js  = 1.0 + O*c_c;
    double phi = 1.0/Js;
    return phi;
}




double fMechMassMP::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 fMechMassMaterial::waveVelocity() const
{
    return sqrt(2.0*_G/density());
}




void fMechMassMP::updateCurrentState(const double theTime, const itensor& F, const double& mu, const ivector& gradMu)
{
    time_c   = theTime;
    F_c      = F;
    J_c      = F.determinant();
    mu_c     = mu;
    gradMu_c = gradMu;
    c_c      = fConcentration(F_c, mu_c, c_c);
}
