//#####################################################################
// 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/SURFACE_TENSION_FORCE.h>
#include <Core/Matrices/MATRIX.h>
#include <Deformables/Particles/DEFORMABLE_PARTICLES.h>
using namespace PhysBAM;
//#####################################################################
// Constructor
//#####################################################################
template<class TV> SURFACE_TENSION_FORCE<TV>::
SURFACE_TENSION_FORCE(SEGMENTED_CURVE<TV>& surface_input,T surface_tension_coefficient_input)
    :BASE(dynamic_cast<DEFORMABLE_PARTICLES<TV>&>(surface_input.particles)),surface(surface_input),surface_tension_coefficient(surface_tension_coefficient_input),dt(0),apply_explicit_forces(true),apply_implicit_forces(true)
{
}
//#####################################################################
// Destructor
//#####################################################################
template<class TV> SURFACE_TENSION_FORCE<TV>::
~SURFACE_TENSION_FORCE()
{
}
//#####################################################################
// Function Add_Forces
//#####################################################################
template<class TV> void SURFACE_TENSION_FORCE<TV>::
Add_Forces(ARRAY_VIEW<TV> F,const T time) const
{
    if(apply_explicit_forces)
        for(int i=0;i<surface.mesh.elements.m;i++){VECTOR<int,2> k=surface.mesh.elements(i);
            TV f=(surface.particles.X(k.x)-surface.particles.X(k.y))*coefficients(i);
            F(k.x)-=f;
            F(k.y)+=f;}
}
//#####################################################################
// Function Tangential_Helper
//#####################################################################
template<class T> static void Tangential_Helper(MATRIX<T,3,2>& tangential,const VECTOR<T,3>& normal)
{
    tangential.Set_Column(1,normal.Unit_Orthogonal_Vector());
    tangential.Set_Column(2,normal.Cross(tangential.Column(1)));
}
template<class T> static void Tangential_Helper(MATRIX<T,2,1>& tangential,const VECTOR<T,2>& normal)
{
    tangential.Set_Column(1,normal.Orthogonal_Vector());
}
template<class T> static void Tangential_Helper(MATRIX<T,1,0>& tangential,const VECTOR<T,1>& normal)
{
}
//#####################################################################
// Function Update_Position_Based_State
//#####################################################################
template<class TV> void SURFACE_TENSION_FORCE<TV>::
Update_Position_Based_State(const T time,const bool is_position_update,const bool update_hessian)
{
    coefficients.Resize(surface.mesh.elements.m);
    normal.Resize(surface.mesh.elements.m);
    for(int i=0;i<surface.mesh.elements.m;i++){VECTOR<int,2> k=surface.mesh.elements(i);
        normal(i)=(surface.particles.X(k.x)-surface.particles.X(k.y)).Orthogonal_Vector();
        coefficients(i)=surface_tension_coefficient/normal(i).Normalize();}
}
//#####################################################################
// Function Add_Velocity_Dependent_Forces
//#####################################################################
template<class TV> void SURFACE_TENSION_FORCE<TV>::
Add_Velocity_Dependent_Forces(ARRAY_VIEW<const TV> V,ARRAY_VIEW<TV> F,const T time) const
{
}
//#####################################################################
// Function Add_Force_Derivatives
//#####################################################################
template<class TV> void SURFACE_TENSION_FORCE<TV>::
Add_Force_Derivatives(ARRAY_VIEW<const TV> V,ARRAY_VIEW<TV> F,const T time,bool transpose) const
{
    // this is different from Add_Velocity_Dependent_Forces because there is no dt in the formular
    for(int i=0;i<surface.mesh.elements.m;i++){VECTOR<int,2> k=surface.mesh.elements(i);
        TV f=(normal(i).Dot(V(k.x)-V(k.y))*coefficients(i))*normal(i);
        F(k.x)-=f;
        F(k.y)+=f;}
}
//#####################################################################
// Function Enforce_Definiteness
//#####################################################################
template<class TV> void SURFACE_TENSION_FORCE<TV>::
Enforce_Definiteness(const bool enforce_definiteness_input)
{
}
//#####################################################################
// Function Initialize_CFL
//#####################################################################
template<class TV> void SURFACE_TENSION_FORCE<TV>::
Initialize_CFL(ARRAY_VIEW<typename BASE::FREQUENCY_DATA> frequency)
{
}
//#####################################################################
// Function CFL_Strain_Rate
//#####################################################################
template<class TV> typename TV::SCALAR SURFACE_TENSION_FORCE<TV>::
CFL_Strain_Rate() const
{
    return FLT_MAX;
}
//#####################################################################
// Function Add_Dependencies
//#####################################################################
template<class TV> void SURFACE_TENSION_FORCE<TV>::
Add_Dependencies(SEGMENT_MESH& dependency_mesh) const
{
    surface.mesh.Add_Dependencies(dependency_mesh);
}
//#####################################################################
// Function Update_Mpi
//#####################################################################
template<class TV> void SURFACE_TENSION_FORCE<TV>::
Update_Mpi(const ARRAY<bool>& particle_is_simulated)
{
}
//#####################################################################
// Function Potential_Energy
//#####################################################################
template<class TV> typename TV::SCALAR SURFACE_TENSION_FORCE<TV>::
Potential_Energy(const T time) const
{
    T pe=0;
    if(apply_explicit_forces)
        for(int i=0;i<surface.mesh.elements.m;i++){VECTOR<int,2> k=surface.mesh.elements(i);
            pe+=surface_tension_coefficient*(surface.particles.X(k.x)-surface.particles.X(k.y)).Magnitude();}
    return pe;
}
//#####################################################################
// Function Dump_Curvatures
//#####################################################################
template<class TV> void SURFACE_TENSION_FORCE<TV>::
Dump_Curvatures() const
{
    T mn=FLT_MAX,mx=-mn,av=0;
    int n=0;
    ARRAY<TV> F(surface.mesh.elements.m);
    Add_Forces(F,0);
    ARRAY<T> A(coefficients.m);
    for(int i=0;i<surface.mesh.elements.m;i++){VECTOR<int,2> k=surface.mesh.elements(i);
        T len=(surface.particles.X(k.x)-surface.particles.X(k.y)).Magnitude()/2;
        A(k.x)+=len;
        A(k.y)+=len;}
    for(int i=0;i<F.m;i++){
        T K=abs(F(i).Magnitude()/A(i)/surface_tension_coefficient);
        TV dx=surface.particles.X(i)-(T).02;
        av+=K;
        n++;
        if(K>mx) mx=K;
        if(K<mn) mn=K;}
    if(n) LOG::cout<<"cstats "<<mn<<"   "<<mx<<std::endl;
    if(coefficients.m) LOG::cout<<"length estimates  "<<coefficients.Max()/coefficients.Min()<<std::endl;
}
namespace PhysBAM{
template class SURFACE_TENSION_FORCE<VECTOR<float,2> >;
template class SURFACE_TENSION_FORCE<VECTOR<double,2> >;
}
