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: Aug. 2014
% 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)];
% Remove any entries that correspond to boundaries (i.e. <= Nbcs)
J(I <= 0) = []; I(I <= 0) = []; I = I(:); J = J(:);

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

% Break the A matrix into three parts corresponding to Dirichlet boundaries,
% Neumann boundaries 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);

I = Node(idsy, idsx) - Nbcs;
I_Nega = (I <= 0); I(I_Nega) = []; I = I(:);

% Mark the center cells.
J = Node(idsy, idsx);
J(I_Nega) = [];  J = J(:);
Adj_C = sparse(I, J, TRUE, Num_Cell, Num_Cell + Nbcs, Num_Cell);

% Mark the west cells.
J = Node(idsy, idsx - 1);
J(I_Nega) = [];  J = J(:);
Adj_W = sparse(I, J, TRUE, Num_Cell, Num_Cell + Nbcs, Num_Cell);

% Mark the east cells.
J = Node(idsy, idsx + 1);
J(I_Nega) = [];  J = J(:);
Adj_E = sparse(I, J, TRUE, Num_Cell, Num_Cell + Nbcs, Num_Cell);

% Mark the south cells.
J = Node(idsy + 1, idsx);
J(I_Nega) = [];  J = J(:);
Adj_S = sparse(I, J, TRUE, Num_Cell, Num_Cell + Nbcs, Num_Cell);

