//#####################################################################
// Copyright 2005-2007, Eran Guendelman, Geoffrey Irving, Nipun Kwatra, Avi Robinson-Mosher, 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.
//#####################################################################
// Class FACE_ARRAYS
//#####################################################################
#ifndef __FACE_ARRAYS__
#define __FACE_ARRAYS__

#include <Core/Arrays/ARRAYS_FORWARD.h>
#include <Core/Arrays/ARRAYS_ND.h>
#include <Core/Arrays/ARRAYS_ND_VIEW.h>
#include <Core/Grids/FACE_INDEX.h>
#include <Core/Vectors/SCALAR_POLICY.h>
namespace PhysBAM{

template<class T,int d>
class ARRAY<T,FACE_INDEX<d> >
{
    typedef VECTOR<T,d> TV;typedef VECTOR<int,d> IV;
    struct UNUSABLE{};
public:
    ARRAY_VIEW<T> array; // one-dimensional data storage

    enum {dimension=TV::m};
    template<class T2> struct REBIND{typedef ARRAY<T2,FACE_INDEX<dimension> > TYPE;};
    typedef T ELEMENT;typedef FACE_INDEX<d> INDEX;typedef T& RESULT_TYPE;typedef const T& CONST_RESULT_TYPE;
    typedef SCALAR_TYPE<T> SCALAR;

    RANGE<IV> domain;
    VECTOR<ARRAY_VIEW<T,IV>,dimension> data;

    ARRAY()
    {}

    ARRAY(const RANGE<IV>& domain_input)
        :ARRAY(domain_input,init_new,T())
    {}

    ARRAY(const RANGE<IV>& domain_input,NO_INIT)
    {Initialize(domain_input);}

    ARRAY(const RANGE<IV>& domain_input,INIT_NEW,const T& initialization_value=T())
    {Initialize(domain_input);array.Fill(initialization_value);}
    
    ARRAY(const ARRAY& old_array)
    {Initialize(old_array.domain);array=old_array.array;}

    ARRAY(ARRAY&& old_array)
    {Exchange(old_array);}

    virtual ~ARRAY()
    {delete[] array.base_pointer;}

    const RANGE<IV>& Domain_Indices() const
    {return domain;}

    ARRAY& operator=(const ARRAY& source)
    {if(source.Domain_Indices()!=Domain_Indices()) Resize(source.Domain_Indices(),no_init);Copy(source);return *this;}

    ARRAY& operator=(ARRAY&&) = default;

    void Resize(const RANGE<IV>& domain_input)
    {Resize(domain_input,true,true,T());}

    void Resize(const RANGE<IV>& domain_input,NO_INIT)
    {Resize(domain_input,false,false,T());}

    void Resize(const RANGE<IV>& domain_input,INIT_NEW,const T& initialization_value=T())
    {Resize(domain_input,true,true,initialization_value);}

    void Resize(const RANGE<IV>& domain_input,INIT_ALL,const T& initialization_value=T())
    {Resize(domain_input,true,false,initialization_value);}

    static int Calculate_Raw_Size(const RANGE<IV>& domain_input)
    {
        IV sizes=domain_input.Edge_Lengths();
        int total=0;
        for(int i=0;i<dimension;i++)
            total+=sizes.Add_Axis(i,1).Product();
        return total;
    }

    void Initialize_Using_Array(const RANGE<IV>& domain_input)
    {
        domain=domain_input;
        T* p=array.Get_Array_Pointer();
        for(int i=0;i<dimension;i++)
        {
            RANGE<IV> new_domain(domain_input);
            new_domain.max_corner(i)++;
            data(i).Set(p,new_domain);
            p+=data(i).array.m;
        }
    }

    void Initialize(const RANGE<IV>& domain_input)
    {
        if(domain_input.Empty()) return;
        int new_length=Calculate_Raw_Size(domain_input);
        assert(!array.Get_Array_Pointer());
        array.Set(new T[new_length],new_length);
        Initialize_Using_Array(domain_input);
    }

  protected:

    virtual void Resize(const RANGE<IV>& domain_input,const bool initialize_new_elements,const bool copy_existing_elements,const T& initialization_value)
    {
        if(domain_input==domain)
        {
            if(initialize_new_elements && !copy_existing_elements)
                array.Fill(initialization_value);
            return;
        }
        ARRAY new_array(domain_input);
        if(initialize_new_elements)
            new_array.Fill(initialization_value);
        
        if(copy_existing_elements && !domain.Empty())
            for(int i=0;i<dimension;i++)
                ARRAY_VIEW<T,IV>::Limited_Shifted_Get(new_array.data(i),data(i),IV());

        Exchange(new_array);
    }
  public:

    void Clean_Memory()
    {
        array.Set(0,0);
        domain=RANGE<IV>::Empty_Box();
        for(auto& a:data) a.Set({});
    }

    void Delete_Pointers_And_Clean_Memory()
    {for(int i=0;i<array.m;i++) delete array.base_pointer[i];Clean_Memory();}

    int Number_Of_Ghost_Cells() const
    {return -domain.min_corner.x;}

    T& operator()(const int axis,const IV& face)
    {return Component(axis)(face);}

    const T& operator()(const int axis,const IV& face) const
    {return Component(axis)(face);}

    T& operator()(const FACE_INDEX<dimension>& index)
    {return Component(index.axis)(index.index);}

    const T& operator()(const FACE_INDEX<dimension>& index) const
    {return Component(index.axis)(index.index);}

    bool Valid_Index(const FACE_INDEX<dimension>& index) const
    {return ((unsigned)index.axis<(unsigned)dimension) && Component(index.axis).Valid_Index(index.index);}

    int Standard_Index(const FACE_INDEX<TV::m>& index) const
    {assert((unsigned)index.axis<(unsigned)dimension);
    return data(index.axis).array.base_pointer-array.base_pointer+data(index.axis).Standard_Index(index.index);}

    ARRAY_VIEW<T,IV>& Component(const int axis)
    {assert((unsigned)axis<dimension);return data(axis);}

    const ARRAY_VIEW<T,IV>& Component(const int axis) const
    {assert((unsigned)axis<dimension);return data(axis);}

    T* begin(){return array.begin();}
    const T* begin() const {return array.begin();}
    T* end(){return array.end();}
    const T* end() const {return array.end();}

    TV Cell_Centered_Average(const IV& cell_index) const
    {TV average;
    for(int i=0;i<dimension;i++) average(i)=(T).5*(data(i)(cell_index)+data(i)(cell_index+IV::Axis_Vector(i)));
    return average;}

    void Set_All_Faces(const T& value,const IV& cell_index)
    {for(int i=0;i<dimension;i++) data(i)(cell_index)=data(i)(cell_index+IV::Axis_Vector(i))=value;}

    static void Extract_Dimension(const ARRAY& old_array,ARRAY<typename ELEMENT_OF_VECTOR<T>::TYPE,FACE_INDEX<dimension> >& extracted_array,const IV& dimensions_to_extract)
    {for(int i=0;i<dimension;i++) ARRAY_VIEW<T,IV>::Extract_Dimension(old_array.data(i),extracted_array.data(i),dimensions_to_extract(i));}

    void Fill(const T& constant)
    {for(int i=0;i<dimension;i++) data(i).Fill(constant);}

    void Copy(const ARRAY& old_copy)
    {for(int i=0;i<dimension;i++) data(i).Copy(old_copy.data(i));}

    template<class T2>
    void Copy(const T2 c,const ARRAY& old)
    {for(int i=0;i<dimension;i++) data(i).Copy(c,old.data(i));}

    void Fill(const TV& value)
    {for(int i=0;i<dimension;i++) data(i).Fill(value.x);}

    template<class T2>
    void Copy(const T2 c1,const ARRAY& v1,const T2 c2,const ARRAY& v2)
    {for(int i=0;i<dimension;i++) data(i).Copy(c1,v1.data(i),c2,v2.data(i));}

