% Author: Layla S. Mayboudi
% 2021-05-27
% This code creates a transient thermal PDE model and plots the geometry 
% and heat flux flow diagrams.
% Selected nodes and their temperatures are also displayed.

% Note:
% A transient model uses ThemalConductivity, MassDensity and SpecificHeat. 


clear all; close all; clf; clc;

set(groot,'DefaultAxesColorOrder',[0 0 0],...
    'DefaultSurfaceEdgeColor',[0 0 0],...
    'DefaultLineLineWidth',1,...
    'DefaultAxesLineStyleOrder','-|--|:|-.|+|o|*|.|x|s|d|^|v|>|<|p|h'); %may be added to set the line style and color.
ax.FontSize = 20; ax.FontName = 'Cambria'; ax.FontSmoothing = 'on'; ax.FontWeight = 'normal'; %'bold'; 

  %'DefaultAxesLineStyleOrder','-|--|:|-.|+|o|*|.|x|s|d|^|v|>|<|p|h' %may be added to set the line style and color.
  %set(gcf,'defaultSurfaceEdgeColor','b')

thermalModelT = createpde('thermal','transient'); % create a model object: thermal
width = 50/1000; height1 = 26/1000; height2 = 30/1000;
geom1 = [3;4;0;width;width;0;0;0;height1;height1]; % create geometry: rectangle solid
geom2 = [3;4;0;width;width;0;height1;height1;height2;height2]; % create geometry: rectangle solid

gd = [geom1 geom2];
sf = 'geom1+geom2';
ns = char('geom1','geom2');
ns = ns';
geometry1 = decsg(gd,sf,ns);
geometryFromEdges(thermalModelT,geometry1); % associate geometry and thermal model object
figure(1);
pdegplot(thermalModelT,'EdgeLabels','on','FaceLabels','on');
xlim([0 width]); ylim([0 height2]);
title({'Geometry with Edge and Face Labels Displayed';''},'FontWeight','normal');
grid on; grid minor; xlabel('x (m)'); ylabel('y (m)');
ax = gca; ax.FontSize = 14;

%Water
k1 = 0.6305; % thermal conductivity (W/mK)
rho1 = 988; % density (kg/m3)
cp1 = 4181; % heat capacity (J/kgK)
q1 = 0; % heat generation (W/m^3)
kFunc1 = @(region,state) k1*region.y;
cFunc1 = @(region,state) cp1*region.y;
qFunc1 = @(region,state) q1*region.y;
rhoFunc1 = @(region,state) rho1*region.y;
tp=thermalProperties(thermalModelT,'ThermalConductivity',k1,...
    'MassDensity',rho1,'SpecificHeat',cp1,'Face',1);
ih1=internalHeatSource(thermalModelT,q1,'Face',1);

%Aluminum
k2 = 210; % thermal conductivity (W/mK)
rho2 = 2698.9; % density (kg/m3)
cp2 = 900; % heat capacity (J/kgK)
q2 = 1000; % heat generation (W/m^3)
kFunc2 = @(region,state) k2*region.y;
cFunc2 = @(region,state) cp2*region.y;
qFunc2 = @(region,state) q2*region.y;
rhoFunc2 = @(region,state) rho2*region.y;
thermalProperties(thermalModelT,'ThermalConductivity',k2,...
    'MassDensity',rho2,'SpecificHeat',cp2,'Face',2);
ih2=internalHeatSource(thermalModelT,q2,'Face',2);

% Boundary Conditions
T_ambient = 22; %degC **********
T_water = 50;   %degC **********
outerCC_V = @(region,~) 10*region.y; %horizontal heat convection coefficient (W/m2K)
top_BC_HF = @(region,~) 100*region.y; %heat source distribution (W/m2)
thermalBC(thermalModelT,'Edge',2,'HeatFlux',top_BC_HF); %boundary condition: heat source (W/m2)
thermalBC(thermalModelT,'Edge',[4 5 6 7],'ConvectionCoefficient',...
outerCC_V,'AmbientTemperature',T_ambient); %boundary condition: convection heat transfer

% Initial Condition
T0_1 = 50; %Initial tempeature of Face 1 (water) **********
T0_2 = 25; %Initial tempeature of Face 2 (pipe body) **********
IC1=thermalIC(thermalModelT,T0_1,'Face',1)
IC2=thermalIC(thermalModelT,T0_2,'Face',2)

