//#####################################################################
// Copyright 2002-2009, Robert Bridson, Ronald Fedkiw, Geoffrey Irving, Sergey Koltakov, Craig Schroeder, Andrew Selle, Tamar Shinar, Eftychios Sifakis, Jonathan Su, Joseph Teran, Wen Zheng.
// This file is part of PhysBAM whose distribution is governed by the license contained in the accompanying file PHYSBAM_COPYRIGHT.txt.
//#####################################################################
// Class DEFORMABLE_WIND_DRAG_3D
//#####################################################################
#ifndef __DEFORMABLE_WIND_DRAG_3D__
#define __DEFORMABLE_WIND_DRAG_3D__

#include <Core/Basic_Geometry/BASIC_SIMPLEX_POLICY.h>
#include <Core/Grids/GRID.h>
#include <Core/Interpolation/LINEAR_INTERPOLATION.h>
#include <Deformables/Forces/DEFORMABLES_FORCES.h>
namespace PhysBAM{

template<class T>
class DEFORMABLE_WIND_DRAG_3D:public DEFORMABLES_FORCES<VECTOR<T,3> >
{
    typedef VECTOR<T,3> TV;
public:
    typedef DEFORMABLES_FORCES<TV> BASE;
    using BASE::particles;
    typedef typename BASE::FREQUENCY_DATA DEFORMABLE_FREQUENCY_DATA;
    typedef MESH_OBJECT<TV,TV::m-1> T_SIMPLICIAL_OBJECT;
    typedef typename BASIC_SIMPLEX_POLICY<TV,TV::m-1>::SIMPLEX T_SIMPLEX;

    TRIANGULATED_SURFACE<T>* triangulated_surface;
    bool use_constant_wind;
    T constant_wind_viscosity;
    TV constant_wind;
    bool use_spatially_varying_wind;
    T spatially_varying_wind_viscosity;
    GRID<TV> V_grid;
    ARRAY<TV,VECTOR<int,3> >* spatially_varying_wind;
    T wind_density;
    ARRAY<T,VECTOR<int,3> > *spatially_varying_wind_density,*spatially_varying_wind_pressure;
    T linear_normal_viscosity; // uses vertex normals
private:
    static const LINEAR_INTERPOLATION_NODE_CONST_EXTRAP<TV,T> interpolation;
    static const LINEAR_INTERPOLATION_NODE_CONST_EXTRAP<TV,TV> vector_interpolation;
    struct OPTIMIZATION{
        OPTIMIZATION()
            :one_third_area((T)0)
        {}

        T one_third_area;TV inward_normal,center,wind_velocity;
    };
    ARRAY<OPTIMIZATION> optimization;
    ARRAY<int> force_elements;
    ARRAY<int> force_particles;
    ARRAY<TV> vertex_normals; // TODO: use triangulated_surface.vertex_normals
public:

    template<class T_OBJECT> DEFORMABLE_WIND_DRAG_3D(T_OBJECT& object);
    DEFORMABLE_WIND_DRAG_3D(DEFORMABLE_PARTICLES<TV>& deformable_body_particles_input);
    virtual ~DEFORMABLE_WIND_DRAG_3D();

    void Use_Constant_Wind(const T viscosity_input,const TV& wind_input=TV())
    {use_constant_wind=true;constant_wind_viscosity=viscosity_input;constant_wind=wind_input;}

    void Use_Spatially_Varying_Wind(const T viscosity_input,const GRID<TV>& grid_input,ARRAY<TV,VECTOR<int,3> >& V_input)
    {use_spatially_varying_wind=true;spatially_varying_wind_viscosity=viscosity_input;
    V_grid=grid_input;spatially_varying_wind=&V_input;}

    void Set_Wind_Density(const T wind_density_input)
    {wind_density=wind_density_input;}

    void Set_Wind_Density(ARRAY<T,VECTOR<int,3> >& density_input)
    {spatially_varying_wind_density=&density_input;}

    void Set_Wind_Pressure(ARRAY<T,VECTOR<int,3> >& pressure_input) // only valid for volumetric objects
    {spatially_varying_wind_pressure=&pressure_input;}

    void Use_Linear_Normal_Viscosity(const T viscosity_input)
    {linear_normal_viscosity=viscosity_input;}

private:
    TV Spatially_Varying_Wind_Velocity(const TV& X) const
    {return vector_interpolation.Interpolate_Node(V_grid,*spatially_varying_wind,X);}

    T Spatially_Varying_Wind_Density(const TV& X) const
    {return interpolation.Interpolate_Node(V_grid,*spatially_varying_wind_density,X);}

    T Spatially_Varying_Wind_Pressure(const TV& X) const
    {return interpolation.Interpolate_Node(V_grid,*spatially_varying_wind_pressure,X);}

    void Add_Force_Derivatives(ARRAY_VIEW<const TV> V,ARRAY_VIEW<TV> F,const T time,bool transpose=false) const override
    {}

    void Enforce_Definiteness(const bool enforce_definiteness_input) override
    {}

    T CFL_Strain_Rate() const override
    {return FLT_MAX;}

    void Initialize_CFL(ARRAY_VIEW<DEFORMABLE_FREQUENCY_DATA> frequency) override
    {}

    VECTOR<T,3> Add_Forces_Helper(TV relative_velocity,int t) const;

public:

//#####################################################################
    void Update_Mpi(const ARRAY<bool>& particle_is_simulated) override;
    void Add_Dependencies(SEGMENT_MESH& dependency_mesh) const 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;
//#####################################################################
};
}
#endif
