% Author: Layla S. Mayboudi
% 2021-05-27
% This code creates a stationary thermal PDE model and plots the geometry 
% and heat flux flow diagrams.

% Note:% A warning is issued at the begging as follows:
% Warning: Steady-state model does not use MassDensity and SpecificHeat. 
% > In pde.ThermalModel/thermalProperties (line 107)
% In thermal_model_pipe (line 37) 
% This is interntionally included for the educational purposes.
% The program functions as it is intended to.


clear all; clc; clf; close all;
thermalModelS = createpde('thermal');

%1. For the circle solid, row one contains 1, the second and
        %third row contain the center x- and y-coordinates respectively
        %Row four contains the radius of the circle
%2. For a polygon solid, row one contains 2, and the second row
        %contains the number, N, of line segments in the boundary.
        %The following N rows contain the x-coordinates of the starting points
        %of the edges, and the following N rows contain the y-coordinates of the
        %starting points of the edges.
%3. For a rectangle solid, row one contains 3. The format is otherwise
        %identical to the polygon format.
%4. For an ellipse solid, row one contains 4, the second and
        %third row contain the center x- and y-coordinates respectively. Row
        %four and five contain the major and minor axes of the ellipse.
        %The rotational angle of the ellipse is stored in row six.

g1 = decsg([3 4 0 0 60 60 26 30 30 26]') % rectangle solid
geometryFromEdges(thermalModelS,g1)
k = 40 % thermal conductivity (W/mK)
rho = 7800 %density (kg/m3)
cp = 500 % heat capacity (J/kgK)
q = 100 %heat generation (W/m^3)
kFunc = @(region,state) k*region.y
cFunc = @(region,state) cp*region.y
qFunc = @(region,state) q*region.y
rhoFunc = @(region,state) rho*region.y

thermalProperties(thermalModelS,'ThermalConductivity',kFunc,...
                               'MassDensity',rhoFunc,...
                               'SpecificHeat',cFunc)

internalHeatSource(thermalModelS,qFunc)
figure(1);
pdegplot(thermalModelS,'EdgeLabels','on')
axis equal
%xlim([-3 3])
title 'Pipe Geoemtry with Edge Labels Displayed'
T_ambient = 25; %degC
T_edge = 50; %degC
outerCC = @(region,~) 10*region.y
thermalBC(thermalModelS,'Edge',1,'ConvectionCoefficient',outerCC,'AmbientTemperature',T_ambient);
thermalBC(thermalModelS,'Edge',2,'ConvectionCoefficient',outerCC,'AmbientTemperature',T_ambient);
thermalBC(thermalModelS,'Edge',3,'ConvectionCoefficient',outerCC,'AmbientTemperature',T_ambient);
thermalBC(thermalModelS,'Edge',4,'Temperature',T_edge);
%leftHF = @(region,~) 5000*region.y %heat sourtce (W/m2)
%thermalBC(thermalModels,'Edge',4,'HeatFlux',leftHF) %BC heat sourtce (W/m2)
generateMesh(thermalModelS,'Hmax',2);
figure(2);
pdeplot( thermalModelS)
axis equal
title 'Pipe with tirangular Element Mesh'
result = solve(thermalModelS);
T = result.Temperature;
figure(3);
pdeplot(thermalModelS,'XYdata',T,'Contour','on');
axis equal
title 'Steady State Temperature'
