% Naviers Stokes Flow Lock Exchange test case
d0=fileparts([pwd,filesep]);
d0=d0(1:strfind(d0,'trunk')+4);
%Forcing Function
Fu = @(time,p) 0;
Fv = @(time,p) 0;

%Set the domain size
if isfield(app,'a')
    a = app.a;
else
    a = 2;
end
if isfield(app,'b')
    b = app.b;
else
    b      = 2;
end
if ~isfield(app,'nu')
    if isfield(app,'Gr')
        app.nu = 1/sqrt(app.Gr);
        if isfield(app,'Sc')
            app.kappa = 1/sqrt(app.Gr)/app.Sc;
        else
            app.kappa = 1/sqrt(app.Gr);
        end
    end
end

dx     = a./(app.Nx);    dy     = b./(app.Ny); 

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

%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);
Nodev(Nodev==1) = -3;
Nodev(Nodev==2) = -4;
Nodev(Nodev==3) = 1;
Nodev(Nodev==4) = 2;
Nodev(Nodev==-3) = 3;
Nodev(Nodev==-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);
XP = XP - 1;
YP = YP - 1;
XU = XU - 1;
YU = YU - 1;
XV = XV - 1;
YV = YV - 1;
%NOTE USED
rhoinit=@(X,Y) ((X>0) - (X<0)) * 0.5 * app.drho;
%USED
pinit = @(X,Y) zeros(size(X));
uinit = @(X,Y) zeros(size(X));
vinit = @(X,Y) zeros(size(X));

rho = [bcsDrho' ;bcsNrho' ; app.drho / 2 * tanh((XP(idsP)) * 2 / app.var)];
% 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))];

if isfield(app,'noise')
    rho(Nbcs+1:end) = rho(Nbcs+1:end) + app.noise*rand(size(rho(Nbcs+1:end)));
end
