Basic Structural Analysis Functions

FEDEASLab functions implement structural analysis concepts in parallel with their theoretical presentation in the introductory structural analysis course, thus illustrating the generality of the methods and their application to large scale and complex structural models. The intent is not only to reinforce the understanding of the structural response under different types and arrangements of loading, but also to facilitate system and parameter studies.

Function Organization

The structural analysis functions for the introductory course are organized in two folders: Basic and Utilities. The Basic folder contains 20 functions that directly implement the analysis concepts of the course in MATLAB . Many functions in this folder are provided in open source form (.m file format), so that interested readers can relate the implementation details with the theory presented in the course.

The FEDEASLab  functions in the Basic folder can be classified in the following categories:

1. Functions for Structural Matrices and Vectors: 8. These are listed in @tbl:StructMat.

2. Analysis Functions and Scripts: 4. These are listed in [@tbl:AnalFunScr].

3. Structural Model and Loading Functions: 3. These are listed in [@tbl:ModLoad].

4. Model and Element Utility Functions: 6, namely Localize, ElmLenOr, Create_PlasticLimitSurface, DefGeom_2dFrm, DeformShape2dFrm and Large2du2v_Frm which are not invoked directly, but are used either by the Model and Loading Functions or by Graphics Utility Functions to be discussed subsequently.

The folder Utilities contains 28 functions. Almost all functions in this folder are provided in content-obscured form (encoded .p file format), since the implementation details have no bearing on the structural analysis concepts and may only be of interest to readers with programming curiosity and skills. They can be classified in the following categories:

1. General Utility Functions: 2, namely CleanStart and Print_PDFile.

2. Graphics Utility Functions: 13. These are listed in [@tbl:GraphUtil].

3. Supporting Functions: 13. These are not invoked directly, but are used by the analysis functions and scripts and by the graphics utility functions. They are not described further in this text, since they are of little interest to the user, but details about input and output arguments may be found in the online help.

The user defines the geometry and properties of the structural model, and then specifies the structural analysis procedure(s) and the graphical post-processing of the results by selecting suitable functions from [@tbl:StructMat]-[@tbl:GraphUtil] in a MATLAB  script file, as described in the following sections. Scripting affords extensive flexibility for the response determination and the presentation of the results, but puts the burden on the user to assemble the whole by understanding its parts. This supports the open-ended educational mission of FEDEASLab , as opposed to the black-box nature of commercial analysis software.

The following sections describe the functions in [@tbl:StructMat]-[@tbl:GraphUtil] with examples. The notation for function arguments and local variables agrees with the notation of the course reader, as summarized in Appendix Chap. A of the reader¹. The examples in this chapter correspond to examples of the course reader with the MATLAB  scripts for all examples available in the folder Examples of the basic FEDEASLab package.

Data Organization

The initial set of analysis functions for the concepts of an introductory course on structural analysis is built around two data structures: Model  for the geometric description of the structural model, and ElemData  for the element and material properties and element loading, as is evident from the input argument list of the functions in [@tbl:StructMat]-[@tbl:GraphUtil]. The organization of structural model information in two data structures serves to emphasize the concept of the introductory structural analysis course that, under linear conditions, the governing static and kinematic relations depend only on the geometry of the structural model and not on any element properties. Accordingly, the functions B_matrix for setting up the linear equilibrium matrix \(\mathbf{B}\) and A_matrix for setting up the linear kinematic matrix \(\mathbf{A}\) of the structure take as single input argument the data structure Model with model geometry information.

Data Structure Model

In FEDEASLab information about the structural model is collected in the data structure Model , which is set up by the function Create_SimpleModel from information supplied by the user. The function Create_SimpleModel is limited to models consisting of 2-node, 2d or 3d truss and frame elements. The function numbers the degrees of freedom (DOFs) of the structure, as discussed in the course reader, thus facilitating their identification for small structural models. With this degree of freedom (DOF) numbering the analysis functions set up full matrices that do not take advantage of the small number of non-zero terms in large structural models. The general model creation function Create_Model does not have limitations on element type and uses a DOF numbering scheme that minimizes the storage of structure matrices. The use of sparse matrices with minimum storage is essential for the analysis of large structural models with 2d and 3d finite elements.

% clear memory; close open windows
CleanStart
% define model geometry
XYZ(1,:) = [  0   0];
XYZ(2,:) = [  8   0];
XYZ(3,:) = [ 16   0];
XYZ(4,:) = [  8   6];  
% element connectivity array
CON(1,:) = [  1   2];
CON(2,:) = [  2   3];
CON(3,:) = [  1   4];
CON(4,:) = [  2   4];
CON(5,:) = [  3   4];
% boundary conditions
BOUN(1,:) = [ 1 1];
BOUN(3,:) = [ 0 1];
% specify element type
ne = size(CON,1);                      
[ElemName{1:ne}] = deal('Truss');

[@fig:ModelInp] shows the input data for the structure of Ex. 3.1 of the reader in the form of numerical arrays XYZ, BOUN, and CON. The rows of XYZ correspond to the nodes and the \(ndm\) columns to the node coordinates, where \(ndm\) is the dimension of the structural model. The rows of BOUN also correspond to the nodes and the \(ndf\) columns to node restraint switches, where \(ndf\) is the maximum number of DOFs per node for the model. The DOF order for each node is: 1 = force in \(X\), 2 = force in \(Y\), 3 = moment about \(Z\) (for 2d model), or force in \(Z\) (for 3d model), 4,5,6 = moment about \(X\), \(Y\) and \(Z\), respectively (for 3d model). The switch value 0 indicates a free and the value 1 a restrained DOF. CON can be a numerical array for a structural model with only 2-node elements. For finite element models with elements having a different number of nodes, CON needs to be specified as a cell array so that the row contents can vary from element to element. Each row of the numerical CON array contains the node numbers to which the corresponding element connects. Finally, ElemName is a cell array that accommodates element names with a different number of characters. The names 'Truss' and '2dFrm' suffice for tasks that do not require element and material properties. For tasks requiring element and material properties the element names need to be 'LinTruss', 'Lin2dFrm', and 'Lin3dFrm' for the 2d/3d truss, the 2d frame, and the 3d frame element, respectively.

At the start of each input script the command CleanStart issues two MATLAB  commands, one for clearing the computer memory of all stored variable values (clearvars), and the second for closing any open figure windows (close all). It may also be expedient to include the call to the function Create_Units at the start of each script, so as to take advantage of the unit definitions in the return argument Units and obviate the need for unit conversion in the values of variables with different units.

The command

Model = Create_SimpleModel(XYZ,CON,BOUN,ElemName)

generates the data structure Model . Once Model is created its contents can be displayed by executing disp(Model).

disp(Model)
          nn: 4
         ndm: 2
         XYZ: [4x2 double]
          ne: 5
         CON: {5x1 cell}
    ElemName: {'Truss'  'Truss'  'Truss'  'Truss'  'Truss'}
         nen: [2 2 2 2 2]
          nq: [1 1 1 1 1]
         ndf: [2 2 2 2 2]
          nt: 8
        BOUN: [4x2 double]
          nf: 5
         DOF: [4x2 double]

@lst:ModelCont shows the fields and some of their contents for the data structure Model  of the truss in [@fig:ModelInp] : \(nn\) is the number of nodes of the structural model, \(ndm\) the model dimension, XYZ the node coordinate array, \(ne\) the number of elements, CON a cell array of element connectivity with the row number corresponding to the element number (by conversion of the numerical array CON supplied as input argument), \(nq\) the number of basic forces in each element, \(ndf\) the number of local DOFs/node and \(nen\) the number of end nodes of each element, ElemName the cell array of element types, \(nt\) the total number of DOFs, BOUN the array of node restraint switches, \(nf\) the number of free DOFs, and DOF the global DOF numbering array. Typing Model or disp(Model) at the command prompt only shows contents of scalar fields of Model  and numeric or cell arrays that can be displayed in a single line, as [@lst:ModelCont] shows (numeric arrays use [ ] operators and cell arrays use { } operators in MATLAB ). To display the contents of fields of other array types, type the data structure name and the field name separated by a period as input argument to the MATLAB  function disp to avoid echoing the variable name:

: Model.BOUN for node restraint switches

Data Structure ElemData for Element Properties and Loading

ElemData  is the second data structure of FEDEASLab after Model . It is specified entirely by the user. Because the element property data vary from element to element, it is necessary to use a cell array for ElemData  so that the array contents can vary from element to element. ElemData  is a one dimensional cell array with the row number corresponding to the element number. It has the fields in [@tbl:ElDatafields] with optional fields with default value 0 shown in parentheses.

[@tbl:ElDatafields] shows that ElemData includes a field Np for the plastic axial capacity and a field Mp for the plastic flexural capacity. These are only required for the plastic analysis of large structural models. Otherwise, the plastic capacities of the elements may be directly specified with vector \(\boldsymbol{Q}_{pl}\), as discussed later in this chapter.

function [xyz,id] = Localize (Model,el)

XYZ = Model.XYZ;        % node coordinates
DOF = Model.DOF;        % array with dof numbers for all nodes

CON = Model.CON{el};    % extract connectivity array for element
ndf = Model.ndf(el);    % extract no of dofs/node for element

% extract element coordinates into array xyz;
% use CON array to extract appropriate rows of global XYZ array
xyz = XYZ(CON(CON>0),:)';

