/****************************************************************************
*
*                                 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 "finitestrain.h"
using namespace muesli;


finiteStrainMaterial::finiteStrainMaterial(const std::string& name)
: material(name),
bulk(0.0),
cp(0.0), cs(0.0),
damageModelActivated(false), deletion(false)
{}




finiteStrainMaterial::finiteStrainMaterial(const std::string& name,
                                           const materialProperties& cl)
:
material(name, cl),
bulk(0.0),
cp(0.0), cs(0.0),
damageModelActivated(false), deletion(false)
{
}




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




double finiteStrainMaterial::getProperty(const propertyName p) const
{
    double ret=0.0;
    std::ostream& mlog = material::getLogger();


    // scan all the possible data
    switch(p)
    {
        case PR_BULK:   ret = bulk;        break;
        case PR_CP:     ret = cp;        break;
        case PR_CS:     ret = cs;        break;

        default:
            mlog << "property not defined";
    }
    return ret;
}




double finiteStrainMaterial::waveVelocity() const
{
    return cp;
}




finiteStrainMP::finiteStrainMP(const finiteStrainMaterial& m) :
theFiniteStrainMaterial(m),
 tn(0.0), tc(0.0),
 Jn(1.0), Jc(1.0),
 Dn(0.0), Dc(0.0),
 tempn(m.referenceTemperature()), tempc(m.referenceTemperature()),
 fullyDamaged(false)
{
    Fn = Fc = itensor::identity();
}




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




void finiteStrainMP::commitCurrentState()
{
    tn    = tc;
    Fn    = Fc;
    Jn    = Jc;
    tempn = tempc;
    if (theFiniteStrainMaterial.damageModelActivated) Dn = Dc;
}




void finiteStrainMP::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,j) += c(i,j,k,l)*v1(j)*v2(l);
                }
}




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




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





/*
 * Given the elasticity material C_{ABCD} = 2 d S_{AB} / d C_{CD} = 4 d^2 W(C) / d C^2,
 * and two material vectors V, W, obtain the (non-symmetric) tensor
 *
 *   T_{AB} = C_{ABCD} V_B W_D
 *
 *
 *  Although in some cases there might be a faster way, we reuse the function contractWithSpatialTangent
 *  through the expression:
 *
 *     C[ V, W] = J F^{-1} c[F^{-T]V, F^{-T}W] F^{-T}
 *
 */
void finiteStrainMP::convectedTangentMatrix(double C[6][6]) const
{
    itensor4 ct;
    convectedTangent(ct);
    muesli::tensorToMatrix(ct, C);
}




void finiteStrainMP::convectedTangentTimesSymmetricTensor(const istensor &M, istensor &CM) const
{
    itensor4 C;
    convectedTangent(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++)
                {
                    CM(i,j) += C(i,j,k,l)*M(k,l);
                }
}




itensor& finiteStrainMP::deformationGradient()
{
    return Fc;
}




const itensor& finiteStrainMP::deformationGradient() const
{
    return Fc;
}



double finiteStrainMP::effectiveStoredEnergy() const
{
    double dt = tc - tn;
    return storedEnergy() + dt * dissipatedEnergy();
}




void finiteStrainMP::energyMomentumTensor(itensor &T) const
{
    itensor P;
    firstPiolaKirchhoffStress(P);

    T = storedEnergy()*itensor::identity() - Fc.transpose()*P;
}




void finiteStrainMP::firstPiolaKirchhoffStress(itensor &P) const
{
    istensor sigma;
    CauchyStress(sigma);
    P = Jc * sigma * Fc.inverse().transpose();
}




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

    state.theTime = tn;
    state.theDouble.push_back(Jn);
    state.theTensor.push_back(Fn);
    if(theFiniteStrainMaterial.damageModelActivated) state.theDouble.push_back(Dn);

    return state;
}




