//#####################################################################
// 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_FACE
//#####################################################################
#ifndef __EXTRAPOLATION_HIGHER_ORDER_FACE__
#define __EXTRAPOLATION_HIGHER_ORDER_FACE__

#include <Core/Arrays/ARRAYS_ND.h>
#include <Core/Grids/GRID.h>
#include <Core/Level_Sets/EXTRAPOLATION_HIGHER_ORDER.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_FACE
{
    typedef VECTOR<int,TV::m> IV;
    typedef VECTOR<bool,TV::m> BV;
    typedef typename TV::SCALAR T;

    VECTOR<EXTRAPOLATION_HIGHER_ORDER<TV>*,TV::m> eho;
public:

    EXTRAPOLATION_HIGHER_ORDER_FACE(const LEVELSET<TV>& phi,const GRID<TV>& grid,
        std::function<bool(const FACE_INDEX<TV::m>& index)> inside_mask,
        const RANGE<IV>& array_domain,int order,int fill_width,
        const BV& periodic={})
    {
        for(int i=0;i<TV::m;i++)
        {
            auto face_domain=array_domain;
            face_domain.max_corner(i)++;
            eho(i)=new EXTRAPOLATION_HIGHER_ORDER<TV>(phi,
                grid.Get_Face_Grid(i),at_nodes,
                [&](const IV& index){return inside_mask(FACE_INDEX<TV::m>(i,index));},
                face_domain,order,fill_width,periodic,grid.numbers_of_cells);
        }
    }

    EXTRAPOLATION_HIGHER_ORDER_FACE(const EXTRAPOLATION_HIGHER_ORDER_FACE&) = delete;
    void operator=(const EXTRAPOLATION_HIGHER_ORDER_FACE&) = delete;
    ~EXTRAPOLATION_HIGHER_ORDER_FACE()
    {
        for(auto p:eho) delete p;
    }

    template<class T2>
    void Extrapolate(ARRAY<T2,FACE_INDEX<TV::m> >& x,int iterations)
    {
        for(int i=0;i<TV::m;i++)
            eho(i)->Extrapolate(x.Component(i),iterations);
    }

//     // void Extrapolate_Face(const GRID<TV>& grid,std::function<bool(const FACE_INDEX<TV::m>& index)> inside_mask,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);
//     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);
//     void Extrapolate_RK2(const ARRAY<VECTOR<STENCIL,TV::m> >& stencil,ARRAY<T2>& x,const ARRAY<T2>* z,ARRAY<T2>& tmp,int o,T dt);
//     void Fill_un(const GRID<TV>& grid,const TV& one_over_dx,const ARRAY<T2>& x,ARRAY<T2>& xn,int o,int mo);
//#####################################################################
};
}
#endif