% extract dof numbers into array id
% use CON array to extract appropriate rows of DOF array
id = DOF(CON(CON>0),1:ndf)';
% reshape id array into vector   
id = id(:);

function [L,dcx] = ElmLenOr (xyz)

% coordinate differences Dx, Dy, Dz depending on problem dimension
Dxyz = xyz (:,2) - xyz(:,1);   
% element length
L    = sqrt (Dxyz'*Dxyz);
% direction cosines of element orientation
dcx  = Dxyz./L;

General Utility Functions

Once the geometry and other relevant information about the structural model is collected in data structure Model  by the function Create_SimpleModel, information about the geometry and the global DOF correspondence for a particular element can be extracted with the help of two general utility functions: Localize and ElmLenOr. All functions with input arguments Model  and ElemData in [@tbl:StructMat] invoke these utility functions for setting up the corresponding structural matrices or vectors.

Function Localize

This function locates the element el in the Model  and returns its end node coordinates in array xyz, and the local-global DOF correspondence in vector id. The first column of the two dimensional array xyz contains the coordinates of end \(i\) of the element and the second column those of end \(j\).

Function ElmLenOr

This function returns the element length L and the direction cosines of the element \(x\)-axis in vector dcx from the end node coordinates supplied in array xyz.

Structural Equilibrium

The global DOF equilibrium equations of the structural model are given by Eq. 3.77 of the reader. As discussed in Chap. 3 of the reader, the equilibrium equations at the free DOFs of the structural model in Eq. 3.79a of the reader are used to solve for the basic forces \(\boldsymbol{Q}\). For determining the solution of these equations it is necessary to set up the applied nodal forces \(\boldsymbol{P}_f\), the equivalent nodal forces \(\boldsymbol{P}_{wf}\) and the equilibrium matrix \(\mathbf{B}_f\). The latter is extracted from the full equilibrium matrix \(\mathbf{B}\).

Assembly of Equilibrium Matrix \(\mathbf{B}\)

The assembly of the equilibrium matrix \(\mathbf{B}\) of the structural model in Eq. 3.4.1 of the reader involves the determination of the equilibrium matrix \(\mathbf{b}_g\) of each element and the re-assignment of the rows of \(\mathbf{b}_g\) to the rows of \(\mathbf{B}\) based on the local-global DOF number correspondence in accordance with the \(\boldsymbol{id}\) array of the element. Setting up the equilibrium matrix \(\mathbf{b}_g\) of each element according to Eq. 3.57 of the reader, Eq. 3.59 of the reader or Eq. 3.60 of the reader for a 2d frame, a 2d truss and a 3d truss element, respectively, requires the element length \(L\) and the direction cosines of the undeformed element \(x\)-axis. These are set up by the utility function ElmLenOr in conjunction with the utility function Localize, which also sets up the \(\boldsymbol{id}\) array for the element, as described earlier.

The assembly of the \(\mathbf{B}\) matrix according to Eq. 3.76 of the reader is implemented in FEDEASLab  function B_matrix , as shown in @lst:Bmatimpl. The function B_matrix  calls the internal function bg_matrix to set up the equilibrium matrix \(\mathbf{b}_g\) of each element according to Eq. 3.57 of the reader, Eq. 3.59 of the reader or Eq. 3.60 of the reader for a 2d frame, a 2d truss and a 3d truss element, respectively. The function syntax is given in [@lst:bgimpl].

Because the data structure Model  does not have information about basic force releases, it is the user’s responsibility to remove the columns of the \(\mathbf{B}\)-matrix corresponding to basic force releases so that the difference between the number of columns and the number of rows of the \(\mathbf{B}\)-matrix matches the degree of static indeterminacy NOS of the structural model. A convenient way for setting up an index of element release locations is furnished by the utility function H_index with information about the element releases in ElemData , as will be discussed later.

After setting up the complete equilibrium matrix \(\mathbf{B}\) for the structural model, the equilibrium matrix \(\mathbf{B}_f\) for the free DOFs corresponds to the upper \(nf\) rows of \(\mathbf{B}\). It can be extracted with the statement

function B = B_matrix (Model)

% assemble structure equilibrium matrix
ne = Model.ne;          % number of elements in structural model
k  = 0;                 % initialize column index into matrix B
for el=1:ne             % loop over all elements
   % locate element in Model and return end coordinates and id array
   [xyz,id] = Localize(Model,el);
   % form element equilibrium matrix bg
   bg = bg_matrix(Model.ElemName{el},xyz);
   % assemble element equilibrium matrix bg into static matrix B
   nc = size(bg,2);     % determine number of columns of matrix bg
   B(id,k+1:k+nc) = bg; % insert bg into id locations of matrix B
   k  = k+nc;           % augment column index by number of columns
end

function bg = bg_matrix(ElemName,xyz)

% determine element length and orientation
[L,dcx] = ElmLenOr(xyz);
if contains(ElemName,'Truss')
   % truss element
   bg = [-dcx; dcx];
elseif contains(ElemName,'2dFrm')
   % 2d frame element
   dXL = dcx(1);
   dYL = dcx(2);
   bg  = [-dXL  -dYL/L  -dYL/L;
          -dYL   dXL/L   dXL/L;
             0     1       0  ;
           dXL   dYL/L   dYL/L;
           dYL  -dXL/L  -dXL/L;
             0     0       1  ];
else
   error ('only truss and 2d frame elements supported; check syntax')
end

Specification of Nodal Forces

%% specify loading
Pe(2,2) = -5; % node 2, direction Y
Pe(4,1) = 10; % node 4, direction X

% generate applied force vector Pf

The next task for setting up the node equilibrium equations of the structural model is the specification of the applied and the equivalent nodal forces. For this task it is impractical to know the global DOF numbering for directly setting up the applied and the equivalent nodal force vectors \(\boldsymbol{P}_f\) and \(\boldsymbol{P}_{wf}\) at the free DOFs, respectively. Instead, it is convenient to specify the nodal forces in a two dimensional array Pe with row numbers corresponding to node numbers and column numbers corresponding to DOF numbers, where 1 = force in \(X\), 2 = force in \(Y\), 3 = moment about \(Z\) (for 2d model), or force in \(Z\) (for 3d model), 4,5,6 = moment about \(X\), \(Y\) and \(Z\), respectively (for 3d model). With the global DOF numbering available in the data structure Model , the FEDEASLab function Create_NodalForces  in @tbl:ModLoad] sets-up the nodal force vector. The input arguments to the function are the data structure Model and the two dimensional nodal force array Pe, while the output argument is the nodal force vector Pf. @fig:PfInp shows the syntax for the applied nodal forces of the truss from Ex. 3.1 of the reader. The equivalent nodal forces \(\boldsymbol{P}_w\) are set up by the function Create_PwForces in [@tbl:ModLoad] for the element loading \(w\) specified in data structure ElemData  according to [@tbl:ElDatafields].

Solution for Statically Determinate Structures

For a stable, statically determinate structure the equilibrium matrix \(\mathbf{B}_f\) is square with full rank. The MATLAB  function rank can be used to make sure that this is indeed the case. Given the nodal forces \(\boldsymbol{P}_f\) and \(\boldsymbol{P}_{wf}\) at the free DOFs of a stable structural model, a unique solution of the system of linear equations in Eq. 3.79a of the reader exists, resulting in the basic forces \(\boldsymbol{Q}\) of the structure. The product of the basic forces \(\boldsymbol{Q}\) with the full equilibrium matrix \(\mathbf{B}\) gives the complete nodal force vector \(\boldsymbol{P}\) according to Eq. 3.77 of the reader. The entries in the upper \(nf\) rows of \(\boldsymbol{P}\) correspond to the nodal forces \(\boldsymbol{P}_f-\boldsymbol{P}_{wf}\), and the entries in the remaining rows correspond to the support reactions \(\boldsymbol{R}\) according to Eq. 3.78 of the reader.

[@lst:QRdet] shows an excerpt from a MATLAB  script file that follows the generation of Model  in [@fig:ModelInp] and of the applied nodal force vector \(\boldsymbol{P}_f\) in [@fig:PfInp]. The script shows the determination of the equilibrium matrix \(\mathbf{B}\) with function B_matrix  and its display, the solution of the equilibrium equations in Eq. 3.79a of the reader with the MATLAB  \(\backslash\) command, and the resulting basic forces \(\boldsymbol{Q}\). The basic forces \(\boldsymbol{Q}\) are then used to determine the support reactions \(\boldsymbol{R}\) after first multiplying these with the full equilibrium matrix \(\mathbf{B}\) and then extracting the rows corresponding to the restrained DOFs.

Solution for Statically Indeterminate Structures

The general solution of the equilibrium equations in Eq. 3.79a of the reader for a equilibrium matrix \(\mathbf{B}_f\) with full rank and with more columns than rows uses the row-echelon form of the matrix. Conversion of a matrix to row-echelon form uses the Gauss-Jordan elimination procedure. This well established procedure of linear algebra is implemented in the function rref of MATLAB . It can be used to establish the rank of \(\mathbf{B}_f\) and provide the particular and the homogeneous solutions of the equilibrium equations according to Fig. 3.5.2 of the reader.

