//#####################################################################
// 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
//#####################################################################
#ifndef __DIAGONALIZED_ISOTROPIC_STRESS_DERIVATIVE__
#define __DIAGONALIZED_ISOTROPIC_STRESS_DERIVATIVE__

#include <Core/Matrices/MATRIX.h>
#include <Core/Matrices/SYMMETRIC_MATRIX.h>
#include <Core/Vectors/VECTOR.h>
namespace PhysBAM{
template<class TV>
class DIAGONALIZED_ISOTROPIC_STRESS_DERIVATIVE
{
    typedef typename TV::SCALAR T;
public:
    SYMMETRIC_MATRIX<T,TV::m> H; // H(i,k) = x_iikk
    typename TV::SPIN B,C; // B = x_ikik; C = x_ikki; order: (1D: none; 2D: 01; 3D: 12 20 01)

//#####################################################################
    void Enforce_Definiteness();
    void Compute_From_Singular_Value_Derivatives(const DIAGONAL_MATRIX<T,TV::m>& F,const TV& dE_ds,const SYMMETRIC_MATRIX<T,TV::m>& dE_dsds); // Not robust
    // div_diff=(dE/ds0-dE/ds1)/(s0-s1), div_sum=(dE/ds0+dE/ds1)/(s0+s1), same ordering as B and C.
    void Compute_From_Derivatives_And_Differences(const SYMMETRIC_MATRIX<T,TV::m>& hess,const typename TV::SPIN& div_diff,const typename TV::SPIN& div_sum); // Robust

    MATRIX<T,TV::m> Contract_03(const MATRIX<T,TV::m>& M) const; // solving Cjk = Mijkl * Mil
    MATRIX<T,TV::m> Contract_13(const MATRIX<T,TV::m>& M) const; // solving Cik = Mijkl * Mjl
    MATRIX<T,TV::m> Contract_23(const MATRIX<T,TV::m>& M) const; // solving Cij = Mijkl * Mkl

    MATRIX<T,TV::m> Differential(const MATRIX<T,TV::m>& M) const // solving Cij = Mijkl * Mkl
    {return Contract_23(M);}
//#####################################################################
};
// Compute Sigma, given SVD of F (U*F*Transpose(V))
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);
// Compute dSigma_dF:M, given SVD of F (U*F*Transpose(V))
template<class TV,class T=typename TV::SCALAR,int d=TV::m> 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);
}
#endif
