//#####################################################################
// Copyright 2002-2007, Ronald Fedkiw, Geoffrey Irving, Nipun Kwatra, Frank Losasso.
// This file is part of PhysBAM whose distribution is governed by the license contained in the accompanying file PHYSBAM_COPYRIGHT.txt.
//#####################################################################
#include <Core/Boundaries/PERIODIC_BOUNDARY_HELPER.h>
#include <Core/Math_Tools/RANGE_ITERATOR.h>
#include <Core/Matrices/MATRIX.h>
namespace PhysBAM{
template<class T2,class IV> void
Fill_Ghost_Periodic(const IV& numbers_of_cells,
    ARRAYS_ND_BASE<T2,VECTOR<int,IV::m> >& u,
    const VECTOR<bool,IV::m>& is_periodic)
{
    VECTOR<RANGE<IV>,2*IV::m> ghost_domains;
    Valid_Periodic_Regions(u.domain,{IV(),numbers_of_cells},is_periodic,ghost_domains);
    for(int a=IV::m-1;a>=0;a--)
    {
        int shift=numbers_of_cells(a);
        for(int s=0;s<2;s++)
        {
            for(RANGE_ITERATOR<IV::m> it(ghost_domains(2*a+s));it.Valid();it.Next())
                u(it.index)=u(it.index.Add_Axis(a,shift));
            shift=-shift;
        }
    }
}
//#####################################################################
// Function Fill_Ghost_Cells_Periodic
//#####################################################################
template<class OP,class T2,class IV> void
Fill_Ghost_Periodic_Accum(const IV& numbers_of_cells,
    ARRAYS_ND_BASE<T2,VECTOR<int,IV::m> >& u,
    const VECTOR<bool,IV::m>& is_periodic,const OP& op)
{
    VECTOR<RANGE<IV>,2*IV::m> ghost_domains;
    Valid_Periodic_Regions(u.domain,{IV(),numbers_of_cells},is_periodic,ghost_domains);
    for(int a=0;a<IV::m;a++)
    {
        int shift=numbers_of_cells(a);
        for(int s=0;s<2;s++)
        {
            for(RANGE_ITERATOR<IV::m> it(ghost_domains(2*a+s));it.Valid();it.Next())
            {
                auto& v=u(it.index.Add_Axis(a,shift));
                v=op(v,u(it.index));
            }
            shift=-shift;
        }
    }
    for(int a=IV::m-1;a>=0;a--)
    {
        int shift=numbers_of_cells(a);
        for(int s=0;s<2;s++)
        {
            for(RANGE_ITERATOR<IV::m> it(ghost_domains(2*a+s));it.Valid();it.Next())
                u(it.index)=u(it.index.Add_Axis(a,shift));
            shift=-shift;
        }
    }
}
template<int d> void
Valid_Periodic_Regions(const RANGE<VECTOR<int,d> >& array_domain,
    const RANGE<VECTOR<int,d> >& grid_domain,
    const VECTOR<bool,d>& is_periodic,
    VECTOR<RANGE<VECTOR<int,d> >,2*d>& ghost_domains)
{
    typedef VECTOR<int,d> IV;
    RANGE<IV> domain(array_domain);
    IV sizes=grid_domain.Edge_Lengths();
    for(int a=0;a<d;a++)
    {
        if(!is_periodic(a)) continue;
        RANGE<IV> lo(domain),hi(domain);
        lo.max_corner(a)=min(lo.max_corner(a),grid_domain.min_corner(a),
            domain.max_corner(a)-sizes(a));
        hi.min_corner(a)=max(hi.min_corner(a),grid_domain.max_corner(a),
            domain.min_corner(a)+sizes(a));
        ghost_domains(2*a)=lo;
        ghost_domains(2*a+1)=hi;

        domain.min_corner(a)=max(domain.min_corner(a),lo.max_corner(a));
        domain.max_corner(a)=min(domain.max_corner(a),hi.min_corner(a));
    }
}
template void Fill_Ghost_Periodic<MATRIX<double,2,2>,VECTOR<int,2> >(
    VECTOR<int,2> const&,ARRAYS_ND_BASE<MATRIX<double,2,2>,VECTOR<int,VECTOR<int,2>::m> >&,
    VECTOR<bool,VECTOR<int,2>::m> const&);
template void Fill_Ghost_Periodic<MATRIX<double,3,3>,VECTOR<int,3> >(
    VECTOR<int,3> const&,ARRAYS_ND_BASE<MATRIX<double,3,3>,VECTOR<int,VECTOR<int,3>::m> >&,
    VECTOR<bool,VECTOR<int,3>::m> const&);
template void Fill_Ghost_Periodic<MATRIX<float,2,2>,VECTOR<int,2> >(
    VECTOR<int,2> const&,ARRAYS_ND_BASE<MATRIX<float,2,2>,VECTOR<int,VECTOR<int,2>::m> >&,
    VECTOR<bool,VECTOR<int,2>::m> const&);
template void Fill_Ghost_Periodic<MATRIX<float,3,3>,VECTOR<int,3> >(
    VECTOR<int,3> const&,ARRAYS_ND_BASE<MATRIX<float,3,3>,VECTOR<int,VECTOR<int,3>::m> >&,
    VECTOR<bool,VECTOR<int,3>::m> const&);
template void Fill_Ghost_Periodic<SYMMETRIC_MATRIX<double,2>,VECTOR<int,2> >(
    VECTOR<int,2> const&,ARRAYS_ND_BASE<SYMMETRIC_MATRIX<double,2>,VECTOR<int,VECTOR<int,2>::m> >&,
    VECTOR<bool,VECTOR<int,2>::m> const&);
template void Fill_Ghost_Periodic<SYMMETRIC_MATRIX<double,3>,VECTOR<int,3> >(
    VECTOR<int,3> const&,ARRAYS_ND_BASE<SYMMETRIC_MATRIX<double,3>,VECTOR<int,VECTOR<int,3>::m> >&,
    VECTOR<bool,VECTOR<int,3>::m> const&);
template void Fill_Ghost_Periodic<SYMMETRIC_MATRIX<float,2>,VECTOR<int,2> >(
    VECTOR<int,2> const&,ARRAYS_ND_BASE<SYMMETRIC_MATRIX<float,2>,VECTOR<int,VECTOR<int,2>::m> >&,
    VECTOR<bool,VECTOR<int,2>::m> const&);
template void Fill_Ghost_Periodic<SYMMETRIC_MATRIX<float,3>,VECTOR<int,3> >(
    VECTOR<int,3> const&,ARRAYS_ND_BASE<SYMMETRIC_MATRIX<float,3>,VECTOR<int,VECTOR<int,3>::m> >&,
    VECTOR<bool,VECTOR<int,3>::m> const&);
template void Fill_Ghost_Periodic<VECTOR<double,3>,VECTOR<int,2> >(
    VECTOR<int,2> const&,ARRAYS_ND_BASE<VECTOR<double,3>,VECTOR<int,VECTOR<int,2>::m> >&,
    VECTOR<bool,VECTOR<int,2>::m> const&);
template void Fill_Ghost_Periodic<VECTOR<float,3>,VECTOR<int,2> >(
    VECTOR<int,2> const&,ARRAYS_ND_BASE<VECTOR<float,3>,VECTOR<int,VECTOR<int,2>::m> >&,
    VECTOR<bool,VECTOR<int,2>::m> const&);
template void Fill_Ghost_Periodic<bool,VECTOR<int,2> >(
    VECTOR<int,2> const&,ARRAYS_ND_BASE<bool,VECTOR<int,VECTOR<int,2>::m> >&,
    VECTOR<bool,VECTOR<int,2>::m> const&);
template void Fill_Ghost_Periodic<bool,VECTOR<int,3> >(
    VECTOR<int,3> const&,ARRAYS_ND_BASE<bool,VECTOR<int,VECTOR<int,3>::m> >&,
    VECTOR<bool,VECTOR<int,3>::m> const&);
template void Fill_Ghost_Periodic<double,VECTOR<int,2> >(
    VECTOR<int,2> const&,ARRAYS_ND_BASE<double,VECTOR<int,VECTOR<int,2>::m> >&,
    VECTOR<bool,VECTOR<int,2>::m> const&);
template void Fill_Ghost_Periodic<double,VECTOR<int,3> >(
    VECTOR<int,3> const&,ARRAYS_ND_BASE<double,VECTOR<int,VECTOR<int,3>::m> >&,
    VECTOR<bool,VECTOR<int,3>::m> const&);
template void Fill_Ghost_Periodic<float,VECTOR<int,2> >(
    VECTOR<int,2> const&,ARRAYS_ND_BASE<float,VECTOR<int,VECTOR<int,2>::m> >&,
    VECTOR<bool,VECTOR<int,2>::m> const&);
template void Fill_Ghost_Periodic<float,VECTOR<int,3> >(
    VECTOR<int,3> const&,ARRAYS_ND_BASE<float,VECTOR<int,VECTOR<int,3>::m> >&,
    VECTOR<bool,VECTOR<int,3>::m> const&);
template void Fill_Ghost_Periodic<int,VECTOR<int,2> >(
    VECTOR<int,2> const&,ARRAYS_ND_BASE<int,VECTOR<int,VECTOR<int,2>::m> >&,
    VECTOR<bool,VECTOR<int,2>::m> const&);
template void Fill_Ghost_Periodic<int,VECTOR<int,3> >(
    VECTOR<int,3> const&,ARRAYS_ND_BASE<int,VECTOR<int,VECTOR<int,3>::m> >&,
    VECTOR<bool,VECTOR<int,3>::m> const&);
template void Fill_Ghost_Periodic<VECTOR<double,2>,VECTOR<int,2> >(
    VECTOR<int,2> const&,ARRAYS_ND_BASE<VECTOR<double,2>,VECTOR<int,VECTOR<int,2>::m> >&,
    VECTOR<bool,VECTOR<int,2>::m> const&);
template void Fill_Ghost_Periodic<VECTOR<double,3>,VECTOR<int,3> >(
    VECTOR<int,3> const&,ARRAYS_ND_BASE<VECTOR<double,3>,VECTOR<int,VECTOR<int,3>::m> >&,
    VECTOR<bool,VECTOR<int,3>::m> const&);
template void Fill_Ghost_Periodic<VECTOR<float,2>,VECTOR<int,2> >(
    VECTOR<int,2> const&,ARRAYS_ND_BASE<VECTOR<float,2>,VECTOR<int,VECTOR<int,2>::m> >&,
    VECTOR<bool,VECTOR<int,2>::m> const&);
template void Fill_Ghost_Periodic<VECTOR<float,3>,VECTOR<int,3> >(
    VECTOR<int,3> const&,ARRAYS_ND_BASE<VECTOR<float,3>,VECTOR<int,VECTOR<int,3>::m> >&,
    VECTOR<bool,VECTOR<int,3>::m> const&);
template void Fill_Ghost_Periodic_Accum<PBH::ADD,VECTOR<double,2>,VECTOR<int,2> >(
    VECTOR<int,2> const&,ARRAYS_ND_BASE<VECTOR<double,2>,VECTOR<int,VECTOR<int,2>::m> >&,
    VECTOR<bool,VECTOR<int,2>::m> const&,PBH::ADD const&);
template void Fill_Ghost_Periodic_Accum<PBH::ADD,VECTOR<double,3>,VECTOR<int,3> >(
    VECTOR<int,3> const&,ARRAYS_ND_BASE<VECTOR<double,3>,VECTOR<int,VECTOR<int,3>::m> >&,
    VECTOR<bool,VECTOR<int,3>::m> const&,PBH::ADD const&);
template void Fill_Ghost_Periodic_Accum<PBH::ADD,VECTOR<float,2>,VECTOR<int,2> >(
    VECTOR<int,2> const&,ARRAYS_ND_BASE<VECTOR<float,2>,VECTOR<int,VECTOR<int,2>::m> >&,
    VECTOR<bool,VECTOR<int,2>::m> const&,PBH::ADD const&);
template void Fill_Ghost_Periodic_Accum<PBH::ADD,VECTOR<float,3>,VECTOR<int,3> >(
    VECTOR<int,3> const&,ARRAYS_ND_BASE<VECTOR<float,3>,VECTOR<int,VECTOR<int,3>::m> >&,
    VECTOR<bool,VECTOR<int,3>::m> const&,PBH::ADD const&);
template void Fill_Ghost_Periodic_Accum<PBH::ADD,double,VECTOR<int,2> >(
    VECTOR<int,2> const&,ARRAYS_ND_BASE<double,VECTOR<int,VECTOR<int,2>::m> >&,
    VECTOR<bool,VECTOR<int,2>::m> const&,PBH::ADD const&);
template void Fill_Ghost_Periodic_Accum<PBH::ADD,double,VECTOR<int,3> >(
    VECTOR<int,3> const&,ARRAYS_ND_BASE<double,VECTOR<int,VECTOR<int,3>::m> >&,
    VECTOR<bool,VECTOR<int,3>::m> const&,PBH::ADD const&);
template void Fill_Ghost_Periodic_Accum<PBH::ADD,float,VECTOR<int,2> >(
    VECTOR<int,2> const&,ARRAYS_ND_BASE<float,VECTOR<int,VECTOR<int,2>::m> >&,
    VECTOR<bool,VECTOR<int,2>::m> const&,PBH::ADD const&);
template void Fill_Ghost_Periodic_Accum<PBH::ADD,float,VECTOR<int,3> >(
    VECTOR<int,3> const&,ARRAYS_ND_BASE<float,VECTOR<int,VECTOR<int,3>::m> >&,
    VECTOR<bool,VECTOR<int,3>::m> const&,PBH::ADD const&);
template void Fill_Ghost_Periodic_Accum<PBH::OR,bool,VECTOR<int,2> >(
    VECTOR<int,2> const&,ARRAYS_ND_BASE<bool,VECTOR<int,VECTOR<int,2>::m> >&,
    VECTOR<bool,VECTOR<int,2>::m> const&,PBH::OR const&);
template void Fill_Ghost_Periodic_Accum<PBH::OR,bool,VECTOR<int,3> >(
    VECTOR<int,3> const&,ARRAYS_ND_BASE<bool,VECTOR<int,VECTOR<int,3>::m> >&,
    VECTOR<bool,VECTOR<int,3>::m> const&,PBH::OR const&);
template void Fill_Ghost_Periodic_Accum<PBH::OR,int,VECTOR<int,2> >(
    VECTOR<int,2> const&,ARRAYS_ND_BASE<int,VECTOR<int,VECTOR<int,2>::m> >&,
    VECTOR<bool,VECTOR<int,2>::m> const&,PBH::OR const&);
template void Fill_Ghost_Periodic_Accum<PBH::OR,int,VECTOR<int,3> >(
    VECTOR<int,3> const&,ARRAYS_ND_BASE<int,VECTOR<int,VECTOR<int,3>::m> >&,
    VECTOR<bool,VECTOR<int,3>::m> const&,PBH::OR const&);
template void Fill_Ghost_Periodic<double,VECTOR<int,1> >(
    VECTOR<int,1> const&,ARRAYS_ND_BASE<double,VECTOR<int,VECTOR<int,1>::m> >&,
    VECTOR<bool,VECTOR<int,1>::m> const&);
template void Fill_Ghost_Periodic<float,VECTOR<int,1> >(
    VECTOR<int,1> const&,ARRAYS_ND_BASE<float,VECTOR<int,VECTOR<int,1>::m> >&,
    VECTOR<bool,VECTOR<int,1>::m> const&);
}
