% Author: Layla S. Mayboudi
% 2021-05-27

% Note:
% This code should be run after running the code Pump_Heating_MATLAB_App.m
% The geometry data, set formula, and labels (gd sf ns) are to be imported herein.
% (gd sf ns) are to be imported before running the program by selecting
% the "Export Geometry Description" under the "Draw" item from
% the PDE Modeler Application toolbar.

% This code creates a PDE thermal transient model for a rectangular object, 2D pipe,
% assigns thermophysical properties, boundary and initial conditions, and internal heat sources, 
% meshes and solves the model, and generates the contour plots.
% It also chooses nodes at certain coordinates and plots their temperature
% distribution as a function of time. 
% Animated temperature and heat flux contours are also plotted.
% The solution time is extended.


clc; close all; %clear all;

set(groot,'DefaultAxesColorOrder',[0 0 0],...
    'DefaultSurfaceEdgeColor',[0 0 0],...
    'DefaultLineLineWidth',1,...
    'DefaultAxesLineStyleOrder','-|--|:|-.|+|o|*|.|x|s|d|^|v|>|<|p|h');

thermalModelTA = createpde('thermal','transient'); % create a model object: thermal
geometryA = decsg(gd,sf,ns);
geo = geometryFromEdges(thermalModelTA,geometryA); % associate geometry and thermal model object
%subplot(2,1,1)
figure(1);
plott = pdegplot(thermalModelTA,'EdgeLabels','on');
grid on; grid minor; xlabel('x (m)'); ylabel('y (m)'); 
title({'Geometry with Edge Labels Displayed';''},'FontWeight','normal');
xlim([-0.7 0.7]); ylim([-0.6 0.8]); ax = gca; ax.FontSize = 14;

%subplot(2,1,2)
figure(2);
pdegplot(thermalModelTA,'FaceLabels','on');
grid on; grid minor; xlabel('x (m)'); ylabel('y (m)'); 
title({'Geometry with Face Labels Displayed';''},'FontWeight','normal');
xlim([-0.7 0.7]); ylim([-0.6 0.8]); ax = gca; ax.FontSize = 14;
                       
%Water
k15 = 0.6305; % thermal conductivity (W/mK)
rho15 = 988; %density (kg/m3)
cp15 = 4181; % heat capacity (J/kgK)
tp15 = thermalProperties(thermalModelTA,'ThermalConductivity',k15,...
    'MassDensity',rho15,'SpecificHeat',cp15,'Face',[2,8,9,11,12,13]);
%Aluminum
k25 = 210; % thermal conductivity (W/mK)
rho25 = 2698.9; %density (kg/m3)
cp25 = 900; % heat capacity (J/kgK)
tp25 = thermalProperties(thermalModelTA,'ThermalConductivity',k25,...
    'MassDensity',rho25,'SpecificHeat',cp25,'Face',[1,3:7,10]);

%Internal heat sources
ih51=internalHeatSource(thermalModelTA,2000,'Face',3);
ih52=internalHeatSource(thermalModelTA,4000,'Face',4);
ih53=internalHeatSource(thermalModelTA,10000,'Face',10);
ih54=internalHeatSource(thermalModelTA,-6000,'Face',[12,13]);

%Boundary conditions
BC51 = thermalBC(thermalModelTA,'Edge',6,'Temperature',35);
BC52 = thermalBC(thermalModelTA,'Edge',2,'Temperature',45);
BC53 = thermalBC(thermalModelTA,'Edge',[21,18],'Temperature',15);

BC54 = thermalBC(thermalModelTA,'Edge',[1,4,22:24,26],...
                       'ConvectionCoefficient',10,...
                       'AmbientTemperature',30);

BC55 = thermalBC(thermalModelTA,'Edge',1,'HeatFlux',500);

%Initial condition
IC56 = thermalIC(thermalModelTA,25);

%internalHeatSource(thermalModelTA,@myHeatSource);

grid_size5 = 0.05;
mesh5 = generateMesh(thermalModelTA,'Hmax',grid_size5,'Hgrad',1.2,...
    'GeometricOrder','quadratic');

thermalModelTA.SolverOptions.ReportStatistics = 'on'
modelD.SolverOptions.RelativeTolerance = 1E-5;
modelD.SolverOptions.AbsoluteTolerance = 1E-9;

