//#####################################################################
// Copyright 2002-2010, Ronald Fedkiw, Eran Guendelman, Michael Lentine, Frank Losasso, Avi Robinson-Mosher, Andrew Selle, Tamar Shinar, Jerry Talton.
// This file is part of PhysBAM whose distribution is governed by the license contained in the accompanying file PHYSBAM_COPYRIGHT.txt.
//#####################################################################
// Class PROJECTION_UNIFORM
//#####################################################################
#ifndef __PROJECTION_UNIFORM__
#define __PROJECTION_UNIFORM__

#include <Core/Interpolation/FACE_LOOKUP_UNIFORM.h>
#include <Core/Poisson/POISSON_UNIFORM.h>
#include <Core/Poisson/PROJECTION.h>
namespace PhysBAM{

template<class TV>
class PROJECTION_UNIFORM:public PROJECTION<typename TV::SCALAR>
{
    typedef typename TV::SCALAR T;typedef VECTOR<int,TV::m> IV;
    typedef ARRAYS_ND_BASE<T,IV> T_ARRAYS_BASE;
public:
    typedef PROJECTION<T> BASE;
    using PROJECTION<T>::use_non_zero_divergence;

    GRID<TV> p_grid; // p_grid is a cell centered MAC grid
    ARRAY<T,IV> p;
    ARRAY<T,IV> p_save_for_projection;
    ARRAY<T,FACE_INDEX<TV::m> > face_velocities_save_for_projection;
    LAPLACE<TV>* elliptic_solver;
    LAPLACE<TV>* laplace;
    POISSON_UNIFORM<TV>* poisson;
    ARRAY<T,IV> divergence; // use this to set up a non-zero divergence
    bool use_divergence_multiplier;
    ARRAY<T,IV> divergence_multiplier;
    VECTOR<bool,TV::m> periodic;

    PROJECTION_UNIFORM(const GRID<TV>& mac_grid,const bool use_variable_beta=false,const bool use_poisson=false);
protected:
    PROJECTION_UNIFORM();
public:
    virtual ~PROJECTION_UNIFORM();

    void Use_Non_Zero_Divergence(const bool use_non_zero_divergence_input=true)
    {use_non_zero_divergence=use_non_zero_divergence_input;
    if(use_non_zero_divergence) divergence.Resize(p_grid.Cell_Indices());else divergence.Clean_Memory();}

    void Use_Divergence_Multiplier(const bool use_divergence_multiplier_input=true)
    {use_divergence_multiplier=use_divergence_multiplier_input;
    if(use_divergence_multiplier) divergence_multiplier.Resize(p_grid.Cell_Indices(3));else divergence_multiplier.Clean_Memory();}

    void Set_Periodic_Boundary(VECTOR<bool,TV::m> periodic_axes)
    {periodic=periodic_axes;elliptic_solver->periodic_boundary=periodic;}

//#####################################################################
    virtual void Initialize_Grid(const GRID<TV>& mac_grid);
    virtual void Make_Divergence_Free(ARRAY<T,FACE_INDEX<TV::m> >& face_velocities,const T dt,const T time);
    void Zero_Out_Neumann_Pocket_Velocities(ARRAY<T,FACE_INDEX<TV::m> >& face_velocities);
    virtual void Apply_Pressure(ARRAY<T,FACE_INDEX<TV::m> >& face_velocities,const T dt,const T time,bool scale_by_dt=false);
    void Enforce_Velocity_Compatibility(ARRAY<T,FACE_INDEX<TV::m> >& face_velocities);
    void Set_Up_For_Projection(ARRAY<T,FACE_INDEX<TV::m> >& face_velocities);
    void Restore_After_Projection(ARRAY<T,FACE_INDEX<TV::m> >& face_velocities);
    void Exchange_Pressures_For_Projection();
    void Compute_Divergence(const ARRAY<T,FACE_INDEX<TV::m> >& face_lookup,LAPLACE<TV>* solver);
//#####################################################################
};
}
#endif
