function Frame_EX
clc; clear;
% Plane frame example - Bhatti Ex 4.11 Page 274
e = 30000; a = 100; inertia = 1000; L = 15*12; q = 1/12;
nodes = [0, 0; L/sqrt(2), L/sqrt(2); L + L/sqrt(2), L/sqrt(2)];
conn=[1,2; 2,3];
lmm=[1,2,3,4,5,6; 4,5,6,7,8,9];
n=3*length(nodes);
debc=[1,2,3,7,8,9]; ebcVals=zeros(length(debc),1);
K=zeros(n); R = zeros(n,1);
% Generate equations for each element and assemble them.
for i=1
    lm=lmm(i,:);
    con=conn(i,:);
    [ke, rq] = PlaneFrameElement(e, inertia, a, 0, -q, nodes(con,:));
    K(lm, lm) = K(lm, lm) + ke;
    R(lm) = R(lm) + rq;
end
for i=2
    lm=lmm(i,:);
    con=conn(i,:);
    [ke, rq] = PlaneFrameElement(e, inertia, a, 0, 0, nodes(con,:));
    K(lm, lm) = K(lm, lm) + ke;
    R(lm) = R(lm) + rq;
end
K
R
% Nodal solution and reactions
d = NodalSoln(K, R, debc, ebcVals)
fa=[]; bma=[]; Va=[];
for i=1
    lm=lmm(i,:);
    con=conn(i,:);
    [f, bm, V]=PlaneFrameResults(e, inertia, a, 0, -q, ...
        nodes(con,:), d(lm));
    fa = [fa; f]; bma = [bma; bm]; Va = [Va; V];
end
for i=2
    lm=lmm(i,:);
    con=conn(i,:);
    [f, bm, V]=PlaneFrameResults(e, inertia, a, 0, 0, ...
        nodes(con,:), d(lm));
    fa = [fa; f]; bma = [bma; bm]; Va = [Va; V];
end
disp('axial force')
fa
disp('Bending moment')
bma
disp('Shear force')
Va


function [ke, rq] = PlaneFrameElement(modulus, inertia, A, qs, qt, coord)
% [ke, rq] = PlaneFrameElement(modulus, inertia, A, qs, qt, coord)
% Generates equations for a plane frame element
% modulus = modulus of elasticity
% inertia = moment of inertia
% A = area of cross-section
% qs = distributed load along the element axis
% qt = distributed load normal to the element axis
% coord = coordinates at the element ends

EI=modulus*inertia; EA = modulus*A;
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;

ke = [(EA*L^2*ls^2 + 12*EI*ms^2)/L^3, ((-12*EI + EA*L^2)*ls*ms)/L^3, ...
        (-6*EI*ms)/L^2, -((EA*L^2*ls^2 + 12*EI*ms^2)/L^3), ...
        ((12*EI - EA*L^2)*ls*ms)/L^3, (-6*EI*ms)/L^2;
    ((-12*EI + EA*L^2)*ls*ms)/L^3, (12*EI*ls^2 + EA*L^2*ms^2)/L^3, ...
        (6*EI*ls)/L^2, ((12*EI - EA*L^2)*ls*ms)/L^3, ...
        -((12*EI*ls^2 + EA*L^2*ms^2)/L^3), (6*EI*ls)/L^2;
    (-6*EI*ms)/L^2, (6*EI*ls)/L^2, (4*EI)/L, ...
        (6*EI*ms)/L^2, (-6*EI*ls)/L^2,(2*EI)/L;
    -((EA*L^2*ls^2 + 12*EI*ms^2)/L^3), ((12*EI -EA*L^2)*ls*ms)/L^3, ...
        (6*EI*ms)/L^2, (EA*L^2*ls^2 + 12*EI*ms^2)/L^3, ...
        ((-12*EI + EA*L^2)*ls*ms)/L^3, (6*EI*ms)/L^2;
    ((12*EI - EA*L^2)*ls*ms)/L^3, -((12*EI*ls^2 + EA*L^2*ms^2)/L^3),...
        (-6*EI*ls)/L^2, ((-12*EI + EA*L^2)*ls*ms)/L^3, ...
        (12*EI*ls^2 + EA*L^2*ms^2)/L^3, (-6*EI*ls)/L^2;
    (-6*EI*ms)/L^2, (6*EI*ls)/L^2, (2*EI)/L, (6*EI*ms)/L^2, ...
        (-6*EI*ls)/L^2, (4*EI)/L];
rq = [(L*(ls*qs - ms*qt))/2; (L*(ms*qs + ls*qt))/2;(L^2*qt)/12; 
    (L*(ls*qs - ms*qt))/2; (L*(ms*qs + ls*qt))/2; -(L^2*qt)/12];



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 [f, bm, V] = PlaneFrameResults(modulus, inertia, A, ...
    qs, qt, coord, dn)
% [f, bm, V] = PlaneFrameResults(modulus, inertia, A, ...
%    qs, qt, coord, dn)
% Generates frame element results
% modulus = modulus of elasticity
% inertia = moment of inertia
% A = area of cross-section
% qs = distributed load along the element axis
% qt = distributed load normal to the element axis
% coord = coordinates at the element ends
% dn = nodal solution
% [f, bm, V] = [axial force, bending moment, shear]

EI=modulus*inertia; EA = modulus*A;
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;

u = dn([1,4]);
v = dn([2,3,5,6]);

f=[]; bm=[]; V=[];
% Change increment to get results at more points
for s=0:L/2:L
    x = x1 + s*ls; y = y1+ s*ms;
    f = [f; [x,y,EA*(u(2)-u(1))/L]];
    dn2 = [(12*s)/L^3 - 6/L^2, (6*s)/L^2 - 4/L, 6/L^2 - (12*s)/L^3, ...
            (6*s)/L^2 - 2/L];
    bm = [bm; [x, y, EI*dn2*v+(qt*(L^2 - 6*s*L + 6*s^2))/(12)]];
    dn3 = [12/L^3, 6/L^2, -(12/L^3), 6/L^2];
    V = [V; [x, y, EI*dn3*v+((qt*(12*s - 6*L))/(12))]];
end