tfinal5 = 5000;
tlist5 = 0:100:tfinal5;
results5 = solve(thermalModelTA,tlist5);
point = thermalModelTA.Mesh.Nodes; % recall node locations from thermal model object
nodesIndex = find(point); % identify far left boundary nodes
xy51=point; x51=xy51(1,:); y51=xy51(2,:);
T51 = results5.Temperature;
Tx51 = results5.XGradients; % identify data gradient along x-coordinate *****
Ty51 = results5.YGradients; % identify data gradient along y-coordinate *****
T_st51 = results5.SolutionTimes; % identify data gradient along y-coordinate *****
[qx51,qy51] = evaluateHeatFlux(results5,x51,y51,1:length(tlist5));

figure(3)
plott = plot(x51,y51,'.')
axis equal;
title({'Node Density';''},'FontWeight','normal');
grid on; grid minor; xlabel('x (m)'); ylabel('y (m)');
xlim([-0.7 0.7]); ylim([-0.6 0.8]); ax = gca; ax.FontSize = 14;

n = length(tlist5)/3;
figure(4)
pdeplot(thermalModelTA,'XYData',T51(:,end),'ColorMap','hot','Contour','on','Levels',24);
hold all; 
pdemesh(thermalModelTA); 
Ef10 = findElements(mesh5,'region','Face',10); 
pdemesh(mesh5.Nodes,mesh5.Elements(:,Ef10),'EdgeColor','blue','ColorMap',hot,'FaceAlpha',0);
Ef12 = findElements(mesh5,'region','Face',12); 
pdemesh(mesh5.Nodes,mesh5.Elements(:,Ef12),'EdgeColor','black','ColorMap',hot,'FaceAlpha',0);
quiver(x51(:),y51(:),qx51(:,end),qy51(:,end),'-','LineWidth',0.5,'Color','black',...
    'AutoScale','on','ShowArrowHead','on');%,'AutoScaleFactor',3);
axis equal; cmap = colormap(gca,gray);
title({['Temperature Contours'];['and Heat Flux Arrows (t = ',...
    num2str(round(tlist5(end),2)),' s)'];''},'FontWeight','normal')
grid on; grid minor; xlabel('x (m)'); ylabel('y (m)');
xlim([-0.7 0.7]); ylim([-0.6 0.8]); ax = gca; ax.FontSize = 20;

figure(5)
%subplot(2,1,1)
pdeplot(thermalModelTA,'XYData',Tx51(:,end),'ColorMap','hot','Contour','on','Levels',24);
hold all; 
%pdemesh(thermalModelTA); 
Ef10 = findElements(mesh5,'region','Face',10); 
pdemesh(mesh5.Nodes,mesh5.Elements(:,Ef10),'EdgeColor','blue','ColorMap',hot,'FaceAlpha',0);
Ef12 = findElements(mesh5,'region','Face',12); 
pdemesh(mesh5.Nodes,mesh5.Elements(:,Ef12),'EdgeColor','red','ColorMap',hot,'FaceAlpha',0);
quiver(x51(:),y51(:),Tx51(:,end),Ty51(:,end),'-','LineWidth',0.5,'Color','black',...
    'AutoScale','on','ShowArrowHead','on','AutoScaleFactor',3);
axis equal; cmap = colormap(gca,turbo);
title({['Temperature Gradient Contours with Respect to'];['x (dT/dx) and their Arrows (t = ',...
    num2str(round(tlist5(end),2)),' s)'];''},'FontWeight','normal')
grid on; grid minor; xlabel('x (m)'); ylabel('y (m)');
xlim([-0.7 0.7]); ylim([-0.6 0.8]); ax = gca; ax.FontSize = 20;

%subplot(2,1,2)
figure(6)
pdeplot(thermalModelTA,'XYData',Ty51(:,end),'ColorMap','hot','Contour','on','Levels',24);
hold all; 
%pdemesh(thermalModelTA); 
Ef10 = findElements(mesh5,'region','Face',10); 
pdemesh(mesh5.Nodes,mesh5.Elements(:,Ef10),'EdgeColor','blue','ColorMap',hot,'FaceAlpha',0);
Ef12 = findElements(mesh5,'region','Face',12); 
pdemesh(mesh5.Nodes,mesh5.Elements(:,Ef12),'EdgeColor','green','ColorMap',hot,'FaceAlpha',0);
quiver(x51(:),y51(:),Tx51(:,end),Ty51(:,end),'-','LineWidth',0.5,'Color','black',...
    'AutoScale','on','ShowArrowHead','on','AutoScaleFactor',3);
