%truss Example
function truss_code
clear all; clc
% Five bar truss example
%element Modulus of Elasticity

e=1.9*10^6;
a=8;
P = 500;
nodes = [0, 0; 36, 0; 0, 36; 36, 36; 72, 36];
conn = [1,2; 2,3; 3,4; 2,4; 2,5; 4,5];
lmm = [1,2,3,4; 3, 4, 5, 6; 5, 6, 7, 8; 3, 4, 7, 8; 3, 4, 9, 10; 7, 8, 9, 10];  
K=zeros(10);
% Generate stiffness matrix for each element and assemble it.
for i=1:6
    lm=lmm(i,:);
    con=conn(i,:);
    k=PlaneTrussElement(e, a, nodes(con,:));
    K(lm, lm) = K(lm, lm) + k;
end


% Define the load vector
R = zeros(10,1); R(8)=-P; R(10)=-P;  

% Nodal solution and reactions
[d, reactions] = NodalSoln(K, R, [1,2,5,6], zeros(4,1))
results=[];
for i=1:5
    results = [results; PlaneTrussResults(e, a, ...
            nodes(conn(i,:),:), d(lmm(i,:)))];
end
format short g
results

   
function k = PlaneTrussElement(e, A, coord)
% PlaneTrussElement(e, A, coord)
% Generates stiffness matrix for a plane truss element
% e = modulus of elasticity
% A = area of cross-section
% coord = coordinates at the element ends

x1=coord(1,1); y1=coord(1,2);
x2=coord(2,1); y2=coord(2,2);
L=sqrt((x2-x1)^2+(y2-y1)^2);
ls=(x2-x1)/L; ms=(y2-y1)/L;
k = e*A/L*[ls^2, ls*ms,-ls^2,-ls*ms;
    ls*ms, ms^2,-ls*ms,-ms^2;
    -ls^2,-ls*ms,ls^2,ls*ms;
    -ls*ms,-ms^2,ls*ms,ms^2];

function [d, rf] = NodalSoln(K, R, debc, ebcVals)
% [nd, rf] = NodalSoln(K, R, debc, ebcVals)
% Computes nodal solution and reactions
% K = global coefficient matrix
% R = global right hand side vector
% debc = list of degrees of freedom with specified values
% ebcVals = specified values
dof = length(R);
df = setdiff(1:dof, debc);
Kf = K(df, df);
Rf = R(df) - K(df, debc)*ebcVals;
dfVals = Kf\Rf;
d = zeros(dof,1);
d(debc) = ebcVals;
d(df) = dfVals;
rf = K(debc,:)*d - R(debc);

function results = PlaneTrussResults(e, A, coord, disps)
% results = PlaneTrussResults(e, A, coord, disps)
% Compute plane truss element results
% e = modulus of elasticity
% A = Area of cross-section
% coord = coordinates at the element ends
% disps = displacements at element ends
% The output quantities are eps = axial strain
% sigma = axial stress and force = axial force.

x1=coord(1,1); y1=coord(1,2);
x2=coord(2,1); y2=coord(2,2);
L=sqrt((x2-x1)^2+(y2-y1)^2);
ls=(x2-x1)/L; ms=(y2-y1)/L;
T=[ls,ms,0,0; 0,0,ls,ms];
d = T*disps;
eps= (d(2)-d(1))/L;
sigma = e.*eps;
force = sigma.*A;
results=[eps, sigma, force];