figure(2);
grid_size = 1.5/1000; %identify grid size
mesh1 = generateMesh(thermalModelT,'Hmax',grid_size); %generate mesh *****
pdeplot(mesh1); %plot PDE thermal model object with node ID's displayed.
xlim([0 width]); ylim([0 height2]);
title({'Triangular Elements';''},'FontWeight','normal');
grid on; grid minor; xlabel('x (m)'); ylabel('y (m)');
ax = gca; ax.FontSize = 14;

figure(3);
pdeplot(mesh1,'NodeLabels','on','ElementLabels','on'); %plot PDE thermal model object with node ID's displayed.
xlim([0 grid_size]); ylim([0 grid_size]); % set axis conditions
title('','FontWeight','normal');
title({'Node and Element IDs';''},'FontWeight','normal');
grid on; grid minor; xlabel('x (m)'); ylabel('y (m)');
ax = gca; ax.FontSize = 14;

%thermalModelT.SolverOptions.PDESolverOptions Properties (or the following)
%AbsoluteTolerance, RelativeTolerance, ResidualTolerance, MaxIterations, MinStep, ResidualNorm, ReportStatistics
%pde.TransientThermalResults.Temperature (or the following)
% Temperature, XGradients, YGradients, and ZGradients 

tfinal = 340*60; % set final time (s)
time_step = 100; % time step (s)
tlist = 0.1:time_step:tfinal; % set time array
t2 = tlist/60;
thermalModelT.SolverOptions.ReportStatistics = 'on';
results = solve(thermalModelT,tlist); % identify PDESolverOptions Properties
T1 = results.Temperature; % identify data to process and plot *****
Tx1 = results.XGradients; % identify data gradient along x-coordinate *****
Ty1 = results.YGradients; % identify data gradient along y-coordinate *****
T_st1 = results.SolutionTimes; % identify data gradient along y-coordinate *****
 
figure(4);
tt=round(tfinal/3600,1);
pdeplot(thermalModelT,'XYData',T1(:,end),'ColorMap','hot','Contour','on','Levels',10); %'ColorMap','hot' ,'Contour','on'
xlim([0 width]); ylim([0 height2]);
title({'Temperature Contours at Last Time Step'; ['(t = ',...
    num2str(tt),' hr)'];''},'FontWeight','normal')
grid on; grid minor; xlabel('x (m)'); ylabel('y (m)');
ax = gca; ax.FontSize = 14;

figure(5);
pdeplot(thermalModelT,'XYData',Tx1(:,end),'ColorMap','hot','Contour','on','Levels',10); %'ColorMap','hot' ,'Contour','on'
xlim([0 width]); ylim([0 height2]); 
title({'Temperature Gradient (dT/dx)'; ['Contours at Last Time step (t = ',...
    num2str(tt),' hr)'];''},'FontWeight','normal')
grid on; grid minor; xlabel('x (m)'); ylabel('y (m)'); xlabel('x (m)'); ylabel('y (m)');
ax = gca; ax.FontSize = 14;

figure(6);
pdeplot(thermalModelT,'XYData',Ty1(:,end),'ColorMap','hot','Contour','on','Levels',10); %'ColorMap','hot' ,'Contour','on'
xlim([0 width]); ylim([0 height2]); 
title({'Temperature Gradient (dT/dy)'; ['Contours at Last Time Step (t = ',...
    num2str(tt),' hr)'];''},'FontWeight','normal')
 grid on; grid minor; xlabel('x (m)'); ylabel('y (m)');
 ax = gca; ax.FontSize = 14;
 
%**********************************************************
%Investigate the node density

point = thermalModelT.Mesh.Nodes; % recall node locations from thermal model object
eps1=eps;
nodesIndex = find(point); % identify far left boundary nodes
xy00=point; x00=xy00(1,:); y00=xy00(2,:);
[qx1,qy1] = evaluateHeatFlux(results,x00,y00,1:length(t2));

figure(7);
plot(x00,'o','MarkerSize',2);
title({'Nodes along x-Coordinate';''},'FontWeight','normal');
grid on; grid minor; xlabel ('Node number'); ylabel ('x (m)'); 
ax = gca; ax.FontSize = 14;

