% Author: Layla S. Mayboudi
% 2021-05-27
% This code creates a transient thermal PDE model for a moving heat source.
% Heat source power and velocity, ambient temperature, and beam size and
% contour path are global variable.
% It plots the geometry and heat flux flow diagrams.
% Selected nodes and their temperatures are also displayed.

% Note:
% The model run time estimated time can exceed 20 minutes, depending on the
% computational capability.

% A transient model uses ThemalConductivity, MassDensity and SpecificHeat. 
% Geometry Pipe_hollow_shortened.stl is imported to this model and then scaled to the desired
% portions.


clear; clc; close all;

%********** % Global data **********
global T_low T_high t_low t_high Velocity q_heat hc T_ambient x_c spot_w y_max grid_size xxx yyy zzz indx
T_low = 35; T_high = 80;
t_low = 5; t_high = 8;
q_heat = 5; T_ambient = 25; Velocity = 10/1000; spot_w = 0.004;

%********** % Model properties - transient thermal analysis **********
disp('Creating PDE...')
thermalModelT = createpde('thermal','transient'); % create a model object: thermal
disp('Importing geometry...')
geom=importGeometry(thermalModelT,'Pipe_hollow_shortened.stl');
scalef = 1000;
scale(geom,[1 1 2]/scalef);

%********** % Geometry input and plot *****
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)'); zlabel('z (m)');
ax = gca; ax.FontSize = 14; axis equal;
%hgexport(figure(1),'-clipboard')

%********** % Pipe material properties (PEX) **********
k1 = 0.41; % thermal conductivity (W/mK)
rho1 = 935; % density (kg/m^3)
cp1 = 2100; % heat capacity (J/kgK)
tp1 = thermalProperties(thermalModelT,'ThermalConductivity',k1,'MassDensity',rho1,'SpecificHeat',cp1,'Cell',1) ;

%********** % Initial conditions **********
tfinal = 10; % final time (s)
time_step = 0.02; % time step (s)
tlist = 0:time_step:tfinal;
t = tlist;
%linspace(0,tfinal,50); % set time array
t_length = length(tlist);
T0_1 = 25; % degC Initial tempeature
IC1 = thermalIC(thermalModelT,T0_1,'Cell',1);

%********** % Mesh properties **********
grid_size = 0.002; %identify grid size (m)
disp('Generating mesh')
mesh = generateMesh(thermalModelT,'Hmax',grid_size,'GeometricOrder','linear');
figure(2);
pdeplot3D(mesh);
title({'Tetrahedral Mesh'},'FontWeight','normal');
ax = gca; ax.FontSize = 14; axis equal; view(175,45)

%********** % Recall node locations from thermal model object
point = thermalModelT.Mesh.Nodes;
eps1=eps;
nodesIndex = find(point); % identify all nodes
xy00=point; x00=xy00(1,:); y00=xy00(2,:); z00=xy00(3,:);
width2 = max(x00)-min(x00); height2 = max(z00)-min(z00); depth2 = max(y00)-min(y00);
nodesLeftEnd_v = find(x00 == min(x00)); % identify far left boundary nodes, x = minx(x)
nodesCenter_v = find(and(0.5*width2 - 0.15*grid_size <= x00,x00 <= 0.5*width2 + 0.15*grid_size)); % identify central nodes, x = width2/2
nodesTop_h = find(and(max(z00) - eps <= z00, z00 <= max(z00) + eps)); % identify top nodes, z = max(z) = height2
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,:);

%********** % Identify center of the cylinder and interior radius points
pp_x = find(and(x3>0.0040,x3<0.0558596));
x33 = x3(pp_x); y33=y3(pp_x); z33=z3(pp_x);
pp_y = find(and(y33>0.0040,y33<0.0558596));
x333 = x33(pp_y); y333=y33(pp_y); z333=z33(pp_y);
min_axis = max(y333)-min(y333);

y0000 = (max(y3)+min(y3))/2;
maj_axis = max(x3)-min(x3);
x0000 = (max(x3)+min(x3))/2;
radius3 = sqrt((x3-x0000).^2+(y3-y0000).^2);
r_index = radius3 <= min_axis*0.5 + grid_size *0.5;
x33 = x3(r_index); y33 = y3(r_index); z33 = z3(r_index);
y_max = max(y00);
x_c = x0000;

figure(3);
hold all;
plot(x33,y33,'o','MarkerSize',5,'MarkerIndices',1:1:length(x33),'Color','black');
title({'Projected Nodes on z-Coordinate';[num2str(y33(:,1)),' m < y < ',...
    num2str(y33(:,end)),' m'];''},'FontWeight','normal');
grid on; grid minor; xlabel ('x (m)'); ylabel ('y (m)'); 
ax = gca; ax.FontSize = 14; axis equal; view(175,45)

pp_y = find(and(y33>0.0040,y33<0.0558596));
x333 = x33(pp_y); y333=y33(pp_y); z333=z33(pp_y);
min_axis = max(y333)-min(y333);

%********** % Identify the moving heat source points

