% 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 used ThemalConductivity, MassDensity and SpecificHeat. 
% Geometry Pipe_hollow_shortened.stl is imported to this model and then scaled to the desired
% portions.


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'; 
  
thermalModelT = createpde('thermal','transient'); % create a model object: thermal
geom=importGeometry(thermalModelT,'Pipe_hollow_shortened.stl')
scalef = 1000;
%scale(geom,[1 1 1]/scalef);
scale(geom,[0.5 1 1.25]/1000);

figure(1);
pdegplot(thermalModelT,'FaceLabels','on','FaceAlpha',0.55,'CellLabels','on')
title({'Geometry with Face and Cell Labels Displayed'},'FontWeight','normal');
grid on; grid minor; xlabel('x (m)'); ylabel('y (m)');
ax = gca; ax.FontSize = 14;
hgexport(figure(1),'-clipboard')

%Aluminum properties
k = 210; % thermal conductivity (W/mK)
rho = 2698.9; %density (kg/m3)
cp = 900; % heat capacity (J/kgK)
q = 1000; %heat generation (W/m^3)
kFunc = @(region,state) k*region.y;
cFunc = @(region,state) cp*region.y;
rhoFunc = @(region,state) rho*region.y;
thermalProperties(thermalModelT,'ThermalConductivity',k,...
    'MassDensity',rho,'SpecificHeat',cp,'Cell',1);
%Optional: Location / Temperature dependent properties
%a, b, c, and d are %constants. u is state vatiable (tempeature)
%kFunc = @(location,state)a+b*state.u+c*state.u^d; 
%cFunc = @(location,state)a*location.x+b*location.y+c*location.z;
%q =@(location,state)state.t*(a*location.x+b)*(state.u)^c;

A = 0./[0 1 2 3 4];
B = 1./[0 1 2 3 4];
isnan(A);
isnan(B);
if(or(isnan(A),isnan(B)))
    qq = NaN;
return
end


% Heat Source (W/m3)
q = 1000;
qFunc = @(region,state) q*region.y;
ih=internalHeatSource(thermalModelT,qFunc,'Cell',1);

% Boundary Conditions (W/m2)
T_ambient = 22; %degC **********
T_water = 50;   %degC **********
qf = 100; %heat flux (W/m^2)
hc = 10; %heat convectyion coefficient (W/m^2K)
outerCC_V = @(region,~) hc*region.y; %horizontal heat convection coefficient (W/m2K)
top_BC_HF = @(region,~) qf*region.y; %heat source distribution (W/m2)
thermalBC(thermalModelT,'Face',[13 14],'HeatFlux',top_BC_HF); %boundary condition: heat source (W/m2)
thermalBC(thermalModelT,'Face',[15:26],'Temperature',T_water); %boundary condition: temperature
thermalBC(thermalModelT,'Face',[1:12],'ConvectionCoefficient',outerCC_V,'AmbientTemperature',T_ambient); %boundary condition: convection heat transfer

% Initial Condition
T0_1 = 25; %Initial tempeature of Face 1 (water) **********
IC1=thermalIC(thermalModelT,T0_1,'Cell',1)

figure(2);
grid_size = 7.5/scalef; %identify grid size (m)
mesh = generateMesh(thermalModelT,'Hmax',grid_size,'GeometricOrder','linear'); %generate mesh *****
pdeplot3D(mesh);
title({'Tetrahedral Elements'},'FontWeight','normal');
grid on; grid minor; xlabel('x (m)'); ylabel('y (m)');
ax = gca; ax.FontSize = 10; view(0,90);

figure(3);
pdeplot3D(mesh,'ElementLabels','on');
title({'Mesh Front View (0,0)';''},'FontWeight','normal');
grid on; grid minor; xlabel('x (m)'); ylabel('y (m)');
ax = gca; ax.FontSize = 14; view(0,0);
xlim([1*grid_size 3.5*grid_size]); zlim([1*grid_size 2.5*grid_size]); ylim([2*grid_size 3*grid_size]); 

