# Blog Archives

## Structural Finite Element Analysis Software Installation

My engineering capstone design team and I recently completed our design and implementation of a structural finite element analysis (SFEA) tool. It allows a user to analyze any two-dimensional structure (beam, frame, truss). The user provides the geometry of the structure along with element parameters (cross-sectional area, modulus of elasticity, etc.) and the loads applied. The program will then calculate the nodal displacements, the reactions at the supports, and the forces acting on each element. The tool is built on the free, open-source platform called Scilab. My desire here is to make the code available to everyone for use or for further modification and development, and to provide instructions on its installation.

In order to run the application it will of course be necessary to have Scilab. Click here and you will be directed to the Scilab site where you can download and install the program. Once you have done this, open Scilab. This will bring up the main console. In the top menu, select ‘Applications->Scinotes’.

This will open the Scinotes editor in a new window. In this window you should paste the SFEA code from the end of this post (sorry for the length)

Once you have pasted the code into the Scinotes editor, save it as ‘SFEA.sci’ in the ‘C:\Program_Files\scilabx.x.x\modules\scinotes\macros’ folder and execute the code from the execute menu or the control button.

Now return to the main console and at the command prompt enter the function name, ‘–>SFEA’ , in all caps. When you hit enter, the program will run and open up two new windows. The left window will provide the controls for user entry while the right window will be used to plot the structure graphically.

To learn how to utilize the program, I have written a user’s guide.

If you have any questions or problems, please let me know in the comments section.

————————————————————————————————————————————————————————–