figure(4);
hold all;
plot(x00,y00,'o','MarkerSize',5,'MarkerIndices',1:25:length(x00),'Color','black');
title({'Projected Nodes on z-Coordinate';''},'FontWeight','normal');
grid on; grid minor; xlabel ('x (m)'); ylabel ('y (m)'); 
xlim([0 max(x00)]); ylim([0 max(y00)])
ax = gca; ax.FontSize = 14; axis equal; view(175,45)

heat_xx = find(and(x0000-2*grid_size<x00,x00<x0000+2*grid_size));
heat_x1 = x00(heat_xx); heat_y1 = y00(heat_xx); heat_z1 = z00(heat_xx);
heat_yy = find(max(y00)-0.5*grid_size < heat_y1);
heat_x2 = heat_x1(heat_yy); heat_y2 = heat_y1(heat_yy); heat_z2 = heat_z1(heat_yy);

figure(5);
hold all;
plot(heat_x2,heat_z2,'o','MarkerSize',5,'MarkerIndices',1:1:length(heat_x2),'Color','black');
title({'Heat Source Nodes';''},'FontWeight','normal');
grid on; grid minor; xlabel ('x (m)'); ylabel ('z (m)'); 
xlim([0 max(x00)]); ylim([0 max(z00)])
ax = gca; ax.FontSize = 14; axis equal; view(175,45)

figure(6);
hold all;
plot(heat_y2,heat_z2,'o','MarkerSize',5,'MarkerIndices',1:1:length(heat_y2),'Color','black');
title({'Heat Source Nodes';''},'FontWeight','normal');
grid on; grid minor; xlabel ('y (m)'); ylabel ('z (m)'); 
xlim([0 max(y00)]); ylim([0 max(z00)])
ax = gca; ax.FontSize = 14; axis equal; view(175,45)

%********** % Boundary Conditions **********
%T_ambient = 25; %degC
hc = 10; %heat convectyion coefficient (W/m^2K)
outerCC_V = hc;
thermalBC(thermalModelT,'Face',7:thermalModelT.Geometry.NumFaces,'ConvectionCoefficient',...
    outerCC_V,'AmbientTemperature',T_ambient,'Vectorized','on'); %boundary condition: convection
%q1 = 0; %W/m^3
%ih1=internalHeatSource(thermalModelT,q1,'Cell',1); %cell condition: heat source

thermalModelT.StefanBoltzmannConstant = 5.670373E-8; 
emis = 0.91;
thermalBC(thermalModelT,'Face',[1,2,3,4,5,6],'ConvectionCoefficient',outerCC_V,'Emissivity',emis,...
    'AmbientTemperature',T_ambient,'Vectorized','on')

thermalBC(thermalModelT,'Face',1,'HeatFlux',@q_func,'Vectorized','on'); %boundary condition: heat flux

%thermalBC(thermalModelT,'Face',[13],'Temperature',@Twater_func); %boundary condition: temperature
bca = findThermalBC(thermalModelT.BoundaryConditions,'Face',[1:26]);

%********** % Solve properties **********
thermalModelT.SolverOptions.ReportStatistics = 'on';
thermalModelT.SolverOptions.AbsoluteTolerance = 1.0000e-05;

tic
timeStart1 = datetime('now');
hSt1 = hour(timeStart1); mSt1 = minute(timeStart1); sSt1 = round(second(timeStart1));
fprintf('**** Solution started at %02d:%02d:%02d ****\n',hSt1,mSt1,sSt1)

results = solve(thermalModelT,tlist);

timeStart2 = datetime('now');
hSt2 = hour(timeStart2); mSt2 = minute(timeStart2); sSt2 = round(second(timeStart2));
fprintf('***** Solution ended at  %02d:%02d:%02d ****\n',hSt2,mSt2,sSt2)

sol_time=toc;
disp(['Solution time = ' num2str(round(sol_time/60,1)) ' min'])

%********** % Model results (tempeature, temperature gradient, and solution time) **********
T1 = results.Temperature; % identify tempeature data for the selected nodes *****
tt = round(tfinal/60,1); % identify final time step in min
Tx1 = results.XGradients; % identify tempeature gradient along x-coordinate *****
Ty1 = results.YGradients; % identify tempeature gradient along y-coordinate *****
Tz1 = results.YGradients; % identify tempeature gradient along y-coordinate *****
T_st1 = results.SolutionTimes; % identify solutine times *****

%********** % Interpolate tempeature data at the identified nodes
Tintrp1 = interpolateTemperature(results,x1,y1,z1,1:t_length);
Tintrp2 = interpolateTemperature(results,x2,y2,z2,1:t_length);
Tintrp3 = interpolateTemperature(results,x3,y3,z3,1:t_length);
Tintrp33 = interpolateTemperature(results,x33,y33,z33,1:t_length);
Tintrp4 = interpolateTemperature(results,x4,y4,z4,1:t_length);

%********** % Results contour plot, tempeature and tempeature gradient *****
figure(7);
pdeplot3D(thermalModelT,'ColorMapData',T1(:,round(t_length*0.5,0)),'FaceAlph',0.5);
title({'Temperature Contours (degC)'; ['(t = ',num2str(tlist(round(t_length*0.5,0))),' s)'];''},'FontWeight','normal')
grid on; grid minor; xlabel('x (m)'); ylabel('y (m)');
ax = gca; ax.FontSize = 14; view(175,45)

