//#####################################################################
// 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_GENERAL.h>
#include <Core/Log/LOG.h>
#include <Core/Math_Tools/cube.h>
#include <Core/Math_Tools/pow.h>
#include <Core/Math_Tools/Robust_Functions.h>
#include <Core/Matrices/MATRIX.h>
#include <Core/Matrices/SYMMETRIC_MATRIX.h>
#include <Deformables/Constitutive_Models/GENERAL_STRESS_DERIVATIVE.h>
using namespace PhysBAM;
//#####################################################################
// Constructor
//#####################################################################
template<class T,int d> NEO_HOOKEAN_GENERAL<T,d>::
NEO_HOOKEAN_GENERAL(const T youngs_modulus,const T poissons_ratio)
{
    assert(poissons_ratio>-1&&poissons_ratio<.5);
    la=youngs_modulus*poissons_ratio/((1+poissons_ratio)*(1-2*poissons_ratio));
    mu=youngs_modulus/(2*(1+poissons_ratio));
}
//#####################################################################
// Destructor
//#####################################################################
template<class T,int d> NEO_HOOKEAN_GENERAL<T,d>::
~NEO_HOOKEAN_GENERAL()
{
}
//#####################################################################
// Function Stress
//#####################################################################
template<class T,int d> MATRIX<T,d> NEO_HOOKEAN_GENERAL<T,d>::
Stress(const MATRIX<T,d>& F) const
{
    T J=F.Determinant();
    return mu*F-(mu-la*log(J))/J*F.Cofactor_Matrix();
}
//#####################################################################
// Function Isotropic_Stress_Derivative
//#####################################################################
template<class T,int d> void NEO_HOOKEAN_GENERAL<T,d>::
Stress_Derivative(const MATRIX<T,d>& F,GENERAL_STRESS_DERIVATIVE<TV>& dP_dF) const
{
    T J=F.Determinant();
    T a=mu-la*log(J);
    auto G=F.Inverse();
    for(int i=0;i<d;i++)
        for(int j=0;j<d;j++)
        {
            MATRIX<T,d> M;
            M(i,j)=1;
            T dJ=G(j,i);
            T da=-la*dJ;
            auto dG=-G*M*G;
            auto dP=mu*M-da*G.Transposed()-a*dG.Transposed();
            dP_dF.A(i,j)=dP;
        }
}
//#####################################################################
// Function Energy_Density
//#####################################################################
template<class T,int d> T NEO_HOOKEAN_GENERAL<T,d>::
Energy_Density(const MATRIX<T,d>& F) const
{
    T I1=F.Frobenius_Norm_Squared(),J=F.Determinant();
    T log_J=log(J);
    return mu*((T).5*(I1-TV::m)-log_J)+(T).5*la*sqr(log_J);
}
namespace PhysBAM{
template class NEO_HOOKEAN_GENERAL<float,1>;
template class NEO_HOOKEAN_GENERAL<float,2>;
template class NEO_HOOKEAN_GENERAL<float,3>;
template class NEO_HOOKEAN_GENERAL<double,1>;
template class NEO_HOOKEAN_GENERAL<double,2>;
template class NEO_HOOKEAN_GENERAL<double,3>;
}