function SFEA() //======================================================================================= //DECLARATIONS OF GLOBAL VARIABLES AND CONSTANTS //======================================================================================= clearglobal(); global main_window_x; //constants global main_window_y; global main_window_w; global main_window_h; global plot_window_x; global plot_window_y; global plot_window_w; global plot_window_h; global margin_x; global margin_y; global spacing_x; global spacing_y; global text_w; global text_h; global button_w; global Truss_Switch; //variables global Beam_Switch; global Node_Entry_Switch; global Member_Entry_Switch; global Reactions_Entry_Switch; global Num_Nodes; global Node_Coord; global Axes_Size; global Node_Type; global Num_Members; global Connect_Mat; global Num_Load_Types; global Loads; global Num_Load_Combos; global Load_Factors; global Ultimate_Loads; global Local_k; global Transf_Mat; global Displacements; global Structure_Forces; //============================================================================================== //GLOBAL CONSTANT DEFINITIONS //============================================================================================== main_window_x =0; main_window_y =0; main_window_w =335; //width main_window_h =600; //height plot_window_x =main_window_x+main_window_w+18; plot_window_y =main_window_y; plot_window_w =640; plot_window_h =main_window_h; margin_x =10; margin_y =15; spacing_x =10; spacing_y =10; text_w =145; text_h =25; button_w =120; //================================================================================================== //CREATE THE PRIMARY WINDOWS AND MENUS //================================================================================================== //create the figure 'Main_Window' Main_Window=figure(1); Main_Window.Figure_name = "Structural Finite Element Analysis"; Main_Window.Position = [main_window_x,main_window_y,main_window_w,main_window_h]; Main_Window.Background = 8; //white //create the figure 'Plot_Window' Plot_Window=figure(2); Plot_Window.Figure_name = "Structure Plot" Plot_Window.Position = [plot_window_x,plot_window_y,plot_window_w,plot_window_h]; Plot_Window.Background = 8; //white //modify the menus of 'Main_Window' delmenu(Main_Window.figure_id, gettext("&File")); delmenu(Main_Window.figure_id, gettext("&Tools")); delmenu(Main_Window.figure_id, gettext("&Edit")); delmenu(Main_Window.figure_id, gettext("&?")); toolbar(Main_Window.figure_id, "off"); // create the "File" menu on the menu bar file_menu = uimenu(Main_Window,'label', '&File'); //create two submenus in the "File" menu new_submenu=uimenu(file_menu,'label', '&New'); close_submenu=uimenu(file_menu,'label', '&Close', 'callback', "Close_Program()"); //create submenus of "New" new_select1=uimenu(new_submenu,'label', '&Beam','callback',"Beam_Analysis()"); new_select2=uimenu(new_submenu,'label', '&Frame','callback',"Frame_Analysis()"); new_select3=uimenu(new_submenu,'label', '&Truss','callback',"Truss_Analysis()"); endfunction //======================================================================================================================== //MENU-ROUTINES //======================================================================================================================== function Close_Program() scf(2); close(); //close the plot window scf(1); close(); //close the main window endfunction //======================================================================================================================== function Beam_Analysis() global Truss_Switch; global Beam_Switch; global Node_Entry_Switch; global Member_Entry_Switch; global Reactions_Entry_Switch; Truss_Switch = 0; Beam_Switch = 1; Node_Entry_Switch = 0; Member_Entry_Switch = 0; Reactions_Entry_Switch = 0; Clear_Figures(); Load_Controls(); endfunction //======================================================================================================================== function Frame_Analysis() global Truss_Switch; global Beam_Switch; global Node_Entry_Switch; global Member_Entry_Switch; global Reactions_Entry_Switch; Truss_Switch = 0; Beam_Switch = 0; Node_Entry_Switch = 0; Member_Entry_Switch = 0; Reactions_Entry_Switch = 0; Clear_Figures(); Load_Controls(); endfunction //======================================================================================================================== function Truss_Analysis() global Truss_Switch; global Beam_Switch; global Node_Entry_Switch; global Member_Entry_Switch; global Reactions_Entry_Switch; Truss_Switch = 1; Beam_Switch = 0; Node_Entry_Switch = 0; Member_Entry_Switch = 0; Reactions_Entry_Switch = 0; Clear_Figures(); Load_Controls(); endfunction //======================================================================================================================== //SUB-ROUTINES //======================================================================================================================== function Clear_Figures() scf(2); clf(); scf(1); f=gcf(); C=f.children; for i=1:length(C) c=C(i); if c.type~="uimenu" delete(c); end end endfunction //======================================================================================================================== //======================================================================================================================== function Load_Controls() global main_window_x; global main_window_y; global main_window_w; global main_window_h; global plot_window_x; global plot_window_y; global plot_window_w; global plot_window_h; global margin_x; global margin_y; global spacing_x; global spacing_y; global text_w; global text_h; global button_w; //================================================================================================ scf(1); //set current figure to 'Main_Window' //================================================================================================ //frames: //inputs frame Frame1=uicontrol(gcf(),"style","frame"); Frame1.Position = [(margin_x/2),435,320,185]; Frame1.HorizontalAlignment = "center"; Frame1.BackgroundColor = [1,1,1]; //white //results frame Frame2=uicontrol(gcf(),"style","frame"); Frame2.Position = [(margin_x/2),300,320,130]; Frame2.HorizontalAlignment = "center"; Frame2.BackgroundColor = [1,1,1]; //white //================================================================================================== //1st Line: number of nodes and their coordinates //creates text in the main window Node_Text=uicontrol(gcf(),"style","text"); Node_Text.Position = [margin_x,main_window_h-margin_y,text_w,text_h]; Node_Text.BackgroundColor = [1,1,1]; //white Node_Text.String = "Number of Nodes"; //creates edit box for entering data in the main window Num_Node_Box=uicontrol(gcf(),"style","edit"); Num_Node_Box.Position = [margin_x+text_w+spacing_x,main_window_h-margin_y,text_h,text_h]; Num_Node_Box.String = "0"; Num_Node_Box.Tag = "NNB" //create a pushbutton in the main window Node_Button=uicontrol(gcf(),"style","pushbutton"); Node_Button.Position = [margin_x+text_w+2*spacing_x+text_h,main_window_h-margin_y,button_w,text_h]; Node_Button.String = "Node Coordinates"; Node_Button.Callback = "Get_Node_Coord()"; //==================================================================================================== //2nd Line: node reactions //creates text in the main window Reactions_Text=uicontrol(gcf(),"style","text"); Reactions_Text.Position = [margin_x,main_window_h-margin_y-text_h-spacing_y,text_w,text_h]; Reactions_Text.BackgroundColor = [1,1,1]; //white Reactions_Text.String = "Node Reactions"; //create a pushbutton in the main window Reactions_Button=uicontrol(gcf(),"style","pushbutton"); Reactions_Button.Position = [margin_x+text_w+2*spacing_x+text_h,main_window_h-margin_y-text_h-spacing_y,button_w,text_h]; Reactions_Button.String = "Select Reactions"; Reactions_Button.Callback = "Get_Node_Reactions()"; //==================================================================================================== //3rd Line: number of members and their connectivity //creates text in the main window Member_Text=uicontrol(gcf(),"style","text"); Member_Text.Position = [margin_x,main_window_h-margin_y-2*text_h-2*spacing_y,text_w,text_h]; Member_Text.BackgroundColor = [1,1,1]; //white Member_Text.String = "Number of Members"; //creates edit box for entering data in the main window Num_Member_Box=uicontrol(gcf(),"style","edit"); Num_Member_Box.Position = [margin_x+text_w+spacing_x,main_window_h-margin_y-2*text_h-2*spacing_y,text_h,text_h]; Num_Member_Box.String = "0"; Num_Member_Box.Tag = "NMB" //create a pushbutton in the main window Member_Button=uicontrol(gcf(),"style","pushbutton"); Member_Button.Position = [margin_x+text_w+2*spacing_x+text_h,main_window_h-margin_y-2*text_h-2*spacing_y,button_w,text_h]; Member_Button.String = "Member Data"; Member_Button.Callback = "Get_Member_Data()"; //==================================================================================================================== //4th Line: applied loads //creates text in the main window Loads_Text=uicontrol(gcf(),"style","text"); Loads_Text.Position = [margin_x,main_window_h-margin_y-3*text_h-3*spacing_y,text_w,text_h]; Loads_Text.BackgroundColor = [1,1,1]; //white Loads_Text.String = "Number of Load Types" //creates edit box for entering data in the main window Num_Loads_Box=uicontrol(gcf(),"style","edit"); Num_Loads_Box.Position = [margin_x+text_w+spacing_x,main_window_h-margin_y-3*text_h-3*spacing_y,text_h,text_h]; Num_Loads_Box.String = "0"; Num_Loads_Box.Tag = "NLTB" //create a pushbutton in the main window Loads_Button=uicontrol(gcf(),"style","pushbutton"); Loads_Button.Position = [margin_x+text_w+2*spacing_x+text_h,main_window_h-margin_y-3*text_h-3*spacing_y,button_w,text_h]; Loads_Button.String = "Loads"; Loads_Button.Callback = "Get_Loads()"; //========================================================================================================================= //5th Line: load factors and load combinations //creates text in the main window Combos_Text=uicontrol(gcf(),"style","text"); Combos_Text.Position = [margin_x,main_window_h-margin_y-4*text_h-4*spacing_y,text_w,text_h]; Combos_Text.BackgroundColor = [1,1,1]; //white Combos_Text.String = "Number of Load Combinations" //creates edit box for entering data in the main window Num_Combos_Box=uicontrol(gcf(),"style","edit"); Num_Combos_Box.Position = [margin_x+text_w+spacing_x,main_window_h-margin_y-4*text_h-4*spacing_y,text_h,text_h]; Num_Combos_Box.String = "0"; Num_Combos_Box.Tag = "NLCB" //create a pushbutton in the main window Combos_Button=uicontrol(gcf(),"style","pushbutton"); Combos_Button.Position = [margin_x+text_w+2*spacing_x+text_h,main_window_h-margin_y-4*text_h-4*spacing_y,button_w,text_h]; Combos_Button.String = "Load Factors"; Combos_Button.Callback = "Get_Load_Combos()"; //========================================================================================================================= //6th Line: find displacements //create a pushbutton in the main window Displacements_Button=uicontrol(gcf(),"style","pushbutton"); Displacements_Button.Position = [margin_x+10*spacing_x,main_window_h-margin_y-5.5*text_h-5.5*spacing_y,button_w,text_h]; Displacements_Button.String = "Find Displacements"; Displacements_Button.Callback = "Find_Displacements()"; //========================================================================================================================= //7th Line: find reactions //create a pushbutton in the main window Reactions_Button=uicontrol(gcf(),"style","pushbutton"); Reactions_Button.Position = [margin_x+10*spacing_x,main_window_h-margin_y-6.5*text_h-6.5*spacing_y,button_w,text_h]; Reactions_Button.String = "Find Reactions"; Reactions_Button.Callback = "Find_Reactions()"; //========================================================================================================================= //8th Line: find member forces //create a pushbutton in the main window Reactions_Button=uicontrol(gcf(),"style","pushbutton"); Reactions_Button.Position = [margin_x+10*spacing_x,main_window_h-margin_y-7.5*text_h-7.5*spacing_y,button_w,text_h]; Reactions_Button.String = "Find Member Forces"; Reactions_Button.Callback = "Find_Forces()"; endfunction //============================================================================================================================= //Data Entry and Callback Functions //============================================================================================================================= function Get_Node_Coord() //Extracts the number of nodes submitted by the user, obtains the node coordinates, and plots the nodes in the plot window global Num_Nodes; global Node_Entry_Switch; global Node_Coord; global Axes_Size; //Extract Num_Nodes from Edit Box in Main_Window Num_Nodes = evstr(get(findobj("Tag","NNB"), "String")); //Used to initialize the node coordinate matrix to zero during the first iteration //but maintain the inputted node coordinate entries during future callbacks if Node_Entry_Switch == 0 then Node_Coord=zeros(Num_Nodes,2); Node_Entry_Switch = 1; end //Let the user select the Node Coordinate data entry method Menu = list('',1,['Singularly','By Sets','From File']); Node_Entry_Method = x_choices('Select Entry Method',list(Menu)); //execute the users selection select Node_Entry_Method case 1 then //====================================================================================================== //Enter Nodes Singularly: //initializes dynamically sized matrix for coordinate labels Row_Labels=zeros(Num_Nodes,1); for i=1:Num_Nodes Row_Labels(i,1) = i; end //displays a matrix dialog box for user inputs and saves the coordinates to the matrix 'Node_Coord' Node_Coord=evstr(x_mdialog('Enter the Coordinates of the Nodes',string(Row_Labels),['X','Y'],string(Node_Coord))); case 2 then //=========================================================================================================== //Enter Nodes via Sets: //Obtain the number of sets to be entered from the user Num_Sets = evstr(x_dialog('How many sets do you want to enter?','0')); //initializes dynamically sized matrix for set labels Row_Labels=zeros(Num_Sets); for i=1:Num_Sets Row_Labels(i) = i; end Set_Data = zeros(Num_Sets,7); //initialize the Set_Data matrix //displays a matrix dialog box for user inputs of the set data Set_Data=evstr(x_mdialog('Enter the Set Data',string(Row_Labels),['Node Number Start','Node Number Increment','X Start','Y Start','X Increment','Y Increment','Number of Nodes in Set'],string(Set_Data))); //take the data entered into the 'Set_Data' matrix and determine the coordinate matrix it produces for i=1:Num_Sets //for each set for j=0:(Set_Data(i,7)-1) //for each node in the set Node_Coord(Set_Data(i,1)+j*Set_Data(i,2),1)=Set_Data(i,3)+j*Set_Data(i,5); Node_Coord(Set_Data(i,1)+j*Set_Data(i,2),2)=Set_Data(i,4)+j*Set_Data(i,6); end end case 3 then //================================================================================================================ //Enter Nodes from a File: File_Directory = uigetfile(["*.txt"],"SCI/modules/scinotes/macros/","Choose a file"); fi=mopen(File_Directory,'r'); Node_Coord=evstr(mgetl(fi,Num_Nodes)); mclose(fi); end //======================================================================================================================== //plots the nodes in the plot window Axes_Size = zeros(4); X=Node_Coord; Y=Node_Coord; X(:,2)=[]; Y(:,1)=[]; scf(2); clf(2); //determine the axes size for plotting //initial estimates Axes_Size(1)=min(Node_Coord(:,1))-((max(Node_Coord(:,1))-min(Node_Coord(:,1)))/10); Axes_Size(2)=min(Node_Coord(:,2))-((max(Node_Coord(:,2))-min(Node_Coord(:,2)))/10); Axes_Size(3)=max(Node_Coord(:,1))+((max(Node_Coord(:,1))-min(Node_Coord(:,1)))/10); Axes_Size(4)=max(Node_Coord(:,2))+((max(Node_Coord(:,2))-min(Node_Coord(:,2)))/10); //check for good proportionality width=Axes_Size(3)-Axes_Size(1); height=Axes_Size(4)-Axes_Size(2); if height==0 then aspect_ratio=3; else aspect_ratio=width/height; end //if not then adjust if aspect_ratio<0.5 then new_width=height/2; Axes_Size(1)=Axes_Size(1)-((new_width-width)/2); Axes_Size(3)=Axes_Size(3)+((new_width-width)/2); elseif aspect_ratio>2 then new_height=width/2; Axes_Size(2)=Axes_Size(2)-((new_height-height)/2); Axes_Size(4)=Axes_Size(4)+((new_height-height)/2); end //plot points plot2d(X,Y,style=-4,rect=Axes_Size); //print out the node coordinate matrix print(6,Node_Coord); endfunction //=========================================================================================================== //=========================================================================================================== function Get_Node_Reactions() //Presents a toggle dialog box for the user to select a type for each node and then //re-plots the nodes in the plot window with representative symbols global Reactions_Entry_Switch; global Num_Nodes; global Node_Type; global Node_Coord; global Axes_Size; //Used to initialize the node type matrix to 1s during the first iteration //but maintain the inputted node type entries during future iterations if Reactions_Entry_Switch == 0 then for j=1:Num_Nodes Node_Type(j) = 1; end Reactions_Entry_Switch = 1; end //generates a toggle selection dialog box for node type selection //current options are: Free, X-Roller, Y-Roller, Pin, Fixed, X-Spring, or Y-String //Free is the default value L=list(); for i=1:Num_Nodes s='Node'+string(i); l = list(s,Node_Type(i),['Free','X-Roller','Y-Roller','Pin','Fixed']); L(i)=l; end Node_Type = x_choices('Select Node Type',L); //re-plot the nodes as their respective node types X=Node_Coord; Y=Node_Coord; X(:,2)=[]; Y(:,1)=[]; scf(2); clf(2); for i=1:Num_Nodes select Node_Type(i) case 1 then s=-4; //Free case 2 then s=-3; //X-Roller case 3 then s=-9; //Y-Roller case 4 then s=-6; //Pin case 5 then s=-11; //Fixed end plot2d(X(i),Y(i),style=s,rect=Axes_Size); end //print out the node type matrix print(6,Node_Type); endfunction //============================================================================================================== //============================================================================================================== function Get_Member_Data() //Extracts the number of members from the user input, presents an x_mdialog box for the //user to input the starting and ending nodes of each member as well as the cross-secional //area, modulus of elasticity and moment of inertia, and plots the members //connecting the nodes in the plot window global Num_Members; global Member_Entry_Switch; global Connect_Mat; global Node_Coord; global Axes_Size; //Extract Num_Members from Edit Box in Main_Window Num_Members = evstr(get(findobj("Tag","NMB"), "String")); //Used to initialize the member data matrix to zero during the first iteration //but maintain the inputted member data entries during future iterations if Member_Entry_Switch == 0 then Connect_Mat=zeros(Num_Members,5); Member_Entry_Switch = 1; end //Let the user select the Member data entry method Menu = list('',1,['Singularly','By Sets','From File']); Member_Entry_Method = x_choices('Select Entry Method',list(Menu)); //execute the users selection select Member_Entry_Method case 1 then //====================================================================================================== //enter Member Data Singularly //initializes dynamically sized matrix for member labels Row_Labels=zeros(Num_Members); for i=1:Num_Members Row_Labels(i) = i; end //displays a matrix dialog box for user inputs and saves the member connections to the matrix 'Connect_Mat' Connect_Mat=evstr(x_mdialog('Enter the Members Information',string(Row_Labels),['Starting Node','Ending Node','Cross-Sectional Area (A)','Modulus of Elasticity (E)','Moment of Inertia (I)'],string(Connect_Mat))); case 2 then //========================================================================================================= //enter Member Data by sets //Obtain the number of sets to be entered from the user Num_Sets = evstr(x_dialog('How many sets do you want to enter?','0')); //initializes dynamically sized matrix for set labels Row_Labels=zeros(Num_Sets); for i=1:Num_Sets Row_Labels(i) = i; end Set_Data = zeros(Num_Sets,9); //initialize the Set_Data matrix //displays a matrix dialog box for user inputs of the set data Set_Data=evstr(x_mdialog('Enter the Set Data',string(Row_Labels),['Member Number Start','Member Number Increment','First Start Node','First End Node','Node Increment','Number of Members in Set','A','E','I'],string(Set_Data))); //take the data entered into the 'Set_Data' matrix and determine the connectivity matrix it produces for i=1:Num_Sets //for each set for j=0:(Set_Data(i,6)-1) //for each member in the set Connect_Mat(Set_Data(i,1)+j*Set_Data(i,2),1)=Set_Data(i,3)+j*Set_Data(i,5); //extracts member starting nodes Connect_Mat(Set_Data(i,1)+j*Set_Data(i,2),2)=Set_Data(i,4)+j*Set_Data(i,5); //extracts member ending nodes Connect_Mat(Set_Data(i,1)+j*Set_Data(i,2),3)=Set_Data(i,7); //extracts member cross-sectional area Connect_Mat(Set_Data(i,1)+j*Set_Data(i,2),4)=Set_Data(i,8); //extracts member modulus of elasticity Connect_Mat(Set_Data(i,1)+j*Set_Data(i,2),5)=Set_Data(i,9); //extracts member moment of inertia end end case 3 then //========================================================================================================= //enter Member Data from a file File_Directory = uigetfile(["*.txt"],"SCI/modules/scinotes/macros/","Choose a file"); fi=mopen(File_Directory,'r'); Connect_Mat=evstr(mgetl(fi,Num_Members)); mclose(fi); end //================================================================================================================= //Draws the members between the nodes with a line for i=1:Num_Members for j=1:2 X(j)=Node_Coord(Connect_Mat(i,j),1); Y(j)=Node_Coord(Connect_Mat(i,j),2); end plot2d(X,Y,style=1,rect=Axes_Size); end //print out the connectivity matrix print(6,Connect_Mat); endfunction //============================================================================================================================ //============================================================================================================================ function Get_Loads() //obtains the number of load types (i.e. dead, live, seismic, etc.) and presents the user //with an x_mdialog box to enter the horizontal and vertical loads as well as applied //moments at each node for each load type global Num_Load_Types; global Loads; global Num_Nodes; //Extract Num_Load_Types from Edit Box in Main_Window Num_Load_Types = evstr(get(findobj("Tag","NLTB"),"String")); Loads = zeros(3*Num_Nodes,Num_Load_Types); for k=1:Num_Load_Types //initializes dynamically sized matrix for loads entry and associated //labels Row_Labels=zeros(Num_Nodes); Loads_Entry=zeros(Num_Nodes,3); for i=1:Num_Nodes Row_Labels(i) = i; end //displays a matrix dialog box for user inputs Loads_Entry =evstr(x_mdialog('Enter the Loads Applied to each Node for Load Type '+string(k),string(Row_Labels),['Horizontal Load','Vertical Load','Applied Moment'],string(Loads_Entry))); for i=1:Num_Nodes //for each node for j=1:3 //for each applied load Loads(3*(i-1)+j,k)=Loads_Entry(i,j); end end end //print out the loads matrix print(6,Loads); endfunction //===================================================================================================================== //===================================================================================================================== function Get_Load_Combos() //obtains the number of load combinations from the user and provides an x_mdialog //box for the user to input the load factors for each load type for each load //combination global Num_Load_Combos; global Load_Factors; global Num_Load_Types; global Loads; global Num_Nodes; global Axes_Size; global Node_Coord; global Ultimate_Loads; global Num_Members; global Connect_Mat; //Extract Num_Load_Combos from Edit Box in Main_Window Num_Load_Combos = evstr(get(findobj("Tag","NLCB"),"String")); //initializes dynamically sized matrix for members and associated //labels Load_Factors = zeros(Num_Load_Types,Num_Load_Combos); for i=1:Num_Load_Types Row_Labels(i)= 'Load Type '+string(i); end for i=1:Num_Load_Combos Column_Labels(i) = 'Combination '+string(i); end //displays a matrix dialog box for user inputs and saves the member connections to the matrix 'Connect_Mat' Load_Factors =evstr(x_mdialog('Enter the Load Factors',string(Row_Labels),string(Column_Labels),string(Load_Factors))); //Apply the load combinations to the loads and find the ultimate load at each node Combined_Loads = Loads*Load_Factors; Ultimate_Loads = zeros(3*Num_Nodes,1); for i=1:(3*Num_Nodes) if max(Combined_Loads(i,:))==max(abs(Combined_Loads(i,:))) then Ultimate_Loads(i,1) = max(Combined_Loads(i,:)); else Ultimate_Loads(i,1) = -max(abs(Combined_Loads(i,:))); end end //define a unit length for graphing the forces L=zeros(Num_Members,1); for i=1:Num_Members Start_Node=Connect_Mat(i,1); //starting node number End_Node=Connect_Mat(i,2); //ending node number X1=Node_Coord(Start_Node,1); //starting node coordinates Y1=Node_Coord(Start_Node,2); X2=Node_Coord(End_Node,1); //ending node coordinates Y2=Node_Coord(End_Node,2); L(i)=(((Y2-Y1)^2)+((X2-X1)^2))^(1/2); //length of the member end unit_size=min(L)/8; //graph the ultimate applied forces for i=1:Num_Nodes for j=1:3 if Ultimate_Loads(3*(i-1)+j,1)~=0 then select j case 1 then //horizontal forces X(1)=Node_Coord(i,1); Y(1)=Node_Coord(i,2); X(2)=Node_Coord(i,1)+unit_size*(Ultimate_Loads(3*(i-1)+j,1)/abs(Ultimate_Loads(3*(i-1)+j,1))); Y(2)=Node_Coord(i,2); case 2 then //vertical forces X(1)=Node_Coord(i,1); Y(1)=Node_Coord(i,2); X(2)=Node_Coord(i,1); Y(2)=Node_Coord(i,2)+unit_size*(Ultimate_Loads(3*(i-1)+j,1)/abs(Ultimate_Loads(3*(i-1)+j,1))); case 3 then //moments x = Node_Coord(i,1)-unit_size/2; y = Node_Coord(i,2)+unit_size/2; w = unit_size; h = unit_size; if Ultimate_Loads(3*(i-1)+j,1)<0 then a1 = -80*64; a2 = 160*64; elseif Ultimate_Loads(3*(i-1)+j,1)>0 then a1 = 100*64; a2 = 160*64; end //plot moment arcs xarcs([x;y;w;h;a1;a2],[5]); //5 - red end //plot force lines plot2d(X,Y,style=5,rect=Axes_Size); //5 - red end end end //print out the load factors matrix and final ultimate loads vector print(6,Load_Factors); print(6,Ultimate_Loads); endfunction //===================================================================================================================== //Find Results Subroutines //===================================================================================================================== function Find_Displacements() //this function calculates the local and global stiffness matrices as well as //the transformation matrices for converting between the two. Finally, it finds //the node displacements and plots the deformed structure global Num_Nodes; global Num_Members; global Connect_Mat; global Node_Coord; global Local_k; global Transf_Mat; global Truss_Switch; global Beam_Switch; global Node_Type; global Ultimate_Loads; global Displacements; global Num_Load_Types; global Axes_Size; global Structure_Forces; //global structure stiffness matrix K=zeros(3*Num_Nodes,3*Num_Nodes); for i=1:Num_Members //obtain the parameters for the member stiffness matrices Start_Node=Connect_Mat(i,1); //starting node number End_Node=Connect_Mat(i,2); //ending node number AE=Connect_Mat(i,3)*Connect_Mat(i,4); //the product of A and E EI=Connect_Mat(i,4)*Connect_Mat(i,5); //the product of E and I X1=Node_Coord(Start_Node,1); //starting node coordinates Y1=Node_Coord(Start_Node,2); X2=Node_Coord(End_Node,1); //ending node coordinates Y2=Node_Coord(End_Node,2); L=(((Y2-Y1)^2)+((X2-X1)^2))^(1/2); //length of the member //find the member directional sines and cosines if X2==X1 then //the member is vertical if Y1<Y2 then //the member points up phi = (%pi/2); //angle of the member with the horizontal elseif Y1>Y2 then //the member points down phi = ((3*%pi)/2); end elseif Y2==Y1 then //the member is horizontal if X1<X2 then //the member points to the right phi = 0; elseif X1>X2 then //the member points to the left phi = %pi; end else //the member is at some intermediate angle phi=atan((Y2-Y1)/(X2-X1)); end C=cos(phi); S=sin(phi); //=================================================================================================== //STIFFNESS MATRIX //=================================================================================================== //compute the stiffness matrix for each member in terms of the local //coordinate axis Local_k(1:3,1:6,i)=[(AE/L),0,0,(-AE/L),0,0;0,(12*EI/L^3),(6*EI/L^2),0,(-12*EI/L^3),(6*EI/L^2);0,(6*EI/L^2),(4*EI/L),0,(-6*EI/L^2),(2*EI/L);]; Local_k(4:6,1:6,i)=[(-AE/L),0,0,(AE/L),0,0;0,(-12*EI/L^3),(-6*EI/L^2),0,(12*EI/L^3),(-6*EI/L^2);0,(6*EI/L^2),(2*EI/L),0,(-6*EI/L^2),(4*EI/L);]; //compute the transformation matrix for converting forces and displacements //from the local coordinate axis to the global structure coordinate system Transf_Mat(1:3,1:6,i)=[C,S,0,0,0,0;-S,C,0,0,0,0;0,0,1,0,0,0;]; Transf_Mat(4:6,1:6,i)=[0,0,0,C,S,0;0,0,0,-S,C,0;0,0,0,0,0,1;]; //stiffness matrix for each member in terms of global coordinate system Global_k(:,:,i)=Transf_Mat(:,:,i)'*Local_k(:,:,i)*Transf_Mat(:,:,i); //k*=T'kT //add the global stiffness matrix for each member into the global structure stiffness matrix,K K((3*Start_Node-2):(3*Start_Node),(3*Start_Node-2):(3*Start_Node))= K((3*Start_Node-2):(3*Start_Node),(3*Start_Node-2):(3*Start_Node))+Global_k(1:3,1:3,i); //adds submatrix global_k11 into K K((3*Start_Node-2):(3*Start_Node),(3*End_Node-2):(3*End_Node))= K((3*Start_Node-2):(3*Start_Node),(3*End_Node-2):(3*End_Node))+Global_k(1:3,4:6,i); //adds submatrix global_k12 into K K((3*End_Node-2):(3*End_Node),(3*Start_Node-2):(3*Start_Node))= K((3*End_Node-2):(3*End_Node),(3*Start_Node-2):(3*Start_Node))+Global_k(4:6,1:3,i); //adds submatrix global_k21 into K K((3*End_Node-2):(3*End_Node),(3*End_Node-2):(3*End_Node))= K((3*End_Node-2):(3*End_Node),(3*End_Node-2):(3*End_Node))+Global_k(4:6,4:6,i); //adds submatrix global_k22 into K //====================================================================================================== end //create an index matrix to keep up with what index of rows are left after deletion due to //boundary conditions Index=zeros(3*Num_Nodes,1); for i=1:(3*Num_Nodes) Index(i,1) = i; end Eliminate = []; Counter = 1; //if the structure is a truss there are no moments if Truss_Switch == 1 then for i = 1:Num_Nodes Eliminate(Counter)=3*i; Counter = Counter+1; end end //if the structure is a beam there are no axial forces if Beam_Switch == 1 then for i = 1:Num_Nodes Eliminate(Counter)=3*i-2; Counter = Counter+1; end end //extract boundary conditions from the Node_Type selections for i=1:Num_Nodes select Node_Type(i) case 1 then u=1;//Free v=1; theta=1; case 2 then u=0; //X-Roller v=1; theta=1; case 3 then u=1; //Y-Roller v=0; theta=1; case 4 then u=0; //Pin v=0; theta=1; case 5 then u=0; //Fixed v=0; theta=0; end //Determine where the boundary conditions give zero displacement //these rows/cols will be deleted to generate the final //invertible matrix if u==0 then Eliminate(Counter)=3*i-2; Counter=Counter+1; end if v==0 then Eliminate(Counter)=3*i-1; Counter=Counter+1; end if theta==0 then Eliminate(Counter)=3*i; Counter=Counter+1; end end //before removing rows and columns, preserve the full vectors and matrices k = K; ultimate_loads = Ultimate_Loads; //delete rows and columns corresponding to given parameters and boundary conditions k(Eliminate,:)=[]; k(:,Eliminate)=[]; Index(Eliminate,:)=[]; ultimate_loads(Eliminate,:)=[]; //calculate the unknown displacements Unknown_Displacements=inv(k)*ultimate_loads; //Develop the full displacement matrix Displacements = zeros(3*Num_Nodes,1); for i=1:size(Index,"r") //for each element in the index //copy the rows from the calculated Unknown_Displacements matrix into the correct rows //of the full Displacements matrix Displacements(Index(i,1),:) = Unknown_Displacements(i,:); end //plot the deflected shape Def_Node_Coord = zeros(Num_Nodes,2); for i=1:Num_Nodes Def_Node_Coord(i,1) = Node_Coord(i,1)+Displacements(3*(i-1)+1,1); //x-coordinates Def_Node_Coord(i,2) = Node_Coord(i,2)+Displacements(3*(i-1)+2,1); //y-coordinates end for i=1:Num_Members for j=1:2 X(j)=Def_Node_Coord(Connect_Mat(i,j),1); Y(j)=Def_Node_Coord(Connect_Mat(i,j),2); end plot2d(X,Y,style=2,rect=Axes_Size); end //Calculate the global force vector for finding the reactions Structure_Forces = K*Displacements; //print out the global structure displacements vector print(6,Displacements); endfunction //===================================================================================================== //===================================================================================================== function Find_Reactions() //this function determines the reactions at the structure supports global Num_Nodes; global Node_Type; global Structure_Forces; //Extract the reactions from the Structure_Forces matrix Reactions = zeros(3*Num_Nodes,1); for i=1:Num_Nodes if Node_Type(i)==1 then Reactions((3*i-2:3*i),1)=0; else Reactions((3*i-2:3*i),1)=Structure_Forces((3*i-2:3*i),1); end end //print out the reactions vector print(6,Reactions); endfunction //===================================================================================================== //===================================================================================================== function Find_Forces() //this function uses the global displacements to determine the forces on each //member of the structure in terms of that member's local coordinate axis global Num_Members; global Connect_Mat; global Displacements; global Local_k; global Transf_Mat; //Adjust the forces to local coordinates Local_Member_Forces = zeros(6,Num_Members); //the i-th column vector represents the forces on the i-th member //in terms of the local coordinate axis Member_Displacements=zeros(6,1); //for extracted each member's displacements in terms of global coordinates for i=1:Num_Members Start_Node=Connect_Mat(i,1); //starting node number End_Node=Connect_Mat(i,2); //ending node number //extract the member's displacements from the structure displacements matrix Member_Displacements(1:3,1)=Displacements(3*Start_Node-2:3*Start_Node); Member_Displacements(4:6,1)=Displacements(3*End_Node-2:3*End_Node); Local_Member_Forces(:,i)=Local_k(:,:,i)*Transf_Mat(:,:,i)*Member_Displacements(:,1); //f=kTD end //print out the local member forces matrix print(6,Local_Member_Forces); endfunction

