function [uh] = vcycle(uh,bh,dcoef_xh,dcoef_yh,level,p0,nu)%  [uh] = vcycle(uh,bh,dcoef_xh,dcoef_yh,level,p0,nu)%%  Apply one Multigrid V-cycle to the system%    Lh*uh = bh.%  Lh is the CCFD discretization of the diffusion operator %    L*v = -div (dcoef grad u).   [nx,ny] = size(uh);  N = nx*ny;  Lh = get_L(dcoef_xh,dcoef_yh);  L = tril(Lh,-1);  D = spdiags(diag(Lh),0,N,N);   U = triu(Lh,1);  bh = Qproject(bh,p0);  %  Perform nu steps of Gauss-Seidel presmoothing.   for gsiter = 1:nu    uh = uh(:) + (D + L)\(bh(:) - Lh*uh(:));    uh = reshape(uh,nx,ny);%    uh = Qproject(reshape(uh,nx,ny),p0);  %  Project onto orthocomp. of VH.  end %  Recursive call to vcycle.  if level > p0 + 1    level = level-1;    rh = reshape( bh(:) - Lh*uh(:) , nx,ny);    [r2h] = restrict(rh);    [dcoef_xH,dcoef_yH] = restrict_dcoef(dcoef_xh,dcoef_yh);    u2h = zeros(size(r2h));    [u2h] = vcycle(u2h,r2h,dcoef_xH,dcoef_yH,level,p0,nu);    uh = uh + prolong(u2h);  else%%%    uh = pinv(full(Lh)) * bh;  end%  Perform nu steps of Gauss-Seidel postsmoothing.  Note that pre- and %  post-smoothing operators are adjoints.   for gsiter = 1:nu    uh(:) = uh(:) + (D + U)\(bh(:) - Lh*uh(:));%    uh = reshape( Qproject(uh,p0) , nx,ny);  %  Project onto orthocomp. of V0.  enduh = reshape( Qproject(uh,p0) , nx,ny);  %  Project onto orthocomp. of V0.