%% Use of PDEtool to solve the EN ISO 10211 Case 1
%% 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 1
%     Analytical solution
%     References
% Use of the PDE Toolbox
%     Solve conduction heat transfer problems
% The MATLAB script
% Verify the calculation
%     Graphical representation of the verification
%% 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 1
% The Case 1, the most simple of the four cases, is particularly interesting, since it admits an analytical solution. It is relative to the heat transfer through a half (for symmetry reasons) square column, with assigned surface temperature. The Standard, in particular, indicate, as depicted in the figure below, the values of the temperature provided by the analytical solution on an equidistant grid of 28 points.
%
%
% The corresponding values are reported in the following table:
%
% The Standard requires that the maximum difference, between the calculated temperature values and those in the table, should not exceed 0.1 °C.
%
%% Analytical solution
% 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
% which, in our case - 2D with constant properties - becames
% \frac{\partial^{2} t}{\partial x^{2}}+\frac{\partial^{2} t}{\partial y^{2}} = 0
% The analytical solution, not provided by the Standard but easily found in the scientific literature [2], is given, in non-dimensional form for the entire square column of side L, by
% \theta\left(x,y\right) = \frac{2}{\pi}\,\sum_{n=1}^{\infty}
%                          \frac{\left(-1\right)^{n+1}+1}{n}
%                          \sin \frac{n\pi x}{L}
%                          \frac{\sinh \left(n\pi y/L\right)}{\sinh\left(n\pi\right)}
% \]
% where \theta is the non-dimensional temperature
% \theta \equiv \frac{T - T_{1}}{T_{2} - T_{1}}
% In this case, we want to solve the Case 1 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 the problem
% Plot the 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.

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

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

% Create the rectangle that represents the half-column:
% For a rectangle, the first row contains 3, and the second row contains 4. 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.

rectHalfColumn = [3; 4; 0.0; 0.5; 0.5; 0.0;...
    0.0; 0.0; 1.0; 1.0];

%
% Create and view the geometry

g = decsg(rectHalfColumn);
geometryFromEdges(HalfColumn,g);
% View the geometry
figure(1);
pdegplot(HalfColumn,'EdgeLabels','on');
axis equal;

%
% Assign thermal conductivity of the material. It is not required in this specific case, but it is provided here for completeness (can be any arbitrary value).

thermalProperties(HalfColumn,'ThermalConductivity',50.);

% Specify boundary conditions

% 1. Dirichlet BC T = 0 [°C] on Edges 1 and 4
thermalBC(HalfColumn,'Edge',1,'Temperature',0);
thermalBC(HalfColumn,'Edge',4,'Temperature',0);
% 2. Dirichlet BC T = 20 [°C] on Edge E3
thermalBC(HalfColumn,'Edge',3,'Temperature',20);
% 3. Neumann (simmetry, e.g. zero heat flux) on Edge 2
thermalBC(HalfColumn,'Edge',2,'HeatFlux',0);

%
% We create the mesh, specifying:
% 'Hmax' -  the maximum mesh edge length
% 'GeometricOrder' - the order (of the shape function) of the element: it can be 'linear' or 'quadratic'

Hmax = 0.08
Order = "Linear"
generateMesh(HalfColumn,'Hmax',Hmax,'GeometricOrder',Order); % Create the mesh
figure(2);
pdeplot(HalfColumn);  % View the mesh
xlim([0.0 0.5]);
ylim([0.0 1.0]);

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

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

%
% We can now solve the problem

result = solve(HalfColumn);

% Visualize a contour plot of the temperature

T = result.Temperature;
figure(3);
pdeplot(HalfColumn,'XYData',T,'Contour','on','ColorMap','Hot');
xlim([0.0 0.5]);
ylim([0.0 1.0]);

%
%% Verify the calculation
% In order to verify the calculation, we have to interpolate our results on the equidistant grid given by the standard, and check the difference between the results of the calculation and the data provided by the Standard.
% X and Y coordinates of the UNI EN ISO 10211 equidistant grid:

