//#####################################################################
// Copyright 2003-2007, Ron Fedkiw, Geoffrey Irving, Frank Losasso, Neil Molino, Igor Neverov, Craig Schroeder, Andrew Selle, Tamar Shinar, 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 FINITE_VOLUME
//#####################################################################
#ifndef __FINITE_VOLUME__
#define __FINITE_VOLUME__

#include <Core/Math_Tools/FACTORIAL.h>
#include <Core/Matrices/MATRIX.h>
#include <Core/Meshes/MESHES_FORWARD.h>
#include <Deformables/Forces/DEFORMABLES_FORCES.h>
namespace PhysBAM{

template<class TV,int d> class DIAGONALIZED_SEMI_IMPLICIT_ELEMENT;
template<class TV,int d> class STRAIN_MEASURE;
template<class T,int d> class CONSTITUTIVE_MODEL;
template<class T,int d> class PLASTICITY_MODEL;
template<class TV> class DIAGONALIZED_ISOTROPIC_STRESS_DERIVATIVE;
template<class TV> class IMPLICIT_OBJECT;
template<class T,int d> class DIAGONALIZED_STRESS_DERIVATIVE;
template<class T,int d> class ISOTROPIC_CONSTITUTIVE_MODEL;
template<class T,int d> class ANISOTROPIC_CONSTITUTIVE_MODEL;

template<class TV,int d>
class FINITE_VOLUME:public DEFORMABLES_FORCES<TV>,public FINITE_VOLUME_TAG
{
    typedef typename TV::SCALAR T;
    typedef MATRIX<T,TV::m,d> T_MATRIX;
    typedef MESH_OBJECT<TV,d> T_OBJECT;
public:
    typedef DEFORMABLES_FORCES<TV> BASE;
    using BASE::particles;using BASE::max_strain_per_time_step;
    using BASE::cfl_number;
    typedef typename BASE::FREQUENCY_DATA FREQUENCY_DATA;

    STRAIN_MEASURE<TV,d>& strain_measure;
    CONSTITUTIVE_MODEL<T,d>& constitutive_model;
    PLASTICITY_MODEL<T,d>* plasticity_model;
    ARRAY<T> Be_scales;
    ARRAY<T_MATRIX> U;
    ARRAY<MATRIX<T,d> > De_inverse_hat;
    ARRAY<DIAGONAL_MATRIX<T,d> > Fe_hat;
    ARRAY<DIAGONALIZED_ISOTROPIC_STRESS_DERIVATIVE<VECTOR<T,d> > >* dPi_dFe;
    ARRAY<DIAGONALIZED_STRESS_DERIVATIVE<T,d> >* dP_dFe;
    ARRAY<T>* Be_scales_save;
    ARRAY<MATRIX<T,d> >* V;
    T twice_max_strain_per_time_step; // for asynchronous
    ARRAY<DIAGONALIZED_SEMI_IMPLICIT_ELEMENT<TV,d> >* semi_implicit_data;
    bool use_uniform_density;
    ARRAY<T>* density_list;
    T density;
    ISOTROPIC_CONSTITUTIVE_MODEL<T,d>* isotropic_model;
    ANISOTROPIC_CONSTITUTIVE_MODEL<T,d>* anisotropic_model;
protected:
    IMPLICIT_OBJECT<TV>* implicit_surface;
    ARRAY<SYMMETRIC_MATRIX<T,TV::m> >* node_stiffness;
    ARRAY<MATRIX<T,TV::m> >* edge_stiffness;
public:
    ARRAY<int> force_elements;
    ARRAY<int>* force_segments;
    ARRAY<int> force_particles;
protected:
    bool destroy_data;
public:

    FINITE_VOLUME(const bool use_uniform_density_input,STRAIN_MEASURE<TV,d>& strain_measure,CONSTITUTIVE_MODEL<T,d>& constitutive_model,PLASTICITY_MODEL<T,d>* plasticity_model=0);
    virtual ~FINITE_VOLUME();

    void Define_Inversion_Via_Implicit_Surface(IMPLICIT_OBJECT<TV>& implicit_surface_input)
    {if(TV::m==d) PHYSBAM_FATAL_ERROR();implicit_surface=&implicit_surface_input;}

//#####################################################################
    static FINITE_VOLUME* Create(T_OBJECT& object,CONSTITUTIVE_MODEL<T,d>* constitutive_model,PLASTICITY_MODEL<T,d>* plasticity_model,const bool verbose,const bool use_uniform_density);
    // initialization
    void Save_V();
    void Save_Stress_Derivative();
    void Use_Quasistatics();
    void Use_Stiffness_Matrix();
    void Update_Be_Scales();
    void Initialize_Be_Scales_Save();
    void Copy_Be_Scales_Save_Into_Be_Scales(const ARRAY<int>& map);
    // solids_forces functions
    void Add_Dependencies(SEGMENT_MESH& dependency_mesh) const override;
    void Update_Mpi(const ARRAY<bool>& particle_is_simulated) override;
    void Update_Position_Based_State(const T time,const bool is_position_update,const bool update_hessian) override;
    void Add_Forces(ARRAY_VIEW<TV> F,const T time) const override;
    void Add_Velocity_Dependent_Forces(ARRAY_VIEW<const TV> V,ARRAY_VIEW<TV> F,const T time) const override;
    void Enforce_Definiteness(const bool enforce_definiteness_input) override;
    void Add_Force_Derivatives(ARRAY_VIEW<const TV> V,ARRAY_VIEW<TV> F,const T time,bool transpose=false) const override;
    void Initialize_CFL(ARRAY_VIEW<FREQUENCY_DATA> frequency) override;
    T CFL_Strain_Rate() const override;
    // asynchronous
    T Potential_Energy(const T time) const override;
    T Compute_Sound_Speed() const;
//#####################################################################
};

template<class T_OBJECT,class T,int d> FINITE_VOLUME<typename T_OBJECT::VECTOR_T,d>*
Create_Finite_Volume(T_OBJECT& object,CONSTITUTIVE_MODEL<T,d>* constitutive_model,const bool limit_time_step_by_strain_rate=true,const T max_strain_per_time_step=.1,
    const bool use_rest_state_for_strain_rate=true,const bool precompute_stiffness_matrix=false,const bool verbose=true,PLASTICITY_MODEL<T,d>* plasticity_model=0,
    const bool use_uniform_density=false);

template<class T_OBJECT,class T,int d> FINITE_VOLUME<typename T_OBJECT::VECTOR_T,d>*
Create_Quasistatic_Finite_Volume(T_OBJECT& object,CONSTITUTIVE_MODEL<T,d>* constitutive_model,const bool precompute_stiffness_matrix=false,const bool verbose=true,
    const bool use_uniform_density=false);

}
#endif
