//#####################################################################
// Copyright 2004-2005, Geoffrey Irving, Eftychios Sifakis.
// This file is part of PhysBAM whose distribution is governed by the license contained in the accompanying file PHYSBAM_COPYRIGHT.txt.
//#####################################################################
#include <Deformables/Constitutive_Models/MOONEY_RIVLIN_3D.h>
#include <Deformables/Constitutive_Models/DIAGONALIZED_ISOTROPIC_STRESS_DERIVATIVE.h>
namespace PhysBAM{
//#####################################################################
// Constructor
//#####################################################################
template<class T> MOONEY_RIVLIN_3D<T>::
MOONEY_RIVLIN_3D(const T mu_10_input,const T mu_01_input,const T kappa_input,const T Rayleigh_coefficient,const T failure_threshold_input)
    :mu_10(mu_10_input),mu_01(mu_01_input),kappa(kappa_input),failure_threshold(failure_threshold_input)
{
    constant_mu=2*(mu_10+mu_01);
    constant_lambda=((T)1/3)*(8*mu_01-4*mu_10)+kappa;
    constant_alpha=Rayleigh_coefficient*constant_lambda;
    constant_beta=Rayleigh_coefficient*constant_mu;
}
//#####################################################################
// Function P_From_Strain
//#####################################################################
template<class T> DIAGONAL_MATRIX<T,3> MOONEY_RIVLIN_3D<T>::
P_From_Strain(const DIAGONAL_MATRIX<T,3>& F,const int id) const
{
    DIAGONAL_MATRIX<T,3> F_threshold(Componentwise_Max(F.x,failure_threshold)),C=F_threshold*F_threshold,F_cube=C*F_threshold,F_inverse=F_threshold.Inverse();
    T I_C=C.Trace(),II_C=(C*C).Trace(),J=F_threshold.Determinant(),Jcc=pow(J,-((T)2/3));
    return (2*Jcc*(mu_10+Jcc*mu_01*I_C))*F_threshold-(2*Jcc*Jcc*mu_01)*F_cube+((kappa*log(J)-((T)2/3)*Jcc*(mu_10*I_C+Jcc*mu_01*(I_C*I_C-II_C))))*F_inverse;
}
//#####################################################################
// Function P_From_Strain_Rate
//#####################################################################
template<class T> MATRIX<T,3> MOONEY_RIVLIN_3D<T>::
P_From_Strain_Rate(const DIAGONAL_MATRIX<T,3>& F,const MATRIX<T,3>& F_dot,const int id) const
{
    T id_alpha=Alpha(id),id_beta=Beta(id);
    SYMMETRIC_MATRIX<T,3> strain_rate=F_dot.Symmetric_Part(); 
    return 2*id_beta*strain_rate+id_alpha*strain_rate.Trace();
}
//#####################################################################
// Function Isotropic_Stress_Derivative
//#####################################################################
template<class T> void MOONEY_RIVLIN_3D<T>::
Isotropic_Stress_Derivative(const DIAGONAL_MATRIX<T,3>& F,DIAGONALIZED_ISOTROPIC_STRESS_DERIVATIVE<VECTOR<T,3> >& dPi_dF,const int id) const
{
    DIAGONAL_MATRIX<T,3> F_threshold(Componentwise_Max(F.x,failure_threshold)),C=F_threshold*F_threshold,F_cube=C*F_threshold,F_inverse=F_threshold.Inverse();
    T I_C=C.Trace(),II_C=(C*C).Trace(),J=F_threshold.Determinant(),Jcc=pow(J,-((T)2/3));
    SYMMETRIC_MATRIX<T,3> alpha;
    alpha(0,0)=2*Jcc*(mu_10+Jcc*mu_01*(I_C-C.x.x-C.x.x));alpha(1,0)=2*Jcc*(mu_10+Jcc*mu_01*(I_C-C.x.x-C.x.x));alpha(2,0)=2*Jcc*(mu_10+Jcc*mu_01*(I_C-C.x.z-C.x.x));
    alpha(1,1)=2*Jcc*(mu_10+Jcc*mu_01*(I_C-C.x.y-C.x.y));alpha(2,1)=2*Jcc*(mu_10+Jcc*mu_01*(I_C-C.x.z-C.x.y));alpha(2,2)=2*Jcc*(mu_10+Jcc*mu_01*(I_C-C.x.z-C.x.z));
    SYMMETRIC_MATRIX<T,3> beta;
    beta(0,0)=(Jcc*((T)2/3)*(mu_10*I_C+Jcc*mu_01*(I_C*I_C-II_C))-kappa*log(J))*F_inverse.x.x*F_inverse.x.x-2*Jcc*Jcc*mu_01*F_threshold.x.x*F_threshold.x.x;
    beta(1,0)=(Jcc*((T)2/3)*(mu_10*I_C+Jcc*mu_01*(I_C*I_C-II_C))-kappa*log(J))*F_inverse.x.y*F_inverse.x.x-2*Jcc*Jcc*mu_01*F_threshold.x.y*F_threshold.x.x;
    beta(2,0)=(Jcc*((T)2/3)*(mu_10*I_C+Jcc*mu_01*(I_C*I_C-II_C))-kappa*log(J))*F_inverse.x.z*F_inverse.x.x-2*Jcc*Jcc*mu_01*F_threshold.x.z*F_threshold.x.x;
    beta(1,1)=(Jcc*((T)2/3)*(mu_10*I_C+Jcc*mu_01*(I_C*I_C-II_C))-kappa*log(J))*F_inverse.x.y*F_inverse.x.y-2*Jcc*Jcc*mu_01*F_threshold.x.y*F_threshold.x.y;
    beta(2,1)=(Jcc*((T)2/3)*(mu_10*I_C+Jcc*mu_01*(I_C*I_C-II_C))-kappa*log(J))*F_inverse.x.z*F_inverse.x.y-2*Jcc*Jcc*mu_01*F_threshold.x.z*F_threshold.x.y;
    beta(2,2)=(Jcc*((T)2/3)*(mu_10*I_C+Jcc*mu_01*(I_C*I_C-II_C))-kappa*log(J))*F_inverse.x.z*F_inverse.x.z-2*Jcc*Jcc*mu_01*F_threshold.x.z*F_threshold.x.z;
    SYMMETRIC_MATRIX<T,3> eta;
    eta(0,0)=4*Jcc*Jcc*mu_01;
    eta(2,0)=-Jcc*((T)1/3)*(4*mu_10+8*Jcc*mu_01*I_C);
    eta(2,1)=8*((T)1/3)*Jcc*Jcc*mu_01;
    eta(2,2)=Jcc*((T)1/9)*(4*mu_10*I_C+Jcc*8*mu_01*(I_C*I_C-II_C))+kappa;
    MATRIX<T,3> F_base{F_threshold.x,F_cube.x,F_inverse.x};
    SYMMETRIC_MATRIX<T,3> gamma=Conjugate(F_base,eta);
    dPi_dF.H(0,0)=alpha(0,0)+beta(0,0)+gamma(0,0);dPi_dF.H(1,1)=alpha(1,1)+beta(1,1)+gamma(1,1);dPi_dF.H(2,2)=alpha(2,2)+beta(2,2)+gamma(2,2);
    dPi_dF.H(1,0)=gamma(1,0);dPi_dF.H(2,0)=gamma(2,0);dPi_dF.H(2,1)=gamma(2,1);
    dPi_dF.B(2)=alpha(1,0);dPi_dF.B(1)=alpha(2,0);dPi_dF.B(0)=alpha(2,1);
    dPi_dF.C(2)=beta(1,0);dPi_dF.C(1)=beta(2,0);dPi_dF.C(0)=beta(2,0);
    if(enforce_definiteness) dPi_dF.Enforce_Definiteness();
}
//#####################################################################
// Function Energy_Density
//#####################################################################
template<class T> T MOONEY_RIVLIN_3D<T>::
Energy_Density(const DIAGONAL_MATRIX<T,3>& F,const int id) const
{
    DIAGONAL_MATRIX<T,3> F_threshold(Componentwise_Max(F.x,failure_threshold)),C=F_threshold*F_threshold;
    T I_C=C.Trace(),II_C=(C*C).Trace(),J=F_threshold.Determinant(),Jcc=pow(J,-((T)2/3));
    return (T)0.5*kappa*sqr(log(J))+(mu_10*I_C+mu_01*II_C*Jcc)*Jcc;
}
template class MOONEY_RIVLIN_3D<float>;
template class MOONEY_RIVLIN_3D<double>;
}