x_iso10211Case1 = [0.125 0.25 0.375 0.5];
y_iso10211Case1 = [0.875 0.75 0.625 0.5 0.375 0.25 0.125];

% Values of the temperature from UNI EN ISO 10211 Case 1:

Tiso10211Case1 = [9.7 13.4 14.7 15.1;
             5.3  8.6 10.3 10.8;
             3.2  5.6  7.0  7.5;
             2.0  3.6  4.7  5.0;
             1.3  2.3  3.0  3.2;
             0.7  1.4  1.8  1.9;
             0.3  0.6  0.8  0.9];

%
% We now use the MATLAB function 'meshgrid' to create the arrays of X and Y coordinates of all nodes of the equidistant grid:
% [Xintrp,Yintrp] = meshgrid(x_iso10211,y_iso10211) returns 2-D grid coordinates based on the coordinates contained in vectors x_iso10211 and y_iso10211.
% Xintrp is a matrix where each row is a copy of x_iso10211, and Yintrp is a matrix where each column is a copy of y_iso10211.

[Xintrp,Yintrp]= meshgrid(x_iso10211Case1,y_iso10211Case1);

% We can now interpolate the temperature on the just defined equidistant grid:

Tintrp = interpolateTemperature(result,Xintrp,Yintrp);

% It is convenient to reshape Tintrp into the equidistant 2D grid in order to facilitate the validation

Tintrp = reshape(Tintrp,size(Xintrp));

%
% Check if the difference is within the limits of the Standard

Tdiff = abs(Tiso10211Case1-Tintrp);
diffT = Tdiff < 0.1;
dTmax = max(max(abs(Tintrp - Tiso10211Case1)));
fprintf('Max DT [K]: %7.3f\n', dTmax)
isOK  = all(all(diffT)); % The difference should be less than 0.1 °C
if(isOK)
    fprintf('Verification OK !\n')
else
    fprintf('Verification failed !\n')
end

%
%% Graphical representation of the verification
% We represent graphically the results of the calculations and the values from the ISO Standard on the equidistant grid:
% Values from the Standard in black;
% Computed values in bleu if they do comply with the Standard, in red if not.

x_box = [0 0.5 0.5 0.0 0.0];
y_box = [0 0.0 1.0 1.0 0.0];
%
lcross = 0.04; % Length of the segment used to plot the small cross at the given grid points
figure('Name','Reference case 1: Temperature values')
hold on;
axis equal;
axis([-0.1 0.6 -0.1 1.1]);
line(x_box,y_box,'Color','r','LineWidth',2)
for igrid = 1:length(x_iso10211Case1)
    for jgrid = 1:length(y_iso10211Case1)
        x_left  = x_iso10211Case1(igrid) - lcross/2;
        x_right = x_iso10211Case1(igrid) + lcross/2;
        line([x_left x_right], ...
            [y_iso10211Case1(jgrid) y_iso10211Case1(jgrid)],...
            'Color','black','LineWidth',2)
        y_down  = y_iso10211Case1(jgrid) - lcross/2;
        y_up    = y_iso10211Case1(jgrid) + lcross/2;
        line([x_iso10211Case1(igrid) x_iso10211Case1(igrid)], ...
            [y_down y_up],...
            'Color','black','LineWidth',2)
        text(x_left-2*lcross/2,y_up,...
            num2str(Tiso10211Case1(jgrid,igrid),'%6.1f'),...
            'FontSize',12,'FontName','Times')
        if(diffT(length(y_iso10211Case1)-jgrid+1,igrid))
            tcolor = 'blu';
        else
            tcolor = 'red';
        end
        text(x_left-2*lcross/2,y_down,...
            num2str(Tintrp(jgrid,igrid),'%6.2f'),...
            'FontSize',12,'FontName','Times','Color',tcolor)
    end
end