%% Heated Plate Demo

% Last modified by Jing Lin in Apr. 2014


%% Set up data and images saving directory


%% Set up problem-specific parameterr.

if ~isfield(app, 'nu'), error('app.nu not specified'); end;


%% Set up mesh and solution domain

% Define size of solution domain
a = 4;
b = 5;

if ~isfield(app, 'Nx'), error('app.Nx not specified'); end;
if ~isfield(app, 'Ny'), error('app.Ny not specified'); end;

% 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]);


%% Set up the 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);

% Mask central plate
Eps = 1e-8;
plate = @(x, y) (x >= 0.25*a*(1-Eps)) & (x <= 0.75*a*(1+Eps)) ...
    & (y >= 0.25*b*(1-Eps)) & (y <= 0.75*b*(1+Eps));
Mask(plate(XP, YP)) = 1; 

% Logical array add_bnd signifies whether to add (add if =true) boundaries 
% for the [west, east, south, north] sides of the solution domain.
% This is almost a standard procedure except when periodic boundaries are used.
add_bnd = [true, true, true, true];

% Create preliminary node numbering matrices using function "NodePad".
% Note: This is a standard procedure and should not be modified.
[NodeP, Nodeu, Nodev, idsP, idsu, idsv] = NodePad(Mask, add_bnd);

% Plot the current boundary IDs to facilitate setting up boundary conditions.
if ~app.NoPlot
    plot_bc_ids(NodeP, idsP);
    fprintf('Paused at boundary id plot. Hit any key to continue.\n\n');
    pause;
end;


%% Set up boundary conditions (BCs)

% User assigns BC types to each boundary ID for u, v, rho in three cell arrays.

% This requires user input. Use the above plot as reference to identify BC IDs.
% e.g. Element 1 in cell array bcid2type_u represents the BC type of
% the boundary with ID 1 for velocity u.

% Note: Only the following four types of BCs can be set for velocities
% and density:
% Single-value Dirichlet BC   (='D')
% Multiple-value Dirichlet BC (='Dmv')
% Single-value Neumann BC     (='N')
% Multiple-value Neumann BC   (='Nmv')

% Warning: Open boundary conditions (='O') are not allowed for velocities
% and density!

bcid2type_u = {'D', 'D', 'D', 'N', 'N'};      
bcid2type_v = {'D', 'D', 'D', 'N', 'N'};
bcid2type_q = {'Dmv', 'Dmv', 'Dmv', 'N', 'N'};

% Set up the final node numbering matrix using function "set_bcs".
% For each single-value BC ('D', 'N', 'O'), all its boundary cells share one ID,
% like in the preliminary node numbering matrices.
% For each multiple-value BC ('Dmv', 'Nmv'), now each boundary cell will
% be assigned a different ID.
% Besides, the IDs are arranged to satisfy: 
%           Dirichlet_BC_IDs < Open_BC_IDs < Neumann_BC_IDs

% Note: This is a standard procedure and should not be modified.

% Without density equation:
[NodeP, nbcq, Nodeu, nbcu, Nodev, nbcv] = set_bcs(...
    bcid2type_q, NodeP, bcid2type_u, Nodeu, bcid2type_v, Nodev);

% With density equation:
% [NodeP, nbcq, Nodeu, nbcu, Nodev, nbcv, NodeRho, nbcrho] = set_bcs(... 
%     bcid2type_q, NodeP, bcid2type_u, Nodeu, bcid2type_v, Nodev, ...
%     bcid2type_rho, NodeP);

% Set preliminary BC values to zero since most BCs are homogeneous.
% Some of them can be modified to be non-homogeneous BCs later.

% Note: This is a standard procedure and should not be modified.

[bcsDq, bcsOq, bcsNq] = init_bcs(nbcq);
[bcsDu, bcsOu, bcsNu] = init_bcs(nbcu);
[bcsDv, bcsOv, bcsNv] = init_bcs(nbcv);

% User sets non-zero boundary conditions.

% The coordinates of the boundary points on each grid can be obtained...
% by function "get_bc_coord".
% Each of the three XY arrays contains 2 columns with x coordinates in
% the first and y coordinates in the second stored in order of BC IDs.
% The three XY arrays correspond to boundary points of type
% Dirichlet ('D' and 'Dmv'), open ('O'), and Neumann ('N' and 'Nmv').