% Mark the north cells.
J = Node(idsy - 1, idsx);
J(I_Nega) = [];  J = J(:);
Adj_N = sparse(I, J, TRUE, Num_Cell, Num_Cell + Nbcs, 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_FW = sparse(I, J, true(length(I), 1), Num_Cell, Num_Cell + Nbcs, Num_Cell);

% 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_FE = sparse(I, J, true(length(I), 1), Num_Cell, Num_Cell + Nbcs, Num_Cell);

% 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_FS = sparse(I, J, true(length(I), 1), Num_Cell, Num_Cell + Nbcs, Num_Cell);

% 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_FN = sparse(I, J, true(length(I), 1), Num_Cell, Num_Cell + Nbcs, Num_Cell);

% Break the Adj matrices into three parts corresponding to Dirichlet boundaries,
% Neumann boundaries and interior cells

Adj_C = Adj_C(:, Int_ID);
Adj_Dbc_W = Adj_W(:,bcD_ID); Adj_Nbc_W = Adj_W(:,bcN_ID); Adj_W = Adj_W(:, Int_ID);
Adj_Dbc_E = Adj_E(:,bcD_ID); Adj_Nbc_E = Adj_E(:,bcN_ID); Adj_E = Adj_E(:, Int_ID);
Adj_Dbc_S = Adj_S(:,bcD_ID); Adj_Nbc_S = Adj_S(:,bcN_ID); Adj_S = Adj_S(:, Int_ID);
Adj_Dbc_N = Adj_N(:,bcD_ID); Adj_Nbc_N = Adj_N(:,bcN_ID); Adj_N = Adj_N(:, Int_ID);
Adj_FW = Adj_FW(:, Int_ID);
Adj_FE = Adj_FE(:, Int_ID); 
Adj_FS = Adj_FS(:, Int_ID); 
Adj_FN = Adj_FN(:, 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_W | Adj_Dbc_E);
        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
        A(Cell_ID,:) = A(Cell_ID,:) + (-2 + 5)*Cx*Adj_C(Cell_ID, :);
        % Modify the inward interior cell coefficients
        Temp_A = (Adj_W | Adj_E);
        A(Cell_ID,:) = A(Cell_ID,:) + (-(-1) + (-2))*Cx*Temp_A(Cell_ID, :);
        % Modify the far inward interior cell coefficients
        Temp_A = (Adj_FW | Adj_FE);
        A(Cell_ID,:) = A(Cell_ID,:) + (-0 + 1/5)*Cx*Temp_A(Cell_ID, :);
    end;
    
    % Dirichlet BC along y-direction
    if (stagr == 'P')||(stagr == 'u')
        Posi = (Adj_Dbc_S | Adj_Dbc_N);
        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
        A(Cell_ID,:) = A(Cell_ID,:) + (-2 + 5)*Cy*Adj_C(Cell_ID, :);
        % Modify the inward interior cell coefficients
        Temp_A = (Adj_S | Adj_N);
        A(Cell_ID,:) = A(Cell_ID,:) + (-(-1) + (-2))*Cy*Temp_A(Cell_ID, :);
        % Modify the far inward interior cell coefficients
        Temp_A = (Adj_FS | Adj_FN);
        A(Cell_ID,:) = A(Cell_ID,:) + (-0 + 1/5)*Cy*Temp_A(Cell_ID, :);
    end;
end;

if ~isempty(bcsN)
    
    % Neumann BC along x-direction
    Posi = (Adj_Nbc_W | Adj_Nbc_E);
    if (stagr == 'P')||(stagr == 'v')
        Cell_ID = find(sum(Posi, 2));
%         % Modify the boundary cell coeffcients
%         A_Nbc(Posi) = A_Nbc(Posi) + (-(-1) + (-24/23)*dx)*Cx;
%         % Modify the center cell coefficients
%         A(Cell_ID,:) = A(Cell_ID,:) + (-2 + 25/23)*Cx*Adj_C(Cell_ID, :);
%         % Modify the inward interior cell coefficients
%         Temp_A = (Adj_W | Adj_E);
%         A(Cell_ID,:) = A(Cell_ID,:) + (-(-1) + (-26/23))*Cx*Temp_A(Cell_ID, :);
%         % Modify the far inward interior cell coefficients
%         Temp_A = (Adj_FW | Adj_FE);
%         A(Cell_ID,:) = A(Cell_ID,:) + (-0 + 1/23)*Cx*Temp_A(Cell_ID, :);
        % Modify the boundary cell coeffcients
        A_Nbc(Posi) = A_Nbc(Posi) + (-(-1) + (-1)*dx)*Cx;
        % Modify the center cell coefficients
        A(Cell_ID,:) = A(Cell_ID,:) + (-2 + 1)*Cx*Adj_C(Cell_ID, :);
    else
        % Remove all interior cells that touch masks and only execute this
        % on cells that correspond to open boundaries
        
%         Diri_in_FreeSlip = true(length(Posi),1);
%         Ends_Cells = [Node(2:end-1,2); Node(2:end-1,end-1)] - Nbcs;
%         Ends_Cells(Ends_Cells<=0) = [];
%         Diri_in_FreeSlip(Ends_Cells) = false;
%         Posi(Diri_in_FreeSlip,:) = false;
        
        Cell_ID = find(sum(Posi, 2));
        % Modify the boundary cell coeffcients
        A_Nbc(Posi) = A_Nbc(Posi) + (-(-1) + (-6/11)*dx)*Cx;
        % Modify the center cell coefficients
        A(Cell_ID,:) = A(Cell_ID,:) + (-2 + 4/11)*Cx*Adj_C(Cell_ID, :);
        % Modify the inward interior cell coefficients
        Temp_A = (Adj_W | Adj_E);
        A(Cell_ID,:) = A(Cell_ID,:) + (-(-1) + (-2/11))*Cx*Temp_A(Cell_ID, :);
        % Modify the far inward interior cell coefficients
        Temp_A = (Adj_FW | Adj_FE);
        A(Cell_ID,:) = A(Cell_ID,:) + (-0 + (-2/11))*Cx*Temp_A(Cell_ID, :);
    end;
    
    % Neumann BC along y-direction
    Posi = (Adj_Nbc_S | Adj_Nbc_N);
    if (stagr == 'P')||(stagr == 'u')
        Cell_ID = find(sum(Posi, 2));
%         % Modify the boundary cell coeffcients
%         A_Nbc(Posi) = A_Nbc(Posi) + (-(-1) + (-24/23)*dy)*Cy;
%         % Modify the center cell coefficients
%         A(Cell_ID,:) = A(Cell_ID,:) + (-2 + 25/23)*Cy*Adj_C(Cell_ID, :);
%         % Modify the inward interior cell coefficients
%         Temp_A = (Adj_S | Adj_N);
%         A(Cell_ID,:) = A(Cell_ID,:) + (-(-1) + (-26/23))*Cy*Temp_A(Cell_ID, :);
%         % Modify the far inward interior cell coefficients
%         Temp_A = (Adj_FS | Adj_FN);
%         A(Cell_ID,:) = A(Cell_ID,:) + (-0 + 1/23)*Cy*Temp_A(Cell_ID, :);
        % Modify the boundary cell coeffcients
        A_Nbc(Posi) = A_Nbc(Posi) + (-(-1) + (-1)*dy)*Cy;
        % Modify the center cell coefficients
        A(Cell_ID,:) = A(Cell_ID,:) + (-2 + 1)*Cy*Adj_C(Cell_ID, :);
    else
        % Remove all interior cells that touch masks and only execute this
        % on cells that correspond to open boundaries
        
%         Diri_in_FreeSlip = true(length(Posi),1);
%         Ends_Cells = [Node(2,2:end-1)'; Node(end-1,2:end-1)'] - Nbcs;
%         Ends_Cells(Ends_Cells<=0) = [];
%         Diri_in_FreeSlip(Ends_Cells) = false;
%         Posi(Diri_in_FreeSlip,:) = false;
        
        Cell_ID = find(sum(Posi, 2));
        % Modify the boundary cell coeffcients
        A_Nbc(Posi) = A_Nbc(Posi) + (-(-1) + (-6/11)*dy)*Cy;
        % Modify the center cell coefficients
        A(Cell_ID,:) = A(Cell_ID,:) + (-2 + 4/11)*Cy*Adj_C(Cell_ID, :);
        % Modify the inward interior cell coefficients
        Temp_A = (Adj_S | Adj_N);
        A(Cell_ID,:) = A(Cell_ID,:) + (-(-1) + (-2/11))*Cy*Temp_A(Cell_ID, :);
        % Modify the far inward interior cell coefficients
        Temp_A = (Adj_FS | Adj_FN);
        A(Cell_ID,:) = A(Cell_ID,:) + (-0 + (-2/11))*Cy*Temp_A(Cell_ID, :);
    end;
end;

if (stagr == 'P')&&(~isempty(bcsO))
    
    % Modify the boundary cell coeffcients
    A_Dbc(:, bcsO) = 0;
    
    % Open BC along x-direction
    Posi = (Adj_Dbc_W(:, bcsO) | Adj_Dbc_E(:, bcsO));
    Cell_ID = find(sum(Posi, 2));
    % Modify the center cell coefficients
    A(Cell_ID,:) = A(Cell_ID,:) + (-5 + (-1/3))*Cx*Adj_C(Cell_ID, :);
    % Modify the inward interior cell coefficients
    Temp_A = (Adj_W | Adj_E);
    A(Cell_ID,:) = A(Cell_ID,:) + (-(-2) + 2/3)*Cx*Temp_A(Cell_ID, :);
    % Modify the far inward interior cell coefficients
    Temp_A = (Adj_FW | Adj_FE);
    A(Cell_ID,:) = A(Cell_ID,:) + (-1/5 + (-1/3))*Cx*Temp_A(Cell_ID, :);
    
    % Open BC along y-direction
    Posi = (Adj_Dbc_S(:, bcsO) | Adj_Dbc_N(:, bcsO));
    Cell_ID = find(sum(Posi, 2));
    % Modify the center cell coefficients
    A(Cell_ID,:) = A(Cell_ID,:) + (-5 + (-1/3))*Cy*Adj_C(Cell_ID, :);
    % Modify the inward interior cell coefficients
    Temp_A = (Adj_S | Adj_N);
    A(Cell_ID,:) = A(Cell_ID,:) + (-(-2) + 2/3)*Cy*Temp_A(Cell_ID, :);
    % Modify the far inward interior cell coefficients
    Temp_A = (Adj_FS | Adj_FN);
    A(Cell_ID,:) = A(Cell_ID,:) + (-1/5 + (-1/3))*Cy*Temp_A(Cell_ID, :);
end;

% Remove the singularity if no Dirichlet BCs are specified.
if (nargin == 10)&&(stagr == 'P')
    switch Rm_Singularity
        case 0
            display('No singularity removal is performed.');
            display(sprintf(...
                'Condition number of the linear system for pressure correction: %0.2e\n',...
                condest(A)));
        case 1
            display('Null space dimension is to be reduced by 1.');
            display(sprintf(...
                'Condition number of the linear system for pressure correction\nBEFORE removing singularity: %0.2e\n',...
                condest(A)));
            A(end, :) = 0; A(end, end) = Cx + Cy;
            A_Dbc(end, :) = 0; A_Nbc(end, :) = 0;
            display(sprintf(...
                'Condition number of the linear system for pressure correction\nAFTER removing singularity: %0.2e\n',...
                condest(A)));
        case 2
            display('Null space dimension is to be reduced by 2.');
            display(sprintf(...
                'Condition number of the linear system for pressure correction\nBEFORE removing singularity: %0.2e\n',...
                condest(A)));
            A(end, :) = 0; A(end, end) = Cx + Cy;
            A_Dbc(end, :) = 0; A_Nbc(end, :) = 0;
            A(1, :) = 0; A(1, 1) = Cx + Cy;
            A_Dbc(1, :) = 0; A_Nbc(1, :) = 0;
            display(sprintf(...
                'Condition number of the linear system for pressure correction\nAFTER removing singularity: %0.2e\n',...
                condest(A)));
    end;
end;


%% Finally, distribute entries of A to outputs

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