%% Program to find (a) displacements, (b) axial force, (c) shear force and
% (d) bending moment of 2D portal frame
clc
clear

%% Initializations [PRE-PROCESSING]

% Geometry
nn=4; % Enter number of nodes
ne=3; % Enter number of elements
nc=[1 2;2 3;3 4]; % Enter element connectivity with nodes. Each row corresponds to each element.
x=[0 0 0 3 3 3]'; % Enter x-coordinates in sequence of nodes as given in 'nc' matrix
y=[0 3 3 3 3 0]'; % Enter y-coordinates in sequence of nodes as given in 'nc' matrix
L=[3 3 3]; % Enter length corresponding to each element

% Material Properties
E=[2.1e11 2.1e11 2.1e11]'; % Enter modulus of elasticity corresponding to each element
Alpha=[11E-6 11E-6 11E-6]'; % Enter coefficient of thermal expansion corresponding to each element

% Section Properties
A=[1e-2 1e-2 1e-2]'; % Enter cross-sectional area corresponding to each element
I=[0.8e-4 0.8e-4 0.8e-4]'; % Enter moment of inertia corresponding to each element

% Boundary Condition
loc=[1 2 3 10 11 12]'; % Enter location of boundary nodes to be initialized as node numbers
bc=[0 0 0 0 0 0]'; % Enter boundary condition displacement values in the order of
                   % node numbers initialized in the location matrix 'p'

% Loading
P=[0 0 0 10000 0 0 10000 -5000 -500 0 0 0]';% Load vector
dt=[0 0 0]'; % % Enter temperature change corresponding to each element

%% [PROCESSING]
%This loop calculates lenght of each element and sin cos of each element
K(3*nn,3*nn)=0;
F(3*nn,1)=0;
%The Loop Calculates local k & f and converts it to Global K& F
for e=1:ne
    X=(x(2*e)-x(2*e-1))^2;
    Y=(y(2*e)-y(2*e-1))^2;
    S(e)=((y(2*e)-y(2*e-1))/L(e));
    C(e)=((x(2*e)-x(2*e-1))/L(e));
    T=[C(e) S(e) 0 0 0 0;
        -S(e) C(e) 0 0 0 0;
        0 0 1 0 0 0;
        0 0 0 C(e) S(e) 0;
        0 0 0 -S(e) C(e) 0;
        0 0 0 0 0 1]
    a=(E(e)*I(e)/(L(e)^3))
    kg=a*[(A(e)*L(e)^2)/I(e) 0 0 -(A(e)*L(e)^2)/I(e) 0 0;
        0 12 6*L(e) 0 -12 6*L(e);
        0 6*L(e) 4*L(e)^2 0 -6*L(e) 2*L(e)^2;
        -(A(e)*L(e)^2)/I(e) 0 0 (A(e)*L(e)^2)/I(e) 0 0;
        0 -12 -6*L(e) 0 12 -6*L(e);
        0 6*L(e) 2*L(e)^2 0 -6*L(e) 4*L(e)^2]
    k=T'*kg*T
    f(1,1:6)=0;
    i=nc(e,1);
    j=nc(e,2);
    p=[3*i-2,3*i-1,3*i, 3*j-2,3*j-1,3*j];
    q=p;
    [c d]=size(p);
    %The local k is added to global K using following for loop
    for l=1:d
        F(p(l))=f(l)+F(p(l));
        for m=1:d
            K(p(l),q(m))=k(l,m)+K(p(l),q(m));
        end
    end
end
% Global Stiffness Matrix
K
% Global Force Vector
F
a=K;
%%The following is the matrix modification of K and R
%so that the size of the matrices remain the same but the midified based
%on the boundry conditions
%Find the Size of the matrix
[n,m]=size(a);
%constant column matrix
R=F+P;
b=R;
%given nodal values in matrix form
k=loc;
p=bc;
%size of the p or k matrix
[c,d]=size(p); %note using c will allow you to have any number of boundry conditions
%Location where the nodal values are known Put them in matrix form %note the size of bc and loc matrix must be the same!!
%size of the bc or loc matrix
%remain the same even after boundry conditions are accounted for.
%Carry out the rule 1 given in the class
for i=1:c
    aa=a(k(i),k(i));
    b(k(i))=aa*p(i);
    for j=1:m
        if j==k(i);
            a(k(i),j)=a(k(i),j);
        else
            a(k(i),j)=0;
        end
    end
end
%carry out rule 2 given in class
for i=1:c
    for j=1:m
        if j~=k(i)
            sum(j)=a(j ,k(i))*p(i);
            b(j)=b(j)-sum(j);
            a(j,k(i))=0; %This will make the first and 5th column zero
            %except i,i and j,j element
        else
        end
    end
end
km=a;
pm=b;
bb=a
cc=b

x=a\b%solution vector U values

xx=bb\cc

for e=1:ne
    X=(x(2*e)-x(2*e-1))^2;
    Y=(y(2*e)-y(2*e-1))^2;
    S(e)=((y(2*e)-y(2*e-1))/L(e));
    C(e)=((x(2*e)-x(2*e-1))/L(e));
    T=[C(e) S(e) 0 0 0 0;
        -S(e) C(e) 0 0 0 0;
        0 0 1 0 0 0;
        0 0 0 C(e) S(e) 0;
        0 0 0 -S(e) C(e) 0;
        0 0 0 0 0 1];
    a=(E(e)*I(e)/(L(e)^3));
    kg=a*[(A(e)*L(e)^2)/I(e) 0 0 -(A(e)*L(e)^2)/I(e) 0 0;
        0 12 6*L(e) 0 -12 6*L(e);
        0 6*L(e) 4*L(e)^2 0 -6*L(e) 2*L(e)^2;
        -(A(e)*L(e)^2)/I(e) 0 0 (A(e)*L(e)^2)/I(e) 0 0;
        0 -12 -6*L(e) 0 12 -6*L(e);
        0 6*L(e) 2*L(e)^2 0 -6*L(e) 4*L(e)^2];
    i=nc(e,1);
    j=nc(e,2);
    p=[xx(3*i-2),xx(3*i-1),xx(3*i), xx(3*j-2),xx(3*j-1),xx(3*j)]';
    SM(:,e)=kg*T*p;
    e;
end
x % Displacement and Slope
SM % Shear Force and Bending Moment