clear all; 
clc; 

%defining materials properties in SI units
E = 210e9;  %Young's modulus
niu = 0.25; %Poisson's ratio
t = 1e-4;     %thickness

%Relationship between stress and strain for plane stress condition
D =E/(1-niu^2)*[1 niu 0; niu 1 0; 0 0 (1-niu)/2];

%Geometry of entire structure; left-right and then bottom top numbering
X = [0 3 0 3];
Y = [0 0 4 4];

%Defining global stiffness and force matrices
KG = zeros(8,8);    %Stiffness matrix 
FG = zeros(8,1);    %Force matrix
FG(2) = (500*4*t)/2;    %Nodal forces on node 2 and 4
FG(4) = FG(2); 

%stiffness matrix calculation
for i = 1:2
    % Vertex of triangular element - anti-clockwise numbering
    x1 = X(i);    x2 = X(i+1);    x3 = X(i+2);
    y1 = Y(i);    y2 = Y(i+1);    y3 = Y(i+2);
    
    %determinant to be used for calculation of interpolation function
    A =0.5*(det([1 x1 y1; 1 x2 y2; 1 x3 y3]));
    
    %defining interpolation function in x,y coordinate system
    syms x y 
    N1 = (0.5/A)*((x2*y3-x3*y2)+(y2-y3)*x+(x3-x2)*y);
    N2 = (0.5/A)*((x3*y1-x1*y3)+(y3-y1)*x+(x1-x3)*y);
    N3 = (0.5/A)*((x1*y2-x2*y1)+(y1-y2)*x+(x2-x1)*y);
    
    %defining [B] matrix - refer notes for symbols
    beta = [diff(N1,x) diff(N2,x) diff(N3,x)]; 
    gamma = [diff(N1,y) diff(N2,y) diff(N3,y)];
    z = [0 0 0];
    B = [beta z; z gamma; gamma beta];
    
    %Defining elemental stiffness matrix - note that area here has to be
    %positive number
    K = (t*(B'*D*B)*abs(A));
    
    %assembly of elemental stiffness matrix into global one
    KG([i:i+2],[i:i+2]) = KG([i:i+2],[i:i+2])+ K([1:3],[1:3]);  %first 3 x 3 entries
    KG([i:i+2],[i+4:i+6]) = KG([i:i+2],[i+4:i+6])+ K([1:3],[4:6]); %first 3 rows and last 3 columns, leaving a column of zeros
    KG([i+4:i+6],[i:i+2]) = KG([i+4:i+6],[i:i+2])+ K([4:6],[1:3]); %the bottom 3 x 3 entries, leaving a row of zeros
    KG([i+4:i+6],[i+4:i+6]) = KG([i+4:i+6],[i+4:i+6])+ K([4:6],[4:6]); %the last 3 x 3 entries, leaving a column & row of zeros
end

%Solution
KA = zeros(5,5); %initialization of active component of stiffness matrix

%Populating the active component of the stiffness matrix
for i = 1:3
    KA(i,:) = [KG(2*i,2) KG(2*i,4) KG(2*i,6) KG(2*i,7) KG(2*i,8)]; %first 3 rows are 2, 4 and 6th rows of Global stifness matrix
end

for i = 4:5
    KA(i,:) = [KG(i+3,2) KG(i+3,4) KG(i+3,6) KG(i+3,7) KG(i+3,8)]; %last 2 rows are 7 and 8 rows of global stifness matrix
end

FA = zeros(5,1);    %5 known forces
FA([1:2]) = [FG(2),FG(4)];
UA = inv(KA)*FA;    %Solution of unknown displacement

%Preparing the entire displacement vector
U = zeros(8,1);
U(2) = UA(1);   U(4) = UA(2); 
for i = 6:8
    U(i) = UA(i-3);
end

%For display
for i = 1:4
    UX(i) = U(i)*1e9;
    UY(i) = U(i+4)*1e9;
end
N = [1 2 3 4];  %node numbers
disp('Displacement'); disp('    Node #   UX (nm)     UY(nm)'); disp([N',UX',UY']);