Following Fig. 3.5.2 of the reader the general solution of the equilibrium equations in Eq. 3.100 of the reader is based on the basic force influence matrix \(\bar{\mathbf{B}}_i\) of the primary structure for the applied and equivalent nodal forces, and the basic force influence matrix \(\bar{\mathbf{B}}_x\) for the redundant basic forces \(\boldsymbol{Q}_x\). FEDEASLab  function BbariBbarx_matrix  takes the equilibrium matrix \(\mathbf{B}_f\) as input argument and uses its reduced-row echelon form to set up the force influence matrices \(\bar{\mathbf{B}}_i\) and \(\bar{\mathbf{B}}_x\). The optional function argument ind_r can be used to select specific basic forces as redundants. If ind_r is not specified, the Gauss-Jordan elimination process selects the redundant basic forces automatically, and the function BbariBbarx_matrix  returns the index array ind_x to the basic force vector \(\boldsymbol{Q}\), so that the redundant basic forces \(\boldsymbol{Q}_x\) can be extracted with Qx = Q(ind_x). The function syntax is

function [Bbari, Bbarx, ind x] = BbariBbarx matrix ( Bf, ind r)

Bf = [ 1/6  1/6  0    0    0  0.8;
        0    1   1    0    0   0 ;
        0    0 -1/4 -1/4  1/4  0 ;
        0    0   0    1    1   0 ];

% select basic force redundants
ind_r = [1,4];
[Bbari,Bbarx,ind_x] = BbariBbarx_matrix(Bf,ind_r);
Pf = [20;0;-30;0];
Qp = Bbari*Pf;

%% Matlab script for force influence matrix generation

% Department of Civil and Environmental Engineering
% University of California, Berkeley
% Professor Filip C. Filippou

% clear memory
CleanStart

% static matrix Bf of braced frame (Ex 2.7)
Bf = [ 1/6  1/6  0    0    0  0.8;
        0    1   1    0    0   0 ;
        0    0 -1/4 -1/4  1/4  0 ;
        0    0   0    1    1   0 ];

% select basic force redundants
ind_r = [1,4];

[@lst:Bbar4brfr] shows the direct specification of the \(\mathbf{B}_f\) matrix for the braced frame of Ex. 3.13 of the reader that corresponds to the equilibrium equations in Eq. 3.122 of the reader. The function BbariBbarx_matrix  then determines the force influence matrices \(\bar{\mathbf{B}}_i\) and \(\bar{\mathbf{B}}_x\) for the selection of \(\boldsymbol{Q}_1\) and \(\boldsymbol{Q}_4\) as redundant basic forces, as specified by ind_r=[1,4]. With the direct specification of the nodal force vector \(\boldsymbol{P}_f\) the particular solution for the given nodal forces can be determined by the product \(\bar{\mathbf{B}}_i\mspace{1.5mu}\boldsymbol{P}_f\). The resulting values agree with those from Ex. 3.13 of the reader (b) and (c).

Structural Kinematics

The structural kinematic relations in Eq. 4.53 of the reader relate the element deformations \(\boldsymbol{V}\) to the free global DOF displacements \(\boldsymbol{U}_f\) through the kinematic matrix \(\mathbf{A}_f\) of the structure. With a given set of compatible element deformations \(\boldsymbol{V}\) the solution of the kinematic relations in Eq. 4.53 of the reader gives the values of the free DOF displacements \(\boldsymbol{U}_f\). This approach to the determination of the free DOF displacements is of particular interest for statically determinate structures, for which the number of element deformations is equal to the number of unknown free DOF displacements, so that the element deformations are not constrained by compatibility conditions, as the examples of Chap. 4 of the reader demonstrate.

Assembly of Kinematic Matrix \(\mathbf{A}\)

Because \(\mathbf{A}_f^T=\mathbf{B}_f\) according to the principle of virtual work, it is not necessary to create separate functions for setting up the equilibrium matrix \(\mathbf{B}\) and the kinematic matrix \(\mathbf{A}\). Nonetheless, FEDEASLab has two separate functions in the interest of completeness.

The assembly of the kinematic matrix \(\mathbf{A}\) is completely analogous to that of the equilibrium matrix \(\mathbf{B}\) and the corresponding functions practically identical. Nonetheless, a brief description of the process is given in the following.

The assembly of the kinematic matrix \(\mathbf{A}\) of the structural model in Fig. 4.4.5 of the reader involves the determination of the kinematic matrix \(\mathbf{a}_g\) of each element and the re-assignment of the columns of \(\mathbf{a}_g\) to the columns of \(\mathbf{A}\) based on the local-global DOF number correspondence in accordance with the \(\boldsymbol{id}\) array of the element. The kinematic matrix \(\mathbf{a}_g\) of each element is set up according to Eq. 4.37 of the reader, Eq. 4.38 of the reader, or Eq. 4.39 of the reader for 2d frame, 2d truss and 3d truss, respectively, and requires the element length \(L\) and the \(x\)-axis orientation in the form of the direction cosines, as discussed for the equilibrium matrix \(\mathbf{b}_g\) in [@sec:FLBfAssembly].

function A = A_matrix (Model)

ne = Model.ne;           % number of elements in structural model                  
k  = 0;                  % initialize row index into matrix A
for el=1:ne              % loop over all elements                   
   % locate element in Model and return end coordinates and id array
   [xyz,id] = Localize(Model,el);
   % form element compatibility matrix ag
   ag = ag_matrix(Model.ElemName{el},xyz);
   % assemble element compatibility matrix ag into structure matrix A
   nc = size(ag,1);      % number of rows of matrix ag
   A(k+1:k+nc,id) = ag;  % insert ag into location of matrix A
   k  = k+nc;            % augment row index by number of rows in ag
end

function ag = ag_matrix(ElemName,xyz)

% determine element length and orientation
[L,dcx] = ElmLenOr(xyz);
if contains(ElemName,'Truss')
   % truss element
   ag = [-dcx' dcx'];
elseif contains(ElemName,'2dFrm')
   % 2d frame element
   dXL = dcx(1);
   dYL = dcx(2);
   ag  = [-dXL     -dYL     0   dXL      dYL      0;
          -dYL/L    dXL/L   1   dYL/L   -dXL/L    0;
          -dYL/L    dXL/L   0   dYL/L   -dXL/L    1];
else
   error ('only truss and 2d frame elements supported')
end

@lst:Amatimpl gives the syntax for function A_matrix . The function A_matrix  calls the internal function ag_matrix to set up the kinematic matrix \(\mathbf{a}_g\) of each element according to Eq. 4.37 of the reader, Eq. 4.38 of the reader, or Eq. 4.39 of the reader for 2d frame, 2d truss and 3d truss, respectively. The function syntax is given in @lst:agimpl.

Because the data structure Model  does not have information about basic force releases, it is the user’s responsibility to remove the rows of the \(\mathbf{A}\)-matrix corresponding to element ends with release deformations \(\boldsymbol{V}_h\). After setting up the complete kinematic matrix \(\mathbf{A}\) for the structural model, the kinematic matrix \(\mathbf{A}_f\) for the free DOFs corresponds to the first \(nf\) columns of \(\mathbf{A}\). For a stable, statically determinate structure the kinematic matrix \(\mathbf{A}_f\) is square with full rank. Given the element deformations \(\boldsymbol{V}\), a unique solution of the system of linear equations in Eq. 4.53 of the reader exists, resulting in the displacements \(\boldsymbol{U}_f\) at the free DOFs of the structural model.

The excerpt from a MATLAB  script file in [@lst:TrussKinem] follows the input data for the plane truss of Ex. 4.2 of the reader and the command Create_SimpleModel. The element deformation vector V0 is populated with zeroes. It has size equal to the total number of basic forces sum(Model.nq) for this model which does not have element releases. After specifying the non-zero values of the element deformations the solution of the kinematic relations in Eq. 4.53 of the reader gives the free DOF displacements Uf. The values are identical with those of Ex. 4.2 of the reader with the corresponding deformed shape in Fig. 4.28 of the reader.

Force-Deformation Relations for Linear Material Response

Chap. 6 of the reader establishes the deformation-force relations of the truss element in Eq. 7.9 of the reader and of the 2d frame element in Eq. 6.39 of the reader. The collection of these relations for all elements of the structural model results in the structure deformation-force relations in Eq. 7.2 of the reader. Consequently, the collection of the element flexibility matrices \(\mathbf{f}\) to form the matrix \(\mathbf{F}_s\) and the collection of the initial element deformation vectors \(\boldsymbol{v}_0\) to form the vector \(\boldsymbol{V}_{0}\) depends on setting up \(\mathbf{f}\) and \(\boldsymbol{v}_0\) for all elements of the structural model. The latter depend on the material properties and the loading of each element which are supplied in ElemData  according to the information in [@tbl:ElDatafields].

Force-Deformation Functions

With the element property and loading information in ElemData , the FEDEASLab function Fs_matrix  collects the element flexibility matrices into the block-diagonal matrix \(\mathbf{F}_s\) according to Eq. 7.2 of the reader. Similarly, the function V0_vector  collects the initial element deformations into vector \(\boldsymbol{V}_{0}\).

function Fs = Fs_matrix (Model,ElemData,Roption)

ne = Model.ne;
f  = cell(ne,1);
% loop over all elements in structural model
for el = 1:ne
  % locate element in Model and return end coordinates xyz
  xyz = Localize(Model,el);
  % determine element length from end coordinates
  L = ElmLenOr(xyz);
  % set up element flexibility matrix
  f{el} = f_matrix (Model.ElemName{el},ElemData{el},L,Roption);
end
% block-diagonal array of element flexibility matrices
Fs = blkdiag (f{:});

