%library for graphic part
classdef FEM
    methods (Static)
    
    

    %% ************************************************************************************************************
    %  ************				  		        lengthElementMesh              	 		               ************
    %  ************************************************************************************************************
    % This function calculates the length of each finite element and adds a field to the structure MESH.
    % MESH : is a structure composed of the following elements:
		% 	MESH.nodes = matrix that contains in the k-th row the coordinates of the node with label k
		% 	MESH.elem = matrix that contains in the j-th row the labels of the nodes of the j-th element 
		% 	MESH.dof =  vector containing in the k-th position the dofs of the k-th node
		% 	MESH.mat =  vector containing in the k-th position the material of the k-th element
		%	MESH.angle = vector containing in the k-th position the angle of rotation of the k -th element
		%	MESH.length = vector containing in the k-th position the length of the k -th element	
    
        function[MESH] = lengthElementMesh(MESH)
        MESH.length = [];
        for k = 1 : size(MESH.elem, 1)
            %coordinates of the nodes of the k-th finite element
            node1 = MESH.nodes(MESH.elem(k, 1), :);
            node2 = MESH.nodes(MESH.elem(k, 2), :);

            %update the vector containing the length of the elements
            MESH.length = [MESH.length ((node2 - node1) * (node2 - node1)')^(1 / 2)];
        end
    end    




    %% ************************************************************************************************************
    %  ************				  		        angleElementMesh              	 		               ************
    %  ************************************************************************************************************
    % This function calculates the angle of rotation of each finite element about global z-axis and adds a field to
    % the structure MESH.
    % MESH : is a structure composed of the following elements:
		% 	MESH.nodes = matrix that contains in the k-th row the coordinates of the node with label k
		% 	MESH.elem = matrix that contains in the j-th row the labels of the nodes of the j-th element 
		% 	MESH.dof =  vector containing in the k-th position the dofs of the k-th node
		% 	MESH.mat =  vector containing in the k-th position the material of the k-th element
		%	MESH.angle = vector containing in the k-th position the angle of rotation of the k -th element
		%	MESH.length = vector containing in the k-th position the length of the k -th element	
    
        function[MESH] = angleElementMesh(MESH)
        MESH.angle = [];
        for k = 1 : size(MESH.elem, 1)
            %coordinates of the nodes of the k-th finite element
            node1 = MESH.nodes(MESH.elem(k, 1), :);
            node2 = MESH.nodes(MESH.elem(k, 2), :);

            %compute the difference of in nodal coordinates
            deltaY = node2(2) - node1(2);
            deltaX = node2(1) - node1(1);

            %update the vector containing the angle of the elements with respect to the global reference frame
            MESH.angle = [MESH.angle atan2(deltaY, deltaX)];
        end
    end



    %% ************************************************************************************************************
    %  ************       						      printStepPlus       					           ************
    %  ************************************************************************************************************
    %This function print to the screen (Matlab prompt) the current step in a for cycle, deleting automatically
    %the previous step. It provides a nicer way of plotting the advance of the code, instead of printing all
    %the cycle indexes. This version shows also the total number of cycle that are required.
    %index : the numerical intex that show the number fo iteration that have been performed
    %maxIndex : the maximum value that the index assumes in the cycle
    %lengthOld : the length of the number printed before, it must be initialized before this first cycle with
    %		a zero value

    function [lengthNew] = printStepPlus(index, maxIndex, lengthOld)
        indexToPrint = num2str(index);
        maxIndexToPrint = num2str(maxIndex);

        strToPrint = strcat(indexToPrint, strcat('/', maxIndexToPrint));
        lengthNew = length(strToPrint);

        strToDelete = [];
            for i = 1 : lengthOld
                strToDelete = strcat(strToDelete, '\b');
            end
        if isempty(strToDelete) == false
            fprintf(strToDelete);
        end
        fprintf(strToPrint);
    end
    
    
    
    
    
    %% ************************************************************************************************************
    %  ************				  		            dofsNode              	 		                   ************
    %  ************************************************************************************************************
    % This function defines the position of the dofs of a given node in the global vector of generalized 
    % displacements U.
    % MESH : structured data containing the information of the mesh 
    % nodeLabel : label of the considered node
    function[dofs] = dofsNode(MESH, nodeLabel)
        dofs = (FEM.sumVectorElements(MESH.dof, nodeLabel - 1) + 1) :1: FEM.sumVectorElements(MESH.dof, nodeLabel);
    end




    %% ************************************************************************************************************
    %  ************				  		        extractDisplacements              	 		           ************
    %  ************************************************************************************************************
    % This extracts the coordinates of the nodes in the deformed configuration (considering only translations form
    % the global vector U).
    % MESH : structured data containing the information of the mesh
    % U : global vector of generalized displacements of the entire structure (solution of the FEM code)
    % xU, yU : global coordinates of the nodes in the deformed configuration (call scatter(xU, yU) to plot the nodes)
    function[xU, yU] = extractDisplacements(U, MESH)
        xU = [];
        yU = [];
        for i = 1 : size(MESH.nodes, 1)
            dofs = FEM.dofsNode(MESH, i);
            xU = [xU (MESH.nodes(i, 1) + U(dofs(1)))];
            yU = [yU (MESH.nodes(i, 2) + U(dofs(2)))];
        end
    end




    %% ************************************************************************************************************
    %  ************				  		        defineStiffnessMatrix          	 		               ************
    %  ************************************************************************************************************
    % This function defines the matlab function (function obtained from a symbolic function, but faster for 
    % numerical computation) of the stiffness matrix, given the formulation of the finite element.
    % codeEF : string defining the formulation of the finite element ("TRUSS" for truss element, "EULER_BERNOULLI", 
    %for complete extensible Euler-Bernoulli beam element)
    % KEF : numeric matlab function of the stiffness matrix of the finite element in the local reference frame,
    % the stiffness matrix can be calculated passing to KEF the geometrical and mechanical parameters 
    % KEF(L, [E, rho], [A, J])
    
    function[KEF] = defineStiffnessMatrix(codeEF)
        %define the symbols used in the matrix formulation
        syms E A J L rho
        
        %define the symbolic matrix
        switch codeEF
            case 'TRUSS'
                %stiffness matrix for a truss element in the local reference system (dofs: [u1, u2])
                KEFs = (E * A / L) * [1, -1;  -1, 1];

            case 'EULER_BERNOULLI'
                %stiffness matrix for an extensible Euler-Bernoulli element in the local reference system 
                %(dofs: [u1, v1, theta1, u2, v2, theta2])
                KEFs = [E * A / L,      0,                       0,                      -E * A / L,     0,                     0;
                        0,              12 * E * J / (L^3) ,     6 * E * J / (L^2),      0,         	 -12 * E * J / (L^3), 	6 * E * J / (L^2);
                        0,              6 * E * J / (L^2), 	     4 * E * J / L,          0,         	 -6 * E * J / (L^2), 	2 * E * J / L;
                        -E * A / L,     0,                       0,                      E * A / L,      0,                     0;
                        0,              -12 * E * J / (L^3), 	 -6 * E * J / (L^2),     0,        		 12 * E * J / (L^3), 	-6 * E * J / (L^2);
                        0,              6 * E * J / (L^2), 	     2 * E * J / L,          0,         	 -6 *E * J / (L^2), 	4 * E * J / L];         
        end

        %definition of the vectors of parameters
        materialParameters = [E, rho];
        sectionParameters = [A, J];

        %define the matlab function of the stiffness matrix
        KEF = matlabFunction(KEFs, 'Vars', {L, materialParameters, sectionParameters});
    end




    %% ************************************************************************************************************
    %  ************				  		        defineMassMatrix          	 		               ************
    %  ************************************************************************************************************
    % This function defines the matlab function (function obtained from a symbolic function, but faster for 
    % numerical computation) of the mass matrix, given the formulation of the finite element.
    % codeEF : string defining the formulation of the finite element ("TRUSS" for truss element, "EULER_BERNOULLI", 
    %for complete extensible Euler-Bernoulli beam element)
    % MEF : numeric matlab function of the stiffness matrix of the finite element in the local reference frame,
    % the stiffness matrix can be calculated passing to MEF the geometrical and mechanical parameters 
    % MEF(L, [E, rho], [A, J])
    
    function[MEF] = defineMassMatrix(codeEF)
        %define the symbols used in the matrix formulation
        syms rho E A J L
        
        %define the symbolic matrix
        switch codeEF
            case 'TRUSS'
                %stiffness matrix for a truss element in the local reference system (dofs: [u1, u2])
                MEFs = (rho * A * L / 6) * [2, 1;  1, 2];

            case 'EULER_BERNOULLI'
                %stiffness matrix for an extensible Euler-Bernoulli element in the local reference system 
                %(dofs: [u1, v1, theta1, u2, v2, theta2])
                MEFs = rho * A * L *...
                        [1 / 3,      0,             0,                  1 / 6,      0,              0;
                        0,           13/35 ,        L * 11/210,         0,         	9/70,           -L * 13/420;
                        0,           L * 11/210, 	(L^2) /105,         0,         	L * 13/420, 	-(L^2) * 1/140;
                        1 / 6,       0,             0,                  1 / 3,      0,              0;
                        0,           9/70,          L * 13/420,         0,        	13/35,          -L * 11/210;
                        0,           -L * 13/420, 	-(L^2) * 1/140,     0,         	-L * 11/210, 	(L^2) * 1/105];         
        end

        %definition of the vectors of parameters
        materialParameters = [E, rho];
        sectionParameters = [A, J];

        %define the matlab function of the stiffness matrix
        MEF = matlabFunction(MEFs, 'Vars', {L, materialParameters, sectionParameters});
    end
    
    
    
    
    
    %% ************************************************************************************************************
    %  ************				  		        defineRotationMatrix          	 		               ************
    %  ************************************************************************************************************
    % This function defines the matlab function (function obtained from a symbolic function, but faster for 
    % numerical computation) of the rotation matrix, given the formulation of the finite element.
    % codeEF : string defining the formulation of the finite element ("TRUSS" for truss element, 
    % "EULER_BERNOULLI", for complete extensible Euler-Bernoulli beam element)
    % T : numeric matlab function of the rotation matrix transforming the stiffness matrix of the finite element  
    % from the local reference frame to the global one, it can be calculated passing to T the angle of rotation of  
    % the finite element T(angle)

    % NOTE: the local reference frame is u (positive when right-oriented), v (positive when down-oriented) and theta
    % (positive when clockwise oriented); the global reference frame is the usual 0xy system with positive 
    % counter-clockwise direction. The global and local reference system are not simply rotated, but there is also
    % an inversion in axis y and in the rotation.
    
    function[T] = defineRotationMatrix(codeEF)
        %define the symbols used in the matrix formulation
        syms angleRot
        
        %define the symbolic matrix performing the rotation from local to global reference frame
        switch codeEF
            case 'TRUSS'
                %rotation matrix for a truss element from local to global reference frame 
                %(local dofs: [u1, u2] --> global dofs: [u1, v1, u2, v2])
                Ts = [	cos(angleRot),     sin(angleRot),         0,              0; 
                        0,			    0,					cos(angleRot), 	sin(angleRot)];

            case 'EULER_BERNOULLI'
                %rotation matrix for an extensible Euler-Bernoulli element from local to global reference frame 
                %(local dofs: [u1, v1, theta1, u2, v2, theta2] --> global dofs: [u1, v1, theta1, u2, v2, theta2])
                Ts = [	cos(angleRot), 	    sin(angleRot),	    0,			0,			        0,			    0;
                        sin(angleRot),     -cos(angleRot), 	    0,			0,			        0,			    0;
                        0,		 	 	    0,			        -1,			0,			        0,			    0;
                        0,				    0,			        0,			cos(angleRot), 	    sin(angleRot), 	0;
                        0,				    0,			        0,			sin(angleRot),     -cos(angleRot), 	0;
                        0,				    0,			        0,			0,		 	        0,			    -1];
                    
        end
        %define the matlab function of the stiffness matrix
        T = matlabFunction(Ts, 'Vars', {angleRot});
    end




    %% ************************************************************************************************************
    %  ************				  		    defineRotationMatrixPostProcessing          	 		   ************
    %  ************************************************************************************************************
    % This function defines the matlab function (function obtained from a symbolic function, but faster for 
    % numerical computation) of the rotation matrix used in the post-processing, given the formulation of the 
    % finite element.
    % codeEF : string defining the formulation of the finite element ("TRUSS" for truss element, 
    % "EULER_BERNOULLI", for complete extensible Euler-Bernoulli beam element)
    % T : numeric matlab function of the rotation matrix transforming the stiffness matrix of the finite element  
    % from the local reference frame to the global one, it can be calculated passing to T the angle of rotation of  
    % the finite element T(angle)
    
    % NOTE : the matrix T here is not the same of the one calulated by FEM.defineRotationMatrix, because when the
    % following function is considered, this output matrix rotates only the axes, without inversion in axis y or 
    % angle theta.

    function[T] = defineRotationMatrixPostProcessing(codeEF)
        %define the symbols used in the matrix formulation
        syms angleRot
        
        %define the symbolic matrix performing the rotation from local to global reference frame
        switch codeEF
            case 'TRUSS'
                %rotation matrix for a truss element from local to global reference frame 
                %(local dofs: [u1, u2] --> global dofs: [u1, v1, u2, v2])
                Ts = [	cos(angleRot), 	sin(angleRot),	    0,			        0;
                        -sin(angleRot), cos(angleRot), 	    0,			        0;
                        0,				0,			        cos(angleRot), 	    sin(angleRot);
                        0,				0,			        -sin(angleRot),     cos(angleRot);];

            case 'EULER_BERNOULLI'
                %rotation matrix for an extensible Euler-Bernoulli element from local to global reference frame 
                %(local dofs: [u1, v1, theta1, u2, v2, theta2] --> global dofs: [u1, v1, theta1, u2, v2, theta2])
                Ts = [	cos(angleRot), 	sin(angleRot),	0,			0,			    0,			    0;
                        -sin(angleRot), cos(angleRot), 	0,			0,			    0,			    0;
                        0,		 	 	0,			    1,			0,			    0,			    0;
                        0,				0,			    0,			cos(angleRot), 	sin(angleRot), 	0;
                        0,				0,			    0,			-sin(angleRot), cos(angleRot), 	0;
                        0,				0,			    0,			0,		 	    0,			    1];
        end

        %define the matlab function of the stiffness matrix
        T = matlabFunction(Ts, 'Vars', {angleRot});
    end




    %% ************************************************************************************************************
    %  ************				  		            defineElementLoad          	 		               ************
    %  ************************************************************************************************************
    % This function defines the matlab function (function obtained from a symbolic function, but faster for 
    % numerical computation) of the rotation matrix used in the post-processing, given the formulation of the 
    % finite element. The distributed load can be linear, p1 and p2 are the values of the axial distributed load
    % in first and second node of the given finite element, q1 and q2 are the values of the transversal distributed
    % load in first and second node of the given finite element.
    % codeEF : string defining the formulation of the finite element ("TRUSS" for truss element, "EULER_BERNOULLI"
    % for complete extensible Euler-Bernoulli beam element)
    % fLin : numeric matlab function of the vector of generalized loads in the local reference frame,
    % it can be calculated passing to fLin the load parameters fLin(L, [p1, p2, q1, q2])

    function[fLin, fTermAx, fTermCurv] = defineElementLoad(codeEF)
        %define the symbols used in the matrix formulation
        syms angleRot p1 p2 q1 q2 L
        syms alphaT h DT
        syms E A J
        
        %define the symbolic matrix performing the rotation from local to global reference frame
        switch codeEF
            case 'TRUSS'
                %vector of generalized loads for a truss element (local dofs: [u1, u2])
                %(p1, p2 positive when directed as u1 and u2)
                fDistrS = L * [(p1 / 3) + (p2 / 6); (p1 / 6) + (p2 / 3)];
                
                %vector of thermal load for a constant dilatation
                fTermAxS = alphaT * DT * E * A * [-1; 1];

                %vector of thermal load in case of constant thermal curvature
                %the sign of the curvature must choose according to the considered reference system
                fTermCurvS = [0; 0] * L;

            case 'EULER_BERNOULLI'
                %vector of generalized loads for an extensible Euler-Bernoulli element 
                %(local dofs: [u1, v1, theta1, u2, v2, theta2])
                %(q1, q2 positive when directed as v1 and v2)
%                 fDistrS = [  (L * p1 / 3) + (L * p2 / 6); 
%                         (7 * L * q1 / 20) + (3 * L * q2 / 20);
%                         -(L^2 * q1 / 20) - (L^2 * q2 / 30);
%                         (L * p1 / 6) + (L * p2 / 3);
%                         (3 * L * q1 / 20) + (7 * L * q2 / 20);
%                         (L^2 * q1 / 30) + (L^2 * q2 / 20) ];
                    
                fDistrS = [ (L * p1 / 3) + (L * p2 / 6); 
                            (7 * L * q1 / 20) + (3 * L * q2 / 20);
                            (L^2 * q1 / 20) + (L^2 * q2 / 30);
                            (L * p1 / 6) + (L * p2 / 3);
                            (3 * L * q1 / 20) + (7 * L * q2 / 20);
                            -(L^2 * q1 / 30) - (L^2 * q2 / 20) ];
                
                %vector of thermal load for a constant dilatation
                fTermAxS = alphaT * DT * E * A * [-1; 0; 0; 1; 0; 0;];
                
                %vector of thermal load in case of constant thermal curvature
                %the sign of the curvature must choose according to the considered reference system
                fTermCurvS = (2 * alphaT * DT * E * J / h) * [0; 0; -1; 0; 0; 1];
        end

        %definition of the vectors of parameters
        loadParametersDistr = [p1, p2, q1, q2];
        loadParametersTerm = [alphaT, h, DT];
        mechanicalParametersTerm = [E, A, J];

        %define the matlab function of the stiffness matrix
        fLin = matlabFunction(fDistrS, 'Vars', {L, loadParametersDistr});
        fTermAx = matlabFunction(fTermAxS, 'Vars', {L, loadParametersTerm, mechanicalParametersTerm});
        fTermCurv = matlabFunction(fTermCurvS, 'Vars', {L, loadParametersTerm, mechanicalParametersTerm});
    end




    %% ************************************************************************************************************
    %  ************       					   	  sumVectorElements       	 			               ************
    %  ************************************************************************************************************
    % this function sums the elements of the vector v, from the first element to the element at the position pos(k)
    %v : vector whose elements have to be summed, v = [v1, v2, v3, ...]
    %pos: vector of the ending positions for the summation, pos = [pos1, pos2,...]

    function [s] = sumVectorElements(v, pos) 
        for i = 1 : length(pos)
            s(i) = sum(v(1 : pos(i)));
        end
    end
           
    


    %% ************************************************************************************************************
    %  ************       				 linkLocalToGlobalDofs		       					           ************
    %  ************************************************************************************************************
    % This function creates the matrix matrixGlobalDofs containing the degrees of freedom related to given 
    % element. It can be convenient to pre-calculate this matrix when the number of the element is not so high, 
    % because less calculation can be done calling the assembly function
    % MESH : structured data containing the information of the mesh 
    % matrixGlobalDofs : matrix defining the dofs that concern a given finite element, i.e. in the k-th row 
    %			there are the dofs that are in the k-th element

    function[matrixGlobalDofs] = linkLocalToGlobalDofs(MESH)
        %extract form the mesh the vector of dofs
        dof = MESH.dof;
        
        %calculate the vector of the global-local dof for each finite element
        for e = 1 : size(MESH.elem, 1)
            elemEF = MESH.elem(e, :);
            %this matrix is used to define the starting and ending dof associated to the considered Finite Element
            %in the global matrix, in the first (second) row there are the starting (ending) dofs for all the nodes
            %in the Finite Element
            matrixPosition = [FEM.sumVectorElements(dof, elemEF - 1) + 1;  FEM.sumVectorElements(dof, elemEF)];
            
            %calculating the vector of the dofs in the global matrix, related to the considered Finite Element
            vGlobalDofs = [];
            for i = 1 : length(elemEF) 
                vGlobalDofs = [vGlobalDofs (matrixPosition(1, i) : matrixPosition(2, i))];
            end
            
            %define the row of the matrix of global dofs
            matrixGlobalDofs(e, :) = vGlobalDofs;
        end
    end


    

    %% ************************************************************************************************************
    %  ************       				 		assemble		       					               ************
    %  ************************************************************************************************************
    % This function assemble the matrix or the vector of the single Finite Element aEF in the global matrix/vector. 
    % This function is generic and can be used with every combination of elements and it is independent of the
    % number of nodes and dofs for each node.
    % aGlobal : global matrix
    % aEF : matrix of the single element
    % matrixGlobalDofs : matrix defining the dofs that concern a given finite element, i.e. in the k-th row 
    %			there are the dofs that are in the k-th element
    % el : label of the considered element

    function[aGlobal] = assemble(aGlobal, aEF, matrixGlobalDofs, el)
        %extract the global dofs of the considered finite element
        vGlobalDofs = matrixGlobalDofs(el, :);
        
        if isvector(aEF) == 1
            %assemble the vector aEF in the global matrix aGlobal
            aGlobal(vGlobalDofs) = aGlobal(vGlobalDofs) + aEF;
        else
            %assemble the matrix aEF in the global matrix aGlobal
            aGlobal(vGlobalDofs, vGlobalDofs) = aGlobal(vGlobalDofs, vGlobalDofs) + aEF;
        end
    end




    %% ************************************************************************************************************
    %  ************       				 finiteElementAnalysis		       					           ************
    %  ************************************************************************************************************
    % This function create the global mass and stiffness matrix for the finite element analysis, it is the "core" 
    % of the code and it recalls the definition of the shape functions for the approximation of the solution.
    % This is the version of the function with pre-defined solution anonymous functions.
    % MESH : structured data containing the information of the mesh 
    % matrixGlobalDofs : matrix defining the dofs that concern a given finite element, i.e. in the k-th row 
    %			there are the dofs that are in the k-th element
    % Kfun : matlab function of the stiffness matrix
    % Mfun : matlab function of the mass matrix
    % Tfun : matlab function of the rotation matrix

    function[KGlobal, MGlobal, KEFGlobal] = finiteElementAnalysis(MESH, matrixGlobalDofs, Kfun, Mfun, Tfun)
        %define the number of elements for the calculation
        nElem = size(MESH.elem, 1);

        %define the global stiffness and mass matrix
        nDofs = sum(MESH.dof);
        KGlobal = sparse(nDofs, nDofs);
        MGlobal = sparse(nDofs, nDofs);

        %define the cell array where the matrix of single EF is stored
        KEFGlobal = {};

        %calculation of the stiffness and mass matrix for each finite element
        for e = 1 : nElem
            %define the parameters of the material
            [materialParameters, sectionParameters] = materialLibrary(MESH.mat(e));
            
            %define the stiffness matrix in the local reference system
            KEFLocal = Kfun(MESH.length(e), materialParameters, sectionParameters);
            MEFLocal = Mfun(MESH.length(e), materialParameters, sectionParameters);
            
            %define the rotation matrix
            T = Tfun(MESH.angle(e));

            %calculate the stiffness matrix of the global reference frame
            KEF = T.' * KEFLocal * T;
            MEF = T.' * MEFLocal * T;
            
            %save the matrix of the single element in a cell array
            KEFGlobal = [KEFGlobal, KEF];

            %assembly the stiffness matrix
            [KGlobal] = FEM.assemble(KGlobal, KEF, matrixGlobalDofs, e);
            [MGlobal] = FEM.assemble(MGlobal, MEF, matrixGlobalDofs, e);
            
            %check the assembly of the matrix in a plot
            %spy(KGlobal)
        end
    end




    %% ************************************************************************************************************
    %  ************       				 addConcentratedLoad		       					           ************
    %  ************************************************************************************************************
    % This function add a nodal load to the global vector of generalized loads.
    % MESH : structured data containing the information of the mesh 
    % fGlobal : matrix defining the dofs that concern a given finite element, i.e. in the k-th row 
    %			there are the dofs that are in the k-th element
    % nodeLoaded : label of the loaded node
    % loadGlobalConc : vector containing the load in the global reference frame [Fh, Fv, M]

    function[fGlobal] = addConcentratedLoad(MESH, fGlobal, nodeLoaded, loadGlobalConc)
        %extract the dofs related to a single node
        vecDofNode = FEM.sumVectorElements(MESH.dof, nodeLoaded - 1) + 1 : FEM.sumVectorElements(MESH.dof, nodeLoaded);
            
        %assign the value in the load vector fGlobal in the right position
        fGlobal(vecDofNode) = fGlobal(vecDofNode) + loadGlobalConc;    
    end




    %% ************************************************************************************************************
    %  ************       				 addDistributedLoad 		       					           ************
    %  ************************************************************************************************************
    % This function add a distributed load on a given element to the global vector of generalized loads.
    % MESH : structured data containing the information of the mesh 
    % matrixGlobalDofs : matrix defining the dofs that concern a given finite element, i.e. in the k-th row 
    %			there are the dofs that are in the k-th element
    % fGlobal : matrix defining the dofs that concern a given finite element, i.e. in the k-th row 
    %			there are the dofs that are in the k-th element
    % nodeLoaded : label of the loaded node
    % ffun : matlab function of the load vector
    % Tfun : matlab function of the rotation matrix
    % elemLoaded : vector containing the label of loaded elements (with the same loadParameters)
    % loadParameters : vector defining the load parameters [p1, p2, q1, q2]

    function[fGlobal] = addDistributedLoad(MESH, matrixGlobalDofs, fGlobal, fFun, Tfun, elemLoaded, loadParameters)
        for k = 1 : length(elemLoaded)
            %compute the load vector of the element in the global reference frame
            fEFLocal = fFun(MESH.length(elemLoaded(k)), loadParameters);

            %define the rotation matrix
            T = Tfun(MESH.angle(elemLoaded(k)));

            %define the load vector of the element in the global reference frame
            fEF = T.' * fEFLocal;

            %assembly the load vector of the element in the global load vector of the entire structure
            [fGlobal] = FEM.assemble(fGlobal, fEF, matrixGlobalDofs, elemLoaded(k));
        end
    end 



    %% ************************************************************************************************************
    %  ************       				     addThermalLoad 		       					           ************
    %  ************************************************************************************************************
    % This function add a thermal load on a given element to the global vector of generalized loads.
    % MESH : structured data containing the information of the mesh 
    % matrixGlobalDofs : matrix defining the dofs that concern a given finite element, i.e. in the k-th row 
    %			there are the dofs that are in the k-th element
    % fGlobal : matrix defining the dofs that concern a given finite element, i.e. in the k-th row 
    %			there are the dofs that are in the k-th element
    % nodeLoaded : label of the loaded node
    % ffun : matlab function of the load vector
    % Tfun : matlab function of the rotation matrix
    % elemLoaded : matlab function of the rotation matrix
    % loadParameters : vector defining the load parameters [p1, p2, q1, q2]

    function[fGlobal] = addThermalLoad(MESH, matrixGlobalDofs, fGlobal, fFun, Tfun, elemLoaded, loadParameters)
        for k = 1 : length(elemLoaded)
            %extract the mechanical parameters
            [materialParameters, sectionParameters] = materialLibrary(MESH.mat(elemLoaded(k)));
            mechanicalParametersTerm = [materialParameters sectionParameters];

            %compute the load vector of the element in the global reference frame
            fEFLocal = fFun(MESH.length(elemLoaded(k)), loadParameters, mechanicalParametersTerm);

            %define the rotation matrix
            T = Tfun(MESH.angle(elemLoaded(k)));

            %define the load vector of the element in the global reference frame
            fEF = T.' * fEFLocal;

            %assembly the load vector of the element in the global load vector of the entire structure
            [fGlobal] = FEM.assemble(fGlobal, fEF, matrixGlobalDofs, elemLoaded(k));
        end
    end 




    %% ************************************************************************************************************
    %  ************       				     defineInternalConstraints    					           ************
    %  ************************************************************************************************************
    % This function defines the matrix of internal constraints A, so that the internal constraints on the global
    % dofs u can be defined as A*u = 0.
    % MESH : structured data containing the information of the mesh 
    % matrixInternalConstraints : matrix containing in each row the dofs that are linked by an internal constraint
    % A : matrix of the internal constraint (A*u = 0)
    
    function[A] = defineInternalConstraints(MESH, matrixInternalConstraints)
        %define the total number of degrees of freedom
        nDofs = sum(MESH.dof);

        %define the matrix A s.t. A*u = 0
        A = zeros(size(matrixInternalConstraints, 1), nDofs);
        for k = 1 : size(matrixInternalConstraints, 1)
            A(k, matrixInternalConstraints(k, 1)) = 1;
            A(k, matrixInternalConstraints(k, 2)) = -1;    
        end
    end




    %% ************************************************************************************************************
    %  ************       				 stiffnessMatrixInternalConstraint 					           ************
    %  ************************************************************************************************************
    % This function allows the transformation between all the global dofs u and the "master" nodes q, i.e. the 
    % nodes that are useful for the solution after the application of internal constraint.
    % A : matrix of the internal constraint (A*u = 0)
    % KGlobal : global stiffness matrix of the entire structure (referring to all dofs u)
    % fGlobal : global load vector (referring to all dofs u)
    % L : matrix linking all the global dofs u and the "master" dofs q, so that u = L * q
    % KGlobalC : global stiffness matrix of the entire structure after the application of the internal constraints
    % (referring to constrained dofs q)
    % fGlobalC : global load vector after the application of the internal constraints
    % (referring to constrained dofs q)
    % nDofsC : dimension of the vector q, the number of "master" dofs

    function[L, KGlobalC, MGlobalC, fGlobalC, nDofsC] = stiffnessMatrixInternalConstraint(A, KGlobal, MGlobal, fGlobal)
        if isempty(A)
            L = sparse(eye(size(KGlobal, 2)));
            KGlobalC = KGlobal;
            MGlobalC = MGlobal;
            fGlobalC = fGlobal;
            nDofsC = size(KGlobal, 2);
        else 
            %compute the matrix that link the global dofs and the master dofs u = L*q form the 
            %constraint equations A*u = 0
            L = null(A, 'r');

            %compute the stiffness matrix and the load vector after the application of the constraints
            KGlobalC = L.' * KGlobal * L;
            MGlobalC = L.' * MGlobal * L;
            fGlobalC = L.' * fGlobal;

            %compute the dimension of the vector q, the number of new dofs
            nDofsC = size(L, 2);
        end
    end
    
    
    
    
    %% ************************************************************************************************************
    %  ************       				 dampingMatrixInternalConstraint 					           ************
    %  ************************************************************************************************************
    % This function allows the transformation between all the global dofs u and the "master" nodes q, i.e. the 
    % nodes that are useful for the solution after the application of internal constraint.
    % A : matrix of the internal constraint (A*u = 0)
    % KGlobal : global stiffness matrix of the entire structure (referring to all dofs u)
    % fGlobal : global load vector (referring to all dofs u)
    % L : matrix linking all the global dofs u and the "master" dofs q, so that u = L * q
    % KGlobalC : global stiffness matrix of the entire structure after the application of the internal constraints
    % (referring to constrained dofs q)
    % fGlobalC : global load vector after the application of the internal constraints
    % (referring to constrained dofs q)
    % nDofsC : dimension of the vector q, the number of "master" dofs

    function[CGlobalC] = dampingMatrixInternalConstraint(A, CGlobal)
        if isempty(A)
            L = sparse(eye(size(CGlobal, 2)));
            CGlobalC = CGlobal;
        else 
            %compute the matrix that link the global dofs and the master dofs u = L*q form the 
            %constraint equations A*u = 0
            L = null(A, 'r');

            %compute the stiffness matrix and the load vector after the application of the constraints
            CGlobalC = L.' * CGlobal * L;
        end
    end
    
    
    
    
    
    %% ************************************************************************************************************
    %  ************       				    defineExternalConstraints   					           ************
    %  ************************************************************************************************************
    % This function allows the creation of the matrix C of the external ground constraints, so that C * U = 0
    % C : matrix of the (ground) external constraint (C*u = 0)
    % MESH : structured data containing the information of the mesh 
    % fixed : vector containing the label of fixed dofs
    function[C] = defineExternalConstraints(MESH, fixed)
        %define the total number of degrees of freedom
        nDofs = sum(MESH.dof);

        %define the matrix C s.t. C*u = g
        C = zeros(size(fixed, 1), nDofs);
        for k = 1 : length(fixed)
            C(k, fixed(k)) = 1; 
        end
    end




    %% ************************************************************************************************************
    %  ************       				            defineUFixed   					                   ************
    %  ************************************************************************************************************
    % This function defines the value of the global vector of degrees of freedom with the given value in the fixed
    % dofs and zero otherwise.
    % MESH : structured data containing the information of the mesh 
    % fixed : vector containing the label of fixed dofs
    % valueFixed : value of the fixed dofs (same size of valueFixed)
    % uFixed : vector of degrees of freedom

    function[uFixed] = defineUFixed(MESH, fixed, valueFixed)
        %define the vector uFixed
        uFixed = sparse(zeros(sum(MESH.dof), 1));
        for k = 1 : length(fixed)
            uFixed(fixed(k)) = valueFixed(k);
        end
    end




    %% ************************************************************************************************************
    %  ************       				     applyExternalConstraints   					           ************
    %  ************************************************************************************************************
    % This function treats the external constraints, it defines the values fixedConst and uFixedConst, the value 
    % of fixed and uFixed when the internal restraints have already been applied (fixed and uFixed refers to u,
    % fixedConst and uFixedConst refers to q).
    % C : matrix of the (ground) external constraint (C*u = 0)
    % L : matrix linking all the global dofs u and the "master" dofs q, so that u = L * q
    % valueFixed : value of the fixed dofs (same size of valueFixed)
    % fixedConst : vector containing the label of fixed dofs (it refers to q)
    % uFixedConst : vector of global "master" degrees of freedom (it refers to q)
    % NOTE : the order is fundamental! Cconst is checked by rows, in each row of Cconst the not-null element is 
    % the constrained degree of freedom q (the vector of unknowns obtained applying the internal release to u). 
    % The vector uFixedConst is computed using the valueFixed node vector (it is the vector g) and is created 
    % spanning Cconst by row and adding the valueFixed element in the right position of uFixedConst

    function[Cconst, fixedConst, uFixedConst] = applyExternalConstraints(C, L, valueFixed)
        %the matrix C must be corrected to consider the internal constrain, 
        %so C*u = g , u = L*q  -->> C*L*q = Cconst*q = g (g is valueFixed)
        Cconst = C * L;

        %define the new fixed dofs and the new uFixedConst
        fixedConst = [];
        uFixedConst = sparse(size(L, 2), 1);
        for k = 1 : size(Cconst, 1)
            %compute the new fixed degree of freedom, considering the displacement vector q
            newFixedConst = find(Cconst(k, :) ~= 0);

            %add the value to the list of fixed dofs (with respect to q)
            fixedConst = [fixedConst newFixedConst];

            %define the new vector of uFixedConst with the value of fixed dofs (with respect to q),
            %(note that k range between 1 and the dimension of valueFixed)
            uFixedConst(newFixedConst) = valueFixed(k);
        end
    end




    %% ************************************************************************************************************
    %  ************       			  		 	        solveStatic	 	    	  	           	       ************
    %  ************************************************************************************************************
    % This function solves the linear static problem, calculating the nodal displacement and reactions
    % KGlobal : global stiffness matrix obtained by assembly the stiffness matrices of the single Finite Elements,
    % the size Ndof x Ndof, where Ndof is the total number of degrees of freedom of the structure
    % fGlobal : vector of (known) external loads, note that the length is equal to Ndof, where Ndof is the total 
    %  number of degrees of freedom of the structure, so where the load is unknown the element of f is zero
    % fixed : vector containing the position (the index) of fixed nodes, so the vector U(fixed) contains all the
    %  nodes whose displacements are known
    % valueFixed : value of the fixed dofs (same size of valueFixed)
    % A : matrix of the internal constraint (A*u = 0)
    % C : matrix of the (ground) external constraint (C*u = 0)
    % U : global vector of nodal displacement, the size is Ndof (it contains both the displacements of fixed and 
    % free nodes)
    % K : matrix KGlobal(free, free), used to solve the displacement problem and find the unknown value of U
    % F : vector fGlobal(free), used to solve the displacement problem and find the unknown value of U
    % R : global vector of reactions, the size is Ndof (it contains the reaction of fixed nodes and the value is
    % zero for free nodes)
    % Rc : global vector of reactions, considering the constraints (it contains in position)

    function [U, K, F, R, Rconst] = solveStatic(KGlobal, fGlobal, fixed, valueFixed, A, C)
        %apply the internal constraint and redefine fixed and uFixed considering the application 
        %of internal constraints (the mass matrix MGlobal is assigned equal to KGlobal, since it is not used)
        [L, KGlobalC, ~, fGlobalC, nDofsC] = FEM.stiffnessMatrixInternalConstraint(A, KGlobal, KGlobal, fGlobal);
        [Cconst, fixedConst, uFixedConst] = FEM.applyExternalConstraints(C, L, valueFixed);

        %sorting the value of fixed, to have an ordered matrix
        fixedConst = sort(fixedConst);
        
        %defining the vector of equivalent forces, due to the imposed nodal displacement
        %(this formulation is correct, because UFixed has zero value in the position UFixed(free), so the 
        % following multiplication gives zero in fD(free) 
        fdispl = - (KGlobalC(:, fixedConst) * uFixedConst(fixedConst));
        fGlobalC = fGlobalC + fdispl; 
        
        %defining the freeConst dofs, taking all the label of the dofs that are not in the vector fixedConst
        freeConst = setdiff(1 : 1 : sum(nDofsC), fixedConst);

        %defining the matrix and vector for global analysis
        K = sparse(KGlobalC(freeConst, freeConst));
        F = sparse(fGlobalC(freeConst, :));
        
        %solution of the system (with respect to q)
        uFreeDofsConst = K \ F; 
        
        %creating the global vector of the solution (with respect to q)
        for j = 1 : size(uFreeDofsConst, 2)
            q(:, j) = uFixedConst;
            q(freeConst, j) = uFreeDofsConst(:, j);
        end
        
        %now the solution must be given in the original unknown u, so the nodal displacement 
        %(not constrained by internal constraints)
        U = L * q;

        %calculation of the reactions, the forces due to imposed displacement (this formulation is correct
        %because it considers that in the position of the reaction, in f the value is zero, this is the definition)
        %of the load vector
        R = sparse(zeros(size(U)));
        R(fixed) = KGlobal(fixed, :) * U;

        Rt = L.' * KGlobal * U;
        Rconst(fixedConst) = Rt(fixedConst);
    end



    %% ************************************************************************************************************
    %  ************       			  		 			solveModal	 	    	  	           		   ************
    %  ************************************************************************************************************
    %This function solves the linear eigenvalue generalized problem (-omega^2 M - K)u = 0, calculating
    % the eigenfrequencies and principal modes of the structure.
    %MGlobal : global mass matrix obtained by assembly the mass matrices of the single Finite Element,
    %	  the size Ndof x Ndof, where Ndof is the total number of degrees of freedom of the structure
    %KGlobal : global stiffness matrix obtained by assembly the stiffness matrices of the single Finite Element,
    %	  the size Ndof x Ndof, where Ndof is the total number of degrees of freedom of the structure
    %valueFixed : vector of the known nodal displacement, so it is an expansion of the vector U(fixed) into a vector
    %		of length Ndof (so where the nodal displacement is known, the value of UFixed is the proper element
    %		of U(fixed), where the nodal displacement is unknown, the value is set to zero)
    % A : matrix of the internal constraint (A*u = 0)
    % C : matrix of the (ground) external constraint (C*u = 0)
    % numEig = number of eigenvalues that have to be found
    % sigmaEig = parameter of the function eigs for solving eigenvalues problem 
    % 		('sm' = find smallest eignevalues / 'lm' = find the greatest eigenvalues / value : the eigenvalues 
    %		near the scalar value )
    %psi : the matrix of eigenmodes, the output of function eigs (the dimension of the matrix is Nfree X Nfree)
    %lambda : diagonal matrix of the eigenvalues (the dimension of the matrix is Nfree X Nfree)
    %omega : vector of the natural frequency in the order given by the function eigs
    %U : global vector of eigenmodes of the structure, the size is Ndof (it contains both the displacements 
    %		of fixed and free nodes)		
    %naturalFrequency : vector of the natural frequency reorder so that the smallest is in the first position

    function [psi, lambda, omega, U, naturalFrequency] = solveModal(MGlobal, KGlobal, valueFixed, A, C, sigmaEig, numEig)
        %apply the internal constraint and redefine fixed and uFixed considering the application 
        %of internal constraints
        [L, KGlobalC, MGlobalC, ~, nDofsC] = FEM.stiffnessMatrixInternalConstraint(A, KGlobal, MGlobal, ones(size(KGlobal, 1), 1));
        [Cconst, fixedConst, uFixedConst] = FEM.applyExternalConstraints(C, L, valueFixed);

        %sorting the value of fixed, to have an ordered matrix
        fixedConst = sort(fixedConst);
        
        %defining the freeConst dofs, taking all the label of the dofs that are not in the vector fixedConst
        freeConst = setdiff(1 : 1 : sum(nDofsC), fixedConst);

        %defining the matrix and vector for global analysis
        K = sparse(KGlobalC(freeConst, freeConst));
        M = sparse(MGlobalC(freeConst, freeConst));   

        %change the number of the eigenvalue if it exceed the size of the matrix
        if numEig > size(K, 1)
            numEig = size(K, 1);
        end
        
        %calculate the eigenvalues and eigenvectors of the constrained mass and stiffness matrix
        [psi, lambda] = eigs(K, M, numEig, sigmaEig);
        omega = sqrt(diag(lambda));

        %find the correct order of the modes and define the naural frequency
        [naturalFrequency, index] = sort(omega);

        %creating the global vector of the solution (with respect to q)
        for j = 1 : size(psi, 2)
            q(:, j) = uFixedConst;
            q(freeConst, j) = psi(:, j);
        end
        
        %now the solution must be given in the original unknown u, so the nodal displacement 
        %(not constrained by internal constraints)
        U = L * q;
        
        %sort the column, so that the order is the same as naturalFrequency
        U = U(:, index);
    end





    %% ************************************************************************************************************
    %  ************       			  		 		defineLoadTimePath			 	    	   		   ************
    %  ************************************************************************************************************
    % This function define the load history as a fuction of time. The symbolic function containing the time
    % variable t is transformed into a symbolic function and then the value of the load in discrete time
    % interval is calculated.
    % symTimeFun : symbolic function of time that is multiplied by the vector of generalized loads F
    % F : vector of generalized load, the discretization of the load in space is done as for the static case
    % deltaT : vector containing the time interval discretization of temporal domain
    % t0 : initial time 
    % Ft = matrix containing in each column the vector of generalized loads at a given time

    function[Ft] = defineLoadTimePath(symTimeFun, F, deltaT, t0)
        %define symbolic function t
        syms t

        %transform the symbolic function in anonymous function
        load = matlabFunction(symTimeFun, 'Vars', {t});

        %define the starting time
        time = t0;

        %define the load at initial time step
        Ft = F * load(time);

        %define the loads at each time step
        for n = 1 : length(deltaT)
            %update the time step
            time = time + deltaT(n);

            %update the matrix of the load
            Ft = [Ft F * load(time)];
        end
    end		


    %% ************************************************************************************************************
    %  ************       			  		 	defineParametersGenAlpha		 	    	   		   ************
    %  ************************************************************************************************************
    % This function defines all the parameters of the Generalized Aplha Method as function of only parameter, the
    % spectral radius at infinity, defining the numerical damping of high frequencies.
    % alphaM, alphaF, theta1, theta2 : parameters of the Generalized alpha methods
    % rhoInfty : spectral radius at infinity, whose value is in the interval [0, 1] 

    function[alphaM, alphaF, theta1, theta2] = defineParametersGenAlpha(rhoInfty)
        alphaF = rhoInfty / (rhoInfty + 1);
        alphaM = (2 * rhoInfty - 1) / (rhoInfty + 1);
        theta1 = 0.5 - alphaM + alphaF;
        theta2 = 0.5 * (1 - alphaM + alphaF)^2;
    end






    %% ************************************************************************************************************
    %  ************       			  		 	timeIntegratorGeneralizedAlpha	 	    	   		   ************
    %  ************************************************************************************************************
    % This function solves the transient problem using the Generalized Alpha Method as time integrator. The
    % Finite Element Method is used for the spatial approximation and for time approximation the Generalized Alpha 
    % Method is used.
    % Kg : global stiffness matrix obtained by assembly the stiffness matrices of the single Finite Elements,
    % 	  the size Ndof x Ndof, where Ndof is the total number of degrees of freedom of the structure
    % Cg : global damping matrix 
    % Mg : global mass matrix obtained by assembly the mass matrices of the single Finite Element,
    %	  the size Ndof x Ndof, where Ndof is the total number of degrees of freedom of the structure
    % Fg : matrix containing in each column the vector of generalized loads at a given time
    % fixed : vector containing the position (the index) of fixed nodes, so the vector U(fixed) contains all the
    % 		nodes whose displacements are known
    % dof : vector of the number of the dofs of the structure 
    % deltaT : vector containing the time interval discretization of temporal domain
    % U0, dUt0 : initial condition for generalized displacements and velocity
    % alphaM, alphaF, theta1, theta2 : parameters of the Generalized alpha methods
    % U : matrix containing in each column the global vector of nodal displacement at each time step, 
    %      the size is Ndof (it contains both the displacements of fixed and free nodes)
    % dUt : matrix containing in each column the global vector of nodal velocity at each time step, 
    %      the size is Ndof (it contains both the displacements of fixed and free nodes)
    % dUtt : matrix containing in each column the global vector of nodal acceleration at each time step, 
    %      the size is Ndof (it contains both the displacements of fixed and free nodes)

    function[U, dUt, dUtt] = timeIntegratorGeneralizedAlpha(KGlobal, CGlobal, MGlobal, fGlobal, valueFixed, A, C, deltaT, ...
                                                                    U0, dUt0, alphaM, alphaF, theta1, theta2)
        
                                                                
        %apply the internal constraint and redefine fixed and uFixed considering the application 
        %of internal constraints (the mass matrix MGlobal is assigned equal to KGlobal, since it is not used)
        [L, KGlobalC, MGlobalC, fGlobalC, nDofsC] = FEM.stiffnessMatrixInternalConstraint(A, KGlobal, MGlobal, fGlobal);
        [CGlobalC] = FEM.dampingMatrixInternalConstraint(A, CGlobal);
        [Cconst, fixedConst, uFixedConst] = FEM.applyExternalConstraints(C, L, valueFixed);
        
        %sorting the value of fixed, to have an ordered matrix
        fixedConst = sort(fixedConst);
        
        %defining the freeConst dofs, taking all the label of the dofs that are not in the vector fixedConst
        freeConst = setdiff(1 : 1 : sum(nDofsC), fixedConst);

        %defining the matrix and vector for global analysis
        K = sparse(KGlobalC(freeConst, freeConst));
        C = sparse(CGlobalC(freeConst, freeConst));
        M = sparse(MGlobalC(freeConst, freeConst));
        F = sparse(fGlobalC(freeConst, :));
                              
        %redefine the initial condition, as the condition only for free dofs
        U0 = U0(freeConst);
        dUt0 = dUt0(freeConst);

        %initialize the displacement, velocity and acceleration as sparse matrices, the number of the columns
        %is equal to the number of load steps stored in F
        UU = sparse(length(freeConst), size(F, 2));
        dUUt = sparse(length(freeConst), size(F, 2));
        dUUtt = sparse(length(freeConst), size(F, 2));

        %assign the initial condition as the solution in the first time interval; the initial acceleration is
        %calculated solving the equation of motion at t = 0
        dUtt0 = M \ (F(:, 1) - (K * U0) - (C * dUt0));
        UU(:, 1) = U0;
        dUUt(:, 1) = dUt0;
        dUUtt(:, 1) = dUtt0;

        %print on the prompt the iteration
        lenPrint = 0;

        %calculation of the stiffness and mass matrix for each finite element
        fprintf('\t\t~ calculation of solution at time interval... ');

        %generalized-alpha algorithm for time analysis
        for n = 1 : (length(deltaT))

            %print the step in the iteration on the prompt
            lenPrint = FEM.printStepPlus(n, length(deltaT), lenPrint);	

            %definition of the predictors for displacement and velocity
            UUU = UU(:, n) + deltaT(n) * dUUt(:, n) + ((1 / 2) * deltaT(n) ^ 2) * (1 - theta2) * dUUtt(:, n);
            dUUU = dUUt(:, n) + deltaT(n) * (1 - theta1) * dUUtt(:, n);

            %define the vector that is used to find the acceleration at time n + 1
            BB = F(:, n + 1) - (C * dUUU * (1 - alphaF)) - (K * UUU * (1 - alphaF)) - (K * alphaF * UU(:, n))...
                                            - (C * alphaF * dUUt(:, n)) - (M * alphaM * dUUtt(:, n));

            %define the matrix that is used to find the acceleration at time n + 1												
            AA = M * (1 - alphaM) + (C * (1 - alphaF) * deltaT(n) * theta1) +...
                ((1 / 2) * K * (1 - alphaF) * (deltaT(n)^2) * theta2);

            %calculation of the acceleration at the time step t + 1
            dUUtt(:,n + 1) = AA \ BB;

            %calculation of the corrector and of the displacement and velocity at time step t + 1
            UU(:,n + 1) = UUU + (1 / 2) * (deltaT(n) ^ 2) * theta2 * dUUtt(:, n + 1);
            dUUt(:,n + 1) = dUUU + deltaT(n) * theta1 * dUUtt(:, n + 1);
        end
        fprintf('\n');

        %creating the global vector of the solution (with respect to q)
        %NB: it is supposed that in the fixed node the value of the velocity is zero
        for j = 1 : size(UU, 2)
            q(:, j) = uFixedConst;
            q(freeConst, j) = UU(:, j);
            
            qt(:, j) = sparse(size(q, 1), 1);
            qt(freeConst, j) = dUUt(:, j);
            
            qtt(:, j) = sparse(size(q, 1), 1);
            qtt(freeConst, j) = dUUtt(:, j);
        end
        
        %now the solution must be given in the original unknown u, so the nodal displacement 
        %(not constrained by internal constraints)
        U = L * q;
        dUt = L * qt;
        dUtt = L * qtt;
    end
        
        
        
        
        
    %% ************************************************************************************************************
    %  ************       			  		 	    computeNodalForces	 	    	  	           	   ************
    %  ************************************************************************************************************
    % This function solves the linear static problem, calculating the nodal displacement and reactions
    % MESH : structured data containing the information of the mesh 
    % U : global vector of nodal displacement, the size is Ndof (it contains both the displacements of fixed and 
    % free nodes)
    % matrixGlobalDofs : matrix defining the dofs that concern a given finite element, i.e. in the k-th row 
    % there are the dofs that are in the k-th element
    % KEFGlobal : cell array containing in the k-th position the stiffness matrix of the k-th element in local 
    % reference frame
    % TPostProc : matlab function containing the rotation matrix for the post-processing defined by the function
    % FEM.defineRotationMatrixPostProcessing(codeEF)
    % nodalForcesGlobal : matrix containing in the k-th column the forces at the nodes of element k-th in global
    % reference system in particular each column is [Fh1, Fv1, W1, Fh2, Fv2, W2].'
    % nodalForcesLocal : matrix containing in the k-th column the forces at the nodes of element k-th in local
    % reference system rotated as the element, in particular each column is [N1, T1, M1, N2, T2, M2].'

    % NOTE : the quantities in the vector nodalForcesLocal are not the internal forces of the element, but it
    % is simply the forces in a reference frame local to the element, obtained by the rotation of the global
    % system of the angle of the considered finite element!

    function[nodalForcesGlobal, nodalForcesLocal] = computeNodalForces(MESH, U, matrixGlobalDofs, ...
                                                                                         KEFGlobal, TPostProc)
        %initialize the outputs
        nodalForcesGlobal = [];
        nodalForcesLocal = [];
        for e = 1 : size(MESH.elem, 1)
            %define the vector of displacements for the considered element
            UEFGlobal = U(matrixGlobalDofs(e, :));
            
            %define the of generalized loads in the global reference frame
            fEFGlobal = KEFGlobal{e} * UEFGlobal;

            %define the rotation matrix for the post processing
            T = TPostProc(MESH.angle(e));
            
            %define the of generalized loads in the local reference frame
            fEFLocal = T * fEFGlobal;
            
            %update the outputs
            nodalForcesGlobal = [nodalForcesGlobal fEFGlobal];
            nodalForcesLocal = [nodalForcesLocal fEFLocal];
        end
    end




    %% ************************************************************************************************************
    %  ************				  		        defineSolutionField          	 		               ************
    %  ************************************************************************************************************
    % This function defines the matlab function (function obtained from a symbolic function, but faster for 
    % numerical computation) of the solution fields, given the formulation of the finite element. 
    % These function are the solutions of axial and transversal displacements (U and V) and internal force (N, M, T)
    % as function of the local coordinates on the element.
    % codeEF : string defining the formulation of the finite element ("TRUSS" for truss element,"EULER_BERNOULLI" for 
    % complete extensible Euler-Bernoulli beam element)
    % solU, solV, solN, solT, solM : matlab function of the solution fields
    
    % NOTE : in order to have a unique function for both truss and beam element, the plot of truss element is more
    % complicated, it is necessary to pass to a beam element with suitable rotation (rigid-body rotation) in the
    % local (!!) reference frame.

    function[solU, solV, solN, solT, solM] = defineSolutionField(codeEF)
        %define the symbols used in the matrix formulation
        syms x E A J L angleRot rho
            
        %shape function for truss and EB beam
        Ntruss = [  1 - (x / L); 
                    (x / L)];
        Nbeam =[    1 - 3 * (x / L)^2 + 2 * (x / L)^3;
                    x - 2 * ((x^2) / L) + ((x^3) / (L^2));
                    3 * (x / L)^2 - 2 * (x / L)^3;
                    -((x^2) / L) + ((x^3) / (L^2));
        ]; 

        %definition of the vectors of parameters
        materialParameters = [E, rho];
        sectionParameters = [A, J];

        %define the symbolic matrix
        switch codeEF
            case 'TRUSS'
                %rotation matrix for a truss element from local to global reference frame 
                %(local dofs: [u1, u2] --> global dofs: [u1, v1, u2, v2])
                TsTruss = [	cos(angleRot),     sin(angleRot),         0,              0; 
                        0,			    0,					cos(angleRot), 	sin(angleRot)];
                
                %rotation matrix for an extensible Euler-Bernoulli element from local to global reference frame 
                %(local dofs: [u1, v1, theta1, u2, v2, theta2] --> global dofs: [u1, v1, theta1, u2, v2, theta2])
                Ts = [	cos(angleRot), 	sin(angleRot),	0,			0,			    0,			    0;
                        sin(angleRot), -cos(angleRot), 	0,			0,			    0,			    0;
                        0,		 	 	0,			    -1,			0,			    0,			    0;
                        0,				0,			    0,			cos(angleRot), 	sin(angleRot), 	0;
                        0,				0,			    0,			sin(angleRot),  -cos(angleRot), 0;
                        0,				0,			    0,			0,		 	    0,			    -1];

                syms u1 u2 u3 u4
                Uglobal = [u1; u2; u3; u4];

                %field solution
                solUs = 0*x + Ntruss.' * TsTruss * Uglobal;
                solNs = 0*x + E * A * diff(Ntruss.', x) * TsTruss * Uglobal;
                solTs = 0 * x;
                solMs = 0 * x;

                
                syms r1 r2
                UglobalBeam = [u1; u2; r1; u3; u4; r2];
                ULoc = Ts * UglobalBeam;
                ULocalBeam = [ULoc(1); ULoc(2); ((ULoc(5) - ULoc(2)) / L); ULoc(4);  ULoc(5);((ULoc(5) - ULoc(2)) / L)];
                
                %define the vector of shape function regarding a beam element
                N = [0, Nbeam(1), Nbeam(2), 0, Nbeam(3), Nbeam(4)];
                
                %field solution (transversal displacement)
                solVs = 0 * x + N * ULocalBeam;


            case 'EULER_BERNOULLI'
                %rotation matrix for an extensible Euler-Bernoulli element from local to global reference frame 
                %(local dofs: [u1, v1, theta1, u2, v2, theta2] --> global dofs: [u1, v1, theta1, u2, v2, theta2])
                Ts = [	cos(angleRot), 	sin(angleRot),	0,			0,			    0,			    0;
                        sin(angleRot), -cos(angleRot), 	0,			0,			    0,			    0;
                        0,		 	 	0,			    -1,			0,			    0,			    0;
                        0,				0,			    0,			cos(angleRot), 	sin(angleRot), 	0;
                        0,				0,			    0,			sin(angleRot),  -cos(angleRot), 0;
                        0,				0,			    0,			0,		 	    0,			    -1];

                syms u1 u2 u3 u4 u5 u6
                Uglobal = [u1; u2; u3; u4; u5; u6];
                
                %define the matrix of shape function
                N = [Ntruss(1),     0,          0,          Ntruss(2),      0,          0;
                    0,             Nbeam(1),   Nbeam(2),   0,              Nbeam(3),   Nbeam(4);];

                %field solution
                solUs = 0*x + N(1, :) * Ts * Uglobal;
                solVs = 0*x + N(2, :) * Ts * Uglobal;
                solNs = 0*x + E * A * diff(N(1, :), x) * Ts * Uglobal;
                solTs = 0*x - E * J * diff(diff(diff(N(2, :), x), x)) * Ts * Uglobal;
                solMs = 0*x - E * J * diff(diff(N(2, :), x), x) * Ts * Uglobal;
        end

        %define the matlab function of the stiffness matrix 
        solU = matlabFunction(solUs, 'Vars', {x, L, angleRot, Uglobal, materialParameters, sectionParameters});
        solV = matlabFunction(solVs, 'Vars', {x, L, angleRot, Uglobal, materialParameters, sectionParameters});
        solN = matlabFunction(solNs, 'Vars', {x, L, angleRot, Uglobal, materialParameters, sectionParameters});
        solT = matlabFunction(solTs, 'Vars', {x, L, angleRot, Uglobal, materialParameters, sectionParameters});
        solM = matlabFunction(solMs, 'Vars', {x, L, angleRot, Uglobal, materialParameters, sectionParameters});
    end

    


    %% ************************************************************************************************************
    %  ************		  		            defineSolutionFieldThermAxial       	 		           ************
    %  ************************************************************************************************************
    % This function defines the matlab function (function obtained from a symbolic function, but faster for 
    % numerical computation) of the solution fields, given the formulation of the finite element. 
    % These function are the solutions of axial and transversal displacements (U and V) and internal force (N, M, T)
    % as function of the local coordinates on the element.
    % codeEF : string defining the formulation of the finite element ("TRUSS" for truss element,"EULER_BERNOULLI" for 
    % complete extensible Euler-Bernoulli beam element)
    % solU, solV, solN, solT, solM : matlab function of the solution fields
    
    function[solN, solT, solM] = defineSolutionFieldThermAxial(codeEF)
        %define the symbols used in the matrix formulation
        syms x E A J L angleRot alphaT DT h rho
            
        %definition of the vectors of parameters
        materialParameters = [E, rho];
        sectionParameters = [A, J];
        loadParametersTerm = [alphaT, h, DT];

        %define the symbolic matrix
        switch codeEF
            case 'TRUSS'
                %field solution
                solNs = 0*x + (E * A * alphaT * DT);
                solTs = 0 * x;
                solMs = 0 * x;

            case 'EULER_BERNOULLI'
                %field solution
                solNs = 0*x + (E * A * alphaT * DT);
                solTs = 0*x;
                solMs = 0*x;
        end

        %define the matlab function of the stiffness matrix
        solN = matlabFunction(solNs, 'Vars', {x, materialParameters, sectionParameters, loadParametersTerm});
        solT = matlabFunction(solTs, 'Vars', {x, materialParameters, sectionParameters, loadParametersTerm});
        solM = matlabFunction(solMs, 'Vars', {x, materialParameters, sectionParameters, loadParametersTerm});
    end




    %% ************************************************************************************************************
    %  ************		  		            defineSolutionFieldThermCurv           	 		           ************
    %  ************************************************************************************************************
    % This function defines the matlab function (function obtained from a symbolic function, but faster for 
    % numerical computation) of the solution fields, given the formulation of the finite element. 
    % These function are the solutions of axial and transversal displacements (U and V) and internal force (N, M, T)
    % as function of the local coordinates on the element.
    % codeEF : string defining the formulation of the finite element ("TRUSS" for truss element,"EULER_BERNOULLI" for 
    % complete extensible Euler-Bernoulli beam element)
    % solU, solV, solN, solT, solM : matlab function of the solution fields

    % NOTE : the solution field have all the same input, so the calculateValueSolution can be called in the same 
    % way to define the values of the solution fields

    function[solN, solT, solM] = defineSolutionFieldThermCurv(codeEF)
        %define the symbols used in the matrix formulation
        syms x E A J L angleRot alphaT DT h rho
            
        %definition of the vectors of parameters
        materialParameters = [E, rho];
        sectionParameters = [A, J];
        loadParametersTerm = [alphaT, h, DT];

        %define the symbolic matrix
        switch codeEF
            case 'TRUSS'
                %field solution
                solNs = 0 * x;
                solTs = 0 * x;
                solMs = 0 * x;

            case 'EULER_BERNOULLI'
                %field solution
                solNs = 0 * x;
                solTs = 0 * x;
                solMs = 0 * x + (2 * alphaT * DT * E * J / h);
        end

        %define the matlab function of the stiffness matrix 
        solN = matlabFunction(solNs, 'Vars', {x, materialParameters, sectionParameters, loadParametersTerm});
        solT = matlabFunction(solTs, 'Vars', {x, materialParameters, sectionParameters, loadParametersTerm});
        solM = matlabFunction(solMs, 'Vars', {x, materialParameters, sectionParameters, loadParametersTerm});
    end



    %% ************************************************************************************************************
    %  ************				  		   calculateAxialForceFunctionEB          	 		           ************
    %  ************************************************************************************************************
    % This function calculate the field of axial force using equilibrium, when the formulation with inextensible
    % beams are used. The solution for inextensible Euler-Bernoulli beam cannot catch the function of axial force, 
    % that must be computed using equilibrium at the node.

    % coordXEF : coordinates of the points in which the function must be calculated (local coordinate x on the beam)
    % nodalForcesLocal : matrix containing in the k-th column the forces at the nodes of element k-th in local
    % reference system in particular each column is [N1, T1, M1, N2, T2, M2].'
    % distributedLoadInfo : structured datum that contains info on the distributed load on the elements; 
    %       distributedLoadInfo.elemLoaded : vector containing the label of the loaded element
    %       distributedLoadInfo.loadParameters  : cell array containing the parameters of the load [p1, p2, q1, q2]
    % elem : label of the considered element

    function[valuesEFN] = calculateAxialForceFunctionEB(coordXEF, nodalForcesLocal, distributedLoadInfo, elem)
        %retrieve the axial force in node 1 (change the sign for convention)
        N1 = - nodalForcesLocal(1, elem);

        %check whether the considered element is loaded by an axial distributed load and retrieve 
        %the value of the load
        p1 = 0;
        p2 = 0;
        pos = find(distributedLoadInfo.elemLoaded == elem);
        if length(pos) > 0
            loadDistr = distributedLoadInfo.loadParameters{pos};
            p1 = loadDistr(1);
            p2 = loadDistr(2);
        end

        %compute the values of the internal axial force
        valuesEFN = N1 + (p1 * coordXEF) + (((p2 - p1) / (2 * L)) * (coordXEF.^2));
        
    end




    %% ************************************************************************************************************
    %  ************				  		    calculateValueSolution          	 		               ************
    %  ************************************************************************************************************
    % This function calculate the numerical value of the solution field at given local coordinate x on each element
    % MESH : structured data containing the information of the mesh 
    % U : global vector of nodal displacement, the size is Ndof (it contains both the displacements of fixed and 
    % free nodes)
    % matrixGlobalDofs : matrix defining the dofs that concern a given finite element, i.e. in the k-th row 
    % nSegment : number of segment in which each element is divided for the computation of the solution
    % nodalForcesLocal : matrix containing in the k-th column the forces at the nodes of element k-th in local
    % distributedLoadInfo : structured datum that contains info on the distributed load on the elements; 
    %       distributedLoadInfo.elemLoaded : vector containing the label of the loaded element
    %       distributedLoadInfo.loadParameters  : cell array containing the parameters of the load [p1, p2, q1, q2] 
    % solU, solV, solN, solT, solM : matlab function of the solution fields
    % coordX : cell array in which the k-th element is the vector of the local coordinate on the finite element
    % valuesU, valuesV, valuesN, valuesT, valuesM : cell array in which the k-th element is the vector of the
    % solution field computed in the given coordinate

    function[coordX, valuesU, valuesV, valuesN, valuesT, valuesM] = calculateValueSolution(MESH, U,...
                                                            matrixGlobalDofs, nSegment, nodalForcesLocal, ...
                                                                distributedLoadInfo, solU, solV, solN, solT, solM)
        %initialize the cell arrays containing the values of the fields
        coordX = {};
        valuesU = {};
        valuesV = {};
        valuesN = {};
        valuesT = {};
        valuesM = {};
        for e = 1 : size(MESH.elem, 1)
            %define the vector of local coordinates in which the unknown field must be calculated
            coordXEF = linspace(0, MESH.length(e), nSegment + 1);

            %define the parameters of the material
            [materialParameters, sectionParameters] = materialLibrary(MESH.mat(e));

            %define the vector of displacements for the considered element
            UEFGlobal = U(matrixGlobalDofs(e, :));

            %values of the fields on the finite element
            valuesEFU = [];
            valuesEFV = [];
            valuesEFT = [];
            valuesEFM = [];
            valuesEFN = [];
            for k = 1 : length(coordXEF)
                valuesEFU = [valuesEFU, solU(coordXEF(k), MESH.length(e), MESH.angle(e), UEFGlobal, materialParameters, sectionParameters)];
                valuesEFV = [valuesEFV, solV(coordXEF(k), MESH.length(e), MESH.angle(e), UEFGlobal, materialParameters, sectionParameters)];
                valuesEFN = [valuesEFN, solN(coordXEF(k), MESH.length(e), MESH.angle(e), UEFGlobal, materialParameters, sectionParameters)];
                valuesEFT = [valuesEFT, solT(coordXEF(k), MESH.length(e), MESH.angle(e), UEFGlobal, materialParameters, sectionParameters)];
                valuesEFM = [valuesEFM, solM(coordXEF(k), MESH.length(e), MESH.angle(e), UEFGlobal, materialParameters, sectionParameters)];
            end

            %update the outputs
            coordX = [coordX, full(coordXEF)];
            valuesU = [valuesU, full(valuesEFU)];
            valuesV = [valuesV, full(valuesEFV)];
            valuesN = [valuesN, full(valuesEFN)];
            valuesT = [valuesT, full(valuesEFT)];
            valuesM = [valuesM, full(valuesEFM)];
        end
    end




    %% ************************************************************************************************************
    %  ************				  		    calculateValueSolutionTherm          	 		           ************
    %  ************************************************************************************************************
    % This function calculate the numerical value of the solution field at given local coordinate x on each element
    % MESH : structured data containing the information of the mesh 
    % U : global vector of nodal displacement, the size is Ndof (it contains both the displacements of fixed and 
    % free nodes)
    % matrixGlobalDofs : matrix defining the dofs that concern a given finite element, i.e. in the k-th row 
    % nSegment : number of segment in which each element is divided for the computation of the solution
    % nodalForcesLocal : matrix containing in the k-th column the forces at the nodes of element k-th in local
    % distributedLoadInfo : structured datum that contains info on the distributed load on the elements; 
    %       distributedLoadInfo.elemLoaded : vector containing the label of the loaded element
    %       distributedLoadInfo.loadParameters  : cell array containing the parameters of the load [p1, p2, q1, q2] 
    % solU, solV, solN, solT, solM : matlab function of the solution fields
    % coordX : cell array in which the k-th element is the vector of the local coordinate on the finite element
    % valuesU, valuesV, valuesN, valuesT, valuesM : cell array in which the k-th element is the vector of the
    % solution field computed in the given coordinate

    function[valuesN, valuesT, valuesM] = calculateValueSolutionTherm(MESH, nSegment, thermalLoadInfo,...
                                                                                                solN, solT, solM)
        %initialize the cell arrays containing the values of the fields
        valuesN = {};
        valuesT = {};
        valuesM = {};
        for e = 1 : size(MESH.elem, 1)
            %define the vector of local coordinates in which the unknown field must be calculated
            coordXEF = linspace(0, MESH.length(e), nSegment + 1);

            %define the parameters of the material
            [materialParameters, sectionParameters] = materialLibrary(MESH.mat(e));

            %find the position of the current element in the structure containing infos on thermal loads
            positionElemTherm = find(thermalLoadInfo.elemLoaded == e);

            %define the parameters for the thermal load
            if isempty(positionElemTherm)
                %no thermal load is applied to the element, the vector of parameter is set to a zero vector with h = 1 to 
                %avoid division by zero
                loadParametersTerm = [0, 1, 0];
            else
                %the correct set of parameter is considered for thermal load
                loadParametersTerm = thermalLoadInfo.loadParameters{positionElemTherm};
            end

            %values of the fields on the finite element
            valuesEFT = [];
            valuesEFM = [];
            valuesEFN = [];
            for k = 1 : length(coordXEF)
                valuesEFN = [valuesEFN, solN(coordXEF(k), materialParameters, sectionParameters, loadParametersTerm)];
                valuesEFT = [valuesEFT, solT(coordXEF(k), materialParameters, sectionParameters, loadParametersTerm)];
                valuesEFM = [valuesEFM, solM(coordXEF(k), materialParameters, sectionParameters, loadParametersTerm)];
            end

            %update the outputs
            valuesN = [valuesN, full(valuesEFN)];
            valuesT = [valuesT, full(valuesEFT)];
            valuesM = [valuesM, full(valuesEFM)];
        end
    end




    %% ************************************************************************************************************
    %  ************				  		            sumCellArraySolution          	 		           ************
    %  ************************************************************************************************************
    % This function sums the vector contained in each element of a cell array. It must be used only when the 
    % vectors have the same size, it is useful to sum the solution of internal forces when also thermal loads are
    % present in the structure, considering that each element of the cell array is the solution computed in a 
    % single finite element.
    % valuesNbase, valuesTbase, valuesMbase : values of the solution of axial force, shear force and bending moment
    % when only forces are present 
    % valuesNta, valuesTta, valuesMta : values of the solution of axial force, shear force and bending moment
    % when only thermal loads with uniform dilatation are present
    % valuesNtc, valuesTtc, valuesMtc : values of the solution of axial force, shear force and bending moment
    % when only thermal loads with constant curvature are present

    function[valuesN, valuesT, valuesM] = sumCellArraySolution(valuesNbase, valuesTbase, valuesMbase, ...
                                                            valuesNta, valuesTta, valuesMta, ... 
                                                                    valuesNtc, valuesTtc, valuesMtc)    
        %initialize the cell array of the values
        valuesN = {};
        valuesT = {};
        valuesM = {};
        for k = 1 : length(valuesNbase)
            valuesN = [valuesN, valuesNbase{k} + valuesNta{k} + valuesNtc{k}];
            valuesT = [valuesT, valuesTbase{k} + valuesTta{k} + valuesTtc{k}];
            valuesM = [valuesM, valuesMbase{k} + valuesMta{k} + valuesMtc{k}]; 
        end
    
    end




    %% ************************************************************************************************************
    %  ************       					  		plotMesh              	 		                   ************
    %  ************************************************************************************************************
    % This function plots on a figure the nodes of the mesh. The plot style is governed by 
    % the variable graphicOptions.
    % MESH : structured data containing the information of the mesh      
    % graphicOptions : structured data that contains info for the plot, for example
    %                   graphicOptions.color = 'k';
    %                   graphicOptions.lineWidth = 2;
    %                   graphicOptions.lineStile = '-';
    %                   graphicOptions.marker = None;
    %                   graphicOptions.markerSize = 2;

    function [] = plotMesh(MESH, graphicOptions)
        %2D plot of the nodes
        xPlot = MESH.nodes(:, 1);
        yPlot = MESH.nodes(:, 2);

        %plot of the nodes
        for k = 1 : size(MESH.elem, 1)
            %plot the nodes						
            line(xPlot(MESH.elem(k, :)), yPlot((MESH.elem(k, :))),...
            'Color', graphicOptions.color,...
            'LineWidth',graphicOptions.lineWidth,...
            'LineStyle', graphicOptions.lineStile,...
            'Marker', graphicOptions.marker,...
            'MarkerSize', graphicOptions.markerSize ...
            );
        end
    end




    %% ************************************************************************************************************
    %  ************       					  		plotDiagram            	 		                   ************
    %  ************************************************************************************************************
    % This function plots on a figure the field solution over the considered element (i.e. classical internal action
    % diagrams). The plot style is governed by the variable graphicOptions.
    % MESH : structured data containing the information of the mesh      
    % coordX : value of the x local coordinate over the element, this is computed calling 
    % FEM.calculateValueSolution(...)
    % valueFun : value of the solution field that must be displayed, this is computed calling 
    % FEM.calculateValueSolution(...)
    % scale : scalar concerning the scale of the plot
    % graphicOptions : structured data that contains info for the plot, for example
    %graphicOptions.faceColor = 'r';
    %graphicOptions.edgeColor = 'r';
    %graphicOptions.faceAlpha = 0.7;
    %graphicOptions.edgeAlpha = 1;

    function[] = plotDiagram(MESH, coordX, valueFun, scale, graphicOptions)
        global toleranceZero
        for e = 1 : size(MESH.elem)
            %define the value of the points of the shape function (in a classical 0xy reference frame)
            xFun = [0];
            yFun = [0];
            xFun = [xFun coordX{e}];
            if max(abs(valueFun{e})) > toleranceZero
                %when the solution is higher than the tolerance (i.e the solution is not simply a numerical error),
                %the solution must be calculated from the valuesFun (the sign is due to the fact that an inversion
                %is needed passing from global 0xy to local system)
                yFun = - scale * [yFun valueFun{e}];
            else
                %when the solution is too small compared to the tolerance (i.e. numerical errors), the plots are
                %forces to be zero
                yFun = - scale * [yFun zeros(size(valueFun{e}))];
            end
            xFun = [xFun MESH.length(e)];
            yFun = [yFun 0];

            %define the matrix of the points
            pointsBase = [xFun; yFun]; 

            %define the vector of the translation (the position vector of the first node of the element)
            vecTr = MESH.nodes(MESH.elem(e, 1), :).';

            %define the rotation of the plot
            angle = MESH.angle(e);
            R = [cos(angle),  -sin(angle); sin(angle), cos(angle)];

            %points of the plot
            pointsPlot = vecTr + R * pointsBase;
            xPlot = pointsPlot(1, :);
            yPlot = pointsPlot(2, :);

            %plot the function
            diagram = patch('XData', xPlot, ...
                    'YData', yPlot, ...
                    'faceColor', graphicOptions.faceColor, ...
                    'edgeColor', graphicOptions.edgeColor, ...
                    'faceAlpha', graphicOptions.faceAlpha, ...
                    'edgeAlpha', graphicOptions.edgeAlpha ...
                    );
        end
    end




    %% ************************************************************************************************************ 
    %  ************       					  		plotDeformed           	 		                   ************
    %  ************************************************************************************************************
    % This function plots on a figure the field solution over the considered element (i.e. classical internal action
    % diagrams). The plot style is governed by the variable graphicOptions.
    % MESH : structured data containing the information of the mesh      
    % coordX : value of the x local coordinate over the element, this is computed calling 
    % FEM.calculateValueSolution(...)
    % valueU : value of the solution field of local axial displacement, this is computed calling 
    % FEM.calculateValueSolution(...)
    % valueV : value of the solution field of local transversal displacement, this is computed calling 
    % FEM.calculateValueSolution(...)
    % scale : scalar concerning the scale of the plot
    % graphicOptions : structured data that contains info for the plot, for example
    %graphicOptions.color = 'b';
    %graphicOptions.lineWidth = 3;
    %graphicOptions.lineStile = '-';
    %graphicOptions.marker = None;
    %graphicOptions.markerSize = 2;

    function[] = plotDeformed(MESH, coordX, valuesU, valuesV, scale, graphicOptions)
       for e = 1 : size(MESH.elem)
           %define the value of the base points subdividing each element
           vecPosElem = [coordX{e}; zeros(size(coordX{e}))];
           
           %define the matrix of the displacement from each base point
           displ = scale * [valuesU{e}; -valuesV{e}]; 

           %define the rotation of the plot
           angleRot = MESH.angle(e);
           R = [cos(angleRot),  -sin(angleRot); sin(angleRot), cos(angleRot)];

           %define the vector of the translation (the position vector of the first node of the element)
           vecTr = MESH.nodes(MESH.elem(e, 1), :).';

           %points of the plot
           pointsPlot = vecTr + R * vecPosElem + R * displ;
           xPlot = pointsPlot(1, :);
           yPlot = pointsPlot(2, :);

           %plot the function
           line(xPlot, yPlot,...
               'Color', graphicOptions.color,...
               'LineWidth',graphicOptions.lineWidth,...
               'LineStyle', graphicOptions.lineStile,...
               'Marker', graphicOptions.marker,...
               'MarkerSize', graphicOptions.markerSize ...
               );
       end    
    end




    %% ************************************************************************************************************ 
    %  ************       					  		plotRefSysLocal           	 		               ************
    %  ************************************************************************************************************
    % This function plot the local reference frame with axis x-u and y-v.
    % MESH : structured data containing the information of the mesh      
    % scaleFraction : ratio between the length of the arrow and the maximum length of the elements
    % graphicOptions : structured data that contains info for the plot, for example
    % graphicOptions.colorX = 'r';
    % graphicOptions.colorY = 'b';
    % graphicOptions.maxHeadSize = 0.3;
     
    function [] = plotRefSysLocal(MESH, scaleFraction, graphicOptions)
		
        hold on    
        %calculate the length of each arrow
        lenArrow = scaleFraction * min(MESH.length);
    
        for e = 1 : size(MESH.elem, 1)
            %define the local x-axis
            axisEF = MESH.nodes(MESH.elem(e, 2), :) - MESH.nodes(MESH.elem(e, 1), :);
            
            %set the scale
            scale = lenArrow / norm(axisEF);
    
            %unit vector to display
            uVector = scale * axisEF;
            vVector = scale * [axisEF(2), -axisEF(1)];
    
            %position of the vector to be displayed
            posX = MESH.nodes(MESH.elem(e, 1), 1);
            posY = MESH.nodes(MESH.elem(e, 1), 2);
    
            %plot the local reference frame (local u vector)
            quiver( posX, ...
                    posY, ...
                    uVector(1), ...
                    uVector(2), ...
                    0, ...
                    'Color', graphicOptions.colorX,...
                   'MaxHeadSize', graphicOptions.maxHeadSize...
                    );
            
            %plot the local reference frame (local v vector)
            quiver( posX, ...
                    posY, ...
                    vVector(1), ...
                    vVector(2), ...
                    0, ...
                    'Color', graphicOptions.colorY,...
                    'MaxHeadSize', graphicOptions.maxHeadSize...
            );
        end    
    
    end

    


    %% ************************************************************************************************************ 
    %  ************       					  		automaticScaling       	 		                   ************
    %  ************************************************************************************************************
    % This function compute an automatic scale factor based on the scaleFactor.
    % MESH : structured data containing the information of the mesh      
    % valueFun : value of the solution field that must be displayed, this is computed calling 
    % scaleFactor : ratio between the maximum value of the plot and the maximum length of the elements
    % scale : scale factor to be used in the plot function for an automatic scale
    
    function[scale] = automaticScaling(MESH, valueFun, scaleFactor)
        %compute the maximum value of the function that must be plotted
        maxValEF = [];
        for e = 1 : size(MESH.elem)
            maxValEF(e) = max(abs(valueFun{e}));
        end
        maxVal =max(maxValEF);

        %retrieve the maximum length of the finite element
        maxLen = max(MESH.length);

        scale = scaleFactor * maxLen / maxVal;

    end




    %% ************************************************************************************************************
    %  ************       					  		libColor		      	 		                   ************
    %  ************************************************************************************************************
    % This function define a palette of colours that are used for all the plots. The colors are taken from the 
    % web sites https://flatuicolors.com/palette/defo. 
    % CodeColor : is the label of the name of the considered color
    % col : a 3-element vector that describes a color in Matlab convention (RGB)
    
    function [col] = libColor(codeColor)
        %from the web site : https://flatuicolors.com/palette/defo
        switch codeColor
            case 'red1'
                %fundamental red
                col = [255, 0, 0];
            case 'red2'
                %alizarin
                col = [231, 76, 60];
            case 'red3'	
                %pomegranate
                col = [192, 57, 43];
            case 'green1'
                %fundamental green
                col = [0, 255, 0];
            case 'green2'
                %turquoise
                col = [26, 188, 156];
            case 'green3'
                %emerald
                col = [46, 204, 113];				
            case 'green4'
                %greenSea
                col = [22, 160, 133];					
            case 'green5'
                %Nephritis
                col = [39, 174, 96];
            case 'blue1'
                %fundamental blue
                col = [0, 0, 255];				
            case 'blue2'
                %peter rivers
                col = [52, 152, 219];
            case 'blue3'
                %belize hole
                col = [41, 128, 185];
            case 'violet1'
                %amethyst
                col = [155, 89, 182];				
            case 'violet2'
                %wisteria
                col = [142, 68, 173];
            case 'black1'
                %wet asphalt
                col = [52, 73, 94];				
            case 'black2'
                %midnight blue
                col = [44, 62, 80];
            case 'yellow1'
                %sunflower
                col = [241, 196, 15];
            case 'orange1'
                %orange
                col = [243, 156, 18];
            case 'orange2'
                %carrot
                col = [230, 126, 34];			
            case 'orange3'
                %pumpkin
                col = [211, 84, 0];
            case 'grey1'
                %cloud
                col = [236, 240, 241];				
            case 'grey2'
                %silver
                col = [189, 195, 199];				
            case 'grey3'
                %concrete
                col = [149, 165, 166];				
            case 'grey4'
                %asbestos
                col = [127, 140, 141];
            case 'timber'
                %color of timber elements	
                col = [255, 192, 72];
                
            case 'timber2'
                %color of timber elements	
                col = [204, 142, 53];
        end
        
        %calculate the proper vector of color for Matlab
        col = col / 255;
    end


    




end % static methods
end % classdef