figure(4);
pdeplot3D(mesh); 
%title({'Mesh: Top View';''},'FontWeight','normal');
grid on; grid minor; xlabel('x (m)'); ylabel('y (m)');
ax = gca; ax.FontSize = 14; view(-37.5,30);
%xlim([0*grid_size 3*grid_size]); ylim([0.5*grid_size 7.5*grid_size]); zlim([2*grid_size 7*grid_size]); % set axis conditions

% Solution settings: identify thermal data to process and plot *****
tfinal = 20*60; % set final time (s)
time_step = 10; % time step (s)
tlist = 0.1:time_step:tfinal; % set time array
t2 = tlist/60;
thermalModelT.SolverOptions.ReportStatistics = 'on'; % identify PDESolverOptions Properties
results = solve(thermalModelT,tlist); 
T1 = results.Temperature; % identify tempeature data for selected nodes *****
Tx1 = results.XGradients; % identify temperature gradients with respect to x-coordinate *****
Ty1 = results.YGradients; % identify temperature  gradients with respect to y-coordinate *****
Tz1 = results.ZGradients; % identify temperature  gradients with respect to y-coordinate *****
T_st1 = results.SolutionTimes; % identify results solution times *****
tt=round(tfinal/60,1); % identify final time step in min

figure(5);
mm = 1;
pdeplot3D(thermalModelT,'ColorMapData',T1(:,mm));
title({'Temperature Contours (degC)'; ['(t = ',num2str(t2(mm)),' s)'];''},'FontWeight','normal')
grid on; grid minor; xlabel('x (m)'); ylabel('y (m)');
ax = gca; ax.FontSize = 14;

figure(6);
nn = length(t2)/5;
pdeplot3D(thermalModelT,'ColorMapData',T1(:,nn));
title({'Temperature Contours (degC)'; ['(t = ',num2str(t2(nn)),' s)'];''},'FontWeight','normal')
grid on; grid minor; xlabel('x (m)'); ylabel('y (m)');
ax = gca; ax.FontSize = 14;

figure(7);
oo = length(t2)/2;
pdeplot3D(thermalModelT,'ColorMapData',Tx1(:,oo))
title({'Temperature Gradient Contours (dT/dx)'; ['(degC/m) (t = ',num2str(t2(oo)),' min)'];''},'FontWeight','normal')
grid on; grid minor; xlabel('x (m)'); ylabel('y (m)');
ax = gca; ax.FontSize = 14;

figure(8);
pdeplot3D(thermalModelT,'ColorMapData',Ty1(:,oo)/2)
title({'Temperature Gradient Contours (dT/dy)'; ['(degC/m) (t = ',num2str(t2(oo)),' min)'];''},'FontWeight','normal')
grid on; grid minor; xlabel('x (m)'); ylabel('y (m)');
ax = gca; ax.FontSize = 14;

%**********************************************************
%Investigate node density
point = thermalModelT.Mesh.Nodes; % recall node locations from thermal model object
eps1=eps;
nodesIndex = find(point); % identify all nodes
xy00=point; x00=xy00(1,:); y00=xy00(2,:); z00=xy00(3,:);
width = max(x00)-min(x00); height = max(z00)-min(z00); length2 = max(y00)-min(y00);
nodesLeftEnd_v = find(x00 == min(x00)); % identify far left boundary nodes, x = minx(x)
nodesCenter_v = find(and(0.5*width - 0.15*grid_size <= x00,x00 <= 0.5*width + 0.15*grid_size)); % identify central nodes, x = width/2
nodesTop_h = find(and(max(z00) - eps <= z00, z00 <= max(z00) + eps)); % identify top nodes, z = height
nodesBottom_h = find(and(min(z00) - eps <= z00, z00 <= min(z00) + eps)); % identify bottom nodes, z = min(z)
xy1=point(:,nodesLeftEnd_v); x1=xy1(1,:); y1=xy1(2,:); z1=xy1(3,:);
xy2=point(:,nodesCenter_v); x2=xy2(1,:); y2=xy2(2,:); z2=xy2(3,:);
xy3=point(:,nodesTop_h); x3=xy3(1,:); y3=xy3(2,:); z3=xy3(3,:);
xy4=point(:,nodesBottom_h); x4=xy4(1,:); y4=xy4(2,:); z4=xy4(3,:);
nodes_Bottom = find(and(z2<=min(z00)+eps,z2>=min(z00)-eps)); 
nodes_Top = find(and(z2<=max(z00)+eps,z2>=max(z00)-eps)); 
xy5=xy2(:,nodes_Bottom); x5=xy5(1,:); y5=xy5(2,:); z5=xy5(3,:);
xy6=xy2(:,nodes_Top); x6=xy6(1,:); y6=xy6(2,:); z6=xy6(3,:);

