% Stokes Flow test case

%Set forcing Functions
Fu = @(time,p) 0;
Fv = @(time,p) 0;

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

%Set boundary conditions
bcsDrho = [];      bcsNrho = [0 0 0 0];        
bcsDP= [];      bcsNP = [0 0 0 0];        
bcsDu= [0 0];      bcsNu= [0 0];        
bcsDv= [0 0 0 0];      bcsNv= [];        

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

T    = 1000;       dt   = T/(Nt);

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

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

%% 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.001*(sqrt((X-1/2).^2+(Y-1).^2) < 0.25) ;
pinit = @(X,Y) zeros(size(X));
uinit = @(X,Y) zeros(size(X));
vinit = @(X,Y) zeros(size(X));
% vinit = @(X,Y) rhoinit(X,Y);
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))];
