%% Program to find (a) deflection and (b) shear force in a beam element
clc
clear

%% Initializations [PRE_PROCESSING]

nn= 3; % Enter number of nodes
ne=2; % Enter number of elements
nc=[1 2; 2 3]; % Node numbers in row-wise corresponding to each element

E=20e6*[1 1];
I=2509*[1 1];
L=500*[1 1];

P=[0 0 0 0 -2000 0]' % Vertical load and moment at each node

loc=[1 2 3 ]'; % Node numbers at which boundary conditions is specified
bc=[0 0 0]'; % Boundary values at specified node numbers in 'loc' matrix

K(2*nn,2*nn)=0; %Size pre-allocation
F(2*nn,1)=0; %Size pre-allocation

%% Following loop create global stiffness matrix (K) and global force
% matrix (F)
for i=1:ne
k(i,1)=E(i)*I(i)/L(i)^3; % Local stiffness for individual elements
end
for e=1:ne
    
    i=nc(e,1);
    j=nc(e,2);
    pp=[2*i-1,2*i,2*j-1,2*j]';
    q=pp;
    ks=k(e)*[12 6*L(e) -12 6*L(e);6*L(e) 4*L(e)^2 -6*L(e) 2*L(e)^2;
        -12 -6*L(e) 12 -6*L(e);6*L(e) 2*L(e)^2 -6*L(e) 4*L(e)^2];
    [c d]=size(pp);
    %Now the local k is added to global K using following for loop
    for l=1:c
        for m=1:c
            l,m;
            K(pp(l),q(m))=ks(l,m)+K(pp(l),q(m));
        end
    end
end

K % Displays global stiffness matrix
P % Displays global force matrix

%% Modifying global stiffness matrix (K) and global force matrix (F) for
% boundary conditions keeping the size same

c=length(bc);

% Step 1
for i=1:c
    aa=K(loc(i),loc(i));
    P(loc(i))=aa*bc(i);
    for j=1:2*nn
        if j==loc(i)
            K(loc(i),j)=K(loc(i),j);
        else
            K(loc(i),j)=0;
        end
    end
end

% Step 2
for i=1:c
    for j=1:2*nn
        if j~=loc(i)
            sum(j)=K(j,loc(i))*bc(i);
            P(j)=P(j)-sum(j);
            K(j,loc(i))=0;
        else
        end
        
    end
end
format long
K % Display modified global stiffness matrix 
P % Display modified global force matrix
x=K\P;


%% Shear force calculations
for e=1:ne
    i=nc(e,1);
    j=nc(e,2);
    pp=[2*i-1,2*i,2*j-1,2*j]';
    q=pp;
    k(e)=E(e)*I(e)/L(e)^3;
    ks=k(e)*[12 6*L(e) -12 6*L(e);6*L(e) 4*L(e)^2 -6*L(e) 2*L(e)^2;...
        -12 -6*L(e) 12 -6*L(e);6*L(e) 2*L(e)^2 -6*L(e) 4*L(e)^2];
    i=nc(e,1);
    j=nc(e,2);
    xx=[x(2*i-1) x(2*i) x(2*j-1) x(2*j)]';
    SM(:,e)=ks*xx;
end

x=K\P % Displays displacements at each node
SM % Displays shear force
