%% Program to find axial (a) forces, (b) stresses and (C) displacements in
% 2D truss systems
clc
clear

%% Initializations [PRE-PROCESSING]
ne=3 % Enter number of elements
nn=4 % Enter number of nodes
nc=[1 2;1 3;1 4]; % Nodes numbers corresponding to each element
dt=[0 0 0 0 0 0]'; % Change in temperature corresponding to each element
Alpha=[0 0 0]'; % Thermal expansion coefficient corresponding to each element
P=[-10000 -20000 0 0 0 0 0 0]'; % Forces corresponding to horizontal and vertical directions at each nodes
x=[0 0 0 2.12 0 3]'; % x-coordinate in the order of each node in 'nc' matrix
y=[0 3 0 2.12 0 0]'; % y-coordinate in the order of each node in 'nc' matrix
A=[4e-4 4e-4 4e-4]'; % Enter cross-sectional area corresponding to each element 'nc'
E=[2.1e11 2.1e11 2.1e11]'; % Enter modulus of elasticity corresponding to each element in 'nc'

kl=[3 4 5 6 7 8]'; % Location matrix containing d.o.f numbers of known boundary values
PP=[ 0 0 0 0 0 0]'; % Matrix containing known boundary values correpsonding to kl matrix

AEOL = zeros(1,3);


%% Finding direction cosines, length of each element [PRE-PROCESSING]
for i=1:ne
    'Following values are for element 1,2 and 3'
    X=(x(2*i)-x(2*i-1))^2;
    Y=(y(2*i)-y(2*i-1))^2;
    L(i)=sqrt(X+Y)
    S(i)=((y(2*i)-y(2*i-1))/L(i))
    C(i)=((x(2*i)-x(2*i-1))/L(i))
    AEOL(i)=A(i)*E(i)/L(i)
end

S=S'; % Arranging as coloumn matrix
C=C'; % Arranging as coloumn matrix
L=L'; % Arranging as coloumn matrix

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

%% To find global stiffness and force matrix [PROCESSING]
for e=1:ne
   
    a=(A(e)*E(e)/L(e)); % Calculates individual element stiffness correponding to index 'e'
   
    % Calculating transformed local stiffness matrix for individual element
    % correponding to index 'e'
    k=a*[C(e)^2 C(e)*S(e) -(C(e))^2 -C(e)*S(e);C(e)*S(e) S(e)^2 -C(e)*S(e) -(S(e)^2);-C(e)^2 -C(e)*S(e) C(e)^2 C(e)*S(e);-C(e)*S(e) -S(e)^2 C(e)*S(e) S(e)^2];
   
    % Calculating temperature force at each node for the element 'e'
    f=[-A(e)*E(e)*Alpha(e)*dt(e)*C(e) -A(e)*E(e)*Alpha(e)*dt(e)*S(e) A(e)*E(e)*Alpha(e)*dt(e)*C(e) A(e)*E(e)*Alpha(e)*dt(e)*S(e)];
   
    i=nc(e,1);
    j=nc(e,2);
    
    p=[2*i-1,2*i,2*j-1,2*j];
    q=p;
    [c d]=size(p);
    
    for l=1:d
        F(p(l))=f(l)+F(p(l)); % Assembling Global force matrix
        for m=1:d
            K(p(l),q(m))=k(l,m)+K(p(l),q(m)); % Assembling Global stiffness matrix
        end
    end
end

K % Display global stiffness matrix
F % Display force matrix

%% Modifying global stiffness and force matrix for boundary known boundary
% conditions keeping the size same

% Rename K as a and R as b for fun
a=K;
[n,m]=size(a)
R=P;
b=R;

[c,d]=size(PP)
for i=1:c
    kl(i);
    aa=a(kl(i),kl(i));
    b(kl(i))=aa*PP(i);
    for j=1:m
        if j==kl(i);
            a(kl(i),kl(i))=a(kl(i),kl(i));
        else
            a(kl(i),j)=0;
        end
    end
end
for i=l:c
    for j=1:c
        if j~=kl(i);
            sum(j)=a(j,kl(i))*PP(i);
            b(j)=b(j)-sum(j);
            a(j,kl(i))=0;
        else
        end
    end
end

a % Modified global stiffness matrix
b % Modified global matrix

%% To calculate force and stresses [PROCESSING]
% Following modification done to avoid MATLAB
%statement of singularity
bb=b(1:2);
aa=a(1:2,1:2);
xx=aa\bb

% U=[];
% for i = 1:length(xx)
%     U=vertcat(U,xx(i));
% end
U=[xx(1) xx(2) 0 0 0 0 0 0]';
for e=1:ne
    i=nc((e),1)
    j=nc((e),2)
    UU=[U(2*i-1) U(2*i) U(2*j-1) U(2*j)]';
    CS=[C(e) S(e) 0 0;0 0 C(e) S(e)];
    Force(:,e)=((A(e)*E(e))/L(e))*[1 -1;-1 1]*CS*UU;
    Stress=(Force/A(e));
end

%% [POST-PROCESSING]
Force
Stress