figure(8);
plot(y00,'o','MarkerSize',2);
title({'Nodes along y-Coordinate';''},'FontWeight','normal');
grid on; grid minor; xlabel ('Node number'); ylabel ('y (m)'); 
ax = gca; ax.FontSize = 14;

figure(9);
plot(x00,y00,'o','MarkerSize',2);
title({'Nodes at Given Coordinates';''},'FontWeight','normal');
grid on; grid minor; xlabel ('x (m)'); ylabel ('y (m)'); 
ax = gca; ax.FontSize = 14;
%**********************************************************

nodesLeftEnd_v = find(x00 == 0) % identify far left boundary nodes
nodesCenter_v = find(and(0.5*width - eps <= x00,x00 <= 0.5*width + eps)) %identify far left boundary nodes
nodesInterior_h = find(and(height1 - eps <= y00, y00 <= height1 + eps)) %identify far left boundary nodes
nodesExterior_h = find(height2 - eps <= y00) % identify far left boundary nodes
xy1=point(:,nodesLeftEnd_v); x1=xy1(1,:); y1=xy1(2,:)
xy2=point(:,nodesCenter_v); x2=xy2(1,:); y2=xy2(2,:)
xy3=point(:,nodesInterior_h); x3=xy3(1,:); y3=xy3(2,:)
xy4=point(:,nodesExterior_h); x4=xy4(1,:); y4=xy4(2,:)
nodesInterior_1 = find(and(- eps <= y2, y2 <= + eps)) % identify far left boundary nodes
nodesInterior_2 = find(and(height1 - eps <= y2, y2 <= height1 + eps)) %identify far left boundary nodes
nodesExterior = find(and(height2 - eps <= y2, y2 <= height2 + eps)) %identify far left boundary nodes
xy5=xy2(:,nodesInterior_1); x5=xy5(1,:); y5=xy5(2,:)
xy6=xy2(:,nodesInterior_2); x6=xy6(1,:); y6=xy6(2,:)
xy7=xy2(:,nodesExterior); x7=xy7(1,:); y7=xy7(2,:)

figure(10);
hold all;
plot(x1,y1,'o','MarkerSize',5,'MarkerIndices',1:1:length(x1),'DisplayName',...
    ['y = ', num2str(y1(:,1))]); %plot central nodes
plot(x2,y2,'d','MarkerSize',5,'MarkerIndices',1:1:length(x2),'DisplayName',...
    ['y = ', num2str(y2(:,1))]); %plot central nodes
legend('Far-left nodes','Central nodes');
xlim([0 width]); ylim([0 height2]); 
title({'Vertical Nodes at Selected Coordinates';''},'FontWeight','normal');
grid on; grid minor; xlabel ('x (m)'); ylabel ('y (m)'); 
ax = gca; ax.FontSize = 14;

figure(11);
hold all;
plot(x3,y3,'o','MarkerSize',5,'MarkerIndices',1:1:length(x3),'DisplayName',...
    ['y = ',num2str(y3(:,1))]); %plot central nodes
plot(x4,y4,'d','MarkerSize',5,'MarkerIndices',1:1:length(x4),'DisplayName',...
    ['y = ',num2str(y4(:,1))]); %plot central nodes
legend('Interior nodes','Exterior nodes','Location','South');
xlim([0 width]); ylim([0 height2]); 
title({'Horizontal Nodes at Selected Coordinates';''},'FontWeight','normal');
grid on; grid minor; xlabel ('x (m)'); ylabel ('y (m)'); 
ax = gca; ax.FontSize = 14;

figure(12);
hold all;
cc = (1:2:length(x2));
for i = 1:2:length(x2)
    p2(i) = plot(x2(i),y2(i),'o','MarkerSize',8,'MarkerFaceColor',[0.5 0.5 0.5],...
        'DisplayName',['y = ',num2str(y2(i)),' m']); %plot central nodes    
    switch i
       case cc(1); p2(i).Marker = 's';  case cc(2); p2(i).Marker = '<';
       case cc(3); p2(i).Marker = '+';  case cc(4); p2(i).Marker = 'p';
       case cc(5); p2(i).Marker = '^';  case cc(6); p2(i).Marker = '>';
       case cc(7); p2(i).Marker = 'X';  case cc(8); p2(i).Marker = 'd';
        otherwise; p2(i).Marker = '*';
    end          
    legend;         
