function [divu] = Div2D(u,v,Nodeu,Nodev,dx,dy,idsP,Neu,Nev)
% function [divu] = Div2D(u,v,Nodeu,Nodev)
% This function takes the divergence of U=[u(Nodeu),v(Nodev)] 
%
% 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 
%       dx:     delta x, grid spacing in x direction
%       dy:     delta y, grid spacing in y direction
%       idsP:   Location of masked Pressure nodes in interior
%       Neu:    Array with two entries giving the ID of the 
%	         lowest and highest numbered Neumann boundary 
%	         conditions for the u-velocity
%       Nev:    Array with two entries giving the ID of the 
%	         lowest and highest numbered Neumann boundary 
%	         conditions for the v-velocity
%       
% OUTPUTS:
%       divu:   du/dx + dv/dy returned at the interior pressure points
%
% Authors:   Matt Ueckermann, Pierre Lermusiaux, Pat Haley for MIT course 2.29


%Find left and right nodes (including boundaries)
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;

%For Neumann boundaries, we need to correct the discrete divergence. 
%First find boundaries that are Neumann
lids1 = logical((Nodeu(2:end-1, 2:end) >= Neu(1)).*...
                (Nodeu(2:end-1, 2:end) <= Neu(2)));
lids2 = logical((Nodeu(2:end-1, 1:end-1) >= Neu(1)).*...
                (Nodeu(2:end-1, 1:end-1) <= Neu(2)));

%Assign the boundary x-derivative
dudx(lids1) = UR(lids1);    
dudx(lids2) = -UL(lids2);    

%Find the top (R) and bottom (L) nodes (including boundaries)
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;

%For Neumann boundaries, we need to correct the discrete divergence. 
%First find boundaries that are Neumann
lids1 = logical((Nodev(1:end-1,2:end-1) >= Nev(1)).*...
                (Nodev(1:end-1,2:end-1) <= Nev(2)));
lids2 = logical((Nodev(2:end,2:end-1) >= Nev(1)).*...
                (Nodev(2:end,2:end-1) <= Nev(2)));
%Assign the boundary y-derivative
dvdy(lids1) = UR(lids1);    
dvdy(lids2) = -UL(lids2);    

%calculate the discrete divergence for interior points
divu = dudx(idsP)+dvdy(idsP);
