%%  This file is a MATLAB script for computing heat flux
%%  and heat flow distributions from a 2D transient conduction
%%  simulation (plate.m).  

%%  The mesh data (p,t,e) and the solution matrix (u) must be
%%  exported from pdetool.  For transient simulations, 
%%  each column of u contains the results for one time step,
%%  i.e.,  u(:,j) is the full solution at the jth time step. 
%%  The function TRI2GRID interpolates u(:,j) onto an xy grid
%%  as uxy(j).  The function uxy can then be plotted, integrated
%%  or differentiated.

%%  Define the grid:  
%%   * dy and dx are the increments of the grid. If they are 
%%     set too small relative to the original mesh, uxy will 
%%     show an uneven, `stair-step' variation.
%%   * xmax, xmin, ymax, ymin are the limits of x and y used
%%     in the simulation domain.  If the domain is drawn with 
%%     the GUI without setting the SNAP-TO-GRID option, 
%%     these values may not correspond precisely to 
%%     the boundaries of the computational domain; if they 
%%     lie outside the domain, TRI2GRID will return `NaN' for
%%     those points.  In that case, the limits of the simulation
%%     should be adjusted in the pdetool .m file (e.g., plate.m).

dy = 0.05;
ymin = -0.6;
ymax = 0.6;

dx = 0.05;
xmin = -1.;
xmax = 1.;

%%  Time interval:
%%    * dt is the time step used in the simulaton
%%    * nt is the total number of time steps, including time zero.

dt = 12500;
nt =  21;   
tmax = dt*(nt-1);

%%  Define the row vectors x and y from which the grid is constructed

x = xmin:dx:xmax;
y = ymin:dy:ymax;

%%  Number of points in x and in y

nx =  ((xmax-xmin)/dx)+1;
ny =  ((ymax-ymin)/dy)+1;

%%  Use a loop to successively process the data for each time step 

for j=1:nt

  uxy = tri2grid(p,t,u(:,j),x,y);  % time dependent case at jth time step
                                     % for time-indep case, replace u(:,j) by u

%%  Differentiate the temperature field in x direction to get the 
%%  x-component of heat flux along the line x=xmin. (This is a 
%%  distribution of flux in the y direction).  The flux is stored
%%  in a matrix q whose columns correspond to the heat flux at each time step.
 
  k = 14;     %% Thermal conductivity of plate

%%  q(:,j)  = -k*(uxy(:,2) - uxy(:,1))/dx;               % first-order forward difference

  q(:,j) = -k*(-3*uxy(:,1)+4*uxy(:,2)-uxy(:,3))/(2*dx);  % second-order forward difference

  qr = q(:,j)' ;  % Transpose the jth column of matrix to get a row vector. 

%% Use a trapezodial-rule integration of qr with respect to y to get
%% the heat flow into the domain through the border at xmin.

  Q(j) = trapz(y,qr); 

  Q(j)  % Display numerical result in MATLAB window while running code

%% Save temperature distribution at xmin for each time-step

  Ty(:,j) = uxy(:,1);

end ;

time = dt:dt:tmax;         % Time vector from t=dt to tmax for plotting.
Qr = Q(2:nt);               % Remove t=0 point (first element of vector) prior to plotting.
qr = q(:,2:nt);              % Remove t=0 point (first column of matrix) prior to plotting.

%%  Plotting commands

figure('NextPlot','replace');

plot(time,Qr);            % Plot heat flow as a function of time 
xlabel('time (sec)');
ylabel('Heat flow (W/m)');
axis([0 tmax 0 2000]);

pause;
figure('NextPlot','replace');

plot(y,qr);             % Plot heat flux versus y for each time-step
axis([ymin ymax 0 1500]);
xlabel('y (m)');
ylabel('Heat flux (W/m^2)');

pause;
figure('NextPlot','replace');

plot(y,Ty);             % Plot temperature versus y for each time-step
xlabel('y (m)');
ylabel('Temperature (C)');