%% Program to mobdify global stiffness matrix to account for homogenous and
% non-homogenous boundary conditions
clc
clear

%% Initialization [PRE-PROCESSING]

% nn: Number of nodes
% ne: Number of Elements
% K: Global stiffness matrix
% k: Local stiffness matrix
% P: Global force matrix
% nc: Nodal matrix (Information on nodes correspinding to each elements)
K=[5.5 -4.6 0 0 0;
    -4.6 14 -4.6 0 0;
     0.4 -4.6 11 -4.6 0.4;
     0 0 -4.6 14.2 -4.6;
     0 0 4.0 -4.6 6.5];
          % Enter global matrix

P=[50 200 100 200 90]; % Enter global force matrix

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

nn=5; % Enter number of nodes
ne=4; % Enter number of elements

%% Initialization- Boundary conditions [PRE-PROCESSING]

loc=[1 5]'; % Enter location of boundary nodes to be initialized as node numbers
% In this case boundaries corresponding to nodes 1 and 5

bc=[15 4]'; % Enter boundary condition displacement values in the order of
% node numbers initialized in the location matrix


% NOTE: 'loc' matrix and 'bc' matrix should be of the same size

%% [PROCESSING]

% For the algorithm to be used for general cases it is required to keep the
% size of K and P to remain the same even after the boundary conditions are
% accounted for.

c=length(bc);

%Step 1
for i=1:c
    aa=K(loc(i),loc(i));
    P(loc(i))=aa*bc(i);
    for j=1: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:nn
        if j~=loc(i)
            sum(j)=K(j,loc(i))*bc(i);
            P(j)=P(j)-sum(j);
            K(j,loc(i))=0 %This line will make 1st
            %and 5 th column=0 (except 1,1 and 5,5 element in our case)
        else
        end
    end
end

%% [POST-PROCESSING]

K % Displays modified global stiffness matrix
P=P' % Displays modified global force matrix

U=inv(K)*P % Displays calculated displacements