//#####################################################################
// Copyright 2024 Jason Goulding
// This file is part of PhysBAM whose distribution is governed by the license contained in the accompanying file PHYSBAM_COPYRIGHT.txt.
//#####################################################################
#ifndef __COMPUTE_LAPLACE_MATRIX__
#define __COMPUTE_LAPLACE_MATRIX__
#include <Core/Arrays/ARRAY.h>
#include <Core/Arrays/ARRAYS_ND.h>
#include <Core/Arrays/ARRAYS_ND_BASE.h>
#include <Core/Arrays/FACE_ARRAYS.h>
#include <Core/Grids/GRID.h>
#include <Core/Matrices/SPARSE_MATRIX_FLAT_MXN.h>
namespace PhysBAM{

template<class TV>
struct BOUNDARY_CONDITIONS
{
    typedef typename TV::SCALAR T;
    typedef VECTOR<int,TV::m> IV;
    enum BC_TYPE {BC_DIRICHLET,BC_NEUMANN,BC_PERIODIC};
    const ARRAY<bool,FACE_INDEX<TV::m>>* psi_N=nullptr;
    const ARRAY<T,FACE_INDEX<TV::m>>* psi_N_values=nullptr;
    const ARRAYS_ND_BASE<bool,IV>* psi_D=nullptr;
    const ARRAYS_ND_BASE<T,IV>* psi_D_values=nullptr;
    VECTOR<BC_TYPE,2*TV::m> bc;
    GRID<TV> grid;
};

template<class TV>
struct LAPLACE_MATRIX
{
    typedef typename TV::SCALAR T;
    typedef VECTOR<int,TV::m> IV;
    ARRAY<SPARSE_MATRIX_FLAT_MXN<T>> A_array;
    ARRAY<int,IV> filled_region_colors;
    ARRAY<ARRAY<IV>> matrix_index_to_cell_index_array;
    ARRAY<int,IV> cell_index_to_matrix_index;
    int number_of_regions=0;
    ARRAY<bool> enforce_compatibility;
};

template<class TV,class T>
void Compute_Laplace_Matrix(LAPLACE_MATRIX<TV>& lm,ARRAY<ARRAY<T>>& rhs,const ARRAYS_ND_BASE<T,VECTOR<int,TV::m>>& f,const BOUNDARY_CONDITIONS<TV>& bc,const ARRAY<T,FACE_INDEX<TV::m>>* beta=nullptr);

template<class TV,class T,class IV>
void Compute_Laplace_Matrix(SPARSE_MATRIX_FLAT_MXN<T>& A,ARRAY<T>& rhs,const ARRAYS_ND_BASE<T,VECTOR<int,TV::m>>& f,const BOUNDARY_CONDITIONS<TV>& bc,ARRAY<IV>& matrix_index_to_cell_index_array,ARRAY<int,IV>& cell_index_to_matrix_index,bool& enforce_compatibility,const ARRAY<T,FACE_INDEX<TV::m>>* beta=nullptr);

}
#endif
