% 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= [20 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);


b2=-0.1; b4=0.0;
rhoinit=@(X,Y) b2*2*cos(pi*X).*sin(pi*Y) + b4*2*cos(2*pi*X).*sin(pi*Y);



pinit = @(X,Y) zeros(size(X));

a1=2.365; a2=3.9266;

uinit = @(X,Y) 0.3*(cos(a1*(X-1/2))./cos(a1/2)-cosh(a1*(X-0.5))./cosh(a1/2)).*(-a1*sin(a1*(Y-1/2))./cos(a1/2)-a1*sinh(a1*(Y-0.5))./cosh(a1/2))...
    +0.2*(sin(a2*(X-1/2))./sin(a2/2)-sinh(a2*(X-0.5))./sinh(a2/2)).*(a2*cos(a2*(Y-1/2))./sin(a2/2)-a2*cosh(a2*(Y-0.5))./sinh(a2/2));
vinit = @(X,Y) -0.3*(cos(a1*(Y-1/2))./cos(a1/2)-cosh(a1*(Y-0.5))./cosh(a1/2)).*(-a1*sin(a1*(X-1/2))./cos(a1/2)-a1*sinh(a1*(X-0.5))./cosh(a1/2))...
    -0.2*(sin(a2*(Y-1/2))./sin(a2/2)-sinh(a2*(Y-0.5))./sinh(a2/2)).*(a2*cos(a2*(X-1/2))./sin(a2/2)-a2*cosh(a2*(X-0.5))./sinh(a2/2));;


rho=[bcsDrho' ;bcsNrho' ;rhoinit(XP(idsP), YP(idsP))];
P = [bcsDP';bcsNP';  pinit(XP(idsP), YP(idsP))];
u = [bcsDu';bcsNu';  0*uinit(XU(idsu), YU(idsu))];
v = [bcsDv';bcsNv';  0*vinit(XV(idsv), YV(idsv))];