end;     
xlim([0 width]); ylim([0 height2]);
title({'Selected Central Nodes at Given Coordinates';''},'FontWeight','normal');
grid on; grid minor; xlabel('x (m)'); ylabel('y (m)');
ax = gca; ax.FontSize = 14;

figure(13);
hold all;
plot(x5,y5,'o','MarkerSize',7,'MarkerFaceColor',[0.5 0.5 0.5],'DisplayName',...
    ['y = ',num2str(y5),' m (Bottom node)']); %plot central nodes
plot(x6,y6,'s','MarkerSize',7,'MarkerFaceColor',[0.5 0.5 0.5],'DisplayName',...
    ['y = ',num2str(y6),' m (Interface node)']); %plot interior nodes
plot(x7,y7,'d','MarkerSize',7,'MarkerFaceColor',[0.5 0.5 0.5],'DisplayName',...
    ['y = ',num2str(y7),' m (Exterior node)']); %plot exterior nodes
xlim([0 width]); ylim([0 height2]); 
title({'Selected Central Nodes';''},'FontWeight','normal');
grid on; grid minor; xlabel ('x (m)'); ylabel ('y (m)'); legend('Location','West');
ax = gca; ax.FontSize = 14;

figure(14);
hold all;
n = 35;
plot(t2(1:n),T1(nodesInterior_1,[1:n]),'--+','MarkerSize',5,'MarkerIndices',...
    1:5:length(t2(1:n)),'MarkerFaceColor',[0.5 0.5 0.5]);
plot(t2(1:n),T1(nodesInterior_2,[1:n]),'--^','MarkerSize',5,'MarkerIndices',...
    1:7:length(t2(1:n)),'MarkerFaceColor',[0.5 0.5 0.5]);
plot(t2(1:n),T1(nodesExterior,[1:n]),'-*','MarkerSize',5,'MarkerIndices',...
    1:9:length(t2(1:n)),'MarkerFaceColor',[0.5 0.5 0.5]);
xlim([0 t2(1,n)]);
title({'Transient Temperature at Selected Nodes';''},'FontWeight','normal');
grid on; grid minor; xlabel 't (min)'; ylabel 'Temperature (degC)'; 
legend('Bottom node','Interface node','Exterior node','Location','southeast');
ax = gca; ax.FontSize = 14;

figure(15);
hold all;
Tintrp1 = interpolateTemperature(results,x1,y1,1:length(t2));
Tintrp2 = interpolateTemperature(results,x2,y2,1:length(t2));
dt1 = 0.5*time_step/60;
for i = 1:4
        timei1(i) = find(and(7.5*i - dt1 <= t2,t2 <= 7.5*i + dt1))
        time1(i) = round( t2(timei1(i)),1)
        %p1(i) = plot(y1,Tintrp1(:,i),'o','MarkerSize',5,'MarkerIndices',1:1:length(x1),...
        %   'MarkerFaceColor',[0.5 0.5 0.5],'DisplayName',['t = ',strcat(' ',num2str(time1(i))), ' min',...
        %  ', x1 = ', num2str(x1(i)), ' m']);
        p2(i) = plot(y2,Tintrp2(:,i),'o','MarkerSize',5,'MarkerIndices',1:1:length(x2),...
           'MarkerFaceColor',[0.5 0.5 0.5],'DisplayName',['t = ',strcat(' ',...
           num2str(time1(i))), ' min',', x2 = ', num2str(x2(i)), ' m']);
        legend('Location','southwest');
         switch i
     case 1;     p2(i).Marker = 's'; %p1(i).Marker = 'd';
     case 2;     p2(i).Marker = '<'; %p1(i).Marker = '>'; 
     case 3;     p2(i).Marker = '+'; %p1(i).Marker = '*'; 
     case 4;     p2(i).Marker = 'X'; %p1(i).Marker = 'h'; 
     case 5;     p2(i).Marker = 'o'; %p1(i).Marker = 'p';
     otherwise;  p2(i).Marker = '*'; %p1(i).Marker = '*';
         end
    hold on
end
title({'Temperature at Selected Interior Nodes along y-Coordinate',''},'FontWeight','normal');
grid on; grid minor; xlabel('y (m)'); ylabel('T (degC)');
ax = gca; ax.FontSize = 14;

