//#####################################################################
// Copyright 2006,2022 Geoffrey Irving, Jason Goulding, Mehrnaz Ayazi
// This file is part of PhysBAM whose distribution is governed by the license contained in the accompanying file PHYSBAM_COPYRIGHT.txt.
//#####################################################################
// Class HELMHOLTZ
//#####################################################################
//
// Class to solve
//   au - div(b*grad(u)) = f + D + N
// D is due to Dirichlet boundary conditions.
// N is due to Neumann boundary conditions.
// a is a non-negative variable coefficient for u that potentially is all zero
// and b is a non-negative variable coefficient for grad(u) that is not all zero
//
//#####################################################################
#ifndef __HELMHOLTZ__
#define __HELMHOLTZ__

#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/Log/DEBUG_UTILITIES.h>
#include <Core/Math_Tools/RANGE.h>
#include <Core/Matrices/SPARSE_MATRIX_FLAT_MXN.h>
namespace PhysBAM{

template<class TV> class GRID;
template<class T> class SPARSE_MATRIX_FLAT_MXN;
template<class TV> struct LAPLACE_MATRIX;

template<class T>
struct PCG_OPTIONS
{
    int max_iterations=1000;
    int restart_iterations=200;
    bool show_results=false;
    bool show_residual=false;
    bool ensure_convergence=true;
    bool enforce_compatibility=true;
    T absolute_tolerance=1e-14;
    T relative_tolerance=1e-7;
};

template<class TV>
class HELMHOLTZ
{
    typedef typename TV::SCALAR T;
    typedef VECTOR<int,TV::m> IV;

    T alpha_c=1;
    T beta_c=1;
    const ARRAYS_ND_BASE<T,IV>* alpha=0;
    const ARRAY<T,FACE_INDEX<TV::m> >* beta=0;
    const ARRAY<bool,FACE_INDEX<TV::m> >* psi_N=0;
    const ARRAY<T,FACE_INDEX<TV::m> >* psi_N_values=0; // User provides value in the normal direction
    const ARRAYS_ND_BASE<bool,IV>* psi_D=0; // move psi_D values out of u (double-check)
    const ARRAYS_ND_BASE<T,IV>* psi_D_values=0;
    VECTOR<bool,TV::m> periodic_boundary;
    ARRAY<int,IV> filled_region_colors;
    ARRAY<int> filled_region_touches_dirichlet;
    GRID<TV> grid;
    bool nodes=false; // true means values on nodes, false on cells
    
public:
    bool print_system=false;
    bool diagonal_scaling=false;
    std::string octave_output; 
    
    HELMHOLTZ(const GRID<TV>& grid_input,bool nodes);
    
    void Set_Beta(T b);
    void Set_Beta(const ARRAY<T,FACE_INDEX<TV::m> >& b);

    void Set_Alpha(T a);
    void Set_Alpha(const ARRAYS_ND_BASE<T,IV>& a);

    void Set_Neumann_Boundaries(const ARRAY<bool,FACE_INDEX<TV::m> >& N, const ARRAY<T,FACE_INDEX<TV::m> >* N_values=0);
    void Set_Dirichlet_Boundaries(const ARRAYS_ND_BASE<bool,IV>& D, const ARRAYS_ND_BASE<T,IV>* D_values);
    void Set_Periodic_Boundaries(const VECTOR<bool,TV::m>& p);
    
    bool Solve(ARRAYS_ND_BASE<T,IV>& u,const ARRAYS_ND_BASE<T,IV>& f,RANGE<IV> grid_domain,PCG_OPTIONS<T>& options=PCG_OPTIONS<T>());
    void Find_A(LAPLACE_MATRIX<TV>& lm,const RANGE<IV>& u_size);
    bool Solve_Subregion(LAPLACE_MATRIX<TV>& lm,SPARSE_MATRIX_FLAT_MXN<T>& A,ARRAY<T>& b,ARRAYS_ND_BASE<T,IV>& u,const PCG_OPTIONS<T>& options,const int color=0);
    
//#####################################################################
};
}
#endif
