/****************************************************************************
*
*                                 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 <math.h>
#include <string.h>
#include <iomanip>
#include <stdexcept>

#include "fisotropic.h"

#define NEWEM 0

using namespace std;
using namespace muesli;


fisotropicMP::fisotropicMP(const f_invariants& m)
:
finiteStrainMP(m)
{
    invar_n[0] = 3.0;
    invar_n[1] = 0.0;
    invar_n[2] = 1.0;
    invar_c    = invar_n;
    dW_n.setZero();
    dW_c  = dW_n;
    ddW_n.setZero();
    ddW_c = ddW_n;
    for (size_t a=0; a<8; a++) G_c[a] = G_n[a];

    Fn     = itensor::identity();
    Fc     = itensor::identity();
}



// Cauchy stress in terms of principal invariants and the partial derivatives of W w/r to the invariants
void fisotropicMP::CauchyStress(istensor& sigma) const
{
    double Inv1 = invar_c[0];
    double Inv3 = invar_c[2];

    double J  = sqrt( Inv3 );
    double c0 =  2.0/J * ( dW_c(0) + Inv1 * dW_c(1) );
    double c1 = -2.0/J *   dW_c(1);
    double c2 =  2.0/J *   Inv3 * dW_c(2);

    sigma.setZero();
    istensor B = istensor::tensorTimesTensorTransposed(Fc);
    if (c0 != 0.0) sigma += c0 * B;
    if (c1 != 0.0) sigma += c1 * istensor::symmetricPartOf(B * B);
    if (c2 != 0.0) sigma += c2 * istensor::identity();
}




void fisotropicMP::commitCurrentState()
{
    finiteStrainMP::commitCurrentState();

    dW_n    = dW_c;
    ddW_n   = ddW_c;
    invar_n = invar_c;
    for (size_t a=0; a<8; a++) G_n[a] = G_c[a];
}




// See Doghri, pg 376
void fisotropicMP::computeGammaCoefficients(const ivector& invariants, const ivector& dW, const istensor& ddW, double G[8])
{
    G[0] =  4.0 * invariants(2) * ( dW(2) + invariants(2) * ddW(2,2) );
    G[1] =  4.0 * invariants(2) * ( ddW(0,2) + invariants(0) * ddW(1,2) );
    G[2] = -4.0 * invariants(2) * ddW(1,2);
    G[3] =  4.0 * (dW(1) + ddW(0,0) + 2.0*invariants(0)*ddW(0,1) + invariants(0)*invariants(0)*ddW(1,1) );
    G[4] = -4.0 * (ddW(0,1) + invariants(0) * ddW(1,1) );
    G[5] =  4.0 * ddW(1,1);
    G[6] = -4.0 * dW(1);
    G[7] = -4.0 * invariants(2) * dW(2);
}



void fisotropicMP::contractWithDeviatoricTangent(const ivector &v1, const ivector& v2, itensor &T) const
{
    istensor B  = istensor::tensorTimesTensorTransposed(Fc);
    istensor BB = B.squared(), Z;
    ivector  Bv1(B*v1), Bv2(B*v2), BBv1(B*Bv1), BBv2(B*Bv2);

    double   Btr = B.trace();
    double   BBtr = BB.trace();

    // construct tangent only with relevant terms
    T.setZero();
    Z.setZero();
    const double tol = theFiniteStrainMaterial.characteristicStiffness()/1e12;
    
    if ( fabs(G_c[0]) > tol )  T += G_c[0] * Z;
    if ( fabs(G_c[1]) > tol )  T += G_c[1] * Z;
    if ( fabs(G_c[2]) > tol )  T += G_c[2] * Z;
    if ( fabs(G_c[3]) > tol )  T += G_c[3] * itensor::dyadic(Bv1, Bv2) - G_c[3] * 1.0/3.0 * Btr * itensor::dyadic(Bv1, v2) - G_c[3] * 1.0/3.0 * Btr * itensor::dyadic(v1, Bv2) + G_c[3] * 1.0/9.0 * pow(Btr,2.0) * itensor::dyadic(v1, v2);
    if ( fabs(G_c[4]) > tol )  T += G_c[4] * itensor::dyadic(Bv1, BBv2) + G_c[4] * itensor::dyadic(BBv1, Bv2) - G_c[4] * 1.0/3.0 * BBtr * itensor::dyadic(Bv1, v2) - G_c[4] * 1.0/3.0 * BBtr * itensor::dyadic(v1, Bv2) - G_c[4] *  1.0/3.0 * Btr * itensor::dyadic(BBv1, v2) - G_c[4] * 1.0/3.0 * Btr * itensor::dyadic(v1, BBv2) + G_c[4] * 2.0/9.0 * BBtr * Btr * itensor::dyadic(v1, v2);
    if ( fabs(G_c[5]) > tol )  T += G_c[5] * itensor::dyadic(BBv1, BBv2) - G_c[5] * 1.0/3.0 * BBtr * itensor::dyadic(BBv1, v2) - G_c[5] * 1.0/3.0 * BBtr * itensor::dyadic(v1, BBv2) + G_c[5] * 1.0/9.0 * pow(BBtr,2.0) * itensor::dyadic(v1, v2);
    if ( fabs(G_c[6]) > tol )  T += G_c[6] * 0.5 * itensor::dyadic(Bv2, Bv1) + G_c[6] * 0.5 * v1.dot(Bv2) * B - G_c[6] * 1.0/3.0 * itensor::dyadic(v1, BBv2) - G_c[6] * 1.0/3.0 * itensor::dyadic(BBv1, v2) + G_c[6] * 1.0/9.0 * BBtr * itensor::dyadic(v1, v2);
    if ( fabs(G_c[7]) > tol )  T += G_c[7] * 0.5 * itensor::dyadic(v2, v1)   + G_c[7] * 0.5 * v1.dot(v2)  * itensor::identity() - G_c[7] * 1.0/3.0 * itensor::dyadic(v1, v2);

    double iJ = 1.0/sqrt(invar_c[2]);
    T *= iJ;
}




void fisotropicMP::contractWithMixedTangent(istensor& CM) const
{
    istensor B = istensor::tensorTimesTensorTransposed(Fc);
    istensor BB = B.squared(), Z;

    // construct tangent only with relevant terms
    CM.setZero();
    Z.setZero();
    CM += G_c[0] * Z;
    CM += G_c[1] * 3.0 * B - G_c[1] * B.trace() * istensor::identity();
    CM += G_c[2] * 3.0 * BB - G_c[2] * BB.trace() * istensor::identity();
    CM += G_c[3] * B.trace() * B - G_c[3] * 1.0/3.0 * pow(B.trace(),2.0) * istensor::identity();
    CM += G_c[4] * BB.trace() * B + G_c[4] * B.trace() * BB - G_c[4] * 2.0/3.0 * BB.trace() * B.trace() * istensor::identity();
    CM += G_c[5] * BB.trace() * BB - G_c[5] * 1.0/3.0 * pow(BB.trace(),2.0) * istensor::identity();
    CM += G_c[6] * BB - G_c[6] * 1.0/3.0 * BB.trace() * istensor::identity();
    CM += G_c[7] * Z;

    double iJ = 1.0/sqrt(invar_c[2]);
    CM *= iJ;
}




/* Given the fourth order tensor of elasticities c, and two vectors v, w
 * compute the second order tensor T with components
 *   T_ij = c_ipjq v_p w_q
 *
 *  Note that the result is not symmetric.
 */
