//#####################################################################
// Copyright 2004-2007, Geoffrey Irving, Igor Neverov, Eftychios Sifakis, Joseph Teran.
// This file is part of PhysBAM whose distribution is governed by the license contained in the accompanying file PHYSBAM_COPYRIGHT.txt.
//#####################################################################
// Class DIAGONALIZED_ISOTROPIC_STRESS_DERIVATIVE
//#####################################################################
#include <Deformables/Constitutive_Models/DIAGONALIZED_ISOTROPIC_STRESS_DERIVATIVE.h>
#include <Core/Matrices/DIAGONAL_MATRIX.h>
#include <Core/Matrices/MATRIX.h>
#include <Core/Matrices/SYMMETRIC_MATRIX.h>
#include <Core/Vectors/VECTOR.h>
namespace PhysBAM{
//#####################################################################
// Function Enforce_Definiteness
//#####################################################################
template<class TV> void DIAGONALIZED_ISOTROPIC_STRESS_DERIVATIVE<TV>::
Enforce_Definiteness()
{
    H=H.Positive_Definite_Part();
    for(int i=0;i<B.m;i++)
    {
        SYMMETRIC_MATRIX<T,2> A(B(i),C(i),B(i));
        A=A.Positive_Definite_Part();
        B(i)=A(0,0);
        C(i)=A(0,1);
    }
}
//#####################################################################
// Function Compute_From_Singular_Value_Derivatives
//#####################################################################
template<class TV> void DIAGONALIZED_ISOTROPIC_STRESS_DERIVATIVE<TV>::
Compute_From_Singular_Value_Derivatives(const DIAGONAL_MATRIX<T,TV::m>& F,const TV& dE_ds,const SYMMETRIC_MATRIX<T,TV::m>& dE_dsds)
{
    H=dE_dsds;
    for(int i=0;i<B.m;i++){
        int j=(i+1)%TV::m,k=(j+1)%TV::m;
        T d=1/(sqr(F.x(j))-sqr(F.x(k)));
        B(i)=(-dE_ds(k)*F.x(k)+dE_ds(j)*F.x(j))*d;
        C(i)=(-dE_ds(k)*F.x(j)+dE_ds(j)*F.x(k))*d;}
}
//#####################################################################
// Function Compute_From_Derivatives_And_Differences
//#####################################################################
template<class TV> void DIAGONALIZED_ISOTROPIC_STRESS_DERIVATIVE<TV>::
Compute_From_Derivatives_And_Differences(const SYMMETRIC_MATRIX<T,TV::m>& hess,const typename TV::SPIN& div_diff,const typename TV::SPIN& div_sum)
{
    H=hess;
    B=(T).5*(div_diff+div_sum);
    C=(T).5*(div_diff-div_sum);
}
//#####################################################################
// Function Contract_03
//#####################################################################
template<class TV> auto DIAGONALIZED_ISOTROPIC_STRESS_DERIVATIVE<TV>::
Contract_03(const MATRIX<T,TV::m>& M) const -> MATRIX<T,TV::m>
{
    if constexpr(TV::m==3)
        return MATRIX<T,3>(H(0,0)*M(0,0)+C(2)*M(1,1)+C(1)*M(2,2),B(2)*M(0,1)+H(1,0)*M(1,0),B(1)*M(0,2)+H(2,0)*M(2,0),
            H(0,1)*M(0,1)+B(2)*M(1,0),C(2)*M(0,0)+H(1,1)*M(1,1)+C(0)*M(2,2),B(0)*M(1,2)+H(2,1)*M(2,1),
            H(0,2)*M(0,2)+B(1)*M(2,0),H(1,2)*M(1,2)+B(0)*M(2,1),C(1)*M(0,0)+C(0)*M(1,1)+H(2,2)*M(2,2));
    else if constexpr(TV::m==2)
        return MATRIX<T,2>(H(0,0)*M(0,0)+C(0)*M(1,1),B(0)*M(0,1)+H(1,0)*M(1,0),H(0,1)*M(0,1)+B(0)*M(1,0),C(0)*M(0,0)+H(1,1)*M(1,1));
    else if constexpr(TV::m==1)
        return MATRIX<T,1>(H(0,0)*M(0,0));
    else return {};
}
//#####################################################################
// Function Contract_13
//#####################################################################
template<class TV> auto DIAGONALIZED_ISOTROPIC_STRESS_DERIVATIVE<TV>::
Contract_13(const MATRIX<T,TV::m>& M) const -> MATRIX<T,TV::m>
{
    if constexpr(TV::m==3)
        return MATRIX<T,3>(H(0,0)*M(0,0)+B(2)*M(1,1)+B(1)*M(2,2),H(1,0)*M(1,0)+C(2)*M(0,1),H(2,0)*M(2,0)+C(1)*M(0,2),
            H(0,1)*M(0,1)+C(2)*M(1,0),H(1,1)*M(1,1)+B(2)*M(0,0)+B(0)*M(2,2),H(2,1)*M(2,1)+C(0)*M(1,2),
            H(0,2)*M(0,2)+C(1)*M(2,0),H(1,2)*M(1,2)+C(0)*M(2,1),H(2,2)*M(2,2)+B(1)*M(0,0)+B(0)*M(1,1));
    else if constexpr(TV::m==2)
        return MATRIX<T,2>(H(0,0)*M(0,0)+B(0)*M(1,1), H(1,0)*M(1,0)+C(0)*M(0,1), H(0,1)*M(0,1)+C(0)*M(1,0), H(1,1)*M(1,1)+B(0)*M(0,0));
    else if constexpr(TV::m==1)
        return MATRIX<T,1>(H(0,0)*M.x[0]);
    else return {};
}
//#####################################################################
// Function Contract_23
//#####################################################################
template<class TV> auto DIAGONALIZED_ISOTROPIC_STRESS_DERIVATIVE<TV>::
Contract_23(const MATRIX<T,TV::m>& M) const -> MATRIX<T,TV::m>
{
    if constexpr(TV::m==3)
        return MATRIX<T,3>(H(0,0)*M.x[0]+H(1,0)*M.x[4]+H(2,0)*M.x[8],B(2)*M.x[1]+C(2)*M.x[3],B(1)*M.x[2]+C(1)*M.x[6],
            C(2)*M.x[1]+B(2)*M.x[3],H(1,0)*M.x[0]+H(1,1)*M.x[4]+H(2,1)*M.x[8],B(0)*M.x[5]+C(0)*M.x[7],
            C(1)*M.x[2]+B(1)*M.x[6],C(0)*M.x[5]+B(0)*M.x[7],H(2,0)*M.x[0]+H(2,1)*M.x[4]+H(2,2)*M.x[8]);
    else if constexpr(TV::m==2)
        return MATRIX<T,2>(H(0,0)*M.x[0]+H(1,0)*M.x[3],B(0)*M.x[1]+C(0)*M.x[2],C(0)*M.x[1]+B(0)*M.x[2],H(1,0)*M.x[0]+H(1,1)*M.x[3]);
    else if constexpr(TV::m==1)
        return MATRIX<T,1>(H(0,0)*M.x[0]);
    else return {};
}
//#####################################################################
// Function Sigma
//#####################################################################
template<class T,int d> SYMMETRIC_MATRIX<T,d>
Sigma(const DIAGONAL_MATRIX<T,d>& P,const DIAGONAL_MATRIX<T,d>& F,const MATRIX<T,d>& U)
{
    return Conjugate(U,P.Times_Transpose(F/F.Determinant()));
}
//#####################################################################
// Function dSigma_dF_Times
//#####################################################################
template<class TV,class T,int d> SYMMETRIC_MATRIX<T,d>
dSigma_dF_Times(const DIAGONALIZED_ISOTROPIC_STRESS_DERIVATIVE<TV>& dP_dF,
    const DIAGONAL_MATRIX<T,TV::m>& P,const DIAGONAL_MATRIX<T,TV::m>& F,
    const MATRIX<T,TV::m>& U,const MATRIX<T,TV::m>& V,const MATRIX<T,TV::m>& M)
{
    auto J=F.Determinant();
    auto G=F.Inverse();
    auto K=U.Transpose_Times(M)*V;
    auto a=(K*G).Trace();
    auto N=K-a*F;
    auto Q=dP_dF.Contract_23(K);
    auto R=Q.Times_Transpose(F)+P.Times_Transpose(N);
    auto Z=(U*R).Times_Transpose(U);
    return Z.Symmetric_Part()/J;
}
template class DIAGONALIZED_ISOTROPIC_STRESS_DERIVATIVE<VECTOR<float,1> >;
template class DIAGONALIZED_ISOTROPIC_STRESS_DERIVATIVE<VECTOR<float,2> >;
template class DIAGONALIZED_ISOTROPIC_STRESS_DERIVATIVE<VECTOR<float,3> >;
template class DIAGONALIZED_ISOTROPIC_STRESS_DERIVATIVE<VECTOR<double,1> >;
template class DIAGONALIZED_ISOTROPIC_STRESS_DERIVATIVE<VECTOR<double,2> >;
template class DIAGONALIZED_ISOTROPIC_STRESS_DERIVATIVE<VECTOR<double,3> >;
template SYMMETRIC_MATRIX<float,2> dSigma_dF_Times<VECTOR<float,2>,float,2>(DIAGONALIZED_ISOTROPIC_STRESS_DERIVATIVE<VECTOR<float,2> > const&,DIAGONAL_MATRIX<float,2> const&,DIAGONAL_MATRIX<float,2> const&,MATRIX<float,2,2> const&,MATRIX<float,2,2> const&,MATRIX<float,2,2> const&);
template SYMMETRIC_MATRIX<float,3> dSigma_dF_Times<VECTOR<float,3>,float,3>(DIAGONALIZED_ISOTROPIC_STRESS_DERIVATIVE<VECTOR<float,3> > const&,DIAGONAL_MATRIX<float,3> const&,DIAGONAL_MATRIX<float,3> const&,MATRIX<float,3,3> const&,MATRIX<float,3,3> const&,MATRIX<float,3,3> const&);
template SYMMETRIC_MATRIX<double,2> dSigma_dF_Times<VECTOR<double,2>,double,2>(DIAGONALIZED_ISOTROPIC_STRESS_DERIVATIVE<VECTOR<double,2> > const&,DIAGONAL_MATRIX<double,2> const&,DIAGONAL_MATRIX<double,2> const&,MATRIX<double,2,2> const&,MATRIX<double,2,2> const&,MATRIX<double,2,2> const&);
template SYMMETRIC_MATRIX<double,3> dSigma_dF_Times<VECTOR<double,3>,double,3>(DIAGONALIZED_ISOTROPIC_STRESS_DERIVATIVE<VECTOR<double,3> > const&,DIAGONAL_MATRIX<double,3> const&,DIAGONAL_MATRIX<double,3> const&,MATRIX<double,3,3> const&,MATRIX<double,3,3> const&,MATRIX<double,3,3> const&);
template SYMMETRIC_MATRIX<double,2> Sigma<double,2>(DIAGONAL_MATRIX<double,2> const&,DIAGONAL_MATRIX<double,2> const&,MATRIX<double,2,2> const&);
template SYMMETRIC_MATRIX<double,3> Sigma<double,3>(DIAGONAL_MATRIX<double,3> const&,DIAGONAL_MATRIX<double,3> const&,MATRIX<double,3,3> const&);
template SYMMETRIC_MATRIX<float,2> Sigma<float,2>(DIAGONAL_MATRIX<float,2> const&,DIAGONAL_MATRIX<float,2> const&,MATRIX<float,2,2> const&);
template SYMMETRIC_MATRIX<float,3> Sigma<float,3>(DIAGONAL_MATRIX<float,3> const&,DIAGONAL_MATRIX<float,3> const&,MATRIX<float,3,3> const&);
template SYMMETRIC_MATRIX<double,1> Sigma<double,1>(DIAGONAL_MATRIX<double,1> const&,DIAGONAL_MATRIX<double,1> const&,MATRIX<double,1,1> const&);
template SYMMETRIC_MATRIX<float,1> Sigma<float,1>(DIAGONAL_MATRIX<float,1> const&,DIAGONAL_MATRIX<float,1> const&,MATRIX<float,1,1> const&);
}
