%% Use of PDEtool to solve the EN ISO 10211 Case 2
%% The European Standard EN ISO 10211: Thermal bridges in building construction - Heat flows and surface temperatures - Detailed calculations (ISO 10211:2017)
% Table of Contents
% The European Standard EN ISO 10211: Thermal bridges in building construction - Heat flows and surface temperatures - Detailed calculations (ISO 10211:2017)
% Introduction
% Validation of calculation methods
%     Case 2
%     References
% Use of the PDE Toolbox
%     Solve conduction heat transfer problems
% The MATLAB script
%     Model
%     Geometry
%     Thermophysical properties
%     Mesh generation
%         Uniform mesh
%         Locally refined mesh
%     Boundary conditions
%     Solution
%     Visualization
% Verify the calculation
%     Heat flow rate
%     Temperature
%     Heat flow and temperatures check
%% Introduction
% The European Standard EN ISO 10211 [1] sets out the specifications for a three-dimensional and two-dimensional geometrical model of a thermal bridge for the numerical calculation of:
% heat flows, in order to estimate the overal heat losses from a building
% minimum surface temperature, in order to assess the potential risk of surface condensation
% These specifications include the geometrical boundaries and subdivisions of the model, the thermal boundary conditions and the thermal values and relationships to be used.
% The document is based on the following assumptions:
% All physical properties are independent of temperature.
% There are not heat sources within the building element.
%
%% Validation of calculation methods
% In order to be classified as a three-dimensional steady-state high precision method, a calculation method shall give results corresponding to those of the test reference cases 1, 2, 3, and 4 reported in Annex C of the Standard.
% In order to be classified as a two-dimensional steady-state high precision method, it should give results corresponding to those of the test reference cases 1 and 2 reported in Annex C of the Standard.
% The Standard specifies also other requirements for the numerical method employed, like e.g. its capoability to provide both temperatures and fluid flow, the way to estimate the minimum number of subdivisions (nodes, finite elements, finite volumes et.), the criterion of convergence if iterative methods are emplyed for the solution of the corresponding system of equations et.
%% Case 2
% The Case 2 is relative to the calculation of temperature distribution and heat flow in a part of a wall (more precisely within a roof construction) depicted in the following figure (not to scale).
%
%
% The dimensions, thermal conductivity and boundary conditions, taken from the Standard, are reported in the following table.
%
% It is interesting to note the large range of thermal conductivity values, close to 4 order of magnitude.
% In addition, the Standard provides the exact temperature in 9 different points (marked in the previous figure) and the exact total heat flow rate per unit width, given in the following table.
%
% The Standard requires that the difference, between the calculated temperature values and those in the table, should not exceed 0.1 °C, AND that the difference between the heat flow calculated by the method being validated and that provided in the previous table, should not exceed 0.1 W/m.
%
% We remind that for this problem - steady state conduction heat transfer without heat generation - the governing equation is
% \nabla \cdot \left(\lambda \nabla t \right) = 0
% In this case, we want to solve the Case 2 of the Standard using the Finite Element based Partial Differential Equation Toolbox (PDE Toolbox) of MATLAB.
%% References
% [1] EN ISO 10211 - Thermal bridges in building construction - Heat flows and surface temperatures - Detailed calculations (ISO 10211:2017).
% [2] F. P. Incropera, D. P. Dewitt, T. L. Bergman, A. S. Lavine, Fundamentals of Heat and Mass Transfer, 6th Ed., Wiley, (2007).
%
%% Use of the PDE Toolbox
% The PDE Toolbox can be used to solve structural mechanics, heat transfer, and general partial differential equations (PDEs) using finite element analysis.
% PDE Toolbox lets to import 2D and 3D geometries from STL or mesh data. It can automatically generate meshes with triangular and tetrahedral elements. And finally it can solve PDEs by using the finite element method, and postprocess results to explore and analyze them.
% The are two ways to use the Toolbox:
% PDE Modeler App Workflow: it lets you solve the problem in an interactive way with the PDE Modeler App.
% Programmatic Workflow: it lets you solve the problem using high-level MATLAB commands programmatically.
% In the following, we will use the latter.
%
%% Solve conduction heat transfer problems
% A typical programmatic workflow for solving a heat transfer problem includes the following steps:
% Create a thermal model container for a steady-state or transient thermal model.
% Define 2-D or 3-D geometry and mesh it.
% Assign thermal properties of the material, such as thermal conductivity k, specific heat c, and mass density ρ.
% Specify internal heat sources q within the geometry, if applicable.
% Specify boundary conditions.
% Set an initial temperature or initial guess.
% Solve and plot results, such as the resulting temperatures, temperature gradients, heat fluxes, and heat rates.
%
%% The MATLAB script
% Preliminaries:

clear all;  % Removes all variables from the workspace.
close all;  % Deletes all figures.
clc;        % Clear Command Window.

%
%% Model
% Create a steady state PDE thermal model container, which we call "RoofConstruction"

RoofConstruction = createpde('thermal','steadystate'); 

%
%% Geometry
% Create the rectangles that represent the materials 1 and 2:
% For a rectangle, the first row contains 3 (indicates a rectangle) and the second row contains 4 (number of line segments). The next four rows contain the x-coordinates of the starting points of the edges, and the four rows after that contain the y-coordinates of the starting points of the edges.

R1 = [3; 4; 0.0; 0.5; 0.5; 0.0;...
    0.0415; 0.0415; 0.0475; 0.0475];
R2 = [3; 4; 0.0; 0.015; 0.015; 0.0;...
    0.0365; 0.0365; 0.0415; 0.0415];

% Create the polygons that represent the materials 3 and 4:
% For a polygon, the first row contains 2 (indicates a polygon) and the second row contains the number of line segments n. The next 3 through 3+n-1 lines contain the x-coordinate of edge starting points, and the following 3+n through 2*n+2 lines contain the y-coordinate of edge starting points.

P3 = [2; 6; 0.015; 0.015; 0.0015; 0.0015; 0.5; 0.5;...
    0.0415; 0.0350; 0.0350; 0.0015; 0.0015; 0.0415];
P4 = [2; 8; 0.015; 0.0; 0.0; 0.5; 0.5; 0.0015; 0.0015; 0.015; ...
    0.0365; 0.0365; 0.0; 0.0; 0.0015; 0.0015; 0.0350; 0.0350];

% Append extra zeros to the rectangles and the polygon with the smaller number of edges, so they have the same number of rows as the polygon with the highest number of edges:

R1 = [R1;zeros(length(P4) - length(R1),1)];
R2 = [R2;zeros(length(P4) - length(R2),1)];
P3 = [P3;zeros(length(P4) - length(P3),1)];

% Combine the two rectangles and the two polygons into one matrix:

gd = [R1,R2,P3,P4];

% Create names for the two rectangles and the two polygons:

ns = char('R1','R2','P3','P4');
ns = ns';

% Define the final geometry as the union of the rectangles and polygons:

sf = '(R1+R2+P3+P4)';

% Create the geometry:

g = decsg(gd,sf,ns);

% Convert the decsg format into a geometry object. Include the geometry in the model.

geometryFromEdges(RoofConstruction,g);

% For convenience and better readability, we plot the geometry in two different figures: in the first we display the face labels, and in the second the edge labels. The face labels will  be used below for defining the values of the thermal conductivity, and the edge labels for setting the boundary conditions.

figure(1);
pdegplot(RoofConstruction,'EdgeLabels','off','FaceLabels','on');
figure(2);
pdegplot(RoofConstruction,'EdgeLabels','on','FaceLabels','off');

%
%% Thermophysical properties
% Assign thermal conductivity to each subdomain:

thermalProperties(RoofConstruction,'Face',1,'ThermalConductivity',0.12);  % Region 2 in ISO Case 2
thermalProperties(RoofConstruction,'Face',2,'ThermalConductivity',0.029); % Region 3 in ISO Case 2
thermalProperties(RoofConstruction,'Face',3,'ThermalConductivity',230.0); % Region 4 in ISO Case 2
thermalProperties(RoofConstruction,'Face',4,'ThermalConductivity',1.15);  % Region 1 in ISO Case 2

%
%% Mesh generation
% We create the mesh, using either one of the two options, uniform or non-uniform:

meshType = "uniform"
%meshType = 'nonuniform';
switch meshType
    case 'uniform'

%% Uniform mesh
% We specify:
% 'GeometricOrder' - the order (of the shape function) of the element: it can be 'linear' or 'quadratic'
% 'Hmax' -  the target maximum element size

generateMesh(RoofConstruction,'Hmax',0.003,'GeometricOrder','Quadratic'); % Create the mesh
    case 'nonuniform'