[XYD, ~, XYN] = get_bc_coord('P', XP, YP, NodeP, bcsDq, bcsOq, bcsNq);
bcsDq(XYD(:,1)<Eps) = 20;
bcsDq((XYD(:,1)>Eps)&(XYD(:,1)<a*(1-Eps))) = 80;
bcsNq(XYN(:,2)<Eps) = -20;
bcsNq(XYN(:,2)>b*(1-Eps)) = 20;

[XYD, ~, XYN] = get_bc_coord('u', XU, YU, Nodeu, bcsDu, bcsOu, bcsNu);
bcsDu(XYD(:,1)<Eps) = 20;
bcsDu((XYD(:,1)>Eps)&(XYD(:,1)<a*(1-Eps))) = 80;
bcsNu(XYN(:,2)<Eps) = -20;
bcsNu(XYN(:,2)>b*(1-Eps)) = 20;

% [XYD, XYO, XYN] = get_bc_coord('v', XV, YV, Nodev, bcsDv, bcsOv, bcsNv);
% bcsDv(XYD(:,1)<Eps) = 20;
% bcsDv((XYD(:,1)>Eps)&(XYD(:,1)<a*(1-Eps))) = 80;
% bcsNv(XYN(:,2)<Eps) = -20;
% bcsNv(XYN(:,2)>b*(1-Eps)) = 20;

% Another way to set boundary conditions when BCs are single-valued
% bcsDq(1:3) = [80, 20, 0]';
% bcsNq(1:2) = [-20, 20]';
% bcsDu(1:3) = [80, 20, 0]';
% bcsNu(1:2) = [-20, 20]';
bcsDv(1:3) = [80, 20, 0]';
bcsNv(1:2) = [-20, 20]';


%% Set up initial conditions


% User specifies the functions that provide initial values at given coordinates.

P_init   = @(X, Y) zeros(size(X));
u_init   = @(X, Y) zeros(size(X));
v_init   = @(X, Y) zeros(size(X));

% Initialize the state variables with values from both boundary and internal
% cells.

% Note: This is a standard procedure and should not be modified.

P   = [bcsDq; bcsNq; P_init(XP(idsP), YP(idsP))];
u   = [bcsDu; bcsNu; u_init(XU(idsu), YU(idsu))];
v   = [bcsDv; bcsNv; v_init(XV(idsv), YV(idsv))];


%% Preliminary setup

% Find number of boundary conditions and number of interior cells
Nbcs = length(bcsDq) + length(bcsNq);

% Get interior ids
pid = Nbcs+1:length(P);
uid = Nbcs+1:length(u);
vid = Nbcs+1:length(v);


%% Constructing A and RHS

display('Creating Operators');

[AP, ~, rhsP] = Diffusion_Operator('P', NodeP, bcsDq, bcsOq, bcsNq, dx, dy, ... 
    app.nu, app.nu);
[Au, ~, rhsu] = Diffusion_Operator('u', Nodeu, bcsDu, bcsOu, bcsNu, dx, dy, ...
    app.nu, app.nu);
[Av, ~, rhsv] = Diffusion_Operator('v', Nodev, bcsDv, bcsOv, bcsNv, dx, dy, ...
    app.nu, app.nu);


%% Solve the linear system

display('Solving the system');
P(pid) = AP\rhsP;
u(uid) = Au\rhsu;
v(vid) = Av\rhsv;


%% Plot the solution

display('Plotting Solution')

figure;
subplot(2,2,1);
surf(XP, YP, P(NodeP(2:app.Ny+1, 2:app.Nx+1))) ;
TT(1)=xlabel('x(m)'); TT(2)=ylabel('y(m)'); TT(3)=zlabel('T(^oC)');
TT(4)=title(sprintf('P-grid Nx=%g, Ny=%g', app.Nx, app.Ny));
formatplot([],gca,TT);

subplot(2,2,2);
surf(XU, YU, u(Nodeu(2:app.Ny+1, 2:app.Nx)));
TT(1)=xlabel('x(m)'); TT(2)=ylabel('y(m)'); TT(3)=zlabel('T(^oC)');
TT(4)=title(sprintf('u-grid Nx=%g, Ny=%g', app.Nx-1, app.Ny));
formatplot([],gca,TT);

subplot(2,2,3);
surf(XV, YV, v(Nodev(2:app.Ny, 2:app.Nx+1)));
TT(1)=xlabel('x(m)'); TT(2)=ylabel('y(m)'); TT(3)=zlabel('T(^oC)');
TT(4)=title(sprintf('v-grid Nx=%g, Ny=%g', app.Nx, app.Ny-1));
formatplot([],gca,TT);