materialState finiteStrainMP::getCurrentState() const
{
    materialState state;

    state.theTime = tc;
    state.theDouble.push_back(Jc);
    state.theTensor.push_back(Fc);
    if(theFiniteStrainMaterial.damageModelActivated) state.theDouble.push_back(Dc);

    return state;
}




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




void finiteStrainMP::KirchhoffStressVector(double tauv[6]) const
{
    istensor tau;
    KirchhoffStress(tau);
    muesli::tensorToVector(tau, tauv);
}




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




void finiteStrainMP::materialTangentMatrix(double C[6][6]) const
{
    itensor4 mt;
    materialTangent(mt);
    muesli::tensorToMatrix(mt, C);
}




void finiteStrainMP::numericalConvectedTangent(itensor4& C) const
{
    itensor4 A;
    itensor dP, Pp1, Pp2, Pm1, Pm2;

    // numerical differentiation stress
    itensor numP;
    numP.setZero();
    itensor Fnp;
    Fnp.setZero();
    Fnp=Fc;
    const double inc = 1.0e-5*Fc.norm();;
    double tn1 = this->tc;
    
    finiteStrainMP& theMPc = const_cast<finiteStrainMP&>(*this);
    
    for (unsigned i=0; i<3; i++)
    {
        for (unsigned j=0; j<3; j++)
        {
            double original = this->Fc(i,j);
            
            Fnp(i,j) = original + inc;
            theMPc.updateCurrentState(tn1, Fnp);
            double Wp1 = theMPc.effectiveStoredEnergy();
            theMPc.firstPiolaKirchhoffStress(Pp1);

            Fnp(i,j) = original + 2.0*inc;
            theMPc.updateCurrentState(tn1, Fnp);
            double Wp2 = theMPc.effectiveStoredEnergy();
            theMPc.firstPiolaKirchhoffStress(Pp2);

            Fnp(i,j) = original - inc;
            theMPc.updateCurrentState(tn1, Fnp);
            double Wm1 = theMPc.effectiveStoredEnergy();
            theMPc.firstPiolaKirchhoffStress(Pm1);

            Fnp(i,j) = original - 2.0*inc;
            theMPc.updateCurrentState(tn1, Fnp);
            double Wm2 = theMPc.effectiveStoredEnergy();
            theMPc.firstPiolaKirchhoffStress(Pm2);

            // fourth order approximation of the derivative
            numP(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++)
                    A(k,l,i,j) = dP(k,l);

            Fnp(i,j) = original;
            theMPc.updateCurrentState(tn1, Fnp);
        }
    }
    
    itensor4 nTg;
    nTg.setZero();

    // transform A to get the convected tangent
    itensor  J  = Fnp.inverse();
    istensor Cn  = istensor::tensorTransposedTimesTensor(Fnp);
    istensor Ci = Cn.inverse();
    istensor S;
    theMPc.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)*A(i,a,j,b)*J(d,j);
                }
    C=nTg;
}




void finiteStrainMP::resetCurrentState()
{
    tc = tn;
    Fc = Fn;
    Jc = Jn;
    tempc = tempn;
    if(theFiniteStrainMaterial.damageModelActivated) Dc = Dn;
}




void finiteStrainMP::secondPiolaKirchhoffStress(istensor& S) const
{
    istensor sigma;
    CauchyStress(sigma);
    S =  Jc * istensor::FSFt(Fc.inverse(), sigma);
}




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




void  finiteStrainMP::setRandom()
{
    tc = tn + muesli::randomUniform(0.1, 1.0);
    Fc.setRandom();
    Jc = Fc.determinant();
    
    if (theFiniteStrainMaterial.damageModelActivated) Dc = muesli::randomUniform(0.0,1.0);

    if (Jc < 0.0)
    {
        Fc *= -1.0;
        Jc  = -Jc;
    }
}




void finiteStrainMP::spatialTangent(itensor4& Cs) const
{
    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 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++)
                                    Cs(i,j,k,l) += Fc(i,A)*Fc(j,B)*Fc(k,C)*Fc(l,D)*Cc(A,B,C,D);

    Cs *= 1.0/Jc;
}




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




