% Stokes Flow Poiseuille flow test case

% Forcing Function
Fu = @(time,p) 0;
Fv = @(time,p) 0;

%Set the size of the domain
a    = 2;          b    = 1;
dx   = a./(app.Nx);    dy   = b./(app.Ny); 


% Set the boundary conditions
    bcsDP= [0];        bcsNP= [0 0 0];        
    bcsOP=[1];
    bcsDu= [1 0 0]  ;  bcsNu= [0];        
    bcsDv= [0 0 0]  ;  bcsNv= [0];               

%Set the number of timesteps
Nt = ceil(app.T/app.dt);

%% create geometry mask
Nbcs = length(bcsDP) + length(bcsNP);
bnd = [1 1 1 1];
Mask = zeros(app.Ny, app.Nx);

%create Node matrices
[NodeP, Nodeu, Nodev, idsP, idsu, idsv] = NodePad(Mask,bnd,1);
%Modify boundaries
NodeP(NodeP==2)=-5;
NodeP(NodeP==1)=2;
NodeP(NodeP==-5)=1;
Nodeu(Nodeu==4)=-5;
Nodeu(Nodeu==2)=4;
Nodeu(Nodeu==-5)=2;
Nodev(Nodev==4)=-5;
Nodev(Nodev==2)=4;
Nodev(Nodev==-5)=2;
Nodev(Nodev==3)=-5;
Nodev(Nodev==1)=3;
Nodev(Nodev==-5)=1;

%% Initialize fields
xp = linspace(dx/2, a-dx/2, app.Nx  )*2-a;yp = fliplr(linspace(dy/2, b-dy/2, app.Ny  ))*2-b;
xu = linspace( dx   ,a-dx , app.Nx-1)*2-a;yu = fliplr(linspace(dy/2, b-dy/2, app.Ny  ))*2-b;
xv = linspace(dx/2, a-dx/2, app.Nx  )*2-a;yv = fliplr(linspace(dy    , b-dy, app.Ny-1))*2-b;
[XP YP] = meshgrid(xp,yp);[XU YU] = meshgrid(xu,yu);[XV YV] = meshgrid(xv,yv);
P = [bcsDP';bcsNP'; -0*XP(:)/a];
u = [bcsDu';bcsNu';1+0.*exp((-a-XU(:)))+idsu*0];
v = [bcsDv';bcsNv'; 0*idsv];
dx=dx*2;dy=dy*2;
