% Heated_Plate Demo
%Define size of plate
a=4;
b=5;

%Calculate x and y divisions
dx  = a./(Nx); 
dy = b./(Ny);  

%Set boundary conditions [First numbers correspond to Dirichlet, and larger
%numbers correspond to Neumann
bcsD = [80 20 0];

%Note, Neumann require multiplication through by dn, that is the spacing in
%the normal direction.
bcsN = [-20*dy 20*dy];

%create geometry mask
Mask = zeros(Ny, Nx);
Mask(round(Ny/4):round(3*Ny/4),round(Nx/4):round(3*Nx/4))=1;
Mask(1:10,1:10) = 0;

%create Node matrices
[NodeP, Nodeu, Nodev] = NodePad(Mask,[1 1 1 1],1);

%% Constructing A and RHS
display('Creating Operators')
[AP rhsP] =  mk_DiffusionOperator(Nx, Ny, dx, dy, bcsD,bcsN, NodeP,'P');
[Au rhsu] =  mk_DiffusionOperator(Nx, Ny, dx, dy, bcsD,bcsN, Nodeu,'u');
[Av rhsv] =  mk_DiffusionOperator(Nx, Ny, dx, dy, bcsD,bcsN, Nodev,'v');
[Au Av rhsu rhsv] = FixDiffuvDbcs(Au,Av,bcsD,bcsD,...
    rhsu,rhsv,Nodeu, Nodev, dx,dy,length(bcsD)+length(bcsN));

%% Solving System
display('Solving the system')
SolP = AP\rhsP; 
SolP=[bcsD'; NaN*bcsN' ;SolP];
Solu = Au\rhsu; 
Solu=[bcsD'; NaN*bcsN' ;Solu];
Solv = Av\rhsv; 
Solv=[bcsD'; NaN*bcsN' ;Solv];

%% Plot the solution
display('Plotting Solution')
xp = linspace(-dx/2, a+dx/2, Nx+2);
yp = fliplr(linspace(-dy/2, b+dy/2, Ny+2));
xu = linspace( 0   , a     , Nx+1);
yu = fliplr(linspace(-dy/2, b+dy/2, Ny+2));
xv = linspace(-dx/2, a+dx/2, Nx+2);
yv = fliplr(linspace(0    , b     , Ny+1));

figure
subplot(221)
surf(xp,yp,SolP(NodeP)) 
xlim([min(xp) max(xp)]);ylim([min(yp) max(yp)]);
TT(1)=xlabel('x(m)');TT(2)=ylabel('y(m)');TT(3)=zlabel('T(^oC)');
TT(4)=title(sprintf('Centrally Heated Plate Nx=%g, Ny=%g',Nx,Ny));
formatplot([],gca,TT);

subplot(222)
surf(xu,yu,Solu(Nodeu)) 
xlim([min(xu) max(xu)]);ylim([min(yu) max(yu)]);
TT(1)=xlabel('x(m)');TT(2)=ylabel('y(m)');TT(3)=zlabel('T(^oC)');
TT(4)=title(sprintf('Centrally Heated Plate',Nx,Ny));
formatplot([],gca,TT);

subplot(223)
surf(xv,yv,Solv(Nodev)) 
xlim([min(xv) max(xv)]);ylim([min(yv) max(yv)]);
TT(1)=xlabel('x(m)');TT(2)=ylabel('y(m)');TT(3)=zlabel('T(^oC)');
TT(4)=title(sprintf('Centrally Heated Plate',Nx,Ny));
formatplot([],gca,TT);
