function [Dfs, Dfs_bc, Dfs_bc_rhs] = Diffusion_Operator(stagr, Node, bcsD, bcsO, bcsN,...
    dx, dy, nu_x, nu_y, Rm_Singularity)
% 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);

NbcsD = length(bcsD);
NbcsO = length(bcsO);
NbcsN = length(bcsN);
Nbcs = NbcsD + NbcsN; % Number of all boundary conditions
bcD_ID = 1 : NbcsD;
bcN_ID = (NbcsD + 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_Nbc = A(:,bcN_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_Nbc = Adj(:,bcN_ID);
Adj = Adj(:, Int_ID);


%% Set up the Dirichlet, Neumann and open boundary conditions

if ~isempty(bcsD)
    
    % Dirichlet BC along x-direction
    if (stagr == 'P')||(stagr == 'v')
        Posi = (Adj_Dbc == 2)|(Adj_Dbc == 3);
        Cell_ID = find(sum(Posi, 2));
        % Modify the boundary cell coeffcients
        A_Dbc(Posi) = A_Dbc(Posi) + (-(-1) + (-16/5))*Cx;
        % Modify the center cell coefficients
        Temp_A = (Adj == 1);
        Temp_A = Temp_A(Cell_ID, :);
        A(Cell_ID,:) = A(Cell_ID,:) + (-2 + 5)*Cx*Temp_A;
        % Modify the inward interior cell coefficients
        Temp_A = (Adj == 2)|(Adj == 3);
        Temp_A = Temp_A(Cell_ID, :);
        A(Cell_ID,:) = A(Cell_ID,:) + (-(-1) + (-2))*Cx*Temp_A;
        % Modify the far inward interior cell coefficients
        Temp_A = (Adj == 6)|(Adj == 7);
        Temp_A = Temp_A(Cell_ID, :);
        A(Cell_ID,:) = A(Cell_ID,:) + (-0 + 1/5)*Cx*Temp_A;
    end;
    
    % Dirichlet BC along y-direction
    if (stagr == 'P')||(stagr == 'u')
        Posi = (Adj_Dbc == 4)|(Adj_Dbc == 5);
        Cell_ID = find(sum(Posi, 2));
        % Modify the boundary cell coeffcients
        A_Dbc(Posi) = A_Dbc(Posi) + (-(-1) + (-16/5))*Cy;
        % Modify the center cell coefficients
        Temp_A = (Adj == 1);
        Temp_A = Temp_A(Cell_ID, :);
        A(Cell_ID,:) = A(Cell_ID,:) + (-2 + 5)*Cy*Temp_A;
        % Modify the inward interior cell coefficients
        Temp_A = (Adj == 4)|(Adj == 5);
        Temp_A = Temp_A(Cell_ID, :);
        A(Cell_ID,:) = A(Cell_ID,:) + (-(-1) + (-2))*Cy*Temp_A;
        % Modify the far inward interior cell coefficients
        Temp_A = (Adj == 8)|(Adj == 9);
        Temp_A = Temp_A(Cell_ID, :);
        A(Cell_ID,:) = A(Cell_ID,:) + (-0 + 1/5)*Cy*Temp_A;
    end;
end;

if ~isempty(bcsN)
    A(Adj == 1) = A(Adj == 1) + sum(A_Nbc, 2);
    
    % Neumann BC along x-direction
    Posi = (Adj_Nbc == 2)|(Adj_Nbc == 3);
    if (stagr == 'P')||(stagr == 'v')
        % Modify the boundary cell coeffcients
        A_Nbc(Posi) = A_Nbc(Posi) * dx;
    else
        Cell_ID = find(sum(Posi, 2));
        % Modify the boundary cell coeffcients
        A_Nbc(Posi) = A_Nbc(Posi) + (-(-1)*Cx + (-6/11)*Cx*dx);
        % Modify the center cell coefficients
        Temp_A = (Adj == 1);
        Temp_A = Temp_A(Cell_ID, :);
        A(Cell_ID,:) = A(Cell_ID,:) + (-(2-1) + 4/11)*Cx*Temp_A;
        % Modify the inward interior cell coefficients
        Temp_A = (Adj == 2)|(Adj == 3);
        Temp_A = Temp_A(Cell_ID, :);
        A(Cell_ID,:) = A(Cell_ID,:) + (-(-1) + (-2/11))*Cx*Temp_A;
        % Modify the far inward interior cell coefficients
        Temp_A = (Adj == 6)|(Adj == 7);
        Temp_A = Temp_A(Cell_ID, :);
        A(Cell_ID,:) = A(Cell_ID,:) + (-0 + (-2/11))*Cx*Temp_A;
    end;
    
    % Neumann BC along y-direction
    Posi = (Adj_Nbc == 4)|(Adj_Nbc == 5);
    if (stagr == 'P')||(stagr == 'u')
        % Modify the boundary cell coeffcients
        A_Nbc(Posi) = A_Nbc(Posi) * dy;
    else
        Cell_ID = find(sum(Posi, 2));
        % Modify the boundary cell coeffcients
        A_Nbc(Posi) = A_Nbc(Posi) + (-(-1)*Cy + (-6/11)*Cy*dy);
        % Modify the center cell coefficients
        Temp_A = (Adj == 1);
        Temp_A = Temp_A(Cell_ID, :);
        A(Cell_ID,:) = A(Cell_ID,:) + (-(2-1) + 4/11)*Cy*Temp_A;
        % Modify the inward interior cell coefficients
        Temp_A = (Adj == 4)|(Adj == 5);
        Temp_A = Temp_A(Cell_ID, :);
        A(Cell_ID,:) = A(Cell_ID,:) + (-(-1) + (-2/11))*Cy*Temp_A;
        % Modify the far inward interior cell coefficients
        Temp_A = (Adj == 8)|(Adj == 9);
        Temp_A = Temp_A(Cell_ID, :);
        A(Cell_ID,:) = A(Cell_ID,:) + (-0 + (-2/11))*Cy*Temp_A;
    end;
end;

if (stagr == 'P')&&(~isempty(bcsO))
    
    % Modify the boundary cell coeffcients
    A_Dbc(:, bcsO) = 0;
    
    Adj_Obc = Adj_Dbc(:, bcsO);
    
    % Open BC along x-direction
    Posi = (Adj_Obc == 2)|(Adj_Obc == 3);
    Cell_ID = find(sum(Posi, 2));
    % Modify the center cell coefficients
    Temp_A = (Adj == 1);
    Temp_A = Temp_A(Cell_ID, :);
    A(Cell_ID,:) = A(Cell_ID,:) + (-5 + (-1/3))*Cx*Temp_A;
    % Modify the inward interior cell coefficients
    Temp_A = (Adj == 2)|(Adj == 3);
    Temp_A = Temp_A(Cell_ID, :);
    A(Cell_ID,:) = A(Cell_ID,:) + (-(-2) + 2/3)*Cx*Temp_A;
    % Modify the far inward interior cell coefficients
    Temp_A = (Adj == 6)|(Adj == 7);
    Temp_A = Temp_A(Cell_ID, :);
    A(Cell_ID,:) = A(Cell_ID,:) + (-1/5 + (-1/3))*Cx*Temp_A;
    
    % Open BC along y-direction
    Posi = (Adj_Obc == 4)|(Adj_Obc == 5);
    Cell_ID = find(sum(Posi, 2));
    % Modify the center cell coefficients
    Temp_A = (Adj == 1);
    Temp_A = Temp_A(Cell_ID, :);
    A(Cell_ID,:) = A(Cell_ID,:) + (-5 + (-1/3))*Cy*Temp_A;
    % Modify the inward interior cell coefficients
    Temp_A = (Adj == 4)|(Adj == 5);
    Temp_A = Temp_A(Cell_ID, :);
    A(Cell_ID,:) = A(Cell_ID,:) + (-(-2) + 2/3)*Cy*Temp_A;
    % Modify the far inward interior cell coefficients
    Temp_A = (Adj == 8)|(Adj == 9);
    Temp_A = Temp_A(Cell_ID, :);
    A(Cell_ID,:) = A(Cell_ID,:) + (-1/5 + (-1/3))*Cy*Temp_A;
end;

% Remove the singularity if no Dirichlet BCs are specified.
if (nargin == 10)&&(Rm_Singularity)&&(NbcsD == NbcsO)&&(stagr == 'P')
    A(end, :) = 0; A(end, end) = 1;
    A_Dbc(end, :) = 0; A_Nbc(end, :) = 0;
end;


%% Finally, distribute entries of A to outputs

Dfs = A;
Dfs_bc = [A_Dbc, A_Nbc];
Dfs_bc_rhs = -Dfs_bc * [bcsD; bcsN];