void fisotropicMP::contractWithSpatialTangent(const ivector& v1, const ivector& v2, itensor& T) const
{
    istensor B = istensor::tensorTimesTensorTransposed(Fc);
    ivector  Bv1(B*v1), Bv2(B*v2), BBv1(B*Bv1), BBv2(B*Bv2);

    // construct tangent only with relevant terms
    T.setZero();
    double tol = theFiniteStrainMaterial.characteristicStiffness()/1e12;
    if ( fabs(G_c[0]) > tol )  T += G_c[0] * itensor::dyadic(v1, v2);
    if ( fabs(G_c[1]) > tol )  T += G_c[1] * (itensor::dyadic(Bv1, v2)  + itensor::dyadic(v1, Bv2) );
    if ( fabs(G_c[2]) > tol )  T += G_c[2] * (itensor::dyadic(BBv1, v2) + itensor::dyadic(v1, BBv2) );
    if ( fabs(G_c[3]) > tol )  T += G_c[3] * itensor::dyadic(Bv1, Bv2);
    if ( fabs(G_c[4]) > tol )  T += G_c[4] * (itensor::dyadic(Bv1, BBv2) + itensor::dyadic(BBv1, Bv2) );
    if ( fabs(G_c[5]) > tol )  T += G_c[5] * itensor::dyadic(BBv1, BBv2);
    if ( fabs(G_c[6]) > tol )  T += G_c[6] * 0.5 * itensor::dyadic(Bv2, Bv1) + G_c[6] * 0.5 * v1.dot(Bv2) * B;
    if ( fabs(G_c[7]) > tol )  T += G_c[7] * 0.5 * itensor::dyadic(v2, v1)  + G_c[7] * 0.5 * v1.dot(v2)  * itensor::identity();

    double iJ = 1.0/sqrt(invar_c[2]);
    T *= iJ;
}