figure(9);
plot(x00,z00,'o','MarkerSize',2);
title({'Projected Nodes at Given Coordinates';''},'FontWeight','normal');
xlim([0 max(x00)]); ylim([0 max(z00)]);
grid on; grid minor; xlabel ('x (m)'); ylabel ('z (m)'); 
ax = gca; ax.FontSize = 14;

figure(10);
plot(x00,y00,'o','MarkerSize',2);
title({'Projected Nodes at Given Coordinates';''},'FontWeight','normal');
xlim([0 max(x00)]); ylim([0 max(y00)]); 
grid on; grid minor; xlabel ('x (m)'); ylabel ('y (m)'); 
ax = gca; ax.FontSize = 14;

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

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

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

%**********************************************************

figure(14);
hold all;
plot(x2,z2,'d','MarkerSize',5,'MarkerFaceColor',[0.5 0.5 0.5],'DisplayName',['y <> ',num2str(y2(:,1)),' m']);
legend;
xlim([0 max(x00)]); ylim([0 max(z00)]);
title({'Vertical Nodes at Selected Coordinates';''},'FontWeight','normal');
grid on; grid minor; xlabel ('x (m)'); ylabel ('z (m)'); 
ax = gca; ax.FontSize = 14;

figure(15);
hold all;
plot(y2,z2,'d','MarkerSize',5,'MarkerFaceColor',[0.5 0.5 0.],'DisplayName',['x <> ',num2str(x2(:,1)),' m']);
legend;
xlim([0 max(x00)]); ylim([0 max(y00)]);
title({'Vertical Nodes at Selected Coordinates';''},'FontWeight','normal');
grid on; grid minor; xlabel ('y (m)'); ylabel ('z (m)'); 
ax = gca; ax.FontSize = 14;

figure(16);
hold all;
plot(x4,z4,'d','MarkerSize',5,'MarkerIndices',1:2:length(x4),'DisplayName',[num2str(y4(:,1)),' < y > ',num2str(y4(:,end)),' m']);
plot(x3,z3,'o','MarkerSize',5,'MarkerIndices',1:2:length(x3),'DisplayName',[num2str(y3(:,1)),' < y > ',num2str(y3(:,end)),' m']);
legend('Location','East');
xlim([0 max(x00)]);
title({'Horizontal Nodes at Selected Coordinates';''},'FontWeight','normal');
grid on; grid minor; xlabel ('x (m)'); ylabel ('z (m)'); 
ax = gca; ax.FontSize = 14;

figure(17);
hold all;
plot(x4,y4,'X','MarkerSize',4,'MarkerIndices',1:1:length(x4),'DisplayName',['z = ',num2str(z4(:,1)),' m']);
plot(x3,y3,'d','MarkerSize',6,'MarkerIndices',1:1:length(x3),'DisplayName',['z = ',num2str(z3(:,1)),' m']);
legend;
xlim([0 max(x00)]);
title({'Horizontal Nodes at Selected Coordinates';''},'FontWeight','normal');
grid on; grid minor; xlabel ('x (m)'); ylabel ('y (m)'); 
ax = gca; ax.FontSize = 14;

figure(18);
plot3(x2,y2,z2,'s','MarkerSize',5,'MarkerIndices',1:1:length(x2));
xlim([0 max(x00)]);
title({'Horizontal Nodes at Selected Coordinates';''},'FontWeight','normal');
grid on; grid minor; xlabel ('x (m)'); ylabel ('y (m)'); zlabel ('z (m)');
ax = gca; ax.FontSize = 14;