## Structural Finite Element Analysis User’s Guide

Now that you have downloaded and installed Scilab, the Structural Finite Element Analysis (SFEA) program, and got it up and running, how do you use it. Well that’s what we want to go through at this point. What you should see on your desktop are two blank screens. The one on the left will provide the user input controls while the one of the right will graphically depict the structure being analyzed. If this is not the case then you should read here first.

**Getting Started: The File Menu**

In the left screen there is a ‘File’ menu. If you click it you will see that you can select ‘New’ or ‘Close’. ‘Close’ is pretty obvious. ‘New’ will give you a sub-menu of three choices: Beam, Frame, or Truss. To make clear the difference, in this program an element of a beam supports transverse and moment forces, a frame supports axial, transverse and moment forces, and a truss element supports axial and transverse forces.

To give an example of why this might be important, consider a simple cantilever beam with a non-orthogonal loading. (In other words, the load does not meet the beam’s longitudinal axis at 90 degrees) Now although it is a ‘beam’ due to its form and positioning, since the load will produce an axial force in the structure, it is not a ‘beam’ with regard to the methodology of this program. If you were to analyze this structure, you would select ‘Frame’.

Once you make this selection, the left window will populate the controls. You will notice that there are two frames. The top frame walks the user through the step-by-step acquisition of the necessary data. The bottom frame contains the controls for obtaining the results of node displacements, support reactions, and elemental or member forces.

