%% Sudden Expansion test case

%% Define Forcing functions
%'XY' is a 2d array. XY(:, 1) contains positions in x, XY(:, 1) contains
%positions in y. Time is just the current time (initial time starts at time=0)
%Fu and Fv are required functions that need to be set for the solver
Fu = @(time, XY) 0;
Fv = @(time, XY) 0;

%Set number of timesteps. Nt is a required variable that needs to be set
%for the solver
Nt = ceil(app.T/app.dt);

%Set size of domain
a = 20;          
b = 1;

% Mesh the domain
% This is a standard call and should not require much modification.
% [0, a] gives the min/max x coordinates, and [0, b] gives the min/max y
% coordinates.
[XP, YP, XU, YU, XV, YV, dx, dy] = meshdomain(app, [0, a], [0, b]);

%% create geometry mask
%This mask is user-created. Zeros are unmasked, and non-zero integers
%indicate a masked region. The domain sides should not be considered here.
%Each interior masked region could have a different number to allow
%different boundary conditions.
%Create mask of correct size. This is a standard call and should not
%require modification.
Mask = zeros(app.Ny, app.Nx);

%User-defined true-false (boolean) functions that define the mask. True
%when masked, false when not. (Note, this is not the only way to do this). 
top_rect = @(x, y) logical((y >= 2 / 3) .* (x <= 4));
bot_rect = @(x, y) logical((y <= 1 / 3) .* (x <= 4));

%Assign masked region, associate 1 with the top domain, and 2 for the
%bottom domain.
Mask(top_rect(XP, YP)) = 1;
Mask(bot_rect(XP, YP)) = 2;

%Additionally, we could add an elliptical obstacle to the flow as follows:
%ellipse = @(x, y) (((x - 10) .^2 / 20 + (y - 0.5) .^ 2) < 0.1);
%Mask(ellipse(XP, YP)) = 2;

%Add_bnd is an array that adds boundaries for the [left, right, bottom, top]
%if these are true. This is almost surely a standard definition.
add_bnd = [true true true true];
%create Node matrices (This is a standard call and should not be modified)
%All these outputs are required by the solver.
[NodeP, Nodeu, Nodev, idsP, idsu, idsv] = NodePad(Mask, add_bnd, 1);

%Plot to determine boundary ids
plot_bc_ids(NodeP, Mask)
display('Paused at boundary id plotting. Hit any key to continue')
pause

%Assign Boundary conditions types to the different ids: This requires user
%input, using the above plot as reference
bcid2type_u = {'D', 'D', 'Dmv', 'O', 'D', 'D'};      
bcid2type_v = {'D', 'D', 'D', 'N', 'D', 'D',};
bcid2type_P = {'N', 'N', 'N', 'O', 'N', 'N'};
%bcid2type_rho = {'N', 'Dmv', 'N', 'N', 'N', 'N'}; %An example for density

%This is a standard call and should not require much modification. In
%addition you could include boundary information for density, as done in
%the commented out section.
 [NodeP, nbcP, Nodeu, nbcu, Nodev, nbcv, NodeRho, nbcrho] = set_bcs(... 
     bcid2type_P, NodeP, bcid2type_u, Nodeu, bcid2type_v, Nodev);
%[NodeP, nbcP, Nodeu, nbcu, Nodev, nbcv, NodeRho, nbcrho] = set_bcs(... 
%    bcid2type_P, NodeP, bcid2type_u, Nodeu, bcid2type_v, Nodev,...
%    bcid2type_rho, NodeP); %Note, NodeRho == NodeP for the inputs only

%Set BC values to zero everywhere (These is a standard call and should not 
% require much modification). Most of the outputs below are required by the
% solver (with the exception of the boundary conditions for density).
[bcsDP, bcsOP, bcsNP] = init_bcs(nbcP);
[bcsDu, bcsOu, bcsNu] = init_bcs(nbcu);
[bcsDv, bcsOv, bcsNv] = init_bcs(nbcv);
[bcsDrho, bcsOrho, bcsNrho] = init_bcs(nbcrho);

%% Modify the boundary conditions
%This is where the user sets any non-zero boundary conditions
%First we obtain the coordinates for the boundary points. Each of the 
%XY arrays contain 2 columns, the first gives the x-coordinate, the second 
% gives the y-coordinate associated with a particular boundary (D, O, or N)
[XYD, XYO, XYN] = get_bc_coord('u', XU, YU, Nodeu, bcsDu, bcsOu, bcsNu);

%When the x-value is 0 (i.e. at the inlet) set the u-velocity dirichlet bc
%equal to 1
ids = (XYD(:, 1) == 0);
bcsDu(ids) = 1;

%On the bottom rectangle, set the right-hand face to have a linear density
%profile
%Note, we don't need XYD, etc from above anymore, so we can overwrite these
%variables
%[XYD, XYO, XYN] = get_bc_coord('P', XP, YP, NodeRho, bcsDrho, bcsOrho, bcsNrho);
%ids = logical(bot_rect(XYD(:, 1), XYD(:, 2)) .* (XYD(:, 1) >= (4 - dx)));
%bcsDrho(ids) = 1/3 - XYD(ids, 2);
%
%And for the ellipse we could add the following, for example:
%ids = ellipse(XYD(:, 1), XYD(:, 2));
%bcsDrho(ids) = (XYD(ids, 1) - mean(XYD(ids, 1))).^2;

%% Initialize fields
if isfield(app,'perturb')
    if isempty(app.perturb)
        perturb = 0.1;
    else
        perturb = app.perturb;
    end
else
    perturb = 0;
end

rhoinit=@(X, Y) perturb * ((X > 0.5) - (X < 0.5)) * 0.5;
pinit = @(X, Y) zeros(size(X));
uinit = @(X, Y) (X <= 4) + 1 / 3 * (X > 4);
vinit = @(X, Y) zeros(size(X));

%Most of the variables below (with the exception of density) are required
%by the solver.
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))];