axis equal; cmap = colormap(gca,turbo);
title({['Temperature Gradient Contours with Respect to'];[' y (dT/dy) and their arrows (t = ',...
    num2str(round(tlist5(end),2)),' s)'];''},'FontWeight','normal')
grid on; grid minor; xlabel('x (m)'); ylabel('y (m)');
xlim([-0.7 0.7]); ylim([-0.6 0.8]); ax = gca; ax.FontSize = 20;

figure(7)
for j = 1:5:length(tlist5)    
    pdeplot(thermalModelTA,'XYData',T51(:,j),'ColorMap','hot','Contour','on','Levels',24);
    axis equal;
    xlim([-0.7 0.7]); ylim([-0.6 0.8]); ax = gca; ax.FontSize = 20;
    hold all;
    quiver(x51(:),y51(:),qx51(:,j),qy51(:,j),'-','LineWidth',0.5,'Color','black',...
    'AutoScale','on','ShowArrowHead','on');%,'AutoScaleFactor',3
    axis equal; cmap = colormap(gca,turbo);
    xlim([-0.7 0.7]); ylim([-0.6 0.8]); ax = gca; ax.FontSize = 20;
    grid on; grid minor; xlabel('x (m)'); ylabel('y (m)');    
    title({['Temperature Contours'];['and Heat Flux Arrows (t = ',...
    num2str(round(tlist5(j),2)),' s)'];''},'FontWeight','normal')
M5(j) = getframe;
end

figure(8)
for k = 1:5:length(tlist5)    
    title({['Temperature Contours'];['and Heat Flux Arrows (t = ',...
    num2str(round(tlist5(k),2)),' s)'];''},'FontWeight','normal')
axis equal; cmap = colormap(gca,turbo);
grid on; grid minor; xlabel('x (m)'); ylabel('y (m)');    
xlim([-0.7 0.7]); ylim([-0.6 0.8]); ax = gca; ax.FontSize = 20;
movie(gca,M5,[7 k],14);
end

figure(9);
hold all;
n = length(tlist5);
data = @(p,x,y)min((x51-x).^2+(y51-y).^2);
[~,nid1] = data( mesh5.Nodes,-0.5,0);
[~,nid2] = data( mesh5.Nodes,0.2,-0.4);
[~,nid3] = data( mesh5.Nodes,-.1,0.2);
A = nid1
B = nid2
C = nid3
plot(tlist5(1:n),T51(A,[1:n]),'-<','MarkerSize',5,'MarkerIndices',...
    1:5:length(tlist5(1:n)),'MarkerFaceColor',[1 1 1]);
plot(tlist5(1:n),T51(B,[1:n]),'-d','MarkerSize',5,'MarkerIndices',...
    1:5:length(tlist5(1:n)),'MarkerFaceColor',[0.5 0.5 0.5]);
plot(tlist5(1:n),T51(C,[1:n]),'->','MarkerSize',5,'MarkerIndices',...
    1:5:length(tlist5(1:n)),'MarkerFaceColor',[0 0 0]);
xlim([0 tlist5(1,n)]); 
legend('(-0.5,0)','(0.2,-0.4)','(-0.1,0.2)','Location','west');
title({'Transient Temperature at the Selected Noded';''},'FontWeight','normal');
grid on; grid minor; xlabel 't (min)'; ylabel 'T (degC)'; 
ax = gca; ax.FontSize = 14;

figure(10);
pdegplot(thermalModelTA);
grid on; grid minor; xlabel('x (m)'); ylabel('y (m)'); 
title({'Geometry with Edge Labels Displayed';''},'FontWeight','normal');
xlim([-0.7 0.7]); ylim([-0.6 0.8]); ax = gca; ax.FontSize = 14;

% Returning a NaN when time equals NaN, indforming the solver that the heat
% source is a function of time.
% Parameters are to be modified to the desired values.
function Q = myHeatSource(location,state)
  Q = zeros(1,numel(location.x));
if(isnan(state.time))
  Q(1,:) = NaN;
  return
end
if state.time < 300
  Q(1,:) = 00*state.time;
end
end