**Entering Node Data**

Now that the controls are up, the user simply starts at the top left and goes across and down just like reading text. Start by entering the number of nodes in the text box at the top and then click the ‘Node Coordinates’ button. This brings up a dialog box that asks the user if they want to enter the nodes ‘Singularly’, ‘By Sets’, or ‘From File’. ‘Singularly’ means manually enter the node coordinates one at a time. ‘By Sets’ means the user will enter some criteria about the nodes that is repetitive and the computer will determine the nodes for you. This is a helpful function when you are entering a large number of nodes that follow some consistent pattern with regard to their geometric arrangement. ‘From File’ means the user can open and import the coordinates from a user created text file that contains all the node coordinates.

**Singularly**

If the user selects this method, a new dialog box will open to let the user input the node coordinates.

The number of rows in the dialog box will dynamically correspond to the number of nodes the user entered would be in the structure.

**By Sets**

If the user selects this method, another dialog will appear to ask the user how many sets they have to enter. This is a simple integer input.

After entering the number of sets, another dialog box will be displayed to obtain the set data from the user.

The rows in the dialog box correspond to the number of sets. For each set the user must enter the node number to start with, how the node number for each subsequent node should increment from the beginning node, the beginning x-y coordinate for the first node, the x and y distance each subsequent node will increment from the previous one, and how many nodes in the set.

