//#####################################################################
// Copyright 2002-2010, Ronald Fedkiw, Eran Guendelman, Geoffrey Irving, Michael Lentine, Andrew Selle, Jerry Talton.
// This file is part of PhysBAM whose distribution is governed by the license contained in the accompanying file PHYSBAM_COPYRIGHT.txt.
//#####################################################################
// Class LAPLACE
//#####################################################################
#ifndef __LAPLACE__
#define __LAPLACE__

#include <Core/Arrays/ARRAY.h>
#include <Core/Arrays/ARRAYS_ND.h>
#include <Core/Grids/GRID.h>
#include <Core/Krylov_Solvers/PCG_SPARSE.h>
#include <Core/Math_Tools/minmax.h>
namespace PhysBAM{

template<class T> class PCG_SPARSE;
template<class T> class SPARSE_MATRIX_FLAT_MXN;

template<class TV>
class LAPLACE
{
    typedef typename TV::SCALAR T;
    typedef VECTOR<int,TV::m> IV;
    typedef ARRAY<int,IV> T_ARRAYS_INT;
public:
    typedef TV VECTOR_T;

    T tolerance;
    bool solve_neumann_regions;
    int number_of_regions=0;
    T relative_tolerance;
    T absolute_tolerance;

    GRID<TV> grid;
    ARRAY<T,IV>& u;
    ARRAY<T,IV> f; // f will be modified and reused as b in Ax=b for PCG
    PCG_SPARSE<T> pcg;
    T_ARRAYS_INT filled_region_colors;
    ARRAY<int> filled_region_touches_dirichlet;
    ARRAY<bool,FACE_INDEX<TV::m> > psi_N;
    const ARRAY<T,FACE_INDEX<TV::m> >* psi_N_values=0;
    ARRAY<T,FACE_INDEX<TV::m> > psi_R;
    ARRAY<bool,IV> psi_D;
    VECTOR<bool,TV::m> periodic_boundary;
    bool enforce_compatibility;
    bool solve_single_cell_neumann_regions;
    bool use_psi_R;
public:

    LAPLACE(const GRID<TV>& grid_input,ARRAY<T,IV>& u_input,const bool initialize_grid,const bool enforce_compatibility_input);
    virtual ~LAPLACE();

    LAPLACE(const LAPLACE&) = delete;
    void operator=(const LAPLACE&) = delete;

    void Set_Relative_Tolerance(const T relative_tolerance_input=1e-7)
    {relative_tolerance=relative_tolerance_input;}

    void Set_Absolute_Tolerance(const T absolute_tolerance_input=1e-14)
    {absolute_tolerance=absolute_tolerance_input;}

    void Solve_Neumann_Regions(const bool solve_neumann_regions_input=true)
    {solve_neumann_regions=solve_neumann_regions_input;}

    virtual void Compute_beta_And_Add_Jumps_To_b(const T dt,const T time)
    {}

    void Find_Tolerance(const ARRAY<T>& b)
    {tolerance=max(relative_tolerance*b.Max_Abs(),absolute_tolerance);}

    bool All_Cell_Faces_Neumann(const IV& cell_index) const
    {for(int axis=0;axis<TV::m;axis++)if(!psi_N.Component(axis)(cell_index) || !psi_N.Component(axis)(cell_index+IV::Axis_Vector(axis))) return false;
    return true;}

    bool Any_Cell_Faces_Neumann(const IV& cell_index) const
    {for(int axis=0;axis<TV::m;axis++)if(psi_N.Component(axis)(cell_index) || psi_N.Component(axis)(cell_index+IV::Axis_Vector(axis))) return true;
    return false;}

    bool Any_Neighbor_Dirichlet(const IV& cell_index) const
    {for(int axis=0;axis<TV::m;axis++){IV offset=IV::Axis_Vector(axis);
        if((!psi_N.Component(axis)(cell_index) && psi_D(cell_index-offset)) || (!psi_N.Component(axis)(cell_index+offset) && psi_D(cell_index+offset))) return true;}
    return false;}

//#####################################################################
    void Set_Neumann_Outer_Boundaries();
    virtual void Initialize_Grid(const GRID<TV>& mac_grid_input);
    virtual void Solve(const T time=0,const bool solution_regions_already_computed=false);
    virtual void Find_A(ARRAY<SPARSE_MATRIX_FLAT_MXN<T> >& A_array,ARRAY<ARRAY<T> >& b_array,const ARRAY<int,VECTOR<int,1> >& filled_region_cell_count,T_ARRAYS_INT& cell_index_to_matrix_index);
    virtual void Find_A_Part_One(T_ARRAYS_INT& cell_index_to_matrix_index,ARRAY<ARRAY<int> >& row_counts);
    virtual void Find_A_Part_Two(ARRAY<SPARSE_MATRIX_FLAT_MXN<T> >& A_array,ARRAY<ARRAY<T> >& b_array,T_ARRAYS_INT& cell_index_to_matrix_index);
    void Find_Solution_Regions();
    void Compute_Matrix_Indices(ARRAY<int,VECTOR<int,1> >& filled_region_cell_count,ARRAY<ARRAY<IV> >& matrix_index_to_cell_index_array,T_ARRAYS_INT& cell_index_to_matrix_index);
    void Compute_Matrix_Indices(const RANGE<IV>& domain,ARRAY<int,VECTOR<int,1> >& filled_region_cell_count,ARRAY<ARRAY<IV> >& matrix_index_to_cell_index_array,T_ARRAYS_INT& cell_index_to_matrix_index);
    void Solve_Subregion(ARRAY<IV>& matrix_index_to_cell_index,SPARSE_MATRIX_FLAT_MXN<T>& A,ARRAY<T>& b,const int color=0);
    void Use_Psi_R();
//#####################################################################
};
}
#endif
