% Plots for stokes problems
Pex = [bcsDq;bcsNq;pexact(T,[XP(idsP) YP(idsP)])];
uex = [bcsDu;bcsNu;uexact(T,[XU(idsu) YU(idsu)])];
vex = [bcsDv;bcsNv;vexact(T,[XV(idsv) YV(idsv)])];

clf
colormap(mk_colmap)

subplot(331)
surf(XP,YP,P(NodeP(2:app.Ny+1,2:app.Nx+1)),'LineStyle','none') 
xlim([min(xp) max(xp)]);ylim([min(yp) max(yp)]);
TT(1)=xlabel('x(m)');TT(2)=ylabel('y(m)');TT(3)=zlabel('T(^oC)');
TT(4)=title(sprintf('Pressure Numerical',app.Nx,app.Ny));
formatplot([],gca,TT);view(2),colorbar,shading interp,caxis([-1 1]*pi)

subplot(334)
surf(xu,yu,u(Nodeu(2:app.Ny+1,2:app.Nx)),'LineStyle','none') 
xlim([min(xu) max(xu)]);ylim([min(yu) max(yu)]);
TT(1)=xlabel('x(m)');TT(2)=ylabel('y(m)');TT(3)=zlabel('T(^oC)');
TT(4)=title(sprintf('U velocity Numerical',app.Nx,app.Ny));
formatplot([],gca,TT);view(2),colorbar,shading interp,caxis([-1 1]*pi)

subplot(337)
surf(xv,yv,v(Nodev(2:app.Ny,2:app.Nx+1)),'LineStyle','none') 
xlim([min(xv) max(xv)]);ylim([min(yv) max(yv)]);
TT(1)=xlabel('x(m)');TT(2)=ylabel('y(m)');TT(3)=zlabel('T(^oC)');
TT(4)=title(sprintf('V velocity Numerical',app.Nx,app.Ny));
formatplot([],gca,TT);view(2),colorbar,shading interp,caxis([-1 1]*pi)

subplot(332)
surf(XP,YP,Pex(NodeP(2:app.Ny+1,2:app.Nx+1)),'LineStyle','none') 
xlim([min(xp) max(xp)]);ylim([min(yp) max(yp)]);
TT(1)=xlabel('x(m)');TT(2)=ylabel('y(m)');TT(3)=zlabel('T(^oC)');
TT(4)=title(sprintf('Pressure Exact',app.Nx,app.Ny));
formatplot([],gca,TT);view(2),colorbar,shading interp,caxis([-1 1]*pi)

subplot(335)
surf(xu,yu,uex(Nodeu(2:app.Ny+1,2:app.Nx)),'LineStyle','none') 
xlim([min(xu) max(xu)]);ylim([min(yu) max(yu)]);
TT(1)=xlabel('x(m)');TT(2)=ylabel('y(m)');TT(3)=zlabel('T(^oC)');
TT(4)=title(sprintf('U velocity Exact',app.Nx,app.Ny));
formatplot([],gca,TT);view(2),colorbar,shading interp,caxis([-1 1]*pi)

subplot(338)
surf(xv,yv,vex(Nodev(2:app.Ny,2:app.Nx+1)),'LineStyle','none') 
xlim([min(xv) max(xv)]);ylim([min(yv) max(yv)]);
TT(1)=xlabel('x(m)');TT(2)=ylabel('y(m)');TT(3)=zlabel('T(^oC)');
TT(4)=title(sprintf('V velocity Exact',app.Nx,app.Ny));
formatplot([],gca,TT);view(2),colorbar,shading interp,caxis([-1 1]*pi)


subplot(333)
du=[Pex-P];
surf(xp,yp,du(NodeP(2:app.Ny+1,2:app.Nx+1)),'LineStyle','none') 
xlim([min(xv) max(xv)]);ylim([min(yv) max(yv)]);
TT(1)=xlabel('x(m)');TT(2)=ylabel('y(m)');TT(3)=zlabel('T(^oC)');
TT(4)=title(sprintf('Pex-P',app.Nx,app.Ny));
formatplot([],gca,TT);view(2),colorbar,shading interp

subplot(336)
tmp=[uex-u];
surf(xu,yu,tmp(Nodeu(2:app.Ny+1,2:app.Nx)),'LineStyle','none') 
xlim([min(xu) max(xu)]);ylim([min(yu) max(yu)]);
TT(1)=xlabel('x(m)');TT(2)=ylabel('y(m)');TT(3)=zlabel('T(^oC)');
TT(4)=title(sprintf('uex-u',app.Nx,app.Ny));
formatplot([],gca,TT);view(2),colorbar,shading interp

subplot(339)
tmp=[vex-v];
surf(xv,yv,tmp(Nodev(2:app.Ny,2:app.Nx+1)),'LineStyle','none') 
xlim([min(xv) max(xv)]);ylim([min(yv) max(yv)]);
TT(1)=xlabel('x(m)');TT(2)=ylabel('y(m)');TT(3)=zlabel('T(^oC)');
TT(4)=title(sprintf('vex-v',app.Nx,app.Ny));
formatplot([],gca,TT);view(2),colorbar,shading interp

ha = axes('Position',[0 0 1 1],'Xlim',[0 1],'Ylim',...
    [0 1],'Box','off','Visible','off','Units','normalized', 'clipping' , 'off');
text(0.5, 1,sprintf('Time = %04.03g, app.Nx=%g,app.Ny=%g',T,app.Nx,app.Ny),'HorizontalAlignment',...
    'center','VerticalAlignment', 'top','fontsize',13,'FontWeight','bold')

drawnow;