void fisotropicMP::contractWithConvectedTangent(const ivector& V1, const ivector& V2, itensor& T) const
{
    const istensor C = istensor::tensorTransposedTimesTensor(Fc);
    istensor Cinv = C.inverse();

    ivector  CiV1(Cinv*V1), CiV2(Cinv*V2);
    ivector  CV1(C*V1), CV2(C*V2);

    // construct tangent only with relevant terms
    T  = G_c[0] * itensor::dyadic(CiV1, CiV2);
    T += G_c[1] * (itensor::dyadic(V1, CiV2)  + itensor::dyadic(CiV1, V2) );
    T += G_c[2] * (itensor::dyadic(CV1, CiV2) + itensor::dyadic(CiV1, CV2) );
    T += G_c[3] * itensor::dyadic(V1, V2);
    T += G_c[4] * (itensor::dyadic(V1, CV2) + itensor::dyadic(CV1, V2) );
    T += G_c[5] * itensor::dyadic(CV1, CV2);
    T += G_c[6] * 0.5 * itensor::dyadic(V2, V1) + G_c[6] * 0.5 * V1.dot(V2) * itensor::identity();
    T += G_c[7] * 0.5 * itensor::dyadic(CiV2, CiV1) + G_c[7] * 0.5 * V1.dot(CiV2) * Cinv;
}




void fisotropicMP::convectedTangent(itensor4& c) const
{
    const istensor id   = istensor::identity();
    const istensor C    = istensor::tensorTransposedTimesTensor(Fc);
    const istensor Cinv = C.inverse();

    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++)
                {
                    c(i,j,k,l) = G_c[0]*Cinv(i,j)*Cinv(k,l)
                    + G_c[1]*(id(i,j)*Cinv(k,l) + Cinv(i,j)*id(k,l))
                    + G_c[2]*(C(i,j)*Cinv(k,l) + Cinv(i,j)*C(k,l))
                    + G_c[3]*id(i,j)*id(k,l)
                    + G_c[4]*(id(i,j)*C(k,l) + C(i,j)*id(k,l))
                    + G_c[5]*C(i,j)*C(k,l)
                    + G_c[6]*0.5*(id(i,k)*id(j,l)+id(i,l)*id(j,k))
                    + G_c[7]*0.5*(Cinv(i,k)*Cinv(j,l)+Cinv(i,l)*Cinv(j,k));
                }
            }
        }
    }
}