For example if I was analyzing a simple cantilever beam and I wanted to divide the beam into 20 separate elements, then I would have 21 nodes. Since I would space the nodes out evenly along the length of the beam, I could enter those nodes by sets without having to enter them one at a time. If the length of the entire beam is 120 inches (10 feet), then I would have 20 elements each 6 inches long. So in entering the data by sets in the dialog above I would input: node number start = 1, node number increment = 1, x start = 0, y start = 0, x increment = 6, y increment 0, number of nodes in set = 21. The program would generate 21 nodes starting at the x-y coordinate (0,0) and at every 6 inches horizontally for the full length of the beam.

**From File**

If the user selects this method, a standard Windows dialog box will open where you can search in the file tree structure to find the text file containing the node coordinates data. The text file should be formatted with one node coordinate per line starting with the x-coordinate followed by a space followed by the y-coordinate.

Which ever method the user selects, once the data has been provided to the program, it will output the node coordinate matrix to the Scilab console. The user can verify the correct coordinates numerically there. Additionally, the nodes will be plotted graphically in the second SFEA window. This allows the user to verify the nodes visually.

**Entering Support Reactions Data**

Once the nodes have been entered, the user can set the boundary conditions for the structure by pushing the ‘Select Reactions’ button. A dialog box will then open and allow the user to select the support type at each node.

