%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