void fisotropicMP::convectedTangentTimesSymmetricTensor(const istensor& T, istensor& CM) const
{
    istensor C        = istensor::tensorTransposedTimesTensor(Fc);
    istensor Cinv     = C.inverse();
    double   CdotT    = C.contract(T);
    double   CinvdotT = Cinv.contract(T);
    double   traceT   = T.trace();

    CM.setZero();
    if (G_c[0] != 0.0) CM += G_c[0] * CinvdotT * Cinv;
    if (G_c[1] != 0.0) CM += G_c[1] * CinvdotT * istensor::identity() + G_c[1] * traceT * Cinv;
    if (G_c[2] != 0.0) CM += G_c[2] * CinvdotT * C + G_c[2] * CdotT * Cinv;
    if (G_c[3] != 0.0) CM += G_c[3] * traceT * istensor::identity();
    if (G_c[4] != 0.0) CM += G_c[4] * CdotT * istensor::identity() + G_c[4] * traceT * C;
    if (G_c[5] != 0.0) CM += G_c[5] * CdotT * C;
    if (G_c[6] != 0.0) CM += G_c[6] * T;
    if (G_c[7] != 0.0) CM += G_c[7] * istensor::FSFt(Cinv, T);
}




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




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




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




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

    for (unsigned a=0; a<8; a++) state.theDouble.push_back(G_n[a]);
    state.theVector.push_back(invar_n);
    state.theVector.push_back(dW_n);
    state.theStensor.push_back(ddW_n);

    return state;
}




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

    for (unsigned a=0; a<8; a++) state.theDouble.push_back(G_c[a]);
    state.theVector.push_back(invar_c);
    state.theVector.push_back(dW_c);
    state.theStensor.push_back(ddW_c);
    
    return state;
}




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




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




void fisotropicMP::secondPiolaKirchhoffStress(istensor &S) const
{
    const istensor C = istensor::tensorTransposedTimesTensor(Fc);
    S = 2.0*(dW_c(0)+invar_c[0]*dW_c(1))*istensor::identity()
        - 2.0*dW_c(1)*C
        + 2.0*invar_c[2]*dW_c(2)*C.inverse();
}




void fisotropicMP::resetCurrentState()
{
    finiteStrainMP::resetCurrentState();

    invar_c = invar_n;
    dW_c    = dW_n;
    ddW_c   = ddW_n;
    for (size_t a=0; a<8; a++) G_c[a] = G_n[a];
}




void fisotropicMP::setRandom()
{
    finiteStrainMP::setRandom();
}




double fisotropicMP::volumetricStiffness() const
{
    const istensor B  = istensor::tensorTimesTensorTransposed(Fc);
    const istensor BB = istensor::symmetricPartOf(B*B);

    const double Btr  = B.trace();
    const double BBtr = BB.trace();
    double bulk = 0.0;

    // construct tangent only with relevant terms

    double tol = theFiniteStrainMaterial.characteristicStiffness()/1e12;
    if ( fabs(G_c[0]) > tol )  bulk += G_c[0] * 9.0;
    if ( fabs(G_c[1]) > tol )  bulk += G_c[1] * 6.0 * Btr;
    if ( fabs(G_c[2]) > tol )  bulk += G_c[2] * 6.0 * BBtr;
    if ( fabs(G_c[3]) > tol )  bulk += G_c[3] * Btr*Btr;
    if ( fabs(G_c[4]) > tol )  bulk += G_c[4] * 2.0 * BBtr * Btr;
    if ( fabs(G_c[5]) > tol )  bulk += G_c[5] * BBtr * BBtr;
    if ( fabs(G_c[6]) > tol )  bulk += G_c[6] * BBtr;
    if ( fabs(G_c[7]) > tol )  bulk += G_c[7] * 3.0;

    const double iJ = 1.0/sqrt(invar_c[2]);
    bulk *= iJ/9.0;

    return bulk;
}
