% Naviers Stokes Flow Lock Exchange test case

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

%Set the domain size
if isfield(app,'a')
    a = app.a;
else
    a = 1;
end
if isfield(app,'b')
    b = app.b;
else
    b      = 1;
end
dx     = a./(app.Nx);    dy     = b./(app.Ny); 

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

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

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

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

%% 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.5 * ones(size(X));
pinit = @(X,Y) zeros(size(X));
uinit = @(X,Y) zeros(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))];


