//#####################################################################
// Copyright 2003-2020, Ron Fedkiw, Geoffrey Irving, Igor Neverov, Craig Schroeder, Eftychios Sifakis, Yunxin Sun, Joseph Teran.
// 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/NEO_HOOKEAN.h>
#include <Core/Math_Tools/cube.h>
#include <Core/Math_Tools/pow.h>
#include <Core/Math_Tools/Robust_Functions.h>
#include <Core/Matrices/DIAGONAL_MATRIX.h>
#include <Core/Matrices/MATRIX.h>
#include <Core/Matrices/SYMMETRIC_MATRIX.h>
#include <Deformables/Constitutive_Models/DIAGONALIZED_ISOTROPIC_STRESS_DERIVATIVE.h>
using namespace PhysBAM;
//#####################################################################
// Constructor
//#####################################################################
template<class T,int d> NEO_HOOKEAN<T,d>::
NEO_HOOKEAN(const T youngs_modulus_input,const T poissons_ratio_input,const T Rayleigh_coefficient,const T failure_threshold_input)
    :youngs_modulus(youngs_modulus_input),poissons_ratio(poissons_ratio_input),failure_threshold(failure_threshold_input),use_constant_ife(false)
{
    assert(poissons_ratio>-1&&poissons_ratio<.5);
    constant_lambda=youngs_modulus*poissons_ratio/((1+poissons_ratio)*(1-2*poissons_ratio));
    constant_mu=youngs_modulus/(2*(1+poissons_ratio));
    constant_alpha=Rayleigh_coefficient*constant_lambda;
    constant_beta=Rayleigh_coefficient*constant_mu;
    dth_root_failure_threshold=pow<1,d>(failure_threshold);
}
//#####################################################################
// Destructor
//#####################################################################
template<class T,int d> NEO_HOOKEAN<T,d>::
~NEO_HOOKEAN()
{
}
//#####################################################################
// Function Clamp_To_Hyperbola
//#####################################################################
template<class T,int d> DIAGONAL_MATRIX<T,1> NEO_HOOKEAN<T,d>::
Clamp_To_Hyperbola(const DIAGONAL_MATRIX<T,1>& F) const
{
    if(sqr(F.x.x)>failure_threshold) return DIAGONAL_MATRIX<T,1>(failure_threshold);
    else return DIAGONAL_MATRIX<T,1>(dth_root_failure_threshold);
}
//#####################################################################
// Function Clamp_To_Hyperbola
//#####################################################################
template<class T,int d> DIAGONAL_MATRIX<T,2> NEO_HOOKEAN<T,d>::
Clamp_To_Hyperbola(const DIAGONAL_MATRIX<T,2>& F) const
{
    if(sqr(F.x.x)>failure_threshold) return DIAGONAL_MATRIX<T,2>(F.x.x,failure_threshold/F.x.x);
    else return DIAGONAL_MATRIX<T,2>(dth_root_failure_threshold,dth_root_failure_threshold);
}
//#####################################################################
// Function Clamp_To_Hyperbola
//#####################################################################
template<class T,int d> DIAGONAL_MATRIX<T,3> NEO_HOOKEAN<T,d>::
Clamp_To_Hyperbola(const DIAGONAL_MATRIX<T,3>& F) const
{
    if(F.x.x*F.x.y*F.x.y>failure_threshold) return DIAGONAL_MATRIX<T,3>(F.x.x,F.x.y,failure_threshold/(F.x.x*F.x.y));
    else if(cube(F.x.x)>failure_threshold){
        T clamped=sqrt(failure_threshold/F.x.x);
        return DIAGONAL_MATRIX<T,3>(F.x.x,clamped,clamped);}
    else return DIAGONAL_MATRIX<T,3>(dth_root_failure_threshold,dth_root_failure_threshold,dth_root_failure_threshold);
}
//#####################################################################
// Function P_From_Strain
//#####################################################################
// clamp to hyperbola to avoid indefiniteness "automatically"
template<class T,int d> DIAGONAL_MATRIX<T,d> NEO_HOOKEAN<T,d>::
P_From_Strain(const DIAGONAL_MATRIX<T,d>& F,const int id) const
{
    T id_mu=Mu(id),id_lambda=Lambda(id),J=F.Determinant();
    if(J>=failure_threshold) return id_mu*F-(id_mu-id_lambda*log(J))*F.Inverse();
    DIAGONAL_MATRIX<T,d> F_clamp=Clamp_To_Hyperbola(F),dF,F_inverse=F_clamp.Inverse();
    if(!use_constant_ife) dF=F-F_clamp;
    T id_mu_minus_lambda_log_J=id_mu-id_lambda*log(failure_threshold);
    return id_mu*F+id_mu_minus_lambda_log_J*(sqr(F_inverse)*dF-F_inverse)+id_lambda*F_inverse.Double_Contract(dF)*F_inverse;
}
//#####################################################################
// Function P_From_Strain_Rate
//#####################################################################
template<class T,int d> MATRIX<T,d> NEO_HOOKEAN<T,d>::
P_From_Strain_Rate(const DIAGONAL_MATRIX<T,d>& F,const MATRIX<T,d>& F_dot,const int id) const
{
    T id_alpha=Alpha(id),id_beta=Beta(id);
    SYMMETRIC_MATRIX<T,d> strain_rate=F_dot.Symmetric_Part(); // use linear damping because of problems with inverting elements...
    return 2*id_beta*strain_rate+id_alpha*strain_rate.Trace();
}
//#####################################################################
// Function Isotropic_Stress_Derivative_Helper
//#####################################################################
template<class T,int d> void NEO_HOOKEAN<T,d>::
Isotropic_Stress_Derivative_Helper(const DIAGONAL_MATRIX<T,1>& F,DIAGONALIZED_ISOTROPIC_STRESS_DERIVATIVE<VECTOR<T,1> >& dP_dF,const int id) const
{
    T id_mu=Mu(id),id_lambda=Lambda(id);
    DIAGONAL_MATRIX<T,1> F_inverse=(F.Determinant()>=failure_threshold?F:Clamp_To_Hyperbola(F)).Inverse();
    T mu_minus_lambda_logJ=id_mu+id_lambda*log(F_inverse.Determinant());
    SYMMETRIC_MATRIX<T,1> F_inverse_outer=Outer_Product(F_inverse.Vector());
    dP_dF.H(0,0)=id_mu+(id_lambda+mu_minus_lambda_logJ)*F_inverse_outer(0,0);//alpha+beta+gamma
    if(enforce_definiteness) dP_dF.Enforce_Definiteness();
}
//#####################################################################
// Function Isotropic_Stress_Derivative_Helper
//#####################################################################
template<class T,int d> void NEO_HOOKEAN<T,d>::
Isotropic_Stress_Derivative_Helper(const DIAGONAL_MATRIX<T,2>& F,DIAGONALIZED_ISOTROPIC_STRESS_DERIVATIVE<VECTOR<T,2> >& dP_dF,const int id) const
{
    T id_mu=Mu(id),id_lambda=Lambda(id);
    DIAGONAL_MATRIX<T,2> F_inverse=(F.Determinant()>=failure_threshold?F:Clamp_To_Hyperbola(F)).Inverse();
    T mu_minus_lambda_logJ=id_mu+id_lambda*log(F_inverse.Determinant());
    SYMMETRIC_MATRIX<T,2> F_inverse_outer=Outer_Product(F_inverse.Vector());
    dP_dF.H(0,0)=id_mu+(id_lambda+mu_minus_lambda_logJ)*F_inverse_outer(0,0);//alpha+beta+gamma
    dP_dF.H(1,1)=id_mu+(id_lambda+mu_minus_lambda_logJ)*F_inverse_outer(1,1);
    dP_dF.H(1,0)=id_lambda*F_inverse_outer(1,0);//gamma
    dP_dF.B(0)=id_mu;//alpha
    dP_dF.C(0)=mu_minus_lambda_logJ*F_inverse_outer(1,0);//beta
    if(enforce_definiteness) dP_dF.Enforce_Definiteness();
}
//#####################################################################
// Function Isotropic_Stress_Derivative_Helper
//#####################################################################
template<class T,int d> void NEO_HOOKEAN<T,d>::
Isotropic_Stress_Derivative_Helper(const DIAGONAL_MATRIX<T,3>& F,DIAGONALIZED_ISOTROPIC_STRESS_DERIVATIVE<VECTOR<T,3> >& dPi_dF,const int id) const
{
    T id_mu=Mu(id),id_lambda=Lambda(id);
    DIAGONAL_MATRIX<T,3> F_inverse=(F.Determinant()>=failure_threshold?F:Clamp_To_Hyperbola(F)).Inverse();
    T mu_minus_lambda_logJ=id_mu+id_lambda*log(F_inverse.Determinant());
    SYMMETRIC_MATRIX<T,3> F_inverse_outer=Outer_Product(F_inverse.Vector());
    dPi_dF.H(0,0)=id_mu+(id_lambda+mu_minus_lambda_logJ)*F_inverse_outer(0,0);
    dPi_dF.H(1,1)=id_mu+(id_lambda+mu_minus_lambda_logJ)*F_inverse_outer(1,1);
    dPi_dF.H(2,2)=id_mu+(id_lambda+mu_minus_lambda_logJ)*F_inverse_outer(2,2);
    dPi_dF.H(1,0)=id_lambda*F_inverse_outer(1,0);
    dPi_dF.H(2,0)=id_lambda*F_inverse_outer(2,0);
    dPi_dF.H(2,1)=id_lambda*F_inverse_outer(2,1);
    dPi_dF.B(2)=id_mu;dPi_dF.B(1)=id_mu;dPi_dF.B(0)=id_mu;
    dPi_dF.C(2)=mu_minus_lambda_logJ*F_inverse_outer(1,0);
    dPi_dF.C(1)=mu_minus_lambda_logJ*F_inverse_outer(2,0);
    dPi_dF.C(0)=mu_minus_lambda_logJ*F_inverse_outer(2,1);
    if(enforce_definiteness) dPi_dF.Enforce_Definiteness();
}
//#####################################################################
// Function Isotropic_Stress_Derivative
//#####################################################################
template<class T,int d> void NEO_HOOKEAN<T,d>::
Isotropic_Stress_Derivative(const DIAGONAL_MATRIX<T,d>& F,DIAGONALIZED_ISOTROPIC_STRESS_DERIVATIVE<TV>& dPi_dF,const int id) const
{
    Isotropic_Stress_Derivative_Helper(F,dPi_dF,id);
}
//#####################################################################
// Function Energy_Density
//#####################################################################
template<class T,int d> T NEO_HOOKEAN<T,d>::
Energy_Density(const DIAGONAL_MATRIX<T,d>& F,const int id) const
{
    T id_mu=Mu(id),id_lambda=Lambda(id);
    T I1=(F*F.Transposed()).Trace(),J=F.Determinant();
    if(J<=0) return (T)1e16; // TODO: Do something smarter here.
    T log_J=log(J);
    return id_mu*((T).5*(I1-TV::m)-log_J)+(T).5*id_lambda*sqr(log_J);
}
//#####################################################################
// Function Robust_Divided_Pressure
//#####################################################################
template<class T,int d> T NEO_HOOKEAN<T,d>::
Robust_Divided_Pressure(T J,const int id) const
{
    T id_mu=Mu(id),id_lambda=Lambda(id);
    PHYSBAM_ASSERT(!id_mu);
    return id_lambda*log1p_x_over_x(J-1)/J;
}
//#####################################################################
// Function Pressure_Bound
//#####################################################################
template<class T,int d> T NEO_HOOKEAN<T,d>::
Pressure_Bound(T J,const int id) const
{
    T id_mu=Mu(id),id_lambda=Lambda(id);
    PHYSBAM_ASSERT(!id_mu);
    return id_lambda;
}
namespace PhysBAM{
template class NEO_HOOKEAN<float,1>;
template class NEO_HOOKEAN<float,2>;
template class NEO_HOOKEAN<float,3>;
template class NEO_HOOKEAN<double,1>;
template class NEO_HOOKEAN<double,2>;
template class NEO_HOOKEAN<double,3>;
}
