//#####################################################################
// Copyright 2002-2007, Robert Bridson, Ronald Fedkiw, Geoffrey Irving, Neil Molino, 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 LINEAR_SPRINGS
//#####################################################################
#ifndef __LINEAR_SPRINGS__
#define __LINEAR_SPRINGS__

#include <Core/Vectors/VECTOR.h>
#include <Deformables/Forces/DEFORMABLES_FORCES.h>
namespace PhysBAM{

template<class TV>
class LINEAR_SPRINGS:public DEFORMABLES_FORCES<TV>,public SPRINGS_TAG
{
    typedef typename TV::SCALAR T;
public:
    typedef DEFORMABLES_FORCES<TV> BASE;
    using BASE::particles;using BASE::Invalidate_CFL;using BASE::cfl_number;
    using BASE::max_strain_per_time_step;
    using BASE::use_rest_state_for_strain_rate;
    typedef typename BASE::FREQUENCY_DATA FREQUENCY_DATA;

public:
    SEGMENT_MESH& segment_mesh;
    ARRAY<T> youngs_modulus; // units of force (i.e. force per unit strain)
    T constant_youngs_modulus; // units of force (i.e. force per unit strain)
    ARRAY<T> restlength,visual_restlength; // visual restlength corresponds to length between particles; restlength may be larger than this to avoid zero/small restlength
    ARRAY<T> damping; // units of force*time (i.e. force per unit strain rate)
    T constant_damping; // units of force*time (i.e. force per unit strain rate)
    bool use_plasticity;
    ARRAY<T> plastic_yield_strain;
    ARRAY<T> plastic_hardening;
    ARRAY<T> plastic_visual_restlength;
    T plasticity_clamp_ratio;
    bool cache_strain;
    mutable ARRAY<VECTOR<T,2> > strains_of_segment; // VECTOR<T,2>(strain_rate, strain)

    bool verbose;

protected:
    struct STATE{
        STATE()
            :coefficient(0)
        {}

        VECTOR<int,2> nodes; // copy of nodes to reduce associativity needed in cache
        T coefficient;
        TV direction;
        TV dX;
    };
    ARRAY<STATE> states;
    ARRAY<T> current_lengths;

public:
    ARRAY<int> force_segments;
    LINEAR_SPRINGS(DEFORMABLE_PARTICLES<TV>& particles,SEGMENT_MESH& segment_mesh_input);

    virtual ~LINEAR_SPRINGS();

//#####################################################################
    void Enforce_Definiteness(const bool enforce_definiteness_input) override;
    virtual void Set_Stiffness(const T youngs_modulus_input);
    virtual void Set_Stiffness(ARRAY_VIEW<const T> youngs_modulus_input);
    virtual void Set_Restlength(ARRAY_VIEW<const T> restlength_input);
    virtual void Set_Damping(const T damping_input);
    virtual void Set_Damping(ARRAY_VIEW<const T> damping_input);
    template<class T_FIELD> void Enable_Plasticity(const T_FIELD& plastic_yield_strain_input,const T_FIELD& plastic_hardening_input,const T plasticity_clamp_ratio_input=4);
    virtual void Set_Restlength_From_Particles();
    virtual void Set_Restlength_From_Material_Coordinates(ARRAY_VIEW<TV> material_coordinates);
    virtual void Clamp_Restlength(const T clamped_restlength);
    void Add_Dependencies(SEGMENT_MESH& dependency_mesh) const override;
    void Update_Mpi(const ARRAY<bool>& particle_is_simulated) override;
    virtual void Set_Overdamping_Fraction(const T overdamping_fraction=1); // 1 is critically damped
    virtual void Set_Overdamping_Fraction(ARRAY_VIEW<const T> overdamping_fraction); // 1 is critically damped
    void Clamp_Restlength_With_Fraction_Of_Springs(const T fraction=.01);
    void Enlarge_Restlength_By(const T enlargement_ratio);
    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 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;
    T Average_Restlength() const;
    void Print_Restlength_Statistics() const;
    void Print_Deformation_Statistics() const;
    T Maximum_Compression_Or_Expansion_Fraction(int* index=0) const;
    T Potential_Energy(const int s,const T time) const;
    T Potential_Energy(const T time) const override;
//#####################################################################
};

template<class TV> LINEAR_SPRINGS<TV>*
Create_Edge_Springs(DEFORMABLE_PARTICLES<TV>& particles,SEGMENT_MESH& segment_mesh,const typename TV::SCALAR stiffness=2e3,
    const typename TV::SCALAR overdamping_fraction=1,const bool limit_time_step_by_strain_rate=true,const typename TV::SCALAR max_strain_per_time_step=.1,
    const bool use_rest_state_for_strain_rate=true,const typename TV::SCALAR restlength_enlargement_fraction=0,const bool verbose=true,const typename TV::SCALAR rest_length_material_length_enlargement_ratio=0);

template<class T_OBJECT> LINEAR_SPRINGS<typename T_OBJECT::VECTOR_T>*
Create_Edge_Springs(T_OBJECT& object,
    const typename T_OBJECT::SCALAR stiffness=2e3,const typename T_OBJECT::SCALAR overdamping_fraction=1,const bool limit_time_step_by_strain_rate=true,
    const typename T_OBJECT::SCALAR max_strain_per_time_step=.1,const bool use_rest_state_for_strain_rate=true,const typename T_OBJECT::SCALAR restlength_enlargement_fraction=0,
    const bool verbose=true,const typename T_OBJECT::SCALAR rest_length_material_length_enlargement_ratio=0);

}
#endif
