//#####################################################################
// 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.
//#####################################################################
#include <Deformables/Constitutive_Models/GENERAL_STRESS_DERIVATIVE.h>
#include <Core/Log/LOG.h>
#include <Core/Matrices/MATRIX.h>
#include <Core/Matrices/SYMMETRIC_MATRIX.h>
#include <Core/Vectors/VECTOR.h>
namespace PhysBAM{
//#####################################################################
// Function Contract_03
//#####################################################################
template<class TV> auto GENERAL_STRESS_DERIVATIVE<TV>::
Contract_03(const MATRIX<T,TV::m>& M) const -> MATRIX<T,TV::m> // solving Cjk = Mijkl * Mil
{
    MATRIX<T,TV::m> R;
    for(int i=0;i<TV::m;i++)
        for(int j=0;j<TV::m;j++)
            R.Add_Row(j,A(i,j)*M.Get_Row(i));
    return R;
}
//#####################################################################
// Function Contract_13
//#####################################################################
template<class TV> auto GENERAL_STRESS_DERIVATIVE<TV>::
Contract_13(const MATRIX<T,TV::m>& M) const -> MATRIX<T,TV::m> // solving Cik = Mijkl * Mjl
{
    MATRIX<T,TV::m> R;
    for(int i=0;i<TV::m;i++)
        for(int j=0;j<TV::m;j++)
            R.Add_Row(i,A(i,j)*M.Get_Row(j));
    return R;
}
//#####################################################################
// Function Contract_23
//#####################################################################
template<class TV> auto GENERAL_STRESS_DERIVATIVE<TV>::
Contract_23(const MATRIX<T,TV::m>& M) const -> MATRIX<T,TV::m> // solving Cij = Mijkl * Mkl
{
    MATRIX<T,TV::m> R;
    for(int i=0;i<TV::m;i++)
        for(int j=0;j<TV::m;j++)
            R(i,j)=A(i,j).Double_Contract(M);
    return R;
}
//#####################################################################
// Function Sigma
//#####################################################################
template<class T,int d> MATRIX<T,d>
Sigma(const MATRIX<T,d>& P,const MATRIX<T,d>& F)
{
    return Times_Transpose(P,F)/F.Determinant();
}
//#####################################################################
// Function dSigma_dF_Times
//#####################################################################
template<class TV,class T,int d> MATRIX<T,d>
dSigma_dF_Times(const GENERAL_STRESS_DERIVATIVE<TV>& dP_dF,
    const MATRIX<T,TV::m>& P,const MATRIX<T,TV::m>& F,const MATRIX<T,TV::m>& M)
{
    auto J=F.Determinant();
    auto G=F.Inverse();
    auto a=(M*G).Trace();
    auto N=M-a*F;
    auto Q=dP_dF.Contract_23(M);
    auto R=Q.Times_Transpose(F)+P.Times_Transpose(N);
    return R/J;
}
template class GENERAL_STRESS_DERIVATIVE<VECTOR<double,1> >;
template class GENERAL_STRESS_DERIVATIVE<VECTOR<double,2> >;
template class GENERAL_STRESS_DERIVATIVE<VECTOR<double,3> >;
template class GENERAL_STRESS_DERIVATIVE<VECTOR<float,1> >;
template class GENERAL_STRESS_DERIVATIVE<VECTOR<float,2> >;
template class GENERAL_STRESS_DERIVATIVE<VECTOR<float,3> >;
template MATRIX<double,2> Sigma<double,2>(MATRIX<double,2,2> const&,MATRIX<double,2,2> const&);
template MATRIX<double,2> dSigma_dF_Times<VECTOR<double,2>,double,2>(GENERAL_STRESS_DERIVATIVE<VECTOR<double,2> > const&,MATRIX<double,VECTOR<double,2>::m,VECTOR<double,2>::m> const&,MATRIX<double,VECTOR<double,2>::m,VECTOR<double,2>::m> const&,MATRIX<double,VECTOR<double,2>::m,VECTOR<double,2>::m> const&);
template MATRIX<double,3> Sigma<double,3>(MATRIX<double,3,3> const&,MATRIX<double,3,3> const&);
template MATRIX<double,3> dSigma_dF_Times<VECTOR<double,3>,double,3>(GENERAL_STRESS_DERIVATIVE<VECTOR<double,3> > const&,MATRIX<double,VECTOR<double,3>::m,VECTOR<double,3>::m> const&,MATRIX<double,VECTOR<double,3>::m,VECTOR<double,3>::m> const&,MATRIX<double,VECTOR<double,3>::m,VECTOR<double,3>::m> const&);
template MATRIX<float,2> Sigma<float,2>(MATRIX<float,2,2> const&,MATRIX<float,2,2> const&);
template MATRIX<float,2> dSigma_dF_Times<VECTOR<float,2>,float,2>(GENERAL_STRESS_DERIVATIVE<VECTOR<float,2> > const&,MATRIX<float,VECTOR<float,2>::m,VECTOR<float,2>::m> const&,MATRIX<float,VECTOR<float,2>::m,VECTOR<float,2>::m> const&,MATRIX<float,VECTOR<float,2>::m,VECTOR<float,2>::m> const&);
template MATRIX<float,3> Sigma<float,3>(MATRIX<float,3,3> const&,MATRIX<float,3,3> const&);
template MATRIX<float,3> dSigma_dF_Times<VECTOR<float,3>,float,3>(GENERAL_STRESS_DERIVATIVE<VECTOR<float,3> > const&,MATRIX<float,VECTOR<float,3>::m,VECTOR<float,3>::m> const&,MATRIX<float,VECTOR<float,3>::m,VECTOR<float,3>::m> const&,MATRIX<float,VECTOR<float,3>::m,VECTOR<float,3>::m> const&);
template MATRIX<double,1,1> Sigma<double,1>(MATRIX<double,1,1> const&,MATRIX<double,1,1> const&);
template MATRIX<float,1,1> Sigma<float,1>(MATRIX<float,1,1> const&,MATRIX<float,1,1> const&);
}

