//##################################################################### 
// Copyright 2021, Ritoban Roy-Chowdhury.
// This file is part of PhysBAM whose distribution is governed by the license contained in the accompanying file PHYSBAM_COPYRIGHT.txt.
//#####################################################################
#ifndef __DIV_FREE_BALSARA__
#define __DIV_FREE_BALSARA__

#include <Core/Arrays/FACE_ARRAYS.h>
#include <Core/Interpolation/INTERPOLATION_FACE.h>

namespace PhysBAM {

template<class TV,class T2,class T_FACE_ARRAY=ARRAY<T2,FACE_INDEX<TV::m> > >
class DIV_FREE_BALSARA:public INTERPOLATION_FACE<TV, T2, T_FACE_ARRAY>
{
    // Implements the interpolation Section III.1 from Balsara, D. S. (2001). Divergence-free adaptive mesh refinement for magnetohydrodynamics. Journal of Computational Physics, 175(2), 614-648.
    typedef typename TV::SCALAR T;typedef VECTOR<int, TV::m> IV;
    typedef VECTOR<T2,TV::m> TV2;typedef SYMMETRIC_MATRIX<T2,TV::m> SM2;
    typedef MATRIX<T2, TV::m> TM2;typedef DIAGONAL_MATRIX<T2,TV::m> DM2;
    typedef INTERPOLATION_FACE<TV, T2> BASE;

    constexpr static int stencil_size=3;
    constexpr static int const_extrap=false;

public:
    DIV_FREE_BALSARA() = default;

    void Get_Coeffs(const GRID<TV>& grid,const T_FACE_ARRAY& u,const IV& center, TV& T, TM2& E, TM2& S) const;

    TV2 Interpolate_Face(const GRID<TV>& grid,const T_FACE_ARRAY& u,const TV& X) const;

    T2 Interpolate_Face(const GRID<TV>& grid,const T_FACE_ARRAY& u,const TV& X,int axis) const ;

    VECTOR<VECTOR<T2,2>,TV::m> Extrema_Face(const GRID<TV>& grid,const T_FACE_ARRAY& u0,const T_FACE_ARRAY& u1,const TV& X) const;

    VECTOR<T2,2> Extrema_Face(const GRID<TV>& grid,const T_FACE_ARRAY& u0,const T_FACE_ARRAY& u1,const TV& X,int axis) const;

    VECTOR<TV2,TV::m> Gradient_Face(const GRID<TV>& grid,const T_FACE_ARRAY& u,const TV& X) const ;


    TV2 Gradient_Face(const GRID<TV>& grid,const T_FACE_ARRAY& u,const TV& X,int axis) const;

    VECTOR<SM2,TV::m> Hessian_Face(const GRID<TV>& grid,const T_FACE_ARRAY& u,const TV& X) const;

    SM2 Hessian_Face(const GRID<TV>& grid,const T_FACE_ARRAY& u,const TV& X,int axis) const;

protected:
    virtual TV2 Interpolate_Face_V(const GRID<TV>& grid,const T_FACE_ARRAY& u,const TV& X) const;
    virtual T2 Interpolate_Face_V(const GRID<TV>& grid,const T_FACE_ARRAY& u,const TV& X,int axis) const;
    virtual VECTOR<VECTOR<T2,2>,TV::m> Extrema_Face_V(const GRID<TV>& grid,const T_FACE_ARRAY& u0,const T_FACE_ARRAY& u1,const TV& X) const;
    virtual VECTOR<T2,2> Extrema_Face_V(const GRID<TV>& grid,const T_FACE_ARRAY& u0,const T_FACE_ARRAY& u1,const TV& X,int axis) const;
    virtual VECTOR<TV2,TV::m> Gradient_Face_V(const GRID<TV>& grid,const T_FACE_ARRAY& u,const TV& X) const;
    virtual TV2 Gradient_Face_V(const GRID<TV>& grid,const T_FACE_ARRAY& u,const TV& X,int axis) const;
    virtual VECTOR<SM2,TV::m> Hessian_Face_V(const GRID<TV>& grid,const T_FACE_ARRAY& u,const TV& X) const;
    virtual SM2 Hessian_Face_V(const GRID<TV>& grid,const T_FACE_ARRAY& u,const TV& X,int axis) const;
//#####################################################################
};

}

#endif
