//#####################################################################
// Copyright 2013.
// This file is part of PhysBAM whose distribution is governed by the license contained in the accompanying file PHYSBAM_COPYRIGHT.txt.
//#####################################################################
#include <Core/Analytic_Tests/ANALYTIC_SYMMETRIC_MATRIX.h>
#include <Core/Diff_Tools/Diff_Test.h>
#include <Core/Random_Numbers/RANDOM_NUMBERS.h>
#include <Core/Symbolics/PROGRAM_COMPILER.h>

namespace PhysBAM{

//#####################################################################
// Destructor
//#####################################################################
template<class TV,int m> ANALYTIC_SYMMETRIC_MATRIX<TV,m>::
~ANALYTIC_SYMMETRIC_MATRIX()
{
}
//#####################################################################
// Function Test
//#####################################################################
template<class TV,int m> void ANALYTIC_SYMMETRIC_MATRIX<TV,m>::
Test(const TV& X,T t) const
{
    T eps=(T)1e-6,ds;
    RANDOM_NUMBERS<T> rand;
    TV dx;
    rand.Fill_Uniform(dx,-eps,eps);
    rand.Fill_Uniform(ds,-eps,eps);
    Diff_Test("analytic stress x",f(X,t),X,dX(X,t),f(X+dx,t),X+dx,dX(X+dx,t));
    Diff_Test("analytic stress t",f(X,t),t,dt(X,t),f(X,t+ds),t+ds,dt(X,t+ds));
}
//#####################################################################
// Constructor
//#####################################################################
template<class TV,int m> ANALYTIC_SYMMETRIC_MATRIX_PROGRAM<TV,m>::
ANALYTIC_SYMMETRIC_MATRIX_PROGRAM(const std::string& str,const std::string& out_name)
    :h(str,out_name,TV::m,m,m,false)
{
}
//#####################################################################
// Destructor
//#####################################################################
template<class TV,int m> ANALYTIC_SYMMETRIC_MATRIX_PROGRAM<TV,m>::
~ANALYTIC_SYMMETRIC_MATRIX_PROGRAM()
{
}
//#####################################################################
// Function f
//#####################################################################
template<class TV,int m> auto ANALYTIC_SYMMETRIC_MATRIX_PROGRAM<TV,m>::
f(const TV& X,T t) const -> T_MATRIX
{
    h.f(X,t);
    T_MATRIX M;
    for(int j=0;j<TV::m;j++)
        for(int i=j;i<TV::m;i++)
            M(i,j)=h.context.data_out(h.idx_f(m*j+i));
    return M;
}
//#####################################################################
// Function dt
//#####################################################################
template<class TV,int m> auto ANALYTIC_SYMMETRIC_MATRIX_PROGRAM<TV,m>::
dt(const TV& X,T t) const -> T_MATRIX
{
    h.dt(X,t);
    T_MATRIX M;
    for(int j=0;j<m;j++)
        for(int i=j;i<m;i++)
            M(i,j)=h.context_dt.data_out(h.idx_dt(m*j+i));
    return M;
}
//#####################################################################
// Function dX
//#####################################################################
template<class TV,int m> auto ANALYTIC_SYMMETRIC_MATRIX_PROGRAM<TV,m>::
dX(const TV& X,T t) const -> T_TENSOR
{
    h.dX(X,t);
    T_TENSOR M;
    for(int j=0;j<m;j++)
        for(int i=j;i<m;i++)
            for(int k=0;k<TV::m;k++)
                M(i,j,k)=h.context_dx.data_out(h.idx_dx(m*m*k+m*j+i));
    return M;
}
template struct ANALYTIC_SYMMETRIC_MATRIX<VECTOR<double,2>,2>;
template struct ANALYTIC_SYMMETRIC_MATRIX<VECTOR<double,3>,3>;
template struct ANALYTIC_SYMMETRIC_MATRIX<VECTOR<float,2>,2>;
template struct ANALYTIC_SYMMETRIC_MATRIX<VECTOR<float,3>,3>;
template struct ANALYTIC_SYMMETRIC_MATRIX_PROGRAM<VECTOR<double,2>,2>;
template struct ANALYTIC_SYMMETRIC_MATRIX_PROGRAM<VECTOR<double,3>,3>;
template struct ANALYTIC_SYMMETRIC_MATRIX_PROGRAM<VECTOR<float,2>,2>;
template struct ANALYTIC_SYMMETRIC_MATRIX_PROGRAM<VECTOR<float,3>,3>;
}
