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

% Note:
% A steady-state model does not use 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')

thermalModelS = createpde('thermal','steadystate'); % create a model object: thermal
width = 10/1000; height1 = 4/1000;
geom1 = [3;4;0;width;width;0;0;0;height1;height1]; % create geometry: rectangle solid

gd = [geom1];
sf = 'geom1';
ns = char('geom1');
ns = ns';
dl = decsg(gd,sf,ns);
geometryFromEdges(thermalModelS,dl); % associate geometry and thermal model object
figure(1);
pdegplot(thermalModelS,'EdgeLabels','on','FaceLabels','on');
xlim([-2/1000 width+2/1000]); ylim([-2/1000 height1+2/1000]);
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;

%Aluminum
k1 = 210; % thermal conductivity (W/mK)
q1 = 1000; %heat generation (W/m^3)
kFunc2 = @(region,state) k1*region.y;
qFunc2 = @(region,state) q1*region.y;
tp1=thermalProperties(thermalModelS,'ThermalConductivity',k1,'Face',1)
q_internal=internalHeatSource(thermalModelS,q1,'Face',1)

% Boundary Conditions
T_ambient = 25; %degC **********
T_water = 50;   %degC **********
hc = 10; %convection heat transfer coefficient
qs = 500; %side-surface heat flux (W/m^2K)
top_BC_HF = @(region,~) qs*region.y; %heat source distribution (W/m2)
outerCC_V = @(region,~) hc*region.y; %horizontal heat convection coefficient (W/m2K)
thermalBC(thermalModelS,'Edge',3,'HeatFlux',top_BC_HF); %boundary condition: heat source (W/m2)
thermalBC(thermalModelS,'Edge',[2 4],'ConvectionCoefficient',outerCC_V,'AmbientTemperature',T_ambient); %boundary condition: convection heat transfer
thermalBC(thermalModelS,'Edge',1,'Temperature',T_water); %boundary condition: temperature

figure(2);
grid_size = 0.5/1000; %identify grid size
mesh1 = generateMesh(thermalModelS,'Hmax',grid_size); %generate mesh *****
pdeplot(mesh1); %plot PDE thermal model object with node ID's displayed.
xlim([-2/1000 width+2/1000]); ylim([-2/1000 height1+2/1000]);
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'); %plot PDE thermal model object with node ID's displayed.
xlim([0 1.5*grid_size]); ylim([0 1.5*grid_size]); % set axis conditions
title('','FontWeight','normal');
title({'Node IDs';''},'FontWeight','normal');
grid on; grid minor; xlabel('x (m)'); ylabel('y (m)');
ax = gca; ax.FontSize = 14;

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

thermalModelS.SolverOptions.ReportStatistics = 'on';
results = solve(thermalModelS); % 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 *****
 
figure(5);
pdeplot(thermalModelS,'XYData',T1(:,end),'ColorMap','hot','Contour','on','Levels',10); %,'ColorMap','hot' ,'Contour','on'
xlim([0 width]); ylim([0 height1]);
title({'Steady-State Temperature Contours (degC)';''},'FontWeight','normal')
grid on; grid minor; xlabel('x (m)'); ylabel('y (m)');
ax = gca; ax.FontSize = 14;

figure(6);
pdeplot(thermalModelS,'XYData',Tx1(:,end),'ColorMap','hot','Contour','on','Levels',10); %,'ColorMap','hot' ,'Contour','on'
xlim([0 width]); ylim([0 height1]); 
title({'Steady-State Temperature Gradient';' (dT/dx) Contours (degC/m)'},'FontWeight','normal')
grid on; grid minor; xlabel('x (m)'); ylabel('y (m)'); xlabel('x (m)'); ylabel('y (m)');
ax = gca; ax.FontSize = 14;

figure(7);
pdeplot(thermalModelS,'XYData',Ty1(:,end),'ColorMap','hot','Contour','on','Levels',10); %,'ColorMap','hot' ,'Contour','on'
xlim([0 width]); ylim([0 height1]); 
title({'Steady-State Temperature Gradient';' (dT/dy) Contours (degC/m)'},'FontWeight','normal')
 grid on; grid minor; xlabel('x (m)'); ylabel('y (m)');
 ax = gca; ax.FontSize = 14;
 
%**********************************************************
%Investigate the node density

point = thermalModelS.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,:);

nodesLeftEnd_v = find(x00 == 0); % identify far left boundary nodes
nodesCenter_v = find(and(0.5*width - 0.5*grid_size <= x00,x00 <= 0.5*width + 0.5*grid_size)); % identify far left boundary nodes
nodesInterior_h = find(and(height1 - eps <= y00, y00 <= height1 + eps)); % 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,:);

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
xy5=xy2(:,nodesInterior_1); x5=xy5(1,:); y5=xy5(2,:);
xy6=xy2(:,nodesInterior_2); x6=xy6(1,:); y6=xy6(2,:);

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

