function [Dfs, Dfs_bc, BV_Op] = Diffusion_Operator_q...
    (Node, Nbcs, Node2, bcid2type, dx, dy, nu_x, nu_y)
% function [Dfs, Dfs_bc, Dfs_bc_rhs] = Diffusion_Operator(stagr, Node,...
%     bcsD, bcsO, bcsN,dx, dy, nu_x, nu_y, Rm_Singularity)
% This function creates the discrete diffusion operator (matrix), its RHS vector
% and the matrix to generate this RHS vector from boundary conditions.
%
% Warning: Note that the diffusion operator here is meant to be
%                     - nu_x d^2/dx^2 - nu_y d^2/dy^2
% This is not simply a Laplace operator! Two things to keep in mind:
% 1. There is a negative sign.
% 2. The two diffusion coeffcients are included and can be different.
%
% To further clarify the usage, we consider the following Poisson equation about P
%                  (- nu_x d^2/dx^2 - nu_y d^2/dy^2)P = b
% where b is some arbitray forcing.
% With the boundary conditions specified as the input of this Matlab function,
% after the output is obtained, the discrete version of the Poisson equation is
%                       Dfs * P_hat = b_hat + Dfs_bc_rhs
% Solving this linear system gives the numerical solution of P:
%                    P_hat = (Dfs)^(-1) * (b_hat + Dfs_bc_rhs)
%
% INPUTS:
%       stagr:  String indicating which grid to use: 'P', 'u', or 'v'.
%       Node:   Node numbering matrix.
%               NOTE: Node needs to include boundary conditions. In the case for
%               periodic boundaries, the periodicity needs to be included. See NodePad.m
%       bcsD:   Vector containing the values of Dirichlet BCs.
%       bcsO:   Vector containing the IDs of Dirichelet BCs that are actually open BCs.
%       bcsN:   Vector containing the values of Neumann BCs.
%               (eg. bcsN = g means du/dn = g on the boundary)
%       dx:     Mesh spacing in x direction
%       dy:     Mesh sapcing in y direction
%       nu_x:   Diffusion coeffcients in x direction
%       nu_y:   Diffusion coeffcients in y direction
%
%       Rm_Singularity: (Optional) Logical variable indicating whether to remove
%               the singularity of the operator matrix Dfs when no Dirichlet BCs are
%               specified. The singularity is removed by changing the last equation to
%               be one setting the last interior cell equal to 0.
%
% OUTPUTS:
%       Dfs:        The diffusion operator matrix
%       Dfs_bc:     The matrix used to create the right-hand-side vector
%       Dfs_bc_rhs: The right-hand-side vector from BC contributions
%
% See also: NodePad.m
%
% Date: Dec. 2013
% Author: Jing Lin


%% Preliminaries

% Nx: Number of cells in the x-direction (include boundaries)
% Ny: Number of cells in the y-direction (include boundaries)
[Ny, Nx] = size(Node);

bcD_ID = 1 : Nbcs;
Int_ID = (Nbcs + 1) : max(Node(:));
Num_Cell = length(Int_ID); % Number of interior cells

% Define the number clusters appearing in entries of diffusion operator matrix
Cx = nu_x/dx^2;
Cy = nu_y/dy^2;


%% Build a preliminary matrix A of size Num_Cell * (Num_Cell + Nbcs)

%Create the indices for the entries of A
idsx = 2 : (Nx - 1);
idsy = 2 : (Ny - 1);
I = repmat(Node(idsy, idsx) - Nbcs, 1, 5);
J = [Node(idsy, idsx), Node(idsy, idsx - 1), Node(idsy, idsx + 1),...
    Node(idsy + 1, idsx), Node(idsy - 1, idsx)];
J(I <= 0) = []; I(I <= 0) = []; I = I(:); J = J(:);

% Create the A matrix in sparse form
One = ones(Num_Cell, 1);
A = sparse(I, J, kron([2*(Cx+Cy);-Cx;-Cx;-Cy;-Cy], One),...
    Num_Cell, Num_Cell + Nbcs, 5 * Num_Cell);

