function [dudx, dvdy] = Grad2D_uv(u, v, Nodeu, Nodev, idsP, dx, dy, bcN_ID_u, bcN_ID_v)
% function [dudx, dvdy] = Grad2D_uv(u, v, Nodeu, Nodev, idsP, dx, dy, bcN_ID_u, bcN_ID_v)
% This function computes du/dx and dv/dy in the centers of active P cells.
%
% INPUTS:
%       u:        Vector containing values of u
%       v:        Vector containing values of v
%       Nodeu:    Reference node numbering matrix for u-velocity
%       Nodev:    Reference node numbering matrix for v-velocity
%       idsP:     Location of active pressure cells in interior
%       dx:       Mesh spacing in x direction
%       dy:       Mesh spacing in y direction
%       bcN_ID_u: Array with two entries giving the lowest and highest ID of the
%	              Neumann boundary conditions for the u-velocity
%       bcN_ID_v: Array with two entries giving the lowest and highest ID of the
%	              Neumann boundary conditions for the u-velocity
%
% OUTPUTS:
%       dudx:     du/dx returned at the centers of active pressure cells
%       dvdy:     dv/dy returned at the centers of active pressure cells
%
% Date: Jan. 2014
% Authors: Jing Lin


%% Preliminary step for du/dx

% Find the left and right u cells (including boundaries) for each interior P cell
uR = u(Nodeu(2:end-1, 2:end  ));
uL = u(Nodeu(2:end-1, 1:end-1));

% Take discrete x derivative
dudx = (uR - uL)/dx;


%% Correct the discrete derivatives at Neumann boundaries for du/dx

% Find right and left u cells that are actually Neumann boundaries
bcN_R = logical((Nodeu(2:end-1, 2:end) >= bcN_ID_u(1))...
    .*(Nodeu(2:end-1, 2:end) <= bcN_ID_u(2)));
bcN_L = logical((Nodeu(2:end-1, 1:end-1) >= bcN_ID_u(1))...
    .*(Nodeu(2:end-1, 1:end-1) <= bcN_ID_u(2)));
Cell_LL = false(size(bcN_R)); Cell_LL(:, 1:end-1) = bcN_R(:, 2:end);
Cell_RR = false(size(bcN_L)); Cell_RR(:, 2:end) = bcN_L(:, 1:end-1);

% Assign the boundary x-derivative
dudx(bcN_R) = 2/3*uR(bcN_R) + (uL(bcN_R) - uL(Cell_LL))/(3*dx);
dudx(bcN_L) = -2/3*uL(bcN_L) - (uR(bcN_L) - uR(Cell_RR))/(3*dx);

dudx = dudx(idsP);


%% Preliminary step for dv/dy

% Find top (R) and bottom (L) v cells (including boundaries) for each interior P cell
uR = v(Nodev(1:end-1, 2:end-1));
uL = v(Nodev(2:end  , 2:end-1));

% Take discrete y derivative
dvdy = (uR - uL)/dy;


%% Correct the discrete derivatives at Neumann boundaries for dv/dy

% Find top (R) and bottom (L) v cells that are actually Neumann boundaries
bcN_R = logical((Nodev(1:end-1, 2:end-1) >= bcN_ID_v(1))...
    .*(Nodev(1:end-1, 2:end-1) <= bcN_ID_v(2)));
bcN_L = logical((Nodev(2:end,2:end-1) >= bcN_ID_v(1))...
    .*(Nodev(2:end, 2:end-1) <= bcN_ID_v(2)));
Cell_LL = false(size(bcN_R)); Cell_LL(2:end, :) = bcN_R(1:end-1, :);
Cell_RR = false(size(bcN_L)); Cell_RR(1:end-1, :) = bcN_L(2:end, :);

% Assign the boundary y-derivative
dvdy(bcN_R) = 2/3*uR(bcN_R) + (uL(bcN_R) - uL(Cell_LL))/(3*dy);
dvdy(bcN_L) = -2/3*uL(bcN_L) - (uR(bcN_L) - uR(Cell_RR))/(3*dy);

dvdy = dvdy(idsP);