figure(16);
hold all;
Tintrp3 = interpolateTemperature(results,x3,y3,1:length(t2));
Tintrp4 = interpolateTemperature(results,x4,y4,1:length(t2));
dt1 = 0.5*time_step/60;
for i = 1:4
        timei2(i) = find(and(7.5*i - dt1 <= t2,t2 <= 7.5*i + dt1))
        time2(i) = round( t2(timei2(i)),1)
        p5(i) = plot(x3,Tintrp3(:,i),'o','MarkerSize',5,'MarkerIndices',1:2:length(x3),...
            'MarkerFaceColor',[0.5 0.5 0.5],'DisplayName',['t = ',strcat(' ',...
            num2str(time2(i))), ' min',', y3 = ', num2str(y3(i)), ' m',]);
        p6(i) = plot(x4,Tintrp4(:,i),'o','MarkerSize',5,'MarkerIndices',1:2:length(x4),...
            'MarkerFaceColor',[0.5 0.5 0.5],'DisplayName',['t = ',strcat(' ',...
            num2str(time2(i))), ' min',', y4 = ', num2str(y4(i)), ' m',]);
        legend; %('Location','east');
        switch i
    case 1;     p5(i).Marker = 's'; p6(i).Marker = 'd';
    case 2;     p5(i).Marker = '<'; p6(i).Marker = '>'; 
    case 3;     p5(i).Marker = '+'; p6(i).Marker = '*'; 
    case 4;     p5(i).Marker = 'X'; p6(i).Marker = 'h'; 
    case 5;     p5(i).Marker = 'o'; p6(i).Marker = 'p';
    otherwise;  p5(i).Marker = '*'; p6(i).Marker = '*';
        end
    hold on
end
title({'Temperature at Selected Interior Nodes along x-Coordinate',''},'FontWeight','normal');
grid on; grid minor; xlabel('x (m)'); ylabel('T (degC)');
ax = gca; ax.FontSize = 14;

%title({'Arrow Direction and Length at Each Data Point';
%'(x- and y- coordinate) represents the';' associated values for qx and qy';
%'heat flux'});
figure(17)
dd = num2str(t2(end));
hold all; 
[qx2,qy2] = evaluateHeatFlux(results,x2,y2,1:length(t2));
quiver(x2(:),y2(:),qx2(:,end),qy2(:,end),'MarkerSize',2,'LineWidth',1,'Color','black',...
    'AutoScale','on','AutoScaleFactor',2.5,'ShowArrowHead','on');
axis equal; xlim([0 width]); ylim([0.015 height2]);
title({'Heat Flux Vectors at Selected Nodes (W/m^2)'; ['(x = ',num2str(x2(1)),...
    ' m, ','t = ',dd,' min)'];''},'FontWeight','normal')
grid on; grid minor; xlabel 'x (m)'; ylabel 'y (m)';
ax = gca; ax.FontSize = 14;

figure(18)
hold all;
[qx3,qy3] = evaluateHeatFlux(results,x3,y3,1:length(t2));
quiver(x3(:),y3(:),qx3(:,end),qy3(:,end),'MarkerSize',2,'LineWidth',1,'Color','black',...
    'AutoScale','on','AutoScaleFactor',2.5,'ShowArrowHead','on');
axis equal; xlim([0 width]); ylim([0.02 height2]);
title({'Heat Flux Vectors at Selected Nodes (W/m^2)'; ['(y = ',num2str(y3(1)),...
    ' m, ','t = ',dd,' min)'];''},'FontWeight','normal')
grid on; grid minor; xlabel 'x (m)'; ylabel 'y (m)';
ax = gca; ax.FontSize = 14;

figure(19)
[qx4,qy4] = evaluateHeatFlux(results,x4,y4,1:length(t2));
quiver(x4(:),y4(:),qx4(:,end),qy4(:,end),'MarkerSize',2,'LineWidth',1,'Color','black',...
    'AutoScale','on','AutoScaleFactor',2.5,'ShowArrowHead','on');
axis equal; xlim([width/3 2*width/3]); ylim([height1 height2]);
title({'Heat Flux Vectors at Selected Nodes (W/m^2)'; ['(y = ',num2str(y4(1)),...
    ' m, ','t = ',dd,' min)'];''},'FontWeight','normal')
grid on; grid minor; xlabel 'x (m)'; ylabel 'y (m)';
ax = gca; ax.FontSize = 14;

