function divu = Div2D(u,v)

% compute the 2D divergence of the vectorfield (u,v)

Globals2D;

u_r = D_r * u; u_s = D_s * u;
v_r = D_r * v; v_s = D_s * v;
divu = u_r .* r_x + u_s .* s_x + v_r .* r_y + v_s .* s_y

end