figure(19);
plot3(x4,y4,z4,'s','MarkerSize',5,'MarkerIndices',1:1:length(x4));
hold all;
plot3(x3,y3,z3,'d','MarkerSize',5,'MarkerIndices',1:1:length(x3));
xlim([0 max(x00)]);
title({'Horizontal Nodes at Selected Coordinates';''},'FontWeight','normal');
grid on; grid minor; xlabel ('x (m)'); ylabel ('y (m)'); zlabel ('z (m)');
ax = gca; ax.FontSize = 14;

figure(20);
hold all;
n = 10;
plot(t2(1:n),T1(nodes_Bottom(1),1:n),'-X','MarkerSize',4,'MarkerIndices',1:1:length(t2(1:n)),...
    'DisplayName',['z = ', num2str(z5(1)),' m']);
plot(t2(1:n),T1(nodes_Top(1),1:n),'--d','MarkerSize',6,'MarkerIndices',1:1:length(t2(1:n)),...
    'DisplayName',['z = ', num2str(z6(1)),' m']);
legend('Location','NorthEast');
xlim([0 t2(1,n)]);
title({'Transient Temperature at Selected Nodes';['(x = ', num2str(x5(1)),' m, y = ',...
    num2str(y5(1)),' m)'];''},'FontWeight','normal');
grid on; grid minor; xlabel 't (min)'; ylabel 'T (degC)'; 
ax = gca; ax.FontSize = 14;

figure(21);
hold all;
n = 10;
plot(t2(1:n),Tx1(nodes_Bottom(1),1:n),'-X','MarkerSize',4,'MarkerIndices',1:1:length(t2(1:n)),...
    'DisplayName',['z = ', num2str(z5(2)),' m']);
plot(t2(1:n),Ty1(nodes_Top(1),1:n),'--d','MarkerSize',6,'MarkerIndices',1:1:length(t2(1:n)),...
    'DisplayName',['z = ', num2str(z6(2)),' m']);
legend('Location','NorthEast');
xlim([0 t2(1,n)]);
title({'Transient Temperature Gradients at Selected Nodes';['(x = ',num2str(x5(2)),' m, y = ',...
    num2str(y5(1)),' m)'];''},'FontWeight','normal');
grid on; grid minor; xlabel 't (min)'; ylabel 'dT/d(x,y) (degC)'; 
ax = gca; ax.FontSize = 14;

figure(22);
hold all;
Tintrp2 = interpolateTemperature(results,x2,y2,z2,1:length(t2));
dt1 = 0.5*time_step/60;
for j=[1 length(t2)/10]
        i = 1*j;
        time2(j) = round(tlist(i)/60,1)
        p4(j) = plot(z2,Tintrp2(:,i),'o','MarkerSize',5,'MarkerIndices',1:2:length(z2),...
            'MarkerFaceColor',[0.5 0.5 0.5],'DisplayName',['t = ',strcat(' ',num2str(time2(j))), ' min']);
        legend('Location','NorthWest');
k = j;
switch k
    case 1; p4(k).Marker = 's'; 
    case 2; p4(k).Marker = '<'; 
    case 3; p4(k).Marker = 'X'; 
    case 4; p4(k).Marker = '*'; 
    otherwise;  p4(k).Marker = '+';
end
    hold on
end
title({'Temperature at Selected Interior Nodes'; ['along z-Coordinate (x = ',num2str(x2(1)),' m)',...
    ];''},'FontWeight','normal')
grid on; grid minor; xlabel('z (m)'); ylabel('T (degC)'); 
ax = gca; ax.FontSize = 14; xlim([0 length2]);

figure(23);
hold all;
Tintrp3 = interpolateTemperature(results,x3,y3,z3,1:length(t2));
Tintrp4 = interpolateTemperature(results,x4,y4,z4,1:length(t2));
dt1 = 0.5*time_step/60;
for j=[1,34]
        i = j;
        time2(j) = round(tlist(i)/60,1);
        p5(j) = plot(x3,Tintrp3(:,i),'o','MarkerSize',6,'MarkerIndices',1:3:length(x3),...
            'MarkerFaceColor',[0.5 0.5 0.5],'DisplayName',...
            ['z = ',strcat(' ',num2str(z3(:,1))), ' m, t = ',strcat(' ',num2str(time2(j))),' min']);
        p52(j) = plot(x4,Tintrp4(:,i),'>','MarkerSize',6,'MarkerIndices',1:5:length(x4),...
            'MarkerFaceColor',[0.5 0.5 0.5],'DisplayName',...
            ['z = ',strcat(' ',num2str(z4(:,1))), ' m, t = ',strcat(' ',num2str(time2(j))),' min']);
        legend('Location','NorthEast');