void finiteStrainMP::updateCurrentState(const double theTime, const itensor& F)
{
    tc = theTime;
    Fc = F;
    Jc = F.determinant();
}




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




bool finiteStrainMP::testImplementation(std::ostream& of, const bool testDE, const bool testDDE) const
{
    bool isok = true;

    // set a random update in the material
    itensor F;
    F.setRandom();

    if (F.determinant() < 0.0) F *= -1.0;
    finiteStrainMP& theMP = const_cast<finiteStrainMP&>(*this);

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

    double tn1 = muesli::randomUniform(0.1,1.0);
    itensor u;
    u.setRandom();
    u *= 1e-2;
    F += u;

    if (F.determinant() < 0.0) F *= -1.0;
    theMP.updateCurrentState(tn1, F);

    // programmed tangent of elasticities
    itensor4 tg;
    convectedTangent(tg);

    itensor P;
    firstPiolaKirchhoffStress(P);

    // material tangent A_{iAjB} = d (P_iA) / d F_jB
    itensor4 A;
    itensor dP, Pp1, Pp2, Pm1, Pm2;

    // numerical differentiation stress
    itensor numP;
    numP.setZero();
    const double inc = 1.0e-5;

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

            F(i,j) = original + inc;
            theMP.updateCurrentState(tn1, F);
            double Wp1 = effectiveStoredEnergy();
            firstPiolaKirchhoffStress(Pp1);

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

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

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

            // fourth order approximation of the derivative
            numP(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++)
                    A(k,l,i,j) = dP(k,l);

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



    // compare DEnergy with the derivative of Energy
    {
        if (testDE)
        {
            itensor errorP = numP - P;
            isok = (errorP.norm()/P.norm() < 1e-4);
            of << "\n   1. Comparing P with derivative of Weff.";
            if (isok)
            {
                of << " Test passed.";
            }
            else
            {
                of << "\n Relative error in DE computation: " << errorP.norm()/P.norm() << ". Test failed.";
                of << "\n P: \n" << P;
                of << "\n numP: \n" << numP;
            }
        }
        else
        {
            of << "\n   1. Comparing P with derivative of Weff ::: not run for this material";
        }
        of << std::flush;
    }


    // test the consistency of the stress tensors sigma and P
    {
        istensor sigma, S;
        CauchyStress(sigma);
        secondPiolaKirchhoffStress(S);

        itensor P1 = F.determinant() * sigma * F.inverse().transpose();
        itensor P2 = F*S;

        itensor errorP1 = P1 - P;
        itensor errorP2 = P2 - P;
        double  error = errorP1.norm() + errorP2.norm();
        isok = error/P.norm() < 1e-4;
        of << "\n   2. Checking consistency of the 1PK, 2PK, and Cauchy stresses.";
        if (isok)
        {
            of << " Test passed.";
        }
        else
        {
            of << "\n Test failed. Relative error: " << error/P.norm();
            of << "\n P: " << P;
            of << "\n P1: " << P1;
            of << "\n P2: " << P2;
            of << "\n numP: " << numP;
            of << "\n Cauchy: " << sigma;
        }
        of << std::flush;
    }


    // test tangent as derivative of the stress
    if (testDDE)
    {
        // numeric material tangent
        itensor4 nTg;
        nTg.setZero();

        // transform 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)*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   3. 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;
    }


    // tangent contracted functions
    if ((true))
    {
        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();
        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   4. 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   5. 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);
        finiteStrainMP::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   6. 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);
        finiteStrainMP::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   7. 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 = finiteStrainMP::volumetricStiffness();
        double kv_error = fabs(kv-kvref);
        isok = (kv_error/fabs(kvref) < 1e-4);
        of << "\n   8. 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   9. 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;
    }


    {
        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   10. 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;
}
