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