//#####################################################################
// Copyright 2010.
// This file is part of PhysBAM whose distribution is governed by the license contained in the accompanying file PHYSBAM_COPYRIGHT.txt.
//#####################################################################
#include <Deformables/Forces/LINEAR_POINT_ATTRACTION.h>
#include <Core/Data_Structures/HASHTABLE.h>
#include <Core/Matrices/MATRIX.h>
#include <Core/Meshes/TRIANGULATED_AREA.h>
#include <Deformables/Particles/DEFORMABLE_PARTICLES.h>
using namespace PhysBAM;
//#####################################################################
// Constructor
//#####################################################################
template<class TV> LINEAR_POINT_ATTRACTION<TV>::
LINEAR_POINT_ATTRACTION(T_MESH& mesh,const TV& pt,T coefficient_input)
    :BASE(dynamic_cast<DEFORMABLE_PARTICLES<TV>&>(mesh.particles)),surface(mesh),coefficient(coefficient_input),point(pt),dt(0),apply_explicit_forces(true),apply_implicit_forces(true)
{
    referenced_particles=Get_Unique(mesh.mesh.elements.Flattened());
}
//#####################################################################
// Destructor
//#####################################################################
template<class TV> LINEAR_POINT_ATTRACTION<TV>::
~LINEAR_POINT_ATTRACTION()
{
}
//#####################################################################
// Function Add_Forces
//#####################################################################
template<class TV> void LINEAR_POINT_ATTRACTION<TV>::
Add_Forces(ARRAY_VIEW<TV> F,const T time) const
{
    if(apply_explicit_forces)
        for(int i=0;i<referenced_particles.m;i++){int p=referenced_particles(i);
            F(p)+=coefficient*(point-surface.particles.X(p));}
}
//#####################################################################
// Function Update_Position_Based_State
//#####################################################################
template<class TV> void LINEAR_POINT_ATTRACTION<TV>::
Update_Position_Based_State(const T time,const bool is_position_update,const bool update_hessian)
{
}
//#####################################################################
// Function Add_Velocity_Dependent_Forces
//#####################################################################
template<class TV> void LINEAR_POINT_ATTRACTION<TV>::
Add_Velocity_Dependent_Forces(ARRAY_VIEW<const TV> V,ARRAY_VIEW<TV> F,const T time) const
{
    if(apply_implicit_forces)
        for(int i=0;i<referenced_particles.m;i++){int p=referenced_particles(i);
            F(p)-=dt*coefficient*V(p);}
}
//#####################################################################
// Function Add_Force_Derivatives
//#####################################################################
template<class TV> void LINEAR_POINT_ATTRACTION<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 LINEAR_POINT_ATTRACTION<TV>::
Enforce_Definiteness(const bool enforce_definiteness_input)
{
}
//#####################################################################
// Function Initialize_CFL
//#####################################################################
template<class TV> void LINEAR_POINT_ATTRACTION<TV>::
Initialize_CFL(ARRAY_VIEW<typename BASE::FREQUENCY_DATA> frequency)
{
}
//#####################################################################
// Function CFL_Strain_Rate
//#####################################################################
template<class TV> typename TV::SCALAR LINEAR_POINT_ATTRACTION<TV>::
CFL_Strain_Rate() const
{
    return FLT_MAX;
}
//#####################################################################
// Function Add_Dependencies
//#####################################################################
template<class TV> void LINEAR_POINT_ATTRACTION<TV>::
Add_Dependencies(SEGMENT_MESH& dependency_mesh) const
{
    surface.mesh.Add_Dependencies(dependency_mesh);
}
//#####################################################################
// Function Update_Mpi
//#####################################################################
template<class TV> void LINEAR_POINT_ATTRACTION<TV>::
Update_Mpi(const ARRAY<bool>& particle_is_simulated)
{
}
//#####################################################################
// Function Potential_Energy
//#####################################################################
template<class TV> typename TV::SCALAR LINEAR_POINT_ATTRACTION<TV>::
Potential_Energy(const T time) const
{
    T pe=0;
    if(apply_explicit_forces)
        for(int i=0;i<referenced_particles.m;i++){int p=referenced_particles(i);
            pe+=coefficient/2*(surface.particles.X(p)-point).Magnitude_Squared();}
    return pe;
}
namespace PhysBAM{
template class LINEAR_POINT_ATTRACTION<VECTOR<float,2> >;
template class LINEAR_POINT_ATTRACTION<VECTOR<double,2> >;
}