function f = f_matrix (ElemName,ElemData,L,Roption)

E  = ElemData.E;    % elastic modulus
A  = ElemData.A;    % area
EA = E*A;           % axial stiffness

switch ElemName
  case 'LinTruss'
    % linear truss element
    f  = L/EA;
  case 'Lin2dFrm'
    % 2d linear frame element
    I  = ElemData.I;
    EI = E*I;
    f  = [L/EA   0          0   ;
           0    L/(3*EI)  -L/(6*EI) ;
           0   -L/(6*EI)   L/(3*EI)];
    % extract necessary submatrix in the presence of release(s)
    ide = 1:3;
    if Roption && isfield(ElemData,'Release')
      idr = find(ElemData.Release);
      ide = setdiff(ide,idr);
    end
    f = f(ide,ide);
  otherwise
    error ('only linear truss and 2d frame elements supported')
end

The syntax of function Fs_matrix  is shown in @lst:Fsmatimpl ². The function Fs_matrix  calls the internal function f_matrix to set up the flexibility matrix \(\mathbf{f}\) for each element according to Eq. 6.39 of the reader with the flexibility matrices for the truss element in Eq. 6.25 of the reader and for the 2d frame element in Eq. 6.36 of the reader. The function syntax is given in @lst:fimpl . The function V0_vector  has similar syntax and is not listed here.

According to Eq. 6.5 of the reader, the element force-deformation relations are established as the inverse of the deformation-force relations of the element with the basic element stiffness matrix \(\mathbf{k}\) equal to the inverse of the element flexibility matrix \(\mathbf{f}\) and the initial force vector \(\boldsymbol{q}_0\) according to Eq. 6.52 of the reader.

function Ks = Ks_matrix (Model,ElemData)

ne = Model.ne;
k  = cell(ne,1);
for el = 1:ne
   % locate element in Model and return end coordinates xyz
   xyz = Localize(Model,el);
   % determine element length from end coordinates
   L = ElmLenOr(xyz);
   % set up element stiffness matrix
   k{el} = k_matrix (Model.ElemName{el},ElemData{el},L);
end
% block-diagonal array of element stiffness matrices
Ks = blkdiag (k{:});

function k = k_matrix (ElemName,ElemData,L)

E  = ElemData.E;    % elastic modulus
A  = ElemData.A;    % area
EA = E*A;           % axial stiffness

switch ElemName
  case 'LinTruss'
    % linear truss element
    k  = EA/L;
  case 'Lin2dFrm'
    % 2d linear frame element
    I  = ElemData.I;
    EI = E*I;
    f  = [L/EA   0          0       ;
           0    L/(3*EI)  -L/(6*EI) ;
           0   -L/(6*EI)   L/(3*EI)];
    % extract necessary submatrix in the presence of release(s)
    ide = 1:3;
    if isfield(ElemData,'Release')
      idr = find(ElemData.Release);
      ide = setdiff(ide,idr);
    end
    f = f(ide,ide);
    % invert element flexibility to get basic element stiffness
    k = inv(f);
  otherwise
    error ('only linear truss and 2d frame elements supported')
end

function q0 = q0_vector (ElemName,ElemData,L)

E  = ElemData.E;    % elastic modulus
A  = ElemData.A;    % area
EA = E*A;           % axial stiffness

switch ElemName
  case 'LinTruss'
    % linear truss element
    if ~isfield(ElemData,'e0'), ElemData.e0 = 0; end
    if ~isfield(ElemData,'q0'), ElemData.q0 = 0; end
    e0 = ElemData.e0;   % initial axial strain
    q0 = ElemData.q0;   % initial prestressing force
    q0 = q0 - EA*e0;    % fixed-end force vector
  case 'Lin2dFrm'
    % 2d linear frame element
    I  = ElemData.I;
    EI = E*I;
    if ~isfield(ElemData,'e0'), ElemData.e0 = [0;0]; end
    if ~isfield(ElemData,'w'),  ElemData.w  = [0;0]; end
    e0 = ElemData.e0;   % initial section deformations
    w  = ElemData.w;    % uniformly element load in x and y
    % 2d linear frame element
    v0 = [e0(1)*L; -e0(2)*L/2 ; e0(2)*L/2];
    v0 = v0 + [w(1)*L*L/(2*EA) ; w(2)*L^3/(24*EI) ; -w(2)*L^3/(24*EI)];
    % flexibility matrix
    f  = [L/EA   0          0       ;
           0    L/(3*EI)  -L/(6*EI) ;
           0   -L/(6*EI)   L/(3*EI)];
    % index of continuous element deformations
    ide = 1:3;
    if isfield(ElemData,'Release')
      idr = find(ElemData.Release);
      ide = setdiff(ide,idr);
    end
    v0 = v0(ide);
    f  = f(ide,ide);
    q0 = -f\v0;
  otherwise
    error ('only linear truss and 2d frame elements supported')
end

FEDEASLab function Ks_matrix  collects the basic element stiffness matrices into the block-diagonal matrix \(\mathbf{K}_s\) according to Eq. 6.11 of the reader. Similarly, the function Q0_vector  collects the initial element forces into vector \(\boldsymbol{Q}_{0}\).

The function Ks_matrix  calls the internal function k_matrix to set up the stiffness matrix \(\mathbf{k}\) for each element according to Eq. 6.52 of the reader with the stiffness for the truss element in Eq. 6.42 of the reader and the stiffness matrix of the 2d frame element in Eq. 6.45 of the reader or Eq. 6.49 of the reader depending on the presence of an end moment release. The syntax of function Ks_matrix  is shown in [@lst:Ksmatimpl], while the syntax for the internal function k_matrix is given in [@lst:kimpl].

The syntax of function Q0_vector  is very similar except for the collection of vectors instead of matrices. The function calls the internal function q0_vector with the function syntax in [@lst:q0impl] to set up the initial force vector \(\boldsymbol{q}_0\) of each element according to Eq. 6.43 of the reader, Eq. 6.46 of the reader and Eq. 6.50 of the reader.

Force Method of Analysis

The steps of the force method of analysis are summarized in Eq. 8.1 of the reader. This step sequence is implemented in the MATLAB  script file S_ForceMethod.

In contrast to a function, a script file operates on the variables in the workspace memory and stores its output also in the workspace memory. FEDEASLab makes use of script files for two purposes:

1. to define several variables at once and place them in the workspace for easy access or to issue a few simple MATLAB  commands with a single call rather than individually;

2. to collect a series of commands involving FEDEASLab functions into solution procedures that can be invoked with a single call from the workspace; such scripts are identified with prefix S_ and can be easily modified to suit individual needs. The limitation of such a script, however, is that variables serving as input arguments to the functions in the script should be present in the workspace with the same syntax as the function arguments in the script.

All script variables reside in the MATLAB  workspace and are, therefore, of global character. Such variables may lead to name conflicts, if care is not exercised. In contrast, variables in functions are local and can only be accessed if they appear in the input or output argument list. Clearly, this is the preferred approach for data protection in a modular computing environment. For this reason, scripts are used sparingly in FEDEASLab.

% form static matrix B for all dofs
B  = B_matrix (Model);
% get element cell array for continuous element deformations
iced = H_index (Model,ElemData);
% concatenate indices to single index vector ic
ic = [iced{:}];
% form Bf matrix without columns corresponding to releases
Bf = B(1:Model.nf,ic);
% determine force influence matrices Bbari, Bbarx
if ~exist('ind_r','var'), ind_r = []; end
[Bbari,Bbarx] = ForceInfl_matrix (Bf,ind_r);
% set up Fs matrix and V0 vector
Fs  = Fs_matrix (Model,ElemData);
V0  = V0_vector (Model,ElemData);