k = j;
switch k
    case 1; p5(k).Marker = 's'; p52(k).Marker = '>'; 
    case 12; p5(k).Marker = 'X'; p52(k).Marker = '<'; 
    otherwise;  p5(k).Marker = '+'; p52(k).Marker = 'd'; 
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; ylim([T0_1 max(T1(:,1))]);

figure(24);
t22 = [length(t2)]/5;
[Tx3,Ty3,Tz3] = evaluateTemperatureGradient(results,x3,y3,z3,1:length(t2));
quiver(x3(:),y3(:),Tx3(:,t22),Ty3(:,t22),'MarkerSize',2,'LineWidth',1,'Color','black',...
    'AutoScale','on','AutoScaleFactor',2.5,'ShowArrowHead','on');
xlim([0 width]); ylim([0 length2]);
title({'Temperature Gradient Vectors, dT/d(x,y) at Selected'; ['Nodes (degC/m),(z = ',...
    num2str(z3(1)),' m, ','t = ',num2str(t22),' min)'];''},'FontWeight','normal')
grid on; grid minor; xlabel 'x (m)'; ylabel 'y (m)';
ax = gca; ax.FontSize = 14; ylim([-0.001 length2*1.05]); xlim([-0.001 width*1.05]); 

figure(25)
hold all;
[qx3,qy3,qz3] = evaluateHeatFlux(results,x3,y3,z3,1:length(t2));
quiver(x3(:),y3(:),qx3(:,t22),qy3(:,t22),'LineWidth',1,'Color','black',...
    'AutoScale','on','AutoScaleFactor',2.5,'ShowArrowHead','on');
xlim([0 width]); ylim([0 length2]);
title({'Heat Flux Vectors at Selected Nodes (W/m^2)'; ['(z = ',num2str(z3(1)),' m, ','t = ',...
    num2str(t22),' min)'];''},'FontWeight','normal')
grid on; grid minor; xlabel 'x (m)'; ylabel 'y (m)';
ax = gca; ax.FontSize = 14; ylim([-0.001 length2*1.15]); xlim([-0.004 width*1.05]); 

figure(26)
hold all;
n = length(t2)/15;
Qn_interior11 = evaluateHeatRate(results,'Face',1);
Qn_interior21 = evaluateHeatRate(results,'Face',15);
Qn_interior12 = evaluateHeatRate(results,'Face',7);
Qn_interior22 = evaluateHeatRate(results,'Face',21);
Qn_exterior1 = evaluateHeatRate(results,'Face',14);
Qn_exterior2 = evaluateHeatRate(results,'Face',13);
plot(t2(1:n),Qn_interior11(1:n),'-d','MarkerSize',4,'MarkerIndices',1:4:length(t2(1:n)));
plot(t2(1:n),Qn_interior21(1:n),'-s','MarkerSize',5,'MarkerIndices',1:5:length(t2(1:n)));
plot(t2(1:n),Qn_interior12(1:n),'-<','MarkerSize',4,'MarkerIndices',1:3:length(t2(1:n)));
plot(t2(1:n),Qn_interior22(1:n),'-+','MarkerSize',5,'MarkerIndices',1:4:length(t2(1:n)));
plot(t2(1:n),Qn_exterior1(1:n),'-X','MarkerSize',4,'MarkerIndices',1:6:length(t2(1:n)));
plot(t2(1:n),Qn_exterior2(1:n),'->','MarkerSize',5,'MarkerIndices',1:5:length(t2(1:n)));
title({'Heat Flow Rate across the Faces';''},'FontWeight','normal'); 'LineWidth'==1; 
legend('Face 1 (external)','Face 15 (internal)','Face 7 (external-left)','Face 21 (internal-left)',...
    'Face 14 (top)','Face 13 (bottom)','Location','NorthEast');