% Break the A matrix into three parts corresponding to Dirichlet, Neumann boundary
% cells and interior cells
A_Dbc = A(:,bcD_ID);
A = A(:, Int_ID);


%% Construct an adjacency matrix Adj

% Construct an adjacency matrix of the same size as A indicating in which direction
% the adjacent cells corresponding to entries of A lie.
% The mapping follows the rule:
% 1 :    Center cell (i  ,  j)
% 2 :      West cell (i  ,j-1)
% 3 :      East cell (i  ,j+1)
% 4 :     South cell (i+1,  j)
% 5 :     North cell (i-1,  j)
% 6 :  Far west cell (i  ,j-2)
% 7 :  Far east cell (i  ,j+2)
% 8 : Far south cell (i+2,  j)
% 9 : Far north cell (i-2,  j)

idsx2 = 2 : (Nx - 2);
idsy2 = 2 : (Ny - 2);
Adj = sparse(I, J, kron([1;2;3;4;5], One), Num_Cell, Num_Cell + Nbcs, 9 * Num_Cell);

% Mark the far west cells.
I = Node(idsy, idsx2+1) - Nbcs; J = Node(idsy, idsx2-1);
J(I <= 0) = []; I(I <= 0) = []; I = I(:); J = J(:);
Adj = Adj + sparse(I, J, 6*ones(length(I), 1), Num_Cell, Num_Cell + Nbcs);

% Mark the far east cells.
I = Node(idsy, idsx2) - Nbcs; J = Node(idsy, idsx2+2);
J(I <= 0) = []; I(I <= 0) = []; I = I(:); J = J(:);
Adj = Adj + sparse(I, J, 7*ones(length(I), 1), Num_Cell, Num_Cell + Nbcs);

% Mark the far south cells.
I = Node(idsy2, idsx) - Nbcs; J = Node(idsy2+2, idsx);
J(I <= 0) = []; I(I <= 0) = []; I = I(:); J = J(:);
Adj = Adj + sparse(I, J, 8*ones(length(I), 1), Num_Cell, Num_Cell + Nbcs);

% Mark the far north cells.
I = Node(idsy2+1, idsx) - Nbcs; J = Node(idsy2-1, idsx);
J(I <= 0) = []; I(I <= 0) = []; I = I(:); J = J(:);
Adj = Adj + sparse(I, J, 9*ones(length(I), 1), Num_Cell, Num_Cell + Nbcs);

% Break the Adj matrix into three parts corresponding to Dirichlet, Neumann boundary
% cells and interior cells
Adj_Dbc = Adj(:,bcD_ID);
Adj = Adj(:, Int_ID);


%% Set up the Neumann and open boundary conditions

I = zeros(Nbcs*2, 1);
J = zeros(Nbcs*2, 1);
Value = zeros(Nbcs*2, 1);
count = 0;

for i  = 1 : Nbcs
    Adj_Node = find((Adj_Dbc(:,i) >= 2)&(Adj_Dbc(:,i) <= 5), 1);
    if isempty(Adj_Node), continue; end;
    if (bcid2type{Node2(find(Node==i, 1))} == 'N')
        count = count + 1;
        I(count) = i; J(count) = Adj_Node; Value(count) = 1;
    else
        % Find the far adjacent interior cell
        Adj_Node_Far = find(Adj_Dbc(:,i) == (Adj_Dbc(Adj_Node,i)+4), 1);
        count = count + 1;
        I(count) = i; J(count) = Adj_Node; Value(count) = 2;
        count = count + 1;
        I(count) = i; J(count) = Adj_Node_Far; Value(count) = -1;
    end;
end;
I = I(1:count); J = J(1:count); Value = Value(1:count);

BV_Op = sparse(I, J, Value, Nbcs, Num_Cell, count);

%% Finally, distribute entries of A to outputs

Dfs = A;
Dfs_bc = A_Dbc;
