% Stokes Flow test case

%# of interior CV's
%Check of app.Ny is a multiple of 3
if mod(app.Ny,3)
    fprintf('NOTE: app.Ny was not a multiple of 3, change app.Ny from %g to %g\n',...
        app.Ny,app.Ny-mod(app.Ny,3));
    app.Ny = app.Ny-mod(app.Ny,3);
end
if app.Ny/3<=2
   display('app.Ny needs to be at least 9 for this test case, app.Ny changed to 9.')
   app.Ny = 9;
end

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

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

%Set boundary conditions
bcsDrho = [0 -1.0 0 0];bcsNrho = [0];        
bcsDP= [0];         bcsNP = [0 0 0 0];        
bcsOP= [1];
bcsDu= [0 1 0 0 0];  bcsNu= [];        
bcsOu= [5];
bcsDv= [0 0 0 0];  bcsNv=  [0];        

%Set 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);
Mask(1:round(app.Ny/3),1:round(app.Nx/5))=1;
Mask(end-round(app.Ny/3-1):end,1:round(app.Nx/5))=1;

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

%% Initialize fields
xp = linspace(dx/2, a-dx/2, app.Nx  );yp = fliplr(linspace(dy/2, b-dy/2, app.Ny  ));
xu = linspace( dx   ,a-dx , app.Nx-1);yu = fliplr(linspace(dy/2, b-dy/2, app.Ny  ));
xv = linspace(dx/2, a-dx/2, app.Nx  );yv = fliplr(linspace(dy    , b-dy, app.Ny-1));
[XP YP] = meshgrid(xp,yp);[XU YU] = meshgrid(xu,yu);[XV YV] = meshgrid(xv,yv);
rhoinit=@(X,Y) 0*((X>0.5) - (X<0.5))*0.5;
pinit = @(X,Y) zeros(size(X));
uinit = @(X,Y) 1*ones(size(X));
vinit = @(X,Y) zeros(size(X));
rho=[bcsDrho' ;bcsNrho' ;rhoinit(XP(idsP), YP(idsP))];
P = [bcsDP';bcsNP';  pinit(XP(idsP), YP(idsP))];
u = [bcsDu';bcsNu';  uinit(XU(idsu), YU(idsu))];
v = [bcsDv';bcsNv';  vinit(XV(idsv), YV(idsv))];
u(Nodeu(2:end-1,round(app.Nx/5)+2:end-1))=1/3;