figure(20)
hold all;
n = 20;
Qn_interior1 = evaluateHeatRate(results,'Edge',1);
Qn_interior2 = evaluateHeatRate(results,'Edge',3);
Qn_exterior = evaluateHeatRate(results,'Edge',2);
%plot(t2(1:n),Qn_interior1(1:n),'-d','MarkerSize',5,'MarkerIndices',1:4:length(t2(1:n)));
plot(t2(1:n),Qn_interior2(1:n),'-s','MarkerSize',5,'MarkerIndices',1:5:length(t2(1:n)));
%plot(t2(1:n),Qn_exterior(1:n),'-p','MarkerSize',5,'MarkerIndices',1:6:length(t2(1:n)));
xlim([0 t2(:,n)]);
title({'Heat Flow Rate Across Edges';''},'FontWeight','normal'); 'LineWidth'==1; 
legend('Edge 3','Location','SouthEast');
grid on; grid minor; grid minor; xlabel('t (min)'); ylabel('Heat flow rate (W)'); 
ax = gca; ax.FontSize = 14;

figure(21)
hold all;
Qn_EdgeLeftB = evaluateHeatRate(results,'Edge',6);
Qn_EdgeLeftT = evaluateHeatRate(results,'Edge',7);
Qn_EdgeRightB = evaluateHeatRate(results,'Edge',4);
Qn_EdgeRightT = evaluateHeatRate(results,'Edge',5);
plot(t2,Qn_EdgeLeftB,':d','MarkerSize',5,'MarkerIndices',1:2:length(t2));
hold all;
plot(t2,Qn_EdgeLeftT,'--s','MarkerSize',5,'MarkerIndices',1:3:length(t2));
plot(t2,Qn_EdgeRightB,'-p','MarkerSize',5,'MarkerIndices',1:4:length(t2));
plot(t2,Qn_EdgeRightT,':X','MarkerSize',5,'MarkerIndices',1:5:length(t2));
legend('Edge 6','Edge 7','Edge 4','Edge 5','Location','SouthEast');
xlim([0 t2(:,10)]);
title({'Heat Flow Rate across Edges';''},'FontWeight','normal');'LineWidth'==1; 
grid on; grid minor; xlabel('t (min)'); ylabel('Heat flow rate (W)'); 
ax = gca; ax.FontSize = 14;

[Tx2,Ty2] = evaluateTemperatureGradient(results,x2,y2,1:length(t2));
figure(22);
hold all;
cc = [nodesInterior_1 nodesInterior_2 nodesExterior];
for i = [nodesInterior_1 nodesInterior_2 nodesExterior]
p1(i) = plot(t2,Tx2(i,:),'-','MarkerSize',5,'MarkerFaceColor',[0.5 0.5 0.5],...
    'DisplayName',['y = ',strcat(' ',num2str(y2(i))), ' m']);
 legend('Location','East');
 switch i
    case cc(1); p1(i).Marker = 's'; 
    case cc(2); p1(i).Marker = '<'; 
    case cc(3); p1(i).Marker = 'd';     
    otherwise;  p1(i).Marker = '*'; 
 end    
end
xlim([0 t2(:,14)]);
title({'Transient Temperature Gradients with respect to';[' x-Coordinate (dT/dx) (x = ',...
    strcat(' ',num2str(x2(i))), ' m)'];''},'FontWeight','normal');
grid on; grid minor; xlabel('t (min)'); ylabel('dT/dx (degC/m)');
legend('Bottom node','Interface node','Exterior node','Location','east');
ax = gca; ax.FontSize = 14;

figure(23);
hold all;
for i = [nodesInterior_1 nodesInterior_2 nodesExterior]
p1(i) = plot(t2,Ty2(i,:),'-','MarkerSize',5,'MarkerFaceColor',[0.5 0.5 0.5],...
    'DisplayName',['y = ',strcat(' ',num2str(y2(i))), ' m']);
 legend('Location','SouthEast');
 switch i
    case cc(1); p1(i).Marker = 's'; 
    case cc(2); p1(i).Marker = '<'; 
    case cc(3); p1(i).Marker = 'd';     
    otherwise;  p1(i).Marker = '*'; 
  end    
