//#####################################################################
// Copyright 2005-2007, Geoffrey Irving, 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.
//#####################################################################
#include <Core/Projection/BOUNDARY_CONDITION_DOUBLE_FINE.h>
#include <Core/Arrays/FACE_ARRAYS.h>
#include <Core/Data_Structures/HASHTABLE.h>
#include <Core/Data_Structures/TUPLE.h>
#include <Core/Geometry_Particles/DEBUG_PARTICLES.h>
#include <Core/Grids/FACE_ITERATOR.h>
#include <Core/Implicit_Objects/IMPLICIT_OBJECT.h>
#include <Core/Intersections/GRID_SURFACE_INTERSECTION.h>
#include <Core/Log/DEBUG_SUBSTEPS.h>
#include <Core/Math_Tools/int_div.h>
#include <Core/Meshes/POINT_SIMPLICES_1D.h>
#include <Core/Meshes/SEGMENTED_CURVE.h>
#include <Core/Meshes/TETRAHEDRALIZED_VOLUME.h>
#include <Core/Meshes/TRIANGULATED_AREA.h>
namespace PhysBAM{

// fine_node_index = (2*face.index+1).Add_Axis(face.axis,-1);

// fine_node_index = 2*mac_cell+1;

//#####################################################################
// Constructor
//#####################################################################
template<class TV> BOUNDARY_CONDITION_DOUBLE_FINE<TV>::
BOUNDARY_CONDITION_DOUBLE_FINE(const GRID<TV>& mac_grid,int ghost)
    :mac_grid(mac_grid),grid(mac_grid.numbers_of_cells*2,mac_grid.domain),
    ghost(ghost),bc_type(grid.Node_Indices(2*ghost)),bc_type_current(grid.Node_Indices(2*ghost))
{
}
//#####################################################################
// Function Reset
//#####################################################################
template<class TV> void BOUNDARY_CONDITION_DOUBLE_FINE<TV>::
Reset()
{
    bc_type.Fill(0);
    bc_type_current.Fill(0);
    bc_p.Remove_All();
    bc_u.Remove_All();
}
template<class IV>
inline void set_bc(ARRAY<char,IV>& bc_type,ARRAY<char,IV>& bc_type_current,const IV& i,char type)
{
    enum bc_enum {bc_slip=-3,bc_noslip=-2,bc_free=-1};
    bc_type_current(i)=type;
    char& c=bc_type(i);
    if(c>=bc_free || type==bc_noslip) c=type;
}
template<int d> inline pair<FACE_INDEX<d>,bool> Node_To_Face(VECTOR<int,d> p)
{
    int a=0;
    for(int k=0;k<d;k++) a+=((p(k)&1)^1)*k;
    if(d==3 && a==d) return {{},false};
    return {{a,fdiv(p,2)},true};
}
//#####################################################################
// Function Set
//#####################################################################
template<class TV> void BOUNDARY_CONDITION_DOUBLE_FINE<TV>::
Set(const IMPLICIT_OBJECT<TV>* io,char type,std::function<T(const CB_DATA& data)> f,
    bool thin,bool invert,T contour,replace_enum rep)
{
    PHYSBAM_ASSERT(!thin || type==bc_noslip || type==bc_slip);
    RANGE<IV> domain=bc_type.domain;
    for(RANGE_ITERATOR<TV::m> it(domain);it.Valid();it.Next())
    {
        if((it.index.Sum()+TV::m)%2) continue; // don't compute if we don't need
        T phi=io->Extended_Phi(grid.Node_X(it.index))-contour;
        if(invert) phi=-phi;
        if(phi>0) bc_type_current(it.index)=0;
        else set_bc(bc_type,bc_type_current,it.index,type);
    }

    HASHTABLE<IV,pair<T,int> > bc_p_cur;
    HASHTABLE<FACE_INDEX<TV::m>,tuple<T,int,bool> > bc_u_cur;

    for(int axis=0;axis<TV::m;axis++)
    {
        int mn=type==bc_noslip?0:axis;
        int mx=type==bc_noslip?TV::m-1:axis;
        for(int a=mn;a<=mx;a++)
        {
            IV off=IV().Add_Axis(axis,1).Add_Axis(a,1);
            RANGE<IV> range(cdiv(bc_type.domain.min_corner-1+off,2),fdiv(bc_type.domain.max_corner-off,2));
            for(FACE_RANGE_ITERATOR<TV::m> it(range,RF::none,axis);it.Valid();it.Next())
            {
                IV i=(it.face.index*2+1).Add_Axis(it.face.axis,-1);
                IV i0=i.Add_Axis(a,-1),i1=i.Add_Axis(a,1);
                char t0=bc_type_current(i0),t1=bc_type_current(i1);
                if(t0>=0 && t1>=0) continue;
                if(t0==type && t1==type) continue;
                TV X0=grid.Node_X(i0),X1=grid.Node_X(i1);
                T phi0=io->Extended_Phi(grid.Node_X(i0))-contour;
                T phi1=io->Extended_Phi(grid.Node_X(i1))-contour;
                T t=phi0/(phi0-phi1);
                CB_DATA data={it.face,it.face.Cell_Index(t0>=0),X0+(X1-X0)*t};
                if(T* x=Add_Local_BC(bc_p_cur,bc_u_cur,data,type,a==axis,rep))
                    if(f) *x+=f(data);
            }
        }
    }
    Commit_BC(bc_p_cur,bc_u_cur,rep);
}

//#####################################################################
// Function Set
//#####################################################################
template<class TV> void BOUNDARY_CONDITION_DOUBLE_FINE<TV>::
Set(const T_SURFACE& surface,char type,std::function<T(const CB_DATA& data)> f,
    bool thin,replace_enum rep)
{
    PHYSBAM_ASSERT(!thin || type==bc_noslip || type==bc_slip);

    ARRAY<bool,IV> flood_fill_array(bc_type.domain.Edge_Lengths()/2+1);

    HASHTABLE<IV,pair<T,int> > bc_p_cur;
    HASHTABLE<FACE_INDEX<TV::m>,tuple<T,int,bool> > bc_u_cur;

    for(RANGE_ITERATOR<TV::m> it(IV()+2);it.Valid();it.Next())
    {
        if((it.index.Sum()+TV::m)%2) continue; // don't compute if we don't need
        bool normal_face=it.index==IV()+1;
        
        TV X0=grid.Node_X(bc_type.domain.min_corner+it.index);
        TV X1=grid.Node_X(bc_type.domain.max_corner-1-it.index);
        IV counts=(bc_type.domain.Edge_Lengths()+1)/2-it.index;
        IV offset=bc_type.domain.min_corner+it.index;
        GRID<TV> vis_grid(counts-1,RANGE<TV>(X0,X1));
        // fine = 2 * vis + offset;
        HASHTABLE<EDGE_INDEX<TV::m>,GRID_SURFACE_INTERSECTION_DATA<TV> > hash;
        Grid_Surface_Intersection(hash,vis_grid,surface,!thin);

        if(!thin)
        {
            Flood_Fill(flood_fill_array,hash);
            for(RANGE_ITERATOR<TV::m> it2(counts);it2.Valid();it2.Next())
                if(flood_fill_array(it2.index))
                    set_bc(bc_type,bc_type_current,2*it2.index+offset,type);
        }
        for(const auto& h:hash)
        {
            IV i[2]={2*h.first.index+offset};
            i[1]=i[0].Add_Axis(h.first.axis,2);
            TV X0=grid.Node_X(i[0]);
            TV X1=grid.Node_X(i[1]);
            auto face_p=Node_To_Face(i[0].Add_Axis(h.first.axis,1));
            if(!face_p.second) continue;
            CB_DATA data={face_p.first,face_p.first.Cell_Index(h.second.in[0]),TV()};

            IV im=i[0].Add_Axis(h.first.axis,1);
            if(bc_type.Valid_Index(im)) bc_type(im)|=1<<h.first.axis;

            T* x=Add_Local_BC(bc_p_cur,bc_u_cur,data,type,normal_face,rep);
            if(!x || !f) continue;
            T y=0;
            for(const auto& t:h.second.cut_elements)
            {
                data.X=X0+(X1-X0)*get<1>(t);
                data.e=get<0>(t);
                y+=f(data);
            }
            *x+=y/h.second.cut_elements.m;
        }
    }
    Commit_BC(bc_p_cur,bc_u_cur,rep);
}
//#####################################################################
// Function Set
//#####################################################################
template<class TV> void BOUNDARY_CONDITION_DOUBLE_FINE<TV>::
Set(T_OBJECT& object,char type,std::function<T(const CB_DATA& data)> f,bool thin,replace_enum rep)
{
    return Set(object.Get_Boundary_Object(),type,f,thin,rep);
}
//#####################################################################
// Function Set_Domain_Walls
//#####################################################################
template<class TV> void BOUNDARY_CONDITION_DOUBLE_FINE<TV>::
Set_Domain_Walls(int side_mask,char type,std::function<T(const CB_DATA& data)>
    f,replace_enum rep)
{
    for(RANGE_ITERATOR<TV::m> it(bc_type.domain,grid.Node_Indices(),RI::side_mask|RI::ghost,side_mask);it.Valid();it.Next())
    {
        if((it.index.Sum()+TV::m)%2) continue; // don't compute if we don't need
        set_bc(bc_type,bc_type_current,it.index,type);
    }
    for(RANGE_ITERATOR<TV::m> it(grid.Node_Indices(),grid.Node_Indices(-1),RI::side_mask|RI::ghost,side_mask);it.Valid();it.Next())
    {
        if((it.index.Sum()+TV::m)%2==0) continue;
        bc_type(it.index)=(1<<TV::m)-1;
    }
    HASHTABLE<IV,pair<T,int> > bc_p_cur;
    HASHTABLE<FACE_INDEX<TV::m>,tuple<T,int,bool> > bc_u_cur;
    if(type==bc_free)
    {
        for(FACE_RANGE_ITERATOR<TV::m> it(mac_grid.Cell_Indices(),0,0,RF::side_mask|RF::ghost,side_mask);it.Valid();it.Next())
        {
            CB_DATA data={it.face,it.face.Cell_Index(it.side%2),mac_grid.Face_X(it.face)};
            if(T* x=Add_Local_BC(bc_p_cur,bc_u_cur,data,type,false,rep))
                if(f) *x+=f(data);
        }
    }
    else
    {
        bool ns=type==bc_noslip;
        auto flags=RF::side_mask|RF::ghost;
        if(ns) flags|=RF::skip_outer;
        for(FACE_RANGE_ITERATOR<TV::m> it(mac_grid.Cell_Indices(),ns,0,flags,side_mask);it.Valid();it.Next())
        {
            CB_DATA data={it.face,IV(),mac_grid.Face_X(it.face)};
            if(T* x=Add_Local_BC(bc_p_cur,bc_u_cur,data,type,it.face.axis==it.side/2,rep))
                if(f) *x+=f(data);
        }
    }
    Commit_BC(bc_p_cur,bc_u_cur,rep);
}
//#####################################################################
// Function Get_Pressure_Boundary_Conditions
//#####################################################################
template<class TV> void BOUNDARY_CONDITION_DOUBLE_FINE<TV>::
Get_Pressure_Boundary_Conditions(ARRAY<bool,IV>& psi_D,
    ARRAY<bool,FACE_INDEX<TV::m> >& psi_N,ARRAY<T,IV>& p,
    ARRAY<T,FACE_INDEX<TV::m> >& u,int ghost_in) const
{
    PHYSBAM_ASSERT(ghost_in<=ghost);
    psi_N.Fill(false);

    // Set psi_D, clear non-dof pressures
    for(RANGE_ITERATOR<TV::m> it(mac_grid.Cell_Indices(ghost_in));it.Valid();it.Next())
    {
        IV bc_i=it.index*2+1;
        char bc_t=bc_type(bc_i);
        if(bc_t<0)
        {
            psi_D(it.index)=true;
            p(it.index)=0;
        }
        else psi_D(it.index)=false;
    }

    // Set bc pressures where they are used
    RANGE<IV> valid_domain(mac_grid.Cell_Indices(1));
    for(auto z:bc_p)
    {
        if(!valid_domain.Lazy_Inside_Half_Open(z.first)) continue;
        bool has_in=false;
        for(int a=0;a<TV::m;a++)
            for(int s=-1;s<2;s+=2)
                if(!psi_D(z.first.Add_Axis(a,s)))
                {
                    has_in=true;
                    break;
                }
        if(has_in) p(z.first)=z.second.first/z.second.second;
    }

    // bc_a/bc_b 0  bc_free  bc_*slip
    // 0         -     -     u=0, psi_N=true
    // bc_free         u=0   u=0
    // bc_*slip              u=0
    for(FACE_RANGE_ITERATOR<TV::m> it(psi_N.domain,mac_grid.Cell_Indices(),RF::ghost|RF::skip_inner);it.Valid();it.Next())
        u(it.face)=0;
    for(FACE_RANGE_ITERATOR<TV::m> it(mac_grid.Cell_Indices());it.Valid();it.Next())
    {
        char bc_a=bc_type(it.face.First_Cell_Index()*2+1);
        char bc_b=bc_type(it.face.Second_Cell_Index()*2+1);
        if(bc_a<0 && bc_b<0) u(it.face)=0;
        else if(bc_a<bc_free || bc_b<bc_free)
        {
            u(it.face)=0;
            psi_N(it.face)=true;
        }
    }
    for(auto z:bc_u)
    {
//        if(!get<2>(z.second)) continue;
        if(!mac_grid.Valid_Face(z.first,0)) continue;
        char bc_a=bc_type(z.first.First_Cell_Index()*2+1);
        char bc_b=bc_type(z.first.Second_Cell_Index()*2+1);
        if(bc_a>=0 || bc_b>=0)
        {
            psi_N(z.first)=true;
            u(z.first)=get<0>(z.second)/get<1>(z.second);
        }
    }
}
//#####################################################################
// Function Get_Viscosity_Boundary_Conditions
//#####################################################################
template<class TV> void BOUNDARY_CONDITION_DOUBLE_FINE<TV>::
Get_Viscosity_Boundary_Conditions(ARRAY<bool,IV>& psi_D,
    ARRAY<bool,FACE_INDEX<TV::m> >& psi_N,ARRAY<T,IV>& u,int axis,int ghost_in) const
{
    PHYSBAM_ASSERT(ghost_in<=ghost);
    IV offset=(IV()+1).Add_Axis(axis,-1);
    GRID<TV> face_grid(mac_grid.Get_Face_MAC_Grid(axis));
    bool dump=DEBUG_SUBSTEPS::write_substeps_level>=2;

    if(dump)
    {
        for(RANGE_ITERATOR<TV::m> it(bc_type.domain);it.Valid();it.Next())
        {
            if((it.index.Sum()+TV::m)%2) continue; // don't compute if we don't need
            char t=bc_type(it.index);
            if(t<0) Add_Point(grid.Node_X(it.index),VECTOR<T,3>(t==bc_free,t==bc_noslip,t==bc_slip));
        }
        PHYSBAM_DEBUG_WRITE_SUBSTEP("bc_type",2);
    }

    for(RANGE_ITERATOR<TV::m> it(mac_grid.Face_Indices(axis,ghost_in));it.Valid();it.Next())
    {
        IV bc_i=it.index*2+offset;
        for(int a=0;a<TV::m;a++)
        {
            IV i0=bc_i.Add_Axis(a,-1),i1=bc_i.Add_Axis(a,1);
            char bc0=bc_type(i0),bc1=bc_type(i1);
            if(bc0<0 && bc1<0)
            {
                psi_D(it.index)=true;
                if(dump) Add_Point(face_grid.Cell_X(it.index),VECTOR<T,3>(1,0,1));
                continue;
            }
            int side=0;
            if(bc0<bc1){std::swap(bc0,bc1);side=1;}
            if(bc1>=0) continue;
            psi_D(it.index)=true;
            if(dump) Add_Point(face_grid.Cell_X(it.index),VECTOR<T,3>(1,0,1));
            if(bc1==bc_free || (bc1==bc_slip && axis!=a))
            {
                psi_N(a,it.index.Add_Axis(a,side))=true;
                if(dump) Add_Point(face_grid.Face_X({a,it.index.Add_Axis(a,side)}),VECTOR<T,3>(0,1,1));
            }
            else
            {
                tuple<T,int,bool> z;
                if(bc_u.Get({axis,it.index},z))
                    u(it.index)=get<0>(z)/get<1>(z);
            }
        }
    }
    for(auto z:bc_u)
    {
        if(z.first.axis!=axis) continue;
        if(!mac_grid.Valid_Face(z.first,ghost_in)) continue;
        psi_D(z.first.index)=true;
        if(dump) Add_Point(face_grid.Cell_X(z.first.index),VECTOR<T,3>(.5,.5,1));
        u(z.first.index)=get<0>(z.second)/get<1>(z.second);
    }
    if(dump) PHYSBAM_DEBUG_WRITE_SUBSTEP("Get_Viscosity_Boundary_Conditions",2);
}
//#####################################################################
// Function Set_Region
//#####################################################################
template<class TV> void BOUNDARY_CONDITION_DOUBLE_FINE<TV>::
Set_Region(IV mac_cell_seed,char type)
{
    for(RANGE_ITERATOR<TV::m> it(IV()+3);it.Valid();it.Next())
    {
        if((it.index.Sum()+TV::m)%2) continue; // don't compute if we don't need
        IV seed_index=2*mac_cell_seed+it.index;
        if(bc_type(seed_index)<0) continue;
        ARRAY<IV> todo;
        bc_type(seed_index)=type;
        todo.Append(seed_index);
        while(todo.m)
        {
            IV cell=todo.Last();
            todo.Pop();
            for(int a=0;a<TV::m;a++)
            {
                for(int s=-1;s<2;s+=2)
                {
                    IV index=cell.Add_Axis(a,s*2);
                    if(bc_type(index)<0) continue;
                    if(bc_type(cell.Add_Axis(a,s))&(1<<a)) continue;
                    bc_type(index)=type;
                    todo.Append(index);
                }
            }
        }
    }
}
//#####################################################################
// Function Commit_BC
//#####################################################################
template<class TV> void BOUNDARY_CONDITION_DOUBLE_FINE<TV>::
Commit_BC(const HASHTABLE<IV,pair<T,int> >& bc_p_cur,
    const HASHTABLE<FACE_INDEX<TV::m>,tuple<T,int,bool> >& bc_u_cur,
    replace_enum rep)
{
    if(rep==average)
    {
        for(auto& z:bc_p_cur)
        {
            auto& p=bc_p.Get_Or_Insert(z.first);
            p.first+=z.second.first/z.second.second;
            p.second++;
        }
        for(auto& z:bc_u_cur)
        {
            auto& p=bc_u.Get_Or_Insert(z.first);
            get<0>(p)+=get<0>(z.second)/get<1>(z.second);
            get<1>(p)++;
            get<2>(p)|=get<2>(z.second);
        }
    }
    else
    {
        for(auto& z:bc_p_cur)
        {
            auto& p=bc_p.Get_Or_Insert(z.first);
            p.first=z.second.first/z.second.second;
            p.second=1;
        }
        for(auto& z:bc_u_cur)
        {
            auto& p=bc_u.Get_Or_Insert(z.first);
            get<0>(p)=get<0>(z.second)/get<1>(z.second);
            get<1>(p)=1;
            get<2>(p)=get<2>(z.second);
        }
    }
}
//#####################################################################
// Function Add_Local_BC
//#####################################################################
template<class TV> auto BOUNDARY_CONDITION_DOUBLE_FINE<TV>::
Add_Local_BC(HASHTABLE<IV,pair<T,int> >& bc_p_cur,
    HASHTABLE<FACE_INDEX<TV::m>,tuple<T,int,bool> >& bc_u_cur,
    const CB_DATA& data,char type,bool normal_face,replace_enum rep) const -> T*
{
    if(type==bc_free)
    {
        if(rep==skip && bc_p.Contains(data.cell)) return 0;
        auto& z=bc_p_cur.Get_Or_Insert(data.cell);
        z.second++;
        return &z.first;
    }
    if(rep==skip && bc_u.Contains(data.face)) return 0;
    auto& z=bc_u_cur.Get_Or_Insert(data.face);
    get<1>(z)++;
    if(normal_face) get<2>(z)=true;
    return &get<0>(z);
}
template class BOUNDARY_CONDITION_DOUBLE_FINE<VECTOR<float,1> >;
template class BOUNDARY_CONDITION_DOUBLE_FINE<VECTOR<float,2> >;
template class BOUNDARY_CONDITION_DOUBLE_FINE<VECTOR<float,3> >;
template class BOUNDARY_CONDITION_DOUBLE_FINE<VECTOR<double,1> >;
template class BOUNDARY_CONDITION_DOUBLE_FINE<VECTOR<double,2> >;
template class BOUNDARY_CONDITION_DOUBLE_FINE<VECTOR<double,3> >;
}
