//#####################################################################
// Copyright 2010.
// This file is part of PhysBAM whose distribution is governed by the license contained in the accompanying file PHYSBAM_COPYRIGHT.txt.
//#####################################################################
// Class EXTRAPOLATION_HIGHER_ORDER
//#####################################################################
#ifndef __EXTRAPOLATION_HIGHER_ORDER__
#define __EXTRAPOLATION_HIGHER_ORDER__

#include <Core/Arrays/ARRAYS_ND.h>
#include <Core/Grids/GRID.h>
#include <Core/Level_Sets/LEVELSET.h>
#include <Core/Math_Tools/INTERVAL.h>
#include <Core/Vectors/VECTOR.h>
#include <functional>
namespace PhysBAM{

template<class TV>
class EXTRAPOLATION_HIGHER_ORDER
{
    typedef VECTOR<int,TV::m> IV;
    typedef VECTOR<bool,TV::m> BV;
    typedef typename TV::SCALAR T;

    struct STENCIL {VECTOR<int,4> nodes;T scale;};
    const LEVELSET<TV>& phi;
    GRID<TV> grid;
    int order;
    int fill_width;
    ARRAY<int,IV> node_to_index;
    ARRAY<IV> index_to_node;
    VECTOR<INTERVAL<int>,3> fill_indices,solve_indices;
    ARRAY<TV> normal;
    ARRAY<VECTOR<STENCIL,TV::m> > stencil;
    BV periodic;
    IV periodic_size;
public:

    EXTRAPOLATION_HIGHER_ORDER(const LEVELSET<TV>& phi,const GRID<TV>& grid,
        AT_NODES,std::function<bool(const IV& index)> inside_mask,
        const RANGE<IV>& array_domain,int order,int fill_width,
        const BV& periodic={});

    EXTRAPOLATION_HIGHER_ORDER(const LEVELSET<TV>& phi,const GRID<TV>& grid,
        AT_CELLS,std::function<bool(const IV& index)> inside_mask,
        const RANGE<IV>& array_domain,int order,int fill_width,
        const BV& periodic={});

    EXTRAPOLATION_HIGHER_ORDER(const LEVELSET<TV>& phi,const GRID<TV>& grid,
        AT_NODES,std::function<bool(const IV& index)> inside_mask,
        const RANGE<IV>& array_domain,int order,int fill_width,
        const BV& periodic,const IV& periodic_size);

    EXTRAPOLATION_HIGHER_ORDER(const LEVELSET<TV>& phi,const GRID<TV>& grid,
        AT_CELLS,std::function<bool(const IV& index)> inside_mask,
        const RANGE<IV>& array_domain,int order,int fill_width,
        const BV& periodic,const IV& periodic_size);

    EXTRAPOLATION_HIGHER_ORDER(const EXTRAPOLATION_HIGHER_ORDER&) = delete;
    void operator=(const EXTRAPOLATION_HIGHER_ORDER&) = delete;
    ~EXTRAPOLATION_HIGHER_ORDER() = default;

    template<class T2>
    void Extrapolate(ARRAYS_ND_BASE<T2,IV>& x,int iterations) const;

    // void Extrapolate_Face(const GRID<TV>& grid,std::function<bool(const FACE_INDEX<TV::m>& index)> inside_mask,ARRAY<T2,FACE_INDEX<TV::m> >& x,int iterations);

private:
    int Lookup_Index(const GRID<TV>& grid,IV index) const;
    int Register_Index(const GRID<TV>& grid,IV index,int only_neg=0);
    void Periodic_Index(const GRID<TV>& grid,IV& index) const;
    void Add_Neighbors(const GRID<TV>& grid,ARRAY<IV>& next,const ARRAY<IV>& neighbors,const IV& index,int unregistered,int registered);
    void Register_Nodes(const GRID<TV>& grid,const RANGE<IV>& array_domain,std::function<bool(const IV& index)> inside_mask);
    template<class T2>
    void Extrapolate_FE(const ARRAY<VECTOR<STENCIL,TV::m> >& stencil,const ARRAY<T2>& x,ARRAY<T2>& y,const ARRAY<T2>* z,int o,T dt,T alpha) const;
    template<class T2>
    void Extrapolate_RK2(const ARRAY<VECTOR<STENCIL,TV::m> >& stencil,ARRAY<T2>& x,const ARRAY<T2>* z,ARRAY<T2>& tmp,int o,T dt) const;
    template<class T2>
    void Fill_un(const GRID<TV>& grid,const TV& one_over_dx,const ARRAY<T2>& x,ARRAY<T2>& xn,int o,int mo) const;
//#####################################################################
};
}
#endif