figure(8);
pdeplot3D(thermalModelT,'ColorMapData',T1(:,t_length),'FaceAlph',0.5);
title({'Temperature Contours (degC)'; ['(t = ',num2str(tlist(t_length)),' s)'];''},'FontWeight','normal')
grid on; grid minor; xlabel('x (m)'); ylabel('y (m)');
ax = gca; ax.FontSize = 14; view(175,45)

figure(9);
pdeplot3D(thermalModelT,'ColorMapData',Tx1(:,t_length),'FaceAlph',0.5);
title({'Temperature Gradient Contours (dT/dx in degC/m)'; ['(t = ',num2str(tlist(t_length)),' s)'];''},'FontWeight','normal')
grid on; grid minor; xlabel('x (m)'); ylabel('y (m)');
ax = gca; ax.FontSize = 14; view(175,45)

figure(10);
pdeplot3D(thermalModelT,'ColorMapData',Ty1(:,t_length),'FaceAlph',0.5);
title({'Temperature Gradient Contours (dT/dy in degC/m)'; ['(t = ',num2str(tlist(t_length)),' s)'];''},'FontWeight','normal')
grid on; grid minor; xlabel('x (m)'); ylabel('y (m)');
ax = gca; ax.FontSize = 14; view(175,45)

%********** % Plot data *****
test2 = round(3*t_length/4,0);
test3 = test2:t_length;
figure(11);
hold all;
plot(tlist(test3),max(Tintrp33(:,test3)),'-d','MarkerSize',6,'MarkerIndices',1:1:length(test3),...
             'MarkerFaceColor',[0.5 0.5 0.5]);
plot(tlist(test3),min(Tintrp33(:,test3)),'-d','MarkerSize',6,'MarkerIndices',1:1:length(test3),...
             'MarkerFaceColor',[0.5 0.5 0.5]);
title({'Interior Surface Tempeature - BC';''},'FontWeight','normal')
grid on; grid minor; xlabel('t (s)'); ylabel('T (degC)'); 
ax = gca; ax.FontSize = 14; legend('Max temp','Min temp','location','best');

figure(12);
plot(x00,y00,'o','MarkerSize',2,'MarkerIndices',1:25:length(x00),'color','black');
title({'Projected Nodes on z-Coordinate';''},'FontWeight','normal');
grid on; grid minor; xlabel ('x (m)'); ylabel ('y (m)'); 
ax = gca; ax.FontSize = 14;

figure(13);
[qx1,qy1,qz1] = evaluateHeatFlux(results,x00,y00,z00,1:t_length);
pdeplot3D(thermalModelT,'ColorMapData',T1(:,t_length)); hold all;
quiver3(x00(:),y00(:),z00(:),qx1(:,t_length),qy1(:,t_length),qz1(:,t_length),'MarkerSize',2,...
    'LineWidth',1,'Color','black','AutoScale','on','AutoScaleFactor',10,'ShowArrowHead','on','LineStyle','-');    
title({'Temperature Contours and Heat Flow Vectors'; [' at Selected Nodes (t = ',num2str(tlist(t_length)),' s)'];''},'FontWeight','normal')
grid on; grid minor; xlabel 'x (m)'; ylabel 'y (m)';
ax = gca; ax.FontSize = 14; view(175,45)

%********** % Heat source nodes *****
figure(14);
hold all;
plot(x00(indx),z00(indx),'o','MarkerSize',5,'MarkerIndices',1:1:length(x00(indx)),'Color','black');
title({'Heat Source Nodes';''},'FontWeight','normal');
grid on; grid minor; xlabel ('x (m)'); ylabel ('z (m)'); 
xlim([0 max(x00)]); ylim([0 max(z00)])
ax = gca; ax.FontSize = 14; axis equal;

%********** % Heat source nodes *****
figure(15);
hold all;
plot(xxx,zzz,'o','MarkerSize',5,'MarkerIndices',1:1:length(xxx),'Color','black');
title({'Heat Source Nodes';''},'FontWeight','normal');
grid on; grid minor; xlabel ('x (m)'); ylabel ('z (m)'); 
xlim([0 max(x00)]); ylim([0 max(z00)])
ax = gca; ax.FontSize = 14; axis equal;

function q_out = q_func(region,state)

global Velocity q_heat hc T_ambient x_c spot_w y_max grid_size xxx yyy zzz indx

t = state.time;
%z_c = 0.05;
z_c = Velocity*t;

x = region.x;
y = region.y;
z = region.z;

q_out = hc*(T_ambient - state.u); % Convection

%if isnan(time2)
%q_out = nan(1,numel(x));
%end

if isnan(t)
    q_out = nan(1,numel(x));
end

indx = (abs(z-z_c)<spot_w*0.5) & (abs(x-x_c)<spot_w*0.5) & (y>(y_max-0.5*grid_size));
%test = find(indx>=1);
xxx = x(indx); yyy = y(indx); zzz = z(indx);

q_out(1,indx) = q_heat/(spot_w^2);
%end
end