%% Locally refined mesh
% We specify:
% 'GeometricOrder' - the order (of the shape function) of the element: it can be 'linear' or 'quadratic'
% 'Hmax' -  the target maximum element size for the entire domain
% 'Hedge' -the target size on different edges, in particular along the edges of the thin metal plate


generateMesh(RoofConstruction,'Hgrad',1.4,'Hmax',0.02,'Hmin',0.002,...
    'Hedge',{[2 5 6 12 13 14 15],0.001,[7 16],0.002},...
    'GeometricOrder','quadratic');
end

% Display the mesh

figure(3);
pdeplot(RoofConstruction);  % View the mesh
xlim([0.0 0.5]);
ylim([0.0 0.05]);

% It is useful to know the number of nodes ('nnodes') and elements ('mtriangles')

[points,edges,tri] = meshToPet(RoofConstruction.Mesh);
nnode = size(points,2);
elements = size(tri,2);
fprintf('\nNo. Nodes = %d     No. elements = %d', nnode, elements)

% The mesh has to be rather fine in order to comply with the ISO Standard, since using the PDEToolbox there are no options to specify different/varying mesh resolutions for the subdomains and/or adaptive meshing
%
%% Boundary conditions

% 1. Convective BC on Edge 1 (Edge AB in the standard):
thermalBC(RoofConstruction,'Edge',1,'ConvectionCoefficient',1/0.06,'AmbientTemperature',0);
% 2. Convective BC on Edge 2 (Edge HI in the standard):
thermalBC(RoofConstruction,'Edge',2,'ConvectionCoefficient',1/0.11,'AmbientTemperature',20);

% There is no need to specify the insulated BC (zero heat flux) on the other boundaries, since this is the default.
%
%% Solution
% We can now solve the problem

result = solve(RoofConstruction);

%
%% Visualization
% Visualize a contour plot of the temperature

T = result.Temperature;
figure(4);
pdeplot(RoofConstruction,'XYData',T,'Contour','on','ColorMap','Hot');
xlim([0.0 0.5]);
ylim([0.0 0.05]);

%
%% Verify the calculation
% In order to verify the calculation, we have to evaluate the total heat flow rate AND  the temperature in some selected points, and check the difference between the results of the calculation and the data provided by the Standard.
%% Heat flow rate
% Total heat flow rate (W/m) from the Standard:

TotalHeatFlowRateIso10211Case2 = 9.5;

% We evaluate the heat flow rate on both two convective boundaries (all the others are perfectly insulated):

TotalHeatFlowRate1 = evaluateHeatRate(result,'Edge',1);
TotalHeatFlowRate2 = evaluateHeatRate(result,'Edge',2);

% We then take the maximum difference for conservative reasons:

HeatFlowDiff = max([abs(abs(TotalHeatFlowRate1) - TotalHeatFlowRateIso10211Case2)...
    abs(abs(TotalHeatFlowRate2) - TotalHeatFlowRateIso10211Case2)]);
fprintf('Max Heat Flow difference: %7.3f\n', HeatFlowDiff);

%% Temperature
% X and Y coordinates of the UNI EN ISO 10211 Case 2 selected points:

xIso10211Case2 = [0.0 0.5 0.0 0.015 0.5 0.0 0.015 0.0 0.5];
yIso10211Case2 = [0.0475 0.0475 0.0415 0.0415 0.0415 0.0365...
    0.0365 0.0 0.0];

%  We can then define a matrix of query points corresponding to the selected points:

coordIso10211Case2 = [xIso10211Case2(:),yIso10211Case2(:)]';

% Temperature values of the UNI EN ISO 10211  Case 2 selected points:

TempIso10211Case2 =  [7.1 0.8 7.9 6.3 0.8 16.4 16.3 16.8 18.3];

% We can now interpolate the temperature on the just defined points:

Tintrp = interpolateTemperature(result,coordIso10211Case2);

% Temperature difference:

Tdiff = abs(TempIso10211Case2-Tintrp');
dTmax = max(max(abs(Tintrp' - TempIso10211Case2)));
fprintf('Max DT [K]: %7.3f\n', dTmax)

%% Heat flow and temperatures check
% Check if differences in heat flow rate and temperatures are within the limits of the Standard:

HeatFlowIsOK  = HeatFlowDiff < 0.1; % The difference should be less than 0.1 W/m
diffTIsOK = Tdiff < 0.1;            % The difference should be less than 0.1 °C
if(HeatFlowIsOK & diffTIsOK);
    fprintf('Verification OK !\n')
else
    fprintf('Verification failed !\n')
end