end
xlim([0 t2(:,5)]);
title({'Transient Temperature Gradients with respect to';[' y-Coordinate (dT/dy) (x = ',strcat(' ',...
    num2str(x2(i))), ' m)'];''},'FontWeight','normal');
grid on; grid minor; xlabel('t (min)'); ylabel('dT/dy (degC/m)');
legend('Bottom node','Interface node','Exterior node','Location','east');
ax = gca; ax.FontSize = 14;

[Tx3,Ty3] = evaluateTemperatureGradient(results,x3,y3,1:length(t2));
figure(24)
hold all;
t11 = [100:time_step*1:length(t2)];
for i = t11
    pxx(i) = plot(x3,Tx3(:,i),'o','MarkerSize',4,'MarkerFaceColor',[0.5 0.5 0.5],...
        'DisplayName',['t (Tx) = ', num2str(t2(i))]);    
    pxy(i) = plot(x3,Ty3(:,i),'<','MarkerSize',4,'MarkerFaceColor',[0.5 0.5 0.5],...
        'DisplayName',['t (Ty) = ', num2str(t2(i))]);
    legend;
end;
pxx(t11(1)).Marker = 's'; pxx(t11(2)).Marker = 'X'; %pxx(t(3)).Marker = '+';
pxy(t11(1)).Marker = 'd'; pxy(t11(2)).Marker = '>'; %pxy(t(3)).Marker = '*';
title({'Spatial Temperature Gradients (Tx and Ty)';'along x-Coordinate';''},'FontWeight','normal');
grid on; grid minor; xlabel('x (m)'); ylabel('Temperature gradient (degC/m)');
ax = gca; ax.FontSize = 14;

figure(25)
hold all;
t22 = [100:time_step*1:length(t2)];
for i = t22
    pyx(i) = plot(y2,Tx2(:,i),'o','MarkerSize',4,'MarkerFaceColor',[0.5 0.5 0.5],...
        'DisplayName',['t (Tx) = ', num2str(t2(i))]);    
    pyy(i) = plot(y2,Ty2(:,i),'<','MarkerSize',4,'MarkerFaceColor',[0.5 0.5 0.5],...
        'DisplayName',['t (Ty) = ', num2str(t2(i))]);
    legend('Location','East');
end;
pyx(t22(1)).Marker = 's'; pyx(t22(2)).Marker = 'X'; %pyx(t(3)).Marker = '+';
pyy(t22(1)).Marker = 'd'; pyy(t22(2)).Marker = '>'; %pyy(t(3)).Marker = '*';

title({'Spatial Temperature Gradients (Tx and Ty)';'along x-Coordinate';''},'FontWeight','normal');
grid on; grid minor; xlabel('y (m)'); ylabel('Temperature gradient (degC/m)');
ax = gca; ax.FontSize = 14;

tlist = 0.1:time_step:tfinal;
t2 = tlist/60;
length1 = height2 + width;
w = linspace(0,width,100);
h = linspace(0,height2,100);
[Xw,Yh] = meshgrid(w,h);
T_xyz = interpolateTemperature(results,Xw,Yh,1:length(t2));
[Tx,Ty] = evaluateTemperatureGradient(results,Xw,Yh,1:length(t2));
[qx,qy] = evaluateHeatFlux(results,Xw,Yh,1:length(t2));
Qn = evaluateHeatRate(results,'Edge',[1 2]);

figure(26)
plot(Xw,Yh,'MarkerSize',6,'MarkerFaceColor',[0.5 0.5 0.5],'DisplayName','Xw');
title({'Vector Coordinates',''},'FontWeight','normal'); 
grid on; grid minor; xlabel('x (m)'); ylabel('y (m)');
ax = gca; ax.FontSize = 14;

figure(27)
for ii = [nodesInterior_1 nodesInterior_2 nodesExterior]
p5(ii) = quiver(x2(:),y2(:),Tx2(:,end),Ty2(:,end),'MarkerSize',3,'MarkerFaceColor',[0.5 0.5 0.5],...
    'DisplayName',[' y = ',num2str(y2(ii)),' m'],'MarkerSize',1,'LineWidth',1,'Color','black',...
    'AutoScale','on','AutoScaleFactor',2.5,'ShowArrowHead','on');
legend; hold all;
end;
xlim([0 width]); ylim([0 height2]);
title({'Temperature Gradient Vectors at Selected Nodes (degC/m)';['(x = ',num2str(x2(ii)),' m, ',...
    't = ',dd,' min)'];''},'FontWeight','normal');