    static void Put(const ARRAY& old_copy,ARRAY& new_copy)
    {for(int i=0;i<dimension;i++) ARRAY_VIEW<T,IV>::Put(old_copy.data(i),new_copy.data(i));}

    TV Max_Abs() const
    {TV maxabs_vector;
    for(int i=0;i<dimension;i++) maxabs_vector(i)=data(i).Max_Abs();
    return maxabs_vector;}

    static void Exchange(ARRAY& a,ARRAY& b)
    {a.Exchange(b);}

    void Exchange(ARRAY& a)
    {exchange(domain,a.domain);array.Exchange(a.array);
    for(int i=0;i<dimension;i++) ARRAY_VIEW<T,IV>::Exchange(data(i),a.data(i));}

    template<class T_INDICES>
    INDIRECT_ARRAY<ARRAY,T_INDICES&> Subset(const T_INDICES& indices)
    {return INDIRECT_ARRAY<ARRAY,T_INDICES&>(*this,indices);}

    template<class T_INDICES>
    INDIRECT_ARRAY<const ARRAY,T_INDICES&> Subset(const T_INDICES& indices) const
    {return INDIRECT_ARRAY<const ARRAY,T_INDICES&>(*this,indices);}

    template<class RW> void Read(UNTYPED_ISTREAM<RW>& input)
    {
        Get(input,'(');
        RANGE<IV> range;
        Read_Binary(input,range);
        ARRAY new_array(range);
        Read_Binary(input,new_array.array);
        if constexpr(!is_same<RW,char>::value)
            for(int i=0;i<d;i++)
                Read_Binary(input,range);
        Get(input,')');
        Exchange(new_array);
    }

    template<class RW> void Write(UNTYPED_OSTREAM<RW>& output) const
    {
        Beg(output,'(');
        Write_Binary(output,domain,array);
        if constexpr(!is_same<RW,char>::value)
            for(int i=0;i<d;i++)
                Write_Binary(output,data(i).domain);
        End(output,')');
    }
//#####################################################################
};

template<class T> inline std::ostream& operator<<(std::ostream& output_stream,const ARRAY<T,FACE_INDEX<1> >& a)
{
    for(int i=a.domain.min_corner.x;i<a.domain.max_corner.x+1;i++)
        output_stream<<a.Component(0)(i)<<" ";
    output_stream<<std::endl;
    return output_stream;
}

template<class T,int d> inline std::ostream& operator<<(std::ostream& output_stream,const ARRAY<T,FACE_INDEX<d> >& a)
{
    for(int i=0;i<d;i++) output_stream<<a.data(i)<<std::endl;
    return output_stream;
}

// Avoid instantiating common non-trivial functions in large numbers of files.
#define COMMA ,
#define COMMON(T,F,IV)                                                  \
    extern template ARRAY<T,F>::ARRAY(const ARRAY<T,F>&);               \
    extern template ARRAY<T,F>::ARRAY(const RANGE<IV>&);                \
    extern template ARRAY<T,F>::ARRAY(const RANGE<IV>&,INIT_NEW,const T&); \
    extern template ARRAY<T,F>::~ARRAY();                               \
    extern template void ARRAY<T,F>::Initialize(const RANGE<IV>&);      \
    extern template void ARRAY<T,F>::Resize(const RANGE<IV>&);          \
    extern template void ARRAY<T,F>::Resize(const RANGE<IV>&,NO_INIT);  \
    extern template void ARRAY<T,F>::Resize(const RANGE<IV>&,bool,bool,const T&);
#define COMMON1(T,d)                                            \
    COMMON(T,FACE_INDEX<d>,VECTOR<int COMMA d>);                \
    COMMON(VECTOR<T COMMA d>,FACE_INDEX<d>,VECTOR<int COMMA d>)
#define COMMON2(T) COMMON1(T,1);COMMON1(T,2);COMMON1(T,3)
COMMON2(bool);
COMMON2(int);
COMMON2(float);
COMMON2(double);
#undef COMMON2
#undef COMMON1
#undef COMMON
#undef COMMA
}
#endif
