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

clear all; clc; clf; close all;

thermalmodel = createpde('thermal');
R1 = [3,4,-1,1,1,-1,1,1,-1,-1]';
g = decsg(R1,'R1',('R1')');
geometryFromEdges(thermalmodel,g);
figure(1);
pdegplot(thermalmodel,'EdgeLabels','on')
xlim([-1.5 1.5]); axis equal;
title('Geometry with Edge Labels Displayed');
thermalProperties(thermalmodel,'ThermalConductivity',79.5,'Face',1);
thermalBC(thermalmodel,'Edge',3,'Temperature',300);
thermalBC(thermalmodel,'Edge',1,'HeatFlux',0);
thermalBC(thermalmodel,'Edge',[2 4], ...
                       'ConvectionCoefficient',25, ...
                       'AmbientTemperature',50);

                   generateMesh(thermalmodel);
results = solve(thermalmodel)
figure(2);
pdeplot(thermalmodel,'FlowData',[results.XGradients results.YGradients]);
v = linspace(-0.5,0.5,11);
[X,Y] = meshgrid(v);
title('Temperature Gradient Flow Data');

[gradTx,gradTy] = evaluateTemperatureGradient(results,X,Y);
gradTx = reshape(gradTx,size(X));
gradTy = reshape(gradTy,size(Y));
figure(3);
quiver(X,Y,gradTx,gradTy)
querypoints = [X(:) Y(:)]';
[gradTx,gradTy] = evaluateTemperatureGradient(results,querypoints);
title('Temperature Gradient Flow');

gradTx = reshape(gradTx,size(X));
gradTy = reshape(gradTy,size(Y));
figure(4);
quiver(X,Y,gradTx,gradTy)
title('Temperature Gradient Flow');