grid on; grid minor; xlabel('x (m)'); ylabel('y (m)');
legend('Bottom node','Interface node','Exterior node','Location','east');
ax = gca; ax.FontSize = 14;

figure(28)
cc = [nodesInterior_1 nodesInterior_2 nodesExterior];
for ii = [nodesInterior_1 nodesInterior_2 nodesExterior]
p5(ii) = quiver(x3(:),y3(:),Tx3(:,end),Ty3(:,end),'MarkerSize',3,'MarkerFaceColor',[0.5 0.5 0.5],...
    'DisplayName',['y = ',num2str(y3(ii)),' m'],'MarkerSize',1,'LineWidth',1,'Color','black',...
    'AutoScale','on','AutoScaleFactor',2.5,'ShowArrowHead','on');
end;
title({'Temperature Gradient Vectors at Selected Nodes (degC/m)';['(y = ',num2str(y3(ii)),' m, ',...
    't = ',dd,' min)'];''},'FontWeight','normal');
xlim([0 width]); ylim([height1 height2+0.003]);
grid on; grid minor; xlabel('x (m)'); ylabel('y (m)');
ax = gca; ax.FontSize = 14;

figure(29)
[Tx4,Ty4] = evaluateTemperatureGradient(results,x4,y4,1:length(t2));
for ii = [nodesInterior_1 nodesInterior_2 nodesExterior]
p5(ii) = quiver(x4(:),y4(:),Tx4(:,end),Ty4(:,end),'MarkerSize',3,'MarkerFaceColor',[0.5 0.5 0.5],...
    'DisplayName',[' y = ',num2str(y4(ii)),' m'],'MarkerSize',1,'LineWidth',1,'Color','black',...
    'AutoScale','on','AutoScaleFactor',2.5,'ShowArrowHead','on');
end;
title({'Temperature Gradient Vectors at Selected Nodes (degC/m)';['(y = ',num2str(y4(ii)),' m, ',...
    't = ',dd,' min)'];''},'FontWeight','normal');
xlim([0 width]); ylim([0.029 height2+0.003]);
grid on; grid minor; xlabel('x (m)'); ylabel('y (m)');
ax = gca; ax.FontSize = 14;

figure(30);
h = pdeplot(thermalModelT,'XYData',T1,'ZData',T1, ...
                   'FaceAlpha',0, ...
                   'FlowData',[Tx1,Ty1], ...
                   'ColorMap','jet', ...
                   'Mesh','on')
xlim([0 width]); ylim([0 height2]);
title({'Flow Data Contours at Last Tiime Step';['(t = ',num2str(tt),' hr)'];''},'FontWeight','normal')
grid on; grid minor; xlabel('x (m)'); ylabel('y (m)'); zlabel('T (degC)');
ax = gca; ax.FontSize = 14;

figure(31);
ax.XAxis.FontSize=20; ax.YAxis.FontSize=20;
plot(T_st1/3600,'-o','MarkerSize',7,'MarkerIndices',1:7:length(T_st1),...
    'MarkerFaceColor',[0.5 0.5 0.5]); %plot solution time versus iteration
title('Solution Time','FontWeight','normal','FontSize',14);
grid on; grid minor; xlabel('Solution step number','FontSize',14); ylabel('t (hr)','FontSize',14);
ax = gca; ax.FontSize = 14;

bc = findThermalBC(thermalModelT.BoundaryConditions,'Edge',6)
ic = findThermalIC(thermalModelT.InitialConditions,'Face',[1 2])

%Specify the number of colors used in a colormap by passing an integer
%as an input argument to the built-in colormap. Use five colors from the
%parula colormap.
figure(32); 
pdeplot(thermalModelT,'XYData',T1(:,end),'ColorMap',hot(5),'Contour','on','Levels',20,...
                     'FlowData',[qx1(:,end),qy1(:,end)],'FlowStyle','arrow');
xlim([0 width]); ylim([0 height2]);
title({'Flow Data Contours at Last Time Step';['(t = ',num2str(tt),' hr)'];''},'FontWeight','normal');
grid on; grid minor; xlabel('x (m)'); ylabel('y (m)'); zlabel('T (degC)');
ax = gca; ax.FontSize = 14;
