function [Au, Av, rhsu, rhsv, AuBcs, AvBcs] = FixDiffuvDbcs(Au,Av,bcsDu,bcsDv,...
    rhsu,rhsv,Nodeu, Nodev, dx,dy,Nbcs,AuBcs,AvBcs)
% function [Au Av rhsu rhsv AuBcs AvBcs] = FixDiffuvDbcs(Au,Av,bcsDu,bcsDv,...
%     rhsu,rhsv,Nodeu, Nodev, dx,dy,Nbc,AuBcs,AvBcs)
% This function fixes the A matrix and RHS vector for the discretized 2D
% Laplace equation using Finite Volume method on a staggerd uv grid.
%
% INPUTS:
%       Au:     Original A matrix for u-velocity 
%	    Av:  	Original A matrix for v-velocity
%       bcsDu:  Vector containing the values of dirichlet boundary conditions for u
%       bcsDv:  Vector containing the values of dirichlet boundary conditions for v
%       rhsu:   The Right-hand-side vector from BC contributions for u
%       rhsv:   The Right-hand-side vector from BC contributions for v
%       Nodeu:  U-velocity node-numbering matrix (see NodePad.m)
%       Nodev:  V-velocity node-numbering matrix (see NodePad.m)
%       dx:     delta x, grid spacing in x direction
%       dy:     delta y, grid spacing in y direction
%	    Nbc:    Number of boundary conditions
%	    AuBcs:  The matrix used to create the Right-hand-side vector for u-velocity
%	    AvBcs:  The matrix used to create the Right-hand-side vector for v-velocity
% OUTPUTS:
%       Au:     Corrected A matrix for u-velocity
%	    Av:  	Corrected A matrix for v-velocity
%       rhsu:   The corrected Right-hand-side vector from BC contributions for u
%       rhsv:   The corrected Right-hand-side vector from BC contributions for v
%	    AuBcs:  The corrected matrix used to create the Right-hand-side vector for u
%	    AvBcs:  The corrected matrix used to create the Right-hand-side vector for v
%
% Authors:   Matt Ueckermann, Pierre Lermusiaux, Pat Haley for MIT course 2.29

if nargin<12, AuBcs=zeros(length(Au),length(bcsDu)); end;
if nargin<13, AvBcs=zeros(length(Av),length(bcsDv)); end;

%% Process u-velocity

[Ny, Nx] = size(Nodeu);
NbcD = length(bcsDu);
ids = find(Nodeu <= NbcD);

%Find Ids of interior cells with a northern boundary
n = ids + 1;
n(n > Ny*Nx) = [];
n(Nodeu(n)<=Nbcs) = [];

%Find Ids of interior cells with a southern boundary
s = ids - 1;
s(s < 1) = [];
s((Nodeu(s)<=Nbcs)) = [];

% Constructing the Inner part of the 'A' Matrix for north boundaries
Tayold = taylor([-dy, dy], 2);
Taynew = taylor([-2*dy, -dy, dy/2], 2);
for j = 1 : length(n)
    %Start with the middle node, subtract out old value, add new value
    Rowij = Nodeu(n(j)) - Nbcs;
    Au(Rowij, Rowij) =Au(Rowij, Rowij) + Tayold(2) - Taynew(2);
    %Now fix the south part in the computational cell
    colno = Nodeu(n(j)+1) - Nbcs;
    Au(Rowij, colno) = Au(Rowij, colno) + Tayold(3) - Taynew(3);
    colno = Nodeu(n(j)+2) - Nbcs;
    Au(Rowij, colno) = Au(Rowij, colno) - Taynew(4);
    rhsu(Rowij) = rhsu(Rowij) + (Taynew(1) - Tayold(1))*bcsDu(Nodeu(n(j)-1));
    AuBcs(Rowij,Nodeu(n(j)-1)) = AuBcs(Rowij,Nodeu(n(j)-1)) + (Taynew(1) - Tayold(1));
