//#####################################################################
// Copyright 2002-2006, Ronald Fedkiw, Geoffrey Irving, Frank Losasso, Andrew Selle, Tamar Shinar.
// This file is part of PhysBAM whose distribution is governed by the license contained in the accompanying file PHYSBAM_COPYRIGHT.txt.
//#####################################################################
// Class ADVECTION
//#####################################################################
//
// Solves Z_t + u Z_x + v Z_y + w Z_z = 0 for one Euler step. Conservative schemes solve Z_t +(uZ)_x + (vZ)_y + (wZ)_z = 0.
// Z is updated using u, v, w, & Z_ghost. Input Z as (1,m) by (1,n) by (1,mn), and Z_ghost as (-2,m+3) by (-2,n+3) by (-2,mn+3).
// Input u, v, & w as (1,m) by (1,n) by (1,mn) for nonconservative schemes. Input u, v, & w as (-2,m+3) by (-2,n+3) by (-2,mn+3) for conservative schemes.
//
//#####################################################################
#ifndef __ADVECTION__
#define __ADVECTION__

#include <Core/Advection/ADVECTION_FORWARD.h>
#include <Core/Grids/FACE_INDEX.h>
#include <Core/Log/DEBUG_UTILITIES.h>
namespace PhysBAM{

template<class TV> struct GRID_ARRAYS_POLICY;
template<class T,class ID> class ARRAY;
template<class TV> class GRID;

template<class TV,class T2,class T_FACE_LOOKUP> // T_FACE_LOOKUP=typename GRID<TV>::FACE_LOOKUP
class ADVECTION
{
    typedef typename TV::SCALAR T;typedef VECTOR<int,TV::m> IV;
public:

    ADVECTION()
    {}

    ADVECTION(const ADVECTION&) = delete;
    void operator=(const ADVECTION&) = delete;

    virtual ~ADVECTION()
    {}

    void Update_Advection_Equation_Cell(const GRID<TV>& grid,ARRAY<T2,IV>& Z,const ARRAY<T2,IV>& Z_ghost,
        const ARRAY<T,FACE_INDEX<TV::m> >& face_velocities,const T dt,const T time,
        const ARRAY<T2,IV>* Z_min_ghost=0,const ARRAY<T2,IV>* Z_max_ghost=0,ARRAY<T2,IV>* Z_min=0,ARRAY<T2,IV>* Z_max=0,
        ARRAY<T2,FACE_INDEX<TV::m> >* flux=0)
    {Update_Advection_Equation_Cell_Lookup(grid,Z,Z_ghost,face_velocities,dt,time,Z_min_ghost,Z_max_ghost,Z_min,Z_max,flux);}

    void Update_Advection_Equation_Face(const GRID<TV>& grid,ARRAY<T,FACE_INDEX<TV::m> >& Z,const ARRAY<T,FACE_INDEX<TV::m> >& Z_ghost,
        const ARRAY<T,FACE_INDEX<TV::m> >& face_velocities,const T dt,const T time,
        const ARRAY<T,FACE_INDEX<TV::m> >* Z_min_ghost=0,const ARRAY<T,FACE_INDEX<TV::m> >* Z_max_ghost=0,ARRAY<T,FACE_INDEX<TV::m> >* Z_min=0,ARRAY<T,FACE_INDEX<TV::m> >* Z_max=0,
        VECTOR<ARRAY<T2,FACE_INDEX<TV::m> >,TV::m>* flux=0)
    {if(Z_min_ghost && Z_max_ghost){
        T_FACE_LOOKUP Z_min_lookup(*Z_min_ghost),Z_max_lookup(*Z_max_ghost);
        Update_Advection_Equation_Face_Lookup(grid,Z,Z_ghost,face_velocities,dt,time,&Z_min_lookup,&Z_max_lookup,Z_min,Z_max,flux);}
        else Update_Advection_Equation_Face_Lookup(grid,Z,Z_ghost,face_velocities,dt,time,0,0,0,0,flux);}

    void Update_Advection_Equation_Cell(const GRID<TV>& grid,ARRAY<T2,IV>& Z,const ARRAY<T2,IV>& Z_ghost,
        const ARRAY<TV,IV>& V,const T dt,const T time,
        const ARRAY<T2,IV>* Z_min_ghost=0,const ARRAY<T2,IV>* Z_max_ghost=0,ARRAY<T2,IV>* Z_min=0,ARRAY<T2,IV>* Z_max=0,
        ARRAY<T2,FACE_INDEX<TV::m> >* flux=0)
    {Update_Advection_Equation_Node(grid.Get_Regular_Grid_At_MAC_Positions(),Z,Z_ghost,V,dt,time,Z_min_ghost,Z_max_ghost,Z_min,Z_max,flux);}

//#####################################################################
    virtual void Update_Advection_Equation_Node(const GRID<TV>& grid,ARRAY<T2,IV>& Z,const ARRAY<T2,IV>& Z_ghost,
        const ARRAY<TV,IV>& V,const T dt,const T time,
        const ARRAY<T2,IV>* Z_min_ghost=0,const ARRAY<T2,IV>* Z_max_ghost=0,ARRAY<T2,IV>* Z_min=0,ARRAY<T2,IV>* Z_max=0,
        ARRAY<T2,FACE_INDEX<TV::m> >* flux=0)
    {PHYSBAM_FUNCTION_IS_NOT_DEFINED();}
    virtual void Update_Advection_Equation_Cell_Lookup(const GRID<TV>& grid,ARRAY<T2,IV>& Z,const ARRAY<T2,IV>& Z_ghost,
        const T_FACE_LOOKUP& face_velocities,const T dt,const T time,
        const ARRAY<T2,IV>* Z_min_ghost=0,const ARRAY<T2,IV>* Z_max_ghost=0,ARRAY<T2,IV>* Z_min=0,ARRAY<T2,IV>* Z_max=0,
        ARRAY<T2,FACE_INDEX<TV::m> >* flux=0)
    {PHYSBAM_FUNCTION_IS_NOT_DEFINED();}
    virtual void Update_Advection_Equation_Face_Lookup(const GRID<TV>& grid,ARRAY<T,FACE_INDEX<TV::m> >& Z,const T_FACE_LOOKUP& Z_ghost,
        const T_FACE_LOOKUP& face_velocities,const T dt,const T time,
        const T_FACE_LOOKUP* Z_min_ghost=0,const T_FACE_LOOKUP* Z_max_ghost=0,ARRAY<T,FACE_INDEX<TV::m> >* Z_min=0,ARRAY<T,FACE_INDEX<TV::m> >* Z_max=0,
        VECTOR<ARRAY<T2,FACE_INDEX<TV::m> >,TV::m>* flux=0)
    {PHYSBAM_FUNCTION_IS_NOT_DEFINED();}
//#####################################################################
};
}
#endif