grid on; grid minor; grid minor; xlabel('t (min)'); ylabel('hfr (W)'); 
ax = gca; ax.FontSize = 14;

figure(27)
hold all;
for i = [t22]
    pyx(i) = plot(y3,Tx3(:,i),'o','MarkerSize',5,'MarkerFaceColor',[0.5 0.5 0.5],'DisplayName',['Tx']);    
    pyy(i) = plot(y3,Ty3(:,i),'X','MarkerSize',5,'MarkerFaceColor',[0.5 0.5 0.5],'DisplayName',['Ty']);
    legend('Location','NorthWest');
    k = i;
switch k
    case 1; pyx(k).Marker = 's'; pyy(k).Marker = '>'; 
    case t22; pyx(k).Marker = 'X'; pyy(k).Marker = 'd'; 
    otherwise;  pyx(k).Marker = 'o'; pyy(k).Marker = 'p'; 
end
    hold on
end;
title({'Spatial Temperature Gradients';['(t = ',num2str(i),' min)'];''},'FontWeight','normal')
grid on; grid minor; xlabel('y (m)'); ylabel('dT/d(x,y) (degC/m)');
ax = gca; ax.FontSize = 14;

figure(28);
hold all;
for j=[1 length(t2)/10]
        i = 1*j;
        time2(j) = round(tlist(i)/60,1)
        p42(j) = plot(y2,Tintrp2(:,i),'o','MarkerSize',5,'MarkerIndices',1:2:length(x2),...
            'MarkerFaceColor',[0.5 0.5 0.5],'DisplayName',['t = ',strcat(' ',num2str(time2(j))), ' min']);
        legend('Location','NorthWest');
k = j;
switch k
    case 1; p42(k).Marker = 's'; 
    case 2; p42(k).Marker = '<'; 
    case 3; p42(k).Marker = 'X'; 
    case 4; p42(k).Marker = '*'; 
    otherwise;  p42(k).Marker = '+';
end
    hold on
end
title({'Temperature at Selected Interior Nodes'; ['along z-Coordinate (x = ',num2str(x2(1)),' m)',...
    ];''},'FontWeight','normal')
grid on; grid minor; xlabel('x (m)'); ylabel('T (degC)'); 
ax = gca; ax.FontSize = 14; xlim([0 length2]);

figure(29);
pdeplot3D(thermalModelT,'ColorMapData',T1(:,t22)); hold all;
quiver3(x00(:),y00(:),z00(:),Tx1(:,t22),Ty1(:,t22),Tz1(:,t22),'MarkerSize',2,'LineWidth',1,...
    'Color','black','AutoScale','on','AutoScaleFactor',1.5,'ShowArrowHead','on','LineStyle','-');
xlim([0 width]); ylim([0 length2]); zlim([0 height]);
title({'Temperature Gradient Vectors, dT/d(x,y,z)'; [' at Selected Nodes (degC/m) (t = ',...
    num2str(t22),' min)'];''},'FontWeight','normal');
grid on; grid minor; xlabel 'x (m)'; ylabel 'y (m)';
ax = gca; ax.FontSize = 14; ylim([-0.001 length2*1.05]); xlim([-0.001 width*1.05]);

figure(30);
pdeplot3D(thermalModelT,'ColorMapData',T1(:,t22)); hold all;
[qx1,qy1,qz1] = evaluateHeatFlux(results,x00,y00,z00,1:length(t2));
quiver3(x00(:),y00(:),z00(:),qx1(:,t22),qy1(:,t22),qz1(:,t22),'MarkerSize',2,'LineWidth',1,...
    'Color','black','AutoScale','on','AutoScaleFactor',1.5,'ShowArrowHead','on','LineStyle','-');
xlim([0 width]); ylim([0 length2]); zlim([0 height]);
title({'Heat Flux Vectors, at Selected Nodes'; ['(W/m^2) (t = ',...
    num2str(t22),' min)'];''},'FontWeight','normal');
grid on; grid minor; xlabel 'x (m)'; ylabel 'y (m)';
ax = gca; ax.FontSize = 14; ylim([-0.001 length2*1.05]); xlim([-0.001 width*1.05]); 
