//#####################################################################
// 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/GENERAL_CONSTITUTIVE_MODEL_PROGRAM.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 <Core/Symbolics/PROGRAM_COMPILER.h>
#include <Deformables/Constitutive_Models/GENERAL_STRESS_DERIVATIVE.h>
using namespace PhysBAM;
//#####################################################################
// Constructor
//#####################################################################
template<class T,int d> GENERAL_CONSTITUTIVE_MODEL_PROGRAM<T,d>::
GENERAL_CONSTITUTIVE_MODEL_PROGRAM(const std::string& code,bool want_dP_dF)
{
    PROGRAM_COMPILER<T> pc;
    pc.Add_Input("F",TV::m,TV::m);
    pc.Add_Output("psi");
    pc.Parse(code);
    pc.Optimize();

    PROGRAM_COMPILER<T> pc_P(pc);
    ARRAY<MEM_IN> F(TV::m*TV::m);
    for(int i=0;i<TV::m*TV::m;i++) F(i)=MEM_IN(i);
    auto idx_P=pc_P.Make_Output(pc_P.Gradient(MEM_OUT(0)));
    pc_P.Optimize();

    if(want_dP_dF)
    {
        PROGRAM_COMPILER<T> pc_dP_dF(pc_P);
        for(int i=0;i<TV::m*TV::m;i++)
            pc_dP_dF.Make_Output(pc_dP_dF.Diff(idx_P,MEM_IN(i)));
        pc_dP_dF.Optimize();
        pc_dP_dF.Finalize(prog_dP_dF);
        context_dP_dF.Initialize(prog_dP_dF);
    }

    pc.Finalize(prog_psi);
    context_psi.Initialize(prog_psi);

    pc_P.Finalize(prog_P);
    context_P.Initialize(prog_P);
}
//#####################################################################
// Destructor
//#####################################################################
template<class T,int d> GENERAL_CONSTITUTIVE_MODEL_PROGRAM<T,d>::
~GENERAL_CONSTITUTIVE_MODEL_PROGRAM()
{
}
//#####################################################################
// Function Energy_Density
//#####################################################################
template<class T,int d> T GENERAL_CONSTITUTIVE_MODEL_PROGRAM<T,d>::
Energy_Density(const MATRIX<T,d>& F) const
{
    context_psi.Set(MEM_IN(0),F);
    prog_psi.Execute(context_psi);
    return context_psi.template Get<T>(MEM_OUT(0));
}
//#####################################################################
// Function Stress
//#####################################################################
template<class T,int d> MATRIX<T,d> GENERAL_CONSTITUTIVE_MODEL_PROGRAM<T,d>::
Stress(const MATRIX<T,d>& F) const
{
    context_P.Set(MEM_IN(0),F);
    prog_P.Execute(context_P);
    return context_P.template Get<MATRIX<T,d>>(MEM_OUT(1));
}
//#####################################################################
// Function Isotropic_Stress_Derivative
//#####################################################################
template<class T,int d> void GENERAL_CONSTITUTIVE_MODEL_PROGRAM<T,d>::
Stress_Derivative(const MATRIX<T,d>& F,GENERAL_STRESS_DERIVATIVE<TV>& dP_dF) const
{
    context_dP_dF.Set(MEM_IN(0),F);
    prog_dP_dF.Execute(context_dP_dF);
    for(int i=0;i<TV::m;i++)
        for(int j=0;j<TV::m;j++)
            for(int k=0;k<TV::m;k++)
                for(int l=0;l<TV::m;l++)
                    dP_dF.A(i,j)(k,l)=context_dP_dF.data_out(MEM_OUT(1+TV::m*TV::m+((i*TV::m+j)*TV::m+k)*TV::m+l));
}
namespace PhysBAM{
template class GENERAL_CONSTITUTIVE_MODEL_PROGRAM<float,1>;
template class GENERAL_CONSTITUTIVE_MODEL_PROGRAM<float,2>;
template class GENERAL_CONSTITUTIVE_MODEL_PROGRAM<float,3>;
template class GENERAL_CONSTITUTIVE_MODEL_PROGRAM<double,1>;
template class GENERAL_CONSTITUTIVE_MODEL_PROGRAM<double,2>;
template class GENERAL_CONSTITUTIVE_MODEL_PROGRAM<double,3>;
}
