//#####################################################################
// Copyright 2002-2005, Ronald Fedkiw, Geoffrey Irving, Frank Losasso, Andrew Selle.
// This file is part of PhysBAM whose distribution is governed by the license contained in the accompanying file PHYSBAM_COPYRIGHT.txt.
//#####################################################################
// Class BOUNDARY
//
// The ghost cell routines copy u into u_ghost and add a layer of three ghost cells around the boundary of u.
// The boundary condition routines modify u directly.
//
//#####################################################################
#ifndef __BOUNDARY__
#define __BOUNDARY__

#include <Core/Boundaries/BOUNDARY.h>
#include <Core/Arrays/ARRAYS_FORWARD.h>
#include <Core/Arrays/ARRAYS_ND.h>
#include <Core/Arrays/FACE_ARRAYS.h>
#include <Core/Grids/FACE_INDEX.h>
#include <Core/Log/DEBUG_UTILITIES.h>
#include <Core/Vectors/VECTOR.h>
#include <cassert>
namespace PhysBAM{

template<class TV> class GRID;
template<class TV,class T2>
class BOUNDARY
{
    typedef typename TV::SCALAR T;typedef VECTOR<bool,2> TV_BOOL2;typedef VECTOR<TV_BOOL2,TV::m> TV_SIDES;typedef VECTOR<int,TV::m> IV;
public:
    bool use_fixed_boundary,clamp_below,clamp_above;
    T2 fixed_boundary_value,lower_threshold,upper_threshold;
    VECTOR<VECTOR<bool,2>,TV::m> constant_extrapolation;

    BOUNDARY();
    BOUNDARY(const BOUNDARY&) = delete;
    void operator=(const BOUNDARY&) = delete;
    virtual ~BOUNDARY() = default;

    virtual void Set_Constant_Extrapolation(const TV_SIDES& constant_extrapolation_input=TV_SIDES::Constant_Vector(VECTOR<bool,2>::Constant_Vector(true)))
    {constant_extrapolation=constant_extrapolation_input;}

    void Turn_Off_Constant_Extrapolation()
    {Set_Constant_Extrapolation(TV_SIDES());}

    virtual bool Constant_Extrapolation(const int side) const
    {assert((unsigned)side<6);int axis=side/2;return constant_extrapolation(axis)(side&1);}

    virtual void Set_Fixed_Boundary(const bool use_fixed_boundary_input=true,const T2 fixed_boundary_value_input=T2())
    {use_fixed_boundary=use_fixed_boundary_input;fixed_boundary_value=fixed_boundary_value_input;
    if(use_fixed_boundary) clamp_above=clamp_below=false;}

    void Limit_Minimum_Boundary_Value(const bool clamp_below_input=true,const T2 lower_threshold_input=T2())
    {clamp_below=clamp_below_input;lower_threshold=lower_threshold_input;}

    void Limit_Maximum_Boundary_Value(const bool clamp_above_input=true,const T2 upper_threshold_input=T2())
    {clamp_above=clamp_above_input;upper_threshold=upper_threshold_input;}

    int Boundary(const int side,const RANGE<IV>& region) const
    {int axis=side/2;return side&1?region.Minimum_Corner()[axis]-1:region.Maximum_Corner()[axis];}

//#####################################################################
    virtual void Apply_Boundary_Condition_Node(const GRID<TV>& grid,ARRAYS_ND_BASE<T2,IV>& u) const;
    virtual void Apply_Boundary_Condition_Cell(const GRID<TV>& grid,ARRAYS_ND_BASE<T2,IV>& u) const;
    virtual void Apply_Boundary_Condition_Face(const GRID<TV>& grid,ARRAY<T2,FACE_INDEX<TV::m> >& u) const;
    virtual void Fill_Ghost_Nodes(const GRID<TV>& grid,const ARRAYS_ND_BASE<T2,IV>& u,ARRAYS_ND_BASE<T2,IV>& u_ghost,const int number_of_ghost_cells) const;
    virtual void Fill_Ghost_Cells(const GRID<TV>& grid,const ARRAYS_ND_BASE<T2,IV>& u,ARRAYS_ND_BASE<T2,IV>& u_ghost,const int number_of_ghost_cells) const;
    virtual void Fill_Ghost_Faces(const GRID<TV>& grid,const ARRAY<T2,FACE_INDEX<TV::m> >& u,ARRAY<T2,FACE_INDEX<TV::m> >& u_ghost,const int number_of_ghost_cells) const;
    virtual void Fill_Single_Ghost_Region_Node(const GRID<TV>& grid,ARRAYS_ND_BASE<T2,IV>& u_ghost,const RANGE<IV>& region,const int side,const int number_of_ghost_cells) const;
    virtual void Fill_Single_Ghost_Region_Cell(const GRID<TV>& grid,ARRAYS_ND_BASE<T2,IV>& u_ghost,const RANGE<IV>& region,const int side,const int number_of_ghost_cells) const;
    static void Find_Ghost_Regions_Node(const GRID<TV>& grid,VECTOR<RANGE<IV>,2*TV::m>& regions,const int number_of_ghost_cells);
    static void Find_Ghost_Regions_Cell(const GRID<TV>& grid,VECTOR<RANGE<IV>,2*TV::m>& regions,const int number_of_ghost_cells);
    void Fill_Single_Ghost_Region_Node(const GRID<TV>& grid,ARRAYS_ND_BASE<T2,IV>& u_ghost,const int side,const RANGE<IV>& region) const;
    void Fill_Single_Ghost_Region_Cell(const GRID<TV>& grid,ARRAYS_ND_BASE<T2,IV>& u_ghost,const int side,const RANGE<IV>& region) const;
//#####################################################################
};
}
#endif