Notice the default node type is free, so the user only has to define those nodes which are supports. An ‘X-Roller’ allows rotation and movement in the y-direction but no movement in the x-direction. Likewise a ‘Y-Roller’ allows rotation and movement in the x-direction but no movement in the y-direction. A pin allows rotations but no x or y movement. Finally a fixed support allows no movement at all.

Once the user defines the node types, the node type matrix will be displayed in the Scilab console and the structure plot will be modified with different symbology for different node types as seen below. The symbols are from left to right, free, x-roller, y-roller, pin, and fixed.

**Entering Member or Element Data**

Once the nodes and reactions are defined the user can input the members that run between the nodes. In much the same way as the entry of the nodes, the user will enter the number of members in the textbox in the SFEA first window and then push the ‘Member Data’ button. As with the nodes, a dialog box will permit the user to input the data ‘Singularly’, ‘By Sets’, or ‘From File’.

**Singularly**

If the user selects this method, a dialog box will appear to allow the user to input the data for each member. The data required are the start node number, the end node number, the member’s cross-sectional area (A), the member’s modulus of elasticity (E), and the member’s moment of inertia (I).

**By Sets**

In the same way as with node entry, the ‘By Sets’ method of member entry allows the user to systematically enter a large number of members with repetitive characteristics. With regard to member entry, each member of a set must have the same A, E, and I.