figure(9);
hold all;
plot(x3,y3,'o','MarkerSize',5,'MarkerFaceColor',[0.5 0.5 0.5]);
xlim([0 width]); ylim([0 height1]);
title({'Horizontal Nodes at Selected Coordinates';''},'FontWeight','normal');
grid on; grid minor; xlabel('x (m)'); ylabel('y (m)');
ax = gca; ax.FontSize = 14;

figure(10); 
hold all;
Tintrp2 = interpolateTemperature(results,x2,y2);
p2 = plot(y2,Tintrp2,'d','MarkerSize',5,'MarkerFaceColor',[0.5 0.5 0.5]);
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(11);
hold all;
Tintrp3 = interpolateTemperature(results,x3,y3);
p5 = plot(x3,Tintrp3,'o','MarkerSize',5,'MarkerFaceColor',[0.5 0.5 0.5]);
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;

figure(12)
hold all; 
[qx2,qy2] = evaluateHeatFlux(results,x2,y2);
quiver(x2(:),y2(:),qx2(:),qy2(:),'MarkerSize',2,'LineWidth',2);
axis equal; xlim([0 width]); ylim([-1/1000 height1+0.5/1000]);
title({'Heat Flux Vectors at Selected Nodes (W/m^2)';''},'FontWeight','normal');
grid on; grid minor; xlabel 'x (m)'; ylabel 'y (m)';
ax = gca; ax.FontSize = 14;

figure(13)
hold all; 
[qx3,qy3] = evaluateHeatFlux(results,x3,y3);
quiver(x3(:),y3(:),qx3(:),qy3(:),'MarkerSize',2,'LineWidth',2);
axis equal; xlim([0 width]); ylim([-1/1000 height1+1/1000]);
title({'Heat Flux Vectors at Selected Nodes (W/m^2)';''},'FontWeight','normal');
grid on; grid minor; xlabel 'x (m)'; ylabel 'y (m)';
ax = gca; ax.FontSize = 14;

Qn_interior = evaluateHeatRate(results,'Edge',1);
Qn_exterior = evaluateHeatRate(results,'Edge',3);
Qn_EdgeLeftB = evaluateHeatRate(results,'Edge',4);
Qn_EdgeRightB = evaluateHeatRate(results,'Edge',2);
formatSpec1 = 'Qn_interior is %4.2f ; '; fprintf(formatSpec1,Qn_interior);
formatSpec2 = 'Qn_exterior is %4.2f'; fprintf(formatSpec2,Qn_exterior);
formatSpec3 = 'Qn_EdgeLeftB is %4.2f ; '; fprintf(formatSpec3,Qn_EdgeLeftB);
formatSpec4 = 'Qn_EdgeRightB is %4.2f'; fprintf(formatSpec4,Qn_EdgeRightB);

figure(14)
hold all;
[Tx3,Ty3] = evaluateTemperatureGradient(results,x3,y3);
pxx = plot(x3,Tx3,'o','MarkerSize',4,'MarkerFaceColor',[0.5 0.5 0.5]);    
pxy = plot(x3,Ty3,'<','MarkerSize',4,'MarkerFaceColor',[0.5 0.5 0.5]);
legend('Tx','Ty','Location','East');
pxx(1).Marker = 's'; pxy(1).Marker = 'd';
xlim([0 width]); %ylim([-1/1000 height1+1/1000]);
title({'Spatial Temperature Gradients (Tx and Ty)';'along x and y-Coordinates';''},...
    'FontWeight','normal');
grid on; grid minor; xlabel('x (m)'); ylabel('dT/d(x,y) (degC/m)');
ax = gca; ax.FontSize = 14;

figure(15)
hold all;
[Tx2,Ty2] = evaluateTemperatureGradient(results,x2,y2);
pxx = plot(x2,Tx2,'o','MarkerSize',4,'MarkerFaceColor',[0.5 0.5 0.5]);    
pxy = plot(x2,Ty2,'<','MarkerSize',4,'MarkerFaceColor',[0.5 0.5 0.5]);
legend('Tx','Ty','Location','East');
pxx(1).Marker = 's'; pxy(1).Marker = 'd';
xlim([0 width]); %ylim([-1/1000 height1+1/1000]);
title({'Spatial Temperature Gradients (Tx and Ty)';'along x and y-Coordinates';''},...
    'FontWeight','normal');
grid on; grid minor; xlabel('x (m)'); ylabel('dT/d(x,y) (degC/m)');
ax = gca; ax.FontSize = 14;

w = linspace(0,width,10);
h = linspace(0,height1,10);
[Xw,Yh] = meshgrid(w,h);
T_xy = interpolateTemperature(results,Xw,Yh);
[Tx,Ty] = evaluateTemperatureGradient(results,Xw,Yh);
[qx,qy] = evaluateHeatFlux(results,Xw,Yh);
Qn = evaluateHeatRate(results,'Edge',[1 2]);





thermalModelS