% determine redundant basic forces
Qp  = Bbari*(Pf-Pwf);
Fxx = Bbarx'*Fs*Bbarx;
Qx  = -Fxx\(Bbarx'*(Fs*Qp+V0));
% determine final basic forces
Q   = Qp + Bbarx*Qx;
% determine element deformations Veps (Ve for short)
Ve  = Fs*Q + V0;
% determine free dof displacements
Uf  = Bbari'*Ve;
% complete Q vector with release values and return Ve 
[Q,Ve] = Complete_QV(Model,ElemData,Q);

The script S_ForceMethod depends on the definition of the structural model geometry in the data structure Model , and on the specification of the element property and loading information in the data structure ElemData . The specification of an index ind_r for the redundant basic forces is optional, as is the specification of the nodal forces Pf and Pwf. The default value for the nodal forces is zero.

The command sequence for the implementation of the force method of analysis according to Eq. 8.1 of the reader is divided in two parts in [@lst:FMsteps]:

1. the setting up of the necessary matrices \(\mathbf{B}\), \(\bar{\mathbf{B}}_i\), \(\bar{\mathbf{B}}_x\), \(\mathbf{F}_s\) and the vector \(\boldsymbol{V}_{\scriptstyle{0}}\) from information in Model  and ElemData , and (b) the solution steps of the force method. The script generates the basic forces \(\boldsymbol{Q}\) of the structural model, the corresponding element deformations \(\boldsymbol{V}_{\varepsilon}\), and the free DOF displacements \(\boldsymbol{U}_f\). These can be used for response interpretation, as will be discussed in [@sec:Post].

The script S_ForceMethod uses two auxiliary functions: (a) H_index to set up a cell array of indices for continuous element deformations from information in the field Release of ElemData , and (b) Complete_QV which uses the same information to pad the basic forces \(\boldsymbol{Q}\) from the force method of analysis with zero values at the releases, so that it can determine the element deformations \(\boldsymbol{V}_{\varepsilon}\) at the ends with flexural releases and plot the deformed shape of the corresponding element, if required.

Displacement Method of Analysis

The steps of the displacement method of analysis are described in Eq. 9.3 of the reader and Eq. 9.4 of the reader. A summary of these steps for implementation in FEDEASLab is:

1. With the kinematic matrix \(\mathbf{A}_f\) set up the stiffness matrix \(\mathbf{K}_f\) and the initial force vector \(\boldsymbol{P}_{0}\) \[ \mathbf{K}_f= \mathbf{A}_f^T\mspace{1.5mu}\mathbf{K}_s\mspace{1.5mu}\mathbf{A}_f \quad \quad \boldsymbol{P}_{0}= \mathbf{A}_f^T\mspace{1.5mu}\boldsymbol{Q}_{0}+ \boldsymbol{P}_{wf} \]

2. Solve the equilibrium equations for the free DOF displacements \(\boldsymbol{U}_f\) \[\boxed{ \boldsymbol{P}_f= \mathbf{K}_f\mspace{1.5mu}\boldsymbol{U}_f+ \boldsymbol{P}_{0}}\]

3. Determine the element deformations \(\boldsymbol{V}\) from the kinematic relations \[\boxed{ \boldsymbol{V}= \mathbf{A}_f\mspace{1.5mu}\boldsymbol{U}_f}\]

4. Determine the basic forces \(\boldsymbol{Q}\) from the collection of element force-deformation relations \[\boxed{ \boldsymbol{Q}= \mathbf{K}_s\mspace{1.5mu}\boldsymbol{V}+ \boldsymbol{Q}_{0}}\]

[@lst:DMsteps] shows the implementation of these steps in the MATLAB  script file S_DisplMethod.

% form kinematic matrix A for all dofs
A  = A_matrix (Model);
% get element cell array for continuous element deformations
iced = H_index (Model,ElemData);
% concatenate indices to single index vector ic
ic = [iced{:}];
% form Af matrix for free dofs without rows with release deformations
Af = A(ic,1:Model.nf);
% set up Ks matrix and Q0 vector
Ks = Ks_matrix (Model,ElemData);
Q0 = Q0_vector (Model,ElemData);

% set up stiffness matrix Kf and initial force vector P0
Kf = Af'*Ks*Af;
P0 = Pwf + Af'*Q0;
% solve equilibrium equations for free dof displacements Uf
Uf = Kf\(Pf-P0);
% determine continuous element deformations V
V  = Af*Uf;
% determine basic forces Q
Q  = Ks*V + Q0;
% complete Q vector with release values and return Ve
[Q,Ve] = Complete_QV(Model,ElemData,Q);

As is the case with the script for the force method of analysis, the script S_DisplMethod depends on the definition of the structural model geometry in the data structure Model , and on the specification of the element property and loading information in the data structure ElemData . The specification of the nodal forces Pf and Pwf is optional with default values equal to zero. The script uses the same auxiliary functions as the script for the force method to remove the rows of the kinematic matrix \(\mathbf{A}\) with release deformations and pad the vector of basic forces \(\boldsymbol{Q}\) from the displacement method of analysis with zero values at the releases, so that it can determine the element deformations \(\boldsymbol{V}_{\varepsilon}\) at the element ends with flexural releases and plot the deformed shape of the corresponding element, if required.

Direct Stiffness Assembly

Fig. 9.6.3 of the reader demonstrates the significant advantage of the displacement method over the force method of analysis: the direct assembly of the structure stiffness matrix \(\mathbf{K}_f\) and the resisting forces \(\boldsymbol{P}_{r}\). We demonstrate briefly the ease of implementing the direct stiffness assembly in FEDEASLab by taking advantage of the array indexing capabilities of MATLAB  described in Appendix [C:ArrayIndx] and avoiding the multiplication by the Boolean array in Eq. 9.31 of the reader. The latter approach is symbolically compact but computationally inefficient, because the Boolean matrix \(\mathbf{A}_b^{(el)}\) in Eq. 9.31 of the reader contains only a few terms of unit value but mostly zeros.

In the direct assembly process the element \(\boldsymbol{id}\) array serves as index of the element stiffness coefficients into the structure stiffness matrix \(\breve{\mathbf{K}}\). The FEDEASLab  function Kf_matrix.m assembles the complete stiffness matrix \(\breve{\mathbf{K}}\) and then extracts the stiffness \(\mathbf{K}\) at the free dofs (in practice only the free dof stiffness \(\mathbf{K}\) is assembled to save storage).

function Kf = Kf_matrix (Model,ElemData)

K  = zeros(nt,nt);                     
for el=1:ne
   % locate element in Model and return end coordinates and id array
   [xyz,id] = Localize(Model,el);
   % form element stiffness matrix ke in global reference system
   ke = ke_matrix(Model.ElemName{el},ElemData{el},xyz);
   % assemble element stiffness matrix ke into structure matrix K 
   K (id,id) = K (id,id) + ke;
end
Kf = K(1:nf,1:nf); % extract stiffness matrix of free dof's

[@lst:DSImpl] shows the function excerpt with the operations of direct stiffness assembly. K corresponds to \(\breve{\mathbf{K}}\) and Kf corresponds to \(\mathbf{K}\) in the code. The function Kf_matrix.m calls each element by its name, as specified in the field ElemName of Model, gives it its end coordinates and element properties and expects its stiffness matrix in global coordinates as return argument. This operation should be implemented in a separate function for each element, so that new elements can be added with ease. Because we plan to pursue a different implementation in FEDEASLab  as will be discussed in the course of Nonlinear Analysis, we include only the truss and the 2d frame element stiffness determination in a subfunction ke_matrix of Kf_matrix.m.

Eq. 9.27 of the reader gives the element stiffness matrix of the truss and the 2d frame element in global coordinates. The MATLAB  commands are excerpted in [@lst:kaImpl] noting that the stiffness determination of the 2d frame element with basic force releases follows the description in Eq. 9.7 of the reader.

% determine element length and orientation (direction cosines of x-axis)
[L,dcx] = ElmLenOr(xyz);

% extract element properties from ElemData
E  = ElemData.E;
A  = ElemData.A;
EA = E * A;

switch ElemName
   case 'LinTruss'
      % linear truss element
      % transformation matrix from basic to global
      ag = [-dcx' dcx'];
      k  = EA/L;
      ke = ag'*k*ag;
   case 'Lin2dFrm'
      % 2d linear frame element with or w/o release
      % transformation matrix from basic to global
      dXL = dcx(1);  dYL = dcx(2);
      ag  = [-dXL     -dYL     0   dXL      dYL      0;
             -dYL/L    dXL/L   1   dYL/L   -dXL/L    0;
             -dYL/L    dXL/L   0   dYL/L   -dXL/L    1];
      I  = ElemData.I;
      EI = E * I;
      k  = [ EA/L       0       0;
               0    4*EI/L  2*EI/L;
               0    2*EI/L  4*EI/L];
      % introduce release indices MR
      MR = zeros(3,1);
      if isfield(ElemData,'Release'), MR(ElemData.Release==1) = 1; end
      % compatibility matrix in the presence of moment release(s)
      ah = [ 1-MR(1)       0                    0;
              0            1-MR(2)        -0.5*(1-MR(3))*MR(2);
              0      -0.5*(1-MR(2))*MR(3)       1-MR(3)        ];
      % transform basic stiffness matrix for release(s)
      k  = ah'*k*ah;
      % transform basic stiffness to global reference
      ke = ag'*k*ag;

Plastic Analysis

According to Chap. 11 of the reader the determination of the collapse load factor of a structural model is possible either with the lower bound theorem of plastic analysis that makes use of the equilibrium equations and the plastic conditions, or with the upper bound theorem of plastic analysis that makes use of the kinematic relations and the plastic work increment of structural mechanisms. Even though the uniqueness theorem of plastic analysis states that the collapse load factor is the same with these methods, FEDEASLab uses two separate functions to better illustrate the resulting information from each method, which is of equilibrium nature for the lower bound theorem and of kinematic nature for the upper bound theorem. The basic forces at incipient collapse from the lower bound theorem can be used to display the moment distribution at incipient collapse, while the displacement and deformation increments from the upper bound theorem can be used to display the mechanism after the onset of collapse.

Lower Bound Theorem

The formulation of the lower bound theorem of plastic analysis in Eq. 11.13 of the reader of Chap. 11 of the reader is a linear programming problem. Functions for the solution of such problems are available in the Optimization toolbox of MATLAB . In particular, the linprog function of the toolbox is suited for structural problems involving linear equalities and inequalities. To set up the input arguments for this function the FEDEASLab function  takes the equilibrium matrix Bf, the array of plastic capacities Qpl, and the reference force vector Pref at the free DOFs of the structural model as input arguments. The latter is set up directly or with the help of the function Create_NodalForces , as discussed in [@sec:SpecNF]. Optional input arguments to the function  are the constant nodal forces \(\boldsymbol{P}_{cf}\) and the data structure Options with parameters for the linear programming solution. The default value for the constant nodal forces is zero, while the default parameters for the linear programming function that MATLAB  selects work in most cases. They need to be modified only if numerical difficulties arise.

% static matrix Bf for column-girder (Ex 12.2)
Bf = [1/6  1/6    0     0     0 ;
      0     1     1     0     0 ;
      0     0   -1/4  -1/4   1/4;
      0     0     0     1     1 ];  
% specify plastic capacities in vector Qpl
Qpl  = [160 160 120 120 120]';
% specify reference load in vector Pref
Pref = [20 0 -20 0]';

[lambdac,Qc]  = PlasticAnalysis_wLBT (Bf,Qpl,Pref);

%% Script for Example 12.2 in First Course on Matrix Structural Analysis
%  plastic analysis of column-girder assembly

%% clear memory and define global variables

With this information the function  determines the collapse load factor lambdac and the corresponding basic forces Qc of the structure under the collapse load. The syntax of the plastic analysis function is

function [lambdac, Qc] = PlasticAnalysis_wLBT(Bf, Qpl, Pref, Pcf, Options)

[@lst:PlAnal] shows the input arguments to  for the column-girder assembly of Ex. 11.2 of the reader. The output arguments are the collapse load factor lambdac and the corresponding basic forces Qc of the structural model under the collapse load.

Upper Bound Theorem

The formulation of the upper bound theorem of plastic analysis in Eq. 11.53 of the reader, Eq. 11.54 of the reader, and Eq. 11.55 of the reader of Chap. 11 of the reader is also a linear programming problem, in fact it is the dual problem of the lower bound theorem in Eq. 11.13 of the reader [@Strang_1988].

The FEDEASLab function [func:PlasticAnalysis_wUBT] sets up the necessary equality and inequality constraints with input arguments \(\mathbf{A}_f\), \(\boldsymbol{Q}_{pl}\), \(\boldsymbol{P}_{ref}\). It invokes the simplex algorithm for the solution of the linear programming problem and returns the collapse load factor \(\lambda_c\), and the displacement rates \(\dot{\boldsymbol{U}}_f\) as displacement increments DUf and the plastic hinge deformation rates \(\dot{\boldsymbol{V}}_{hp}\) as plastic deformation increments DVhp of the collapse mechanism. The syntax is

%% Script for 
%  plastic analysis of one story portal frame

%  =========================================================================================
%  FEDEASLab - Release 3.2, July 2012
%  Matlab Finite Elements for Design, Evaluation and Analysis of Structures
%  Professor Filip C. Filippou (filippou@ce.berkeley.edu)
%  Department of Civil and Environmental Engineering, UC Berkeley
%  Copyright(c) 1998-2012. The Regents of the University of California. All Rights Reserved.
%  =========================================================================================

%% clear memory and define global variables
CleanStart;

%% define structural model (coordinates, connectivity, boundary conditions, element types)
XYZ(1,:) = [  0   0];  % first node
XYZ(2,:) = [  0   5];  % second node, etc
XYZ(3,:) = [  4   5];  %
XYZ(4,:) = [  8   5];  % 
XYZ(5,:) = [  8   0];  % 
  
% element connectivity array
CON(1,:) = [  1   2];
CON(2,:) = [  2   3];
CON(3,:) = [  3   4];
CON(4,:) = [  4   5];

% boundary conditions (1 = restrained,  0 = free) (specify only restrained dof's)
BOUN(1,:) = [ 1 1 1];
BOUN(5,:) = [ 1 1 1];

% specify element type
ne = length(CON);                       % number of elements
[ElemName{1:ne}] = deal('2dFrm');       % 2d frame element

%% Model creation
Model = Create_SimpleModel (XYZ,CON,BOUN,ElemName);
% plot and label model for checking (optional)
Create_Window (0.80,0.80);         % open figure window
Plot_Model  (Model);               % plot model
Label_Model (Model);               % label model

%% static matrix Af for portal frame
Af = [1/5   1/5    0     0     0     0   1/5   1/5;
       0     1     1     0     0     0    0     0 ;
       0     0   -1/4  -1/4   1/4   1/4   0     0 ;
       0     0     0     1     1     0    0     0 ;
       0     0     0     0     0     1    1     0  ]';   
% specify plastic capacities in vector Qpl
Qpl  = [150 150 120 120 120 120 150 150]';
% specify reference load in vector Pref
Pref = [30 0 -50 0  0]';
%% call function for upper bound plastic analysis in FEDEASLab
[lambdac,DUf,DVhp] = PlasticAnalysis_wUBT (Af,Qpl,Pref);
diary on
disp('the collapse mechanism displacement rates are')
disp(DUf/DUf(1))
disp('the plastic deformation rates are')
disp(DVhp/DUf(1))
diary off

% expand displacements and deformations for plotting
Model.InextEList = 1:ne;
ied    = D_index(Model);
A      = A_matrix(Model);
Af     = A(:,1:Model.nf);
ic     = setdiff(1:size(A,1),[ied{:}]);
Ac     = Ac_matrix(Af,ic,[Model.DOF(2,1),Model.DOF(3,2)]);
Aftild = Af*Ac;

%% plot the collapse mode
Create_Window (0.80,0.80);
Plot_Model (Model);
PlotOpt.MAGF = 50;
Plot_Model (Model,Ac*DUf,PlotOpt);
Plot_PlasticHinges (Model,[],Ac*DUf,DVhp,PlotOpt);
print -painters -dpdf -r600 Ex12_2F1

: Upper bound plastic analysis of the portal frame from Ex. 11.4 of the reader

[@lst:UB4PorFr] shows the input arguments to [func:PlasticAnalysis_wUBT] for the portal frame of Ex. 11.4 of the reader. The function output consists of the collapse load factor lambdac, the corresponding displacement increments DUf, and the plastic hinge deformation increments DVhp of the collapse mechanism. The values match those from Ex. 11.4 of the reader.

Graphics Utility Functions

FEDEASLab  has functions for displaying on the computer screen the geometry of the model and the results of the analysis in the form of internal force distributions, deformed shape, and plastic hinge locations. The graphics utility functions are listed in [@tbl:GraphUtil] and are briefly described in the following. These functions depend on several supporting functions in the folder Utilities/, which are of no further interest for the user.

Specifically:

• The function Create_Window opens an empty window with two input arguments DX and DY for the percentage of the computer screen that the window occupies in the horizontal and in the vertical direction, respectively. Thus, the function call Create_Window(0.80,0.60) creates an empty window centered on the computer display and occupying 80% of the horizontal display and 60% of the vertical display dimension.

• The function Plot_Model displays the nodes of the structural model connected by straight lines (the element chords).

• The presence of the optional input argument Uf in any plotting function ensures that the information is displayed in the deformed configuration.

• The optional input argument Options in a plotting function is a data structure with fields specifying the attributes of graphical elements. Field descriptions of the data structure Options are given in [@tbl:OptPlotModel]-[@tbl:OptPlHng] for the different plotting functions. Because some fields appear in more than one function, the user should make sure to change the field value before invoking the corresponding function, lest the last definition of the field is used, which may not be the intention.

• The function Label_Model displays node and element numbers of the structural model as well as the coordinate axes for the geometry description.

• The functions for plotting the axial force, the bending moment and the curvature distribution support an optional numerical input argument Scale for the magnification of the distribution display without affecting the numerical values of the internal forces and deformations.

• The functions for labeling the axial force and the bending moment distribution support an optional integer input argument Digit for the number of digits after the comma to display. The default is 1.

• The functions for labeling the axial force and the bending moment distribution support an optional character input argument Units, so that the values of the corresponding internal force are scaled to the specified unit of axial force or bending moment. The function Create_Units.m generates different units of force, moment, length and area in the return argument. The user can define additional dependent units in terms of those already present in the fields of the return argument.

• Plotting and labeling can be limited to some elements of the structural model with the optional argument ElemList. The default option is to display the requested information for all elements of the structural model.

• In the absence of the input argument ElemData the corresponding plot function will display the moment distribution and the deformed shape for the homogeneous solution of the element equilibrium. The field w must be specified in ElemData for the function to display the moment distribution and the deformed shape for the complete solution that includes the particular solution under the element loading.

• The variable Ve corresponds to \(\boldsymbol{V}_{\varepsilon}\), the element deformations of continuous nature.

Examples

The examples in this section showcase the main graphics utility functions with further details about features. A separate MATLAB  script file is available for each example in the course reader with the file name starting with the prefix ChmExn corresponding to the example n in Chapter m of the reader. Salient excerpts from some of these script files are used in the following to showcase the FEDEASLab functions in the preceding sections and, in particular, the graphics utility functions.

Plotting and Labeling of Model

Once the function Create_SimpleModel creates the data structure Model   the information about the geometry of the model can be displayed with the functions Plot_Model and Label_Model. The default values for the graphic objects of these functions may result in an unsatisfactory image: disproportionately small size of labels and disproportionately large coordinate axes. Moreover, the default option for the function Plot_Model does not display node and boundary symbols in [@fig:PoorSizeMod].

Create_Window (0.50,0.80);
Plot_Model  (Model);               
Label_Model (Model);               

After turning on the fields PlNod and PlBnd in PlotOpt, adjusting the size of the node and boundary symbols with NodSF, adjusting the size of the coordinate axes with AxsSF, and adjusting the size and offset of labels with FntSz and LOfSz, respectively, the resulting display in [@fig:GoodSizeMod] is improved significantly.

Create_Window (0.50,0.80);                
PlotOpt.PlNod = 'yes';
PlotOpt.PlBnd = 'yes';
PlotOpt.NodSF = 0.6;
Plot_Model  (Model,[],PlotOpt);           
LblOpt.AxsSF = 0.5;
LblOpt.FntSF = 2;
LblOpt.LOfSF = 1.3;
Label_Model (Model,LblOpt);               

The default sizing and labeling values produce satisfactory results for the 3d truss model in [@fig:GoodSz3dTruss], as long as the user requests that the node and boundary symbols be displayed.

Plotting of Nodal Forces and Element Loading

The function Plot_NodalForces displays the nodal forces Pe before these are assigned to the free DOFs of the structural model by the function Create_NodalForces. It is expedient to prevent the display of node symbols in Plot_Model. The nodal force values are displayed by turning on the field Label in PlotOpt. The size of the arrow can be customized with the field FrcSF.

The function Plot_ElemLoading displays the element loads specified in field w and the element deformations in field e0 of the data object ElemData . The size of the arrow for the element loading can be customized with the field FrcSF and the size of the arrow tip with the field TipSF. Calling the function Plot_Model after calling Plot_ElemLoading and not before ensures that the node and boundary symbols are clearly visible in [@fig:ElemLoad].

Plotting of Releases for Internal Forces

The function Plot_Releases displays the internal force releases. End moment releases for truss elements are displayed automatically with default values when invoking the function. The size and offset of the moment release symbol can be adjusted with the fields HngSF and HOfSF.

For frame elements the display of internal force releases depends on the field Release of the data structure ElemData . The entries of the field in [@fig:MomRel] indicate the presence of an end moment release at end \(i\) of element \(c\) for the model of a three hinge portal frame.

Plotting and Labeling of Axial Force Distribution

The script excerpt in [@lst:Qplotcom] displays and labels the axial forces \(\boldsymbol{Q}\) in [@fig:Qplot] for the plane truss of Ex. 3.1 of the reader. The optional input argument ScaleN in Plot_AxialForces controls the scaling of the axial force distribution, while the optional input argument NDigt=2 in Label_AxialForces requests 2 digits after the comma for the label of the axial force value. A tensile axial force distribution is shown in blue and a compressive axial force in red.

ScaleN = 1/3;    
NDigt  = 2;
Create_Window (WinXr,WinYr);       
Plot_Model(Model);
Plot_AxialForces(Model,Q,[],ScaleN);
Label_AxialForces(Model,Q,[],NDigt);

: Command sequence for the display of the axial forces in @fig:Qplot

Plotting and Labeling of Bending Moment Distribution

: Command sequence for display of the bending moment distribution in @fig:MomPlotwow

The script excerpt in [@lst:MPlotwow] uses the basic forces \(\boldsymbol{Q}\) to display and label the bending moments in [@fig:MomPlotwow] for the beam with overhang of Ex. 3.5 of the reader. Without specification of the value for the distributed element loading in the field w of ElemData  the function "is not aware" of the presence of element loads and only displays the homogeneous solution for the bending moment distribution.

The script excerpt in [@lst:MPlotww] includes the input argument ElemData  in the call to the function Plot_2dMomntDistr, so that it can display the complete solution by superposing the particular solution for element \(b\) to the homogeneous solution of [@fig:MomPlotwow].

Plotting of Deformed Shape

The script excerpt in [@fig:DefoTrwDV] uses the free DOF displacements Uf under the initial element deformations V0 in [@lst:TrussKinem] to display the resulting deformed shape for the truss of Ex. 4.2 of the reader with the help of the function Plot_DeformedStructure. The input arguments ElemData and Veps to the function are not required for the case in hand.

[@fig:DefoTrwDV] shows the deformed shape of the truss with default values under a magnification factor MAGF of 30. The function displays the releases of internal forces by default, which may not be desirable for truss models. The function also displays the boundary symbols by default, but the default size value is disproportionately large for the model in [@fig:DefoTrwDV].

[@fig:DefoTrwCV] shows the improvement of the display for the deformed shape of the truss by suppressing the display of internal force releases and scaling the boundary symbols to a quarter of the default value.

[@fig:DefoTrCrd] shows that the deformed shape of a truss model can also be displayed with the function Plot_Model by including the free dof displacements Uf as the second input argument. This instructs the function to connect the nodes of the model in the deformed configuration by straight lines and can, therefore, be used for elements with straight edges to display the deformed shape of the structure. The default option is not to display the boundary conditions. The default values for the line color and line style reflect the color and style selection for element chords of deformed frame elements. These values can be overwritten with the field PlBnd for displaying the the boundary symbols, and with the fields LnClr and LnStl of PlotOpt for setting the line color and style for the deformed shape of the structure, respectively.

[@lst:SpecV04Frm] shows the specification of the element deformations \(\boldsymbol{V}_{\varepsilon}\) for the three hinge portal frame of Ex. 4.4 of the reader. It is important to initialize \(\boldsymbol{V}_{\varepsilon}\) to its full size. The script uses the function Localize to extract the end coordinates of each element, which are then used with the function ElmLenOr to determine the element length. The axial deformations in elements a through d are set equal to zero on account of their inextensible character. The index of continuous element deformations ic is set up with the MATLAB  function setdiff after specification of the flexural release locations at end \(i\) of element a, end \(i\) of element c, and end \(j\) of element d for the three hinge portal frame in Fig. 4.35 of the reader (a).

Veps = zeros(sum(Model.nq),1);
for el=1:Model.ne
   [xyz] = Localize (Model,el);
   L(el) = ElmLenOr (xyz);
end
Veps(2:3) = -0.00072.*[-1;1].*L(1)/2;
Veps(5:6) = -0.0006 .*[-1;1].*L(2)/2;
Veps(8:9) = -0.0006 .*[-1;1].*L(3)/2;
% ih=release index, ic=continuous deformation index
ih = [2 8 12];  
ic = setdiff(1:sum(Model.nq),ih);

The command sequence in [@lst:Uf4Frm] sets up the kinematic matrix \(\mathbf{A}\) and extracts the submatrix \(\mathbf{A}_f\) corresponding to the free DOFs of the model in Fig. 4.35 of the reader (a) without taking advantage of the inextensibility of elements a through d to reduce the size of the matrix, as was done for the hand calculations in Ex. 4.4 of the reader. The matrix \(\mathbf{A}_f\) refers to element deformations at continuous element ends. In this case the kinematic relations in Eq. 4.53 of the reader result in a system of 9 equations with 9 unknowns that gives a unique solution for the free DOF displacements under the given element deformations \(\boldsymbol{V}_{\varepsilon}\). Following the determination of the free DOF displacements \(\boldsymbol{U}_f\) the release deformations can be determined with the relation

\[ \boxed{ \boldsymbol{V}_h= \breve{\mathbf{A}}_f\mspace{1.5mu}\boldsymbol{U}_f- \boldsymbol{V}_{\varepsilon} } \]{#eq:VhVeps}

% form kinematic matrix A
A  = A_matrix(Model);
% submatrix for free dofs and continuous deformations
Af = A(ic,1:Model.nf);
% solve for free dof displacements
Uf = Af\Veps(ic);
% determine release deformation(s)
Vh = A(:,1:Model.nf)*Uf-Veps;
disp('the release deformations are')
disp(Vh(ih))

Create_Window(0.50,0.80);
PlotOpt.NodSF = 0.60;
Plot_Model(Model);
% plot element chords in deformed configuration
PlotOpt.MAGF = 60; 
Plot_Model(Model,Uf,PlotOpt);
PlotOpt.HngSF = 0.60;
ElemData{1}.Release = [0;1;0];
ElemData{3}.Release = [0;1;0];
ElemData{4}.Release = [0;0;1];
Plot_DeformedStructure (Model,ElemData,Uf,Veps,PlotOpt);

The command sequence in [@lst:DefoFrmcom] draws the model in [@fig:DefoFrm] in the original configuration in blue color. The command Plot_Model with the optional argument Uf draws the nodes and the element chords in the deformed configuration. The element chords are displayed with a black, dash-dotted line by default. Finally, the function call Plot_DeformedStructure draws the deformed shape of the model in the deformed configuration with a red solid line by default. The default line color, line style and line width of most plotting functions can be overwritten with specification of the PlotOpt fields LnClr, LnStl, and LnWth, respectively.

The presence of flexural releases at end \(i\) of element a, end \(i\) of element c, and end \(j\) of element d results in concentrated rotations at these locations. For the correct deformed shape it is imperative to specify the element deformations \(\boldsymbol{V}_{\varepsilon}\) as the fourth input argument to Plot_DeformedStructure in [@lst:DefoFrmcom]. Furthermore, the function displays the end moment releases in elements a, c and d, as long as ElemData  appears as the second input argument of the function with the release specification in the field Release, as [@lst:DefoFrmcom] shows.

If the element deformations Veps are not specified as the fourth input argument to the function Plot_DeformedStructure, as in [@lst:DefoWrong], the function draws the deformed shape in [@fig:WrongDefoFrm] under the assumption that the element tangents are continuous at the nodes. Consequently, instead of determining \(\boldsymbol{V}_h\) according to [@eq:VhVeps] and the element deformations \(\boldsymbol{V}\) from \(\boldsymbol{V}_{\varepsilon}+\boldsymbol{V}_h\), the function simply uses \(\boldsymbol{V}= \mathbf{A}_f\mspace{1.5mu}\boldsymbol{U}_f\) resulting in the wrong deformed shape in [@fig:WrongDefoFrm] with fixed ends at the two supports and a continuous tangent at the midspan of the girder.

Create_Window(0.50,0.80);
PlotOpt.NodSF = 0.60;
Plot_Model(Model,[],PlotOpt);
% plot element chords in deformed configuration
PlotOpt.MAGF = 60; 
Plot_Model(Model,Uf,PlotOpt);
Plot_DeformedStructure (Model,[],Uf,[],PlotOpt);

The specification of the element release locations in [@lst:SpecV04Frm] by setting up the global index ic into the element deformation vector \(\boldsymbol{V}\) is inconvenient for large structural models. In such case it is advisable to take advantage of the field Release in ElemData  for specifying the release locations of each element. Release is a 3-component column vector, with component 1 corresponding to the axial force, component 2 to the moment at end \(i\), and component 3 to the moment at end \(j\) of the element. A value of 1 indicates the presence of a release for the corresponding basic element force. [@lst:Icgen] shows the specification of the flexural releases at end \(i\) of element a, end \(i\) of element c, and end \(j\) of element d for the three hinge portal frame in B-F:KineFrEx2F1. This information is used by FEDEASLab function H_index to set up the collection of cell arrays iced in [@lst:Icgen] with the indices for the continuous element deformations. Concatenating these cell arrays produces the global index ic for the continuous element deformations and the use of the MATLAB  function setdiff gives the global index ih for the release locations in [@lst:Icgen].

% specify release locations in ElemData
ElemData{1}.Release = [0;1;0];
ElemData{3}.Release = [0;1;0];
ElemData{4}.Release = [0;0;1];
% cell array of element continuity indices
iced = H_index(Model,ElemData);
% concatenate to ic index for structure
ic = [iced{:}];
% ih release index for structure
ih = setdiff(1:sum(Model.nq),ic);

Linear Analysis

Following the specification of the model geometry and the generation of the data structure Model with the function Create_SimpleModel [@lst:ElProp4DM] shows the specification of the element properties for the portal frame of Ex. 9.3 of the reader in Fig. 9.23 of the reader (a). To account for the inextensible character of elements a, b, c and d a large value for the cross section area \(A\) is specified in ElemData . The value of \(A=10^6\) suffices for the purpose.

% specify element properties
for el=1:Model.ne
   ElemData{el}.E = 1000;
   ElemData{el}.A = 1e6;
   ElemData{el}.I = 60;
end
% insert release at base of element a
ElemData{1}.Release = [0;1;0];

[@lst:Load4DM] shows the determination of the equivalent nodal forces for the uniform element load in elements a and b in Fig. 9.32 of the reader. The uniform element load values are specified in ElemData  and used as input to the function Create_PwForces to generate the equivalent nodal forces at all DOFs. The equivalent nodal forces \(\boldsymbol{P}_{wf}\) at the free DOFs are extracted from Pw. The element load is defined in the local reference system with the \(x\)-axis pointing from end \(i\) to end \(j\). Consequently, the \(y\)-axis of element a points to the left and that of element b upward, with the element load \(w\) in Fig. 9.32 of the reader pointing in the negative \(y\) direction in both in agreement with the sign for the field w of ElemData  in [@lst:Load4DM].

% Load case: uniform load in elements a and b
ElemData{1}.w  = [0; -5];
ElemData{2}.w  = [0;-10];
Pw  = Create_PwForces(Model,ElemData);
Pwf = Pw(1:Model.nf);

Following the element property and loading specification, [@lst:Mdistr4DM] shows the command sequence for the displacement method of analysis of the portal frame and for plotting and labeling the bending moment distribution under the uniform element load in elements a and b with a scale factor of \(0.5\).

% displacement method of analysis
S_DisplMethod
% bending moment distribution
Create_Window(0.50,0.80);
Plot_Model (Model);
Plot_2dMomntDistr (Model,ElemData,Q,[],0.5);
Label_2dMoments(Model,Q)

[@lst:e0Prop4DM] shows the specification of the uniform thermal curvature in elements a and b for the portal frame of Ex. 9.3 of the reader. Note the clearing of the values from the preceding load case of uniform element loading. Finally, [@lst:MCurv4DM] shows the command sequence for producing the result display in [@fig:Ex8_3F3] under the uniform thermal curvature in elements a and b.

% no nodal forces, clear earlier values
clear Pwf
% initial element deformations
ElemData{1}.e0 = [0;0.002];
ElemData{2}.e0 = [0;0.002];
% clear distributed load from ElemData
ElemData{1}.w  = [0; 0];
ElemData{2}.w  = [0; 0];

% displacement method of analysis
S_DisplMethod
% bending moment distribution
Create_Window(0.50,0.80);
Plot_Model (Model);
Plot_2dMomntDistr (Model,ElemData,Q,[],0.6);
Label_2dMoments(Model,Q)

% plot curvature distribution
Create_Window(0.50,0.80);
Plot_Model(Model);
Plot_2dCurvDistr(Model,ElemData,Q,[],0.6);

% plot deformed shape
Create_Window(0.50,0.80);
PlotOpt.MAGF = 50;               
Plot_Model(Model);       
% plot element chords in deformed geometry
Plot_Model (Model,Uf,PlotOpt);
% plot deformed shape
Plot_DeformedStructure(Model,ElemData,Uf,[],PlotOpt);

1. Bending moment distribution

2. Curvature distribution

3. Deformed shape

Plastic Analysis

[@sec:PlAnwFlab] discusses the use of FEDEASLab function  for determining the collapse load factor \(\lambda_c\) and the corresponding basic forces \(\boldsymbol{Q}_c\) of a structural model under a given reference load vector \(\boldsymbol{P}_{ref}\) at the free DOFs. The input arguments to the function are the equilibrium matrix Bf at the free DOFs, the vector of plastic capacities Qpl and the reference load vector Pref. The script in [@lst:PlAnal] shows the set-up of the input arguments to the function  for the column-girder assembly of Ex. 11.2 of the reader.

Model.Qmis{1} = 1;
Model.Qmis{2} = 1;
Model.Qmis{3} = [1 3];

% open new window
Create_Window(0.50,0.80);
Plot_Model (Model);
Plot_2dMomntDistr (Model,[],Qc,[],0.4);
PlotOpt.HngSF = 0.6;
Plot_PlasticHinges (Model,Qpl,[],Qc,PlotOpt)
Label_2dMoments (Model,Qc)

Because the end forces Qc in [@lst:PlAnal] do not include the axial basic forces nor the end moment at the roller support, it is not possible to use the post-processing functions of FEDEASLab to display the moment distribution and the plastic hinge locations at collapse without information about the missing basic element forces. To accomplish this we append the field Qmis to the data structure Model  to supply the index of the missing basic forces. The plotting and labeling functions of FEDEASLab use the entries in this field, if present, to adjust the incoming basic force vector \(\boldsymbol{Q}\), so that its terms end up in the right location of the full basic force vector \(\boldsymbol{Q}\) before performing the display and labeling operations.

[@lst:PlotColGirdCom] shows the specification of the element index Qmis for the missing basic element forces and the command sequence for plotting the moment distribution and the plastic hinge locations at collapse for the girder-column assembly of Ex. 11.2 of the reader in [@fig:PlotColGird]. Because the plastic flexural capacity of elements b and c is the same, a plastic hinge is displayed at end \(j\) of element b as well as at end \(i\) of element c indicating, in fact, that a plastic hinge forms at the girder midspan. To avoid the display of a double hinge the plastic capacity of element b should be slightly different from element c.

The last example deals with the determination and display of the collapse mechanism of the portal frame in Ex. 11.4 of the reader. In [@lst:UB4PorFr] the collapse load factor and the kinematics of the collapse mechanism were determined by specifying explicitly the kinematic matrix \(\mathbf{A}_f\) of the portal frame and invoking the FEDEASLab function for the upper bound theorem. In [@lst:UB4PorFrB] we let FEDEASLab set up the kinematic matrix with the function A_matrix  based on the geometric information of the structural model in data structure Model . Because this function includes also the axial deformations of the frame elements in the entries of the kinematic matrix \(\mathbf{A}_f\) it is important to specify a plastic capacity vector of equal length with the number of rows of \(\mathbf{A}_f\). For this reason the plastic axial capacity is included in the plastic capacity vector Qpl in [@lst:UB4PorFrB]. It is set to a relatively high value to exclude the possibility of an axial plastic hinge forming in the portal frame.

[@lst:UB4PorFrB] shows the required sequence of commands for plotting the collapse mechanism of the portal frame with the plastic hinge locations in [@fig:Mech4PorFrm]. After creating a new window and displaying the structural model in the undeformed configuration the magnification factor MAGF is set and the Plot_Model function is invoked a second time with an additional argument the displacement increments DUf of the collapse mechanism. Finally the function Plot_PlasticHinges with the optional argument DUf displays the plastic hinge locations on the collapse mechanism.

We note that the third input argument of the Plot_PlasticHinges function can be either the basic forces Qc at incipient collapse, as is the case in [@lst:PlotColGirdCom], in which case the plastic capacity at the corresponding location needs to be specified as the second argument to the function, or the plastic hinge deformation increments DVhp, as shown in [@lst:UB4PorFrC].

Summary of Basic FEDEASLab Functions and Scripts

The basic FEDEASLab  functions and script files are summarized below.

Model and Loading Functions

Functions for Structural Matrices and Vectors

Analysis Scripts and Functions

General Utility Functions

Graphics Utility Functions

Auxiliary Functions