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];