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