**From File**

As before, this allows the user to load the member data from a text file on the computer. The format of the file should be one member per line with each line having: start node number-space-end node number-space-A-space-E-space-I

Once the member data is entered, this data is output to the Scilab console as the connectivity matrix. Also the members will be plotted in the SFEA plot window.

**Entering Applied Loads**

There are two aspects to the entry of the loads that affect the structure. To understand this fully, let’s define some of the terminology. The first piece of information to enter in the SFEA first window control is the ‘Number of Load Types’. A load type is typically something like dead loads, live loads, wind loads, seismic loads, etc. These are each different load types. If a structure was going to have each of these load types listed then it would have 4 load types. Another term in the program is ‘Number of Load Combinations’. A load combination is a linear combination of the form:

Where each ‘x’ is a load type and each ‘A’ is a load factor. A load factor is something that is typically used in Load and Resistance Factor Design (LRFD) as a form of scaling factor or safety factor. So for example we might have:

Where D stands for dead load and L stands for live load. In many structural applications and building codes there are sets of multiple load combinations and corresponding load factors. Each of these combinations must be calculated and the worst case scenario is then selected for design. The SFEA program provides the capacity for this type of analysis.

**Load Types and Loads**

Once the user enters the number of load types and pushes the ‘Loads’ button, a dialog box that allows the entry of the different loads for each type will appear.