end;

% Constructing the Inner part of the 'A' Matrix for south boundaries
Tayold = taylor([-dy, dy], 2);
Taynew = taylor([-dy/2, dy, dy*2], 2);
for j = 1 : length(s)
    %Start with the middle node, subtract out old value, add new value
    Rowij = Nodeu(s(j)) - Nbcs;
    Au(Rowij, Rowij) = Au(Rowij, Rowij) + Tayold(2) - Taynew(3);
    %Now fix the north part in the computational cell
    colno = Nodeu(s(j)-1) - Nbcs;
    Au(Rowij, colno) = Au(Rowij, colno) + Tayold(1) - Taynew(2);
    colno = Nodeu(s(j)-2) - Nbcs;
    Au(Rowij, colno) = Au(Rowij, colno) - Taynew(1);
    rhsu(Rowij) = rhsu(Rowij) + (Taynew(4) - Tayold(3))*bcsDu(Nodeu(s(j)+1));
    AuBcs(Rowij,Nodeu(s(j)+1)) = AuBcs(Rowij,Nodeu(s(j)+1)) + (Taynew(4) - Tayold(3));
end;


%% Process v-velocity

[Ny, Nx] = size(Nodev);
NbcD = length(bcsDv);
ids = find(Nodev <= NbcD);

%Find Ids of interior cells with a western boundary
w = ids + Ny;
w(w > Ny*Nx) = [];
w(Nodev(w)<=Nbcs) = [];

%Find Ids of interior cells with a eastern boundary
e = ids - Ny;
e(e < 1) = [];
e(Nodev(e)<=Nbcs) = [];

% Constructing the Inner part of the 'A' Matrix for west boundaries
Tayold = taylor([-dx, dx], 2);
Taynew = taylor([-dx/2, dx, 2*dx], 2);
for j = 1 : length(w)
    %Start with the middle node, subtract out old value, add new value
    Rowij = Nodev(w(j)) - Nbcs;
    Av(Rowij, Rowij ) =Av(Rowij, Rowij) + Tayold(2) - Taynew(3);
    %Now fix the east part in the computational cell
    colno = Nodev(w(j)+Ny) - Nbcs;
    Av(Rowij, colno) = Av(Rowij, colno) + Tayold(1) - Taynew(2);
    colno = Nodev(w(j)+2*Ny) - Nbcs;
    Av(Rowij, colno) = Av(Rowij, colno) - Taynew(1);
    rhsv(Rowij) = rhsv(Rowij) + (Taynew(4) - Tayold(3)) * bcsDv(Nodev(w(j) - Ny));
    AvBcs(Rowij, Nodev(w(j)-Ny)) = AvBcs(Rowij,Nodev(w(j)-Ny)) + (Taynew(4) - Tayold(3));
end;

% Constructing the Inner part of the 'A' Matrix for east boundaries
Tayold = taylor([-dx, dx],2);
Taynew = taylor([-2*dx, -dx, dx/2], 2);
for j = 1 : length(e)
    %Start with the middle node, subtract out old value, add new value
    Rowij = Nodev(e(j)) - Nbcs; 
    Av(Rowij, Rowij ) =Av(Rowij, Rowij ) + Tayold(2) - Taynew(2);
    %Now fix the west part in the computational cell
    colno = Nodev(e(j)-Ny) - Nbcs;  
    Av(Rowij, colno) = Av(Rowij, colno) + Tayold(3) - Taynew(3);
    colno = Nodev(e(j)-2*Ny) - Nbcs;
    Av(Rowij, colno) = Av(Rowij, colno) - Taynew(4);
    rhsv(Rowij) = rhsv(Rowij) + (Taynew(1) - Tayold(1)) * bcsDv(Nodev(e(j)+Ny));
    AvBcs(Rowij,Nodev(e(j)+Ny)) = AvBcs(Rowij,Nodev(e(j)+Ny)) + (Taynew(1) - Tayold(1));
end;
