//#####################################################################
// Copyright 2007-2008, Michael Lentine, Craig Schroeder, Tamar Shinar.
// This file is part of PhysBAM whose distribution is governed by the license contained in the accompanying file PHYSBAM_COPYRIGHT.txt.
//#####################################################################
#include <Deformables/Forces/DEFORMABLES_FORCES.h>
#include <Core/Arrays/ARRAY.h>
#include <Core/Data_Structures/HASHTABLE.h>
#include <Core/Diff_Tools/Diff_Test.h>
#include <Core/Log/LOG.h>
#include <Core/Random_Numbers/RANDOM_NUMBERS.h>
#include <Core/Read_Write/READ_WRITE.h>
#include <Core/Vectors/VECTOR.h>
#include <Deformables/Particles/DEFORMABLE_PARTICLES.h>
namespace PhysBAM{
//#####################################################################
// Constructor
//#####################################################################
template<class TV> DEFORMABLES_FORCES<TV>::
DEFORMABLES_FORCES(DEFORMABLE_PARTICLES<TV>& particles)
    :particles(particles),cfl_number((T)1),allow_external_cfl_number(true),cfl_initialized(false),
    unique_id(Get_Unique_Id())
{
    Use_Rest_State_For_Strain_Rate(false);
    Limit_Time_Step_By_Strain_Rate();
}
//#####################################################################
// Destructor
//#####################################################################
template<class TV> DEFORMABLES_FORCES<TV>::
~DEFORMABLES_FORCES()
{}
//#####################################################################
// Function Use_Rest_State_For_Strain_Rate
//#####################################################################
template<class TV> void DEFORMABLES_FORCES<TV>::
Use_Rest_State_For_Strain_Rate(const bool use_rest_state_for_strain_rate_input)
{
    use_rest_state_for_strain_rate=use_rest_state_for_strain_rate_input;
}
//#####################################################################
// Function Limit_Time_Step_By_Strain_Rate
//#####################################################################
template<class TV> void DEFORMABLES_FORCES<TV>::
Limit_Time_Step_By_Strain_Rate(const bool limit_time_step_by_strain_rate_input,const T max_strain_per_time_step_input)
{
    limit_time_step_by_strain_rate=limit_time_step_by_strain_rate_input;
    assert(max_strain_per_time_step_input>0);max_strain_per_time_step=max_strain_per_time_step_input;
}
//#####################################################################
// Function Add_Force_Derivatives
//#####################################################################
template<class TV> void DEFORMABLES_FORCES<TV>::
Add_Force_Derivatives(ARRAY_VIEW<const TV> V,ARRAY_VIEW<TV> F,const T time,bool transpose) const
{
}
//#####################################################################
// Function Enforce_Definiteness
//#####################################################################
template<class TV> void DEFORMABLES_FORCES<TV>::
Enforce_Definiteness(const bool enforce_definiteness_input)
{
    PHYSBAM_FUNCTION_IS_NOT_DEFINED();
}
//#####################################################################
// Function Potential_Energy
//#####################################################################
template<class TV> typename TV::SCALAR DEFORMABLES_FORCES<TV>::
Potential_Energy(const T time) const
{
    PHYSBAM_WARN_IF_NOT_OVERRIDDEN();
    return 0;
}
//#####################################################################
// Function Update_Position_Based_State
//#####################################################################
template<class TV> void DEFORMABLES_FORCES<TV>::
Update_Position_Based_State(const T time,const bool is_position_update,const bool update_hessian)
{
}
//#####################################################################
// Function Test_Diff
//#####################################################################
template<class TV> void DEFORMABLES_FORCES<TV>::
Test_Diff(const T time) 
{
    PHYSBAM_ASSERT(sizeof(T)==sizeof(double));
    RANDOM_NUMBERS<T> random;
    T e=(T)1e-5;
    ARRAY<TV> dX(particles.X.m);
    random.Fill_Uniform(dX,-e,e);
    ARRAY<TV> X2a(particles.X+dX);
    ARRAY_VIEW<TV> X1(X2a);
    ARRAY<TV> F0(particles.X.m);
    ARRAY<TV> F1(particles.X.m);
    ARRAY<TV> G0(particles.X.m);
    ARRAY<TV> G1(particles.X.m);

    Update_Position_Based_State(time,true,false);
    T PE0=Potential_Energy(time);
    Add_Forces(F0,time);
    Add_Force_Derivatives(dX,G0,time);
    particles.X.Exchange(X1);

    Update_Position_Based_State(time,true,false);
    T PE1=Potential_Energy(time);
    Add_Forces(F1,time);
    Add_Force_Derivatives(dX,G1,time);
    particles.X.Exchange(X1);

    Update_Position_Based_State(time,true,false);

    LOG::printf("testing force %P\n",typeid(*this).name());
    Diff_Test("potential energy",-PE0,particles.X,F0,-PE1,X1,F1);
    Diff_Times_Test("force derivatives",F0,G0,particles.X,F1,G1,X1);
}
template class DEFORMABLES_FORCES<VECTOR<double,1> >;
template class DEFORMABLES_FORCES<VECTOR<double,2> >;
template class DEFORMABLES_FORCES<VECTOR<double,3> >;
template class DEFORMABLES_FORCES<VECTOR<float,1> >;
template class DEFORMABLES_FORCES<VECTOR<float,2> >;
template class DEFORMABLES_FORCES<VECTOR<float,3> >;
}