Each row in the dialog box corresponds to each node in the structure. Consequently you can enter a horizontal, vertical, and moment load at each node. After pressing ‘OK’ another identical dialog box will appear for the next load type.

The loads matrix will be displayed in the Scilab console after entry. Each column of the matrix is a different load type. The rows follow the following pattern. Each node has 3 potential loads therefore the first 3 rows correspond to the 3 loads on node 1. The next 3 rows correspond to node 2 and so on. Therefore there should be 3N rows where N is the number of nodes. Additionally, the order of the loads in each set of 3 rows is the same as it is entered in the dialog box: horizontal, vertical, moment.

**Load Combinations and Factors**

As discussed above, the number of load combinations and load factors are often listed in manuals or code books. Once the number of combinations is entered and the user pushes the ‘Load Factors’ button, the following dialog box appears.

Here the user enters the load factors that correspond to each load type for each combination.

Once this data is entered the load factors will be displayed to the Scilab console as well as the ultimate loads. The ultimate loads are the maximum loads at each node calculated from all the given loads, load combinations, and load factors. Also, the ultimate loads will be depicted in the plot window on the structure with red lines or curves. The line goes from the node in the direction of the force. A semicircular curve around the left side of a node is a positive moment according to the right-hand rule while a curve around the right side of a node is negative.

**Find Displacements**

At this point the structure is fully defined and graphically displayed in the plot window. Now when the user pushes the ‘Find Displacements’ button, the displacements vector will be displayed in the Scilab console and the new displaced shape of the structure will be plotted in blue. The displacements vector follows the same convention as discussed above with the loads matrix. The first 3 rows correspond to the horizontal, vertical, and angular displacements of node 1. The second 3 for node 2 and so on.

**Find Reactions**

When the user presses this button next, the reactions vector is displayed in the Scilab console. Again, as with the loads and displacements you have the horizontal, vertical, and moment reactions for each node starting with node 1 and going down.

**Find Member Forces**

Finally when the user pushes this button, the local member forces matrix is displayed in the Scilab console. Each column corresponds to a member. Column 1 for member 1 and so on. The first 3 rows are the horizontal, vertical, and moment forces on the member at the starting node while rows 4,5, and 6 are the forces on the member at the ending node. It is important to understand that these forces are with reference to the member’s local axis. The local axis for each member has its origin at the starting node with the horizontal or x-axis directed toward the ending node. The vertical or y-axis is then appropriately orthogonal to the local x-axis.

I hope this user’s guide has been helpful and provided sufficient information or answers to your questions. However if you are still stuck or want further information just let me know in the comments.