clear all;
clc;

syms r s     %Variables

E = 210e9;  %Young's modulus
niu = 0.25;   % Poisson's ratio

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

%Defining the interpolation function in parent coordinate system (i.e., in
%terms of r and s)
N1 = 0.25*(1-r)*(1-s);  N2 = 0.25*(1+r)*(1-s);   N3 = 0.25*(1+r)*(1+s);    N4 = 0.25*(1-r)*(1+s);

% Differential of interpolation function with respect to r and s
dN_dr = [diff(N1,r) diff(N2,r) diff(N3,r) diff(N4,r)];
dN_ds = [diff(N1,s) diff(N2,s) diff(N3,s) diff(N4,s)];

%defining the coordinates of the quadrilateral element in global or x-y
%coordinate system
x = [1 3 2.5 1.25];
y = [1 1 2 1.75];

%if you need to check if the quadrilateral element is already in square
%form in the global coordinate system
% x = [-1 1 -1 1];
% y = [-1 -1 1 1];

%Jacobian 
J(1,1) = dN_dr*x';  J(1,2) = dN_dr*y';
J(2,1) = dN_ds*x';  J(2,2) = dN_ds*y';

%G matrix - as per the lecture notes
G = (1/det(J))*[J(2,2) -J(1,2) 0 0; 0 0 -J(2,1) J(1,1); - J(2,1) J(1,1) J(2,2) -J(1,2)];

z = [0 0 0 0];   %zero matrix

P = [dN_dr z; dN_ds z; z dN_dr; z dN_ds];  %P Matrix as per the lecture notes
B = G*P;
K_i = B'*D*B;
K_int = K_i*det(J);       %The thing that we need to integrate

% intergration using Gauss Quadrature
w = [1 1];      %weight function 
rx = [-sqrt(1/3) sqrt(1/3)];    %Gauss points
K = zeros(8,8);   %Stiffness matrix - initialization 

%Integration
for i = 1:2
    for j = 1:2
        r = rx(i);
        s = rx(j);
        K = K+ eval(K_int);
    end
end

