function [stiffness,mass]=formStiffnessMassKirchhoff(GDof,numberElements,elementNodes,Length,Width,...
Thickness,rho,numberElementsX,numberElementsY,D,Poisson)
% This function generates the submatrix for a horizontal plate
% k is the returned stiffness matrix
% m is the returned mass matrix
% E = Young's modulus
% D = E*(Thickness^3) / (12*(1-Poisson*Poisson));
% 2a = length of plate between nodes 1 and 2, or 3 and 4
% 2b = length of plate between nodes 1 and 4, or 2 and 3
a = Length/(2*numberElementsX); b = Width/(2*numberElementsY);
% stiffness matrix
K1=[7/10/b/a*D-1/5/b/a*D*Poisson+b/a^3*D+1/b^3*D*a , zeros(1,11)];
K1=[K1;1/10/a*D+2/5/a*D*Poisson+1/b^2*D*a , 4/15*b/a*D-4/15*b/a*D*Poisson+4/3/b*a*D, zeros(1,10)];
K1=[K1;-D*b/a^2-2/5/b*D*Poisson-1/10/b*D , -D*Poisson , 4/3*b/a*D+4/15/b*a*D-4/15/b*a*D*Poisson, zeros(1,9)];
K1=[K1;-7/10/b/a*D+1/5/b/a*D*Poisson-b/a^3*D+1/2/b^3*D*a, -1/10/a*D-2/5/a*D*Poisson+1/2/b^2*D*a , D*b/a^2+1/10/b*D-1/10/b*D*Poisson , 7/10/b/a*D-1/5/b/a*D*Poisson+b/a^3*D+1/b^3*D*a,zeros(1,8)];
K1=[K1;-1/10/a*D-2/5/a*D*Poisson+1/2/b^2*D*a , -4/15*b/a*D+4/15*b/a*D*Poisson+2/3/b*a*D , 0, 1/10/a*D+2/5/a*D*Poisson+1/b^2*D*a, 4/15*b/a*D-4/15*b/a*D*Poisson+4/3/b*a*D, zeros(1,7)];
K1=[K1; -D*b/a^2-1/10/b*D+1/10/b*D*Poisson, 0,2/3*b/a*D-1/15/b*a*D+1/15/b*a*D*Poisson, 2/5/b*D*Poisson+D*b/a^2+1/10/b*D, D*Poisson, 4/3*b/a*D+4/15/b*a*D-4/15/b*a*D*Poisson,zeros(1,6)];
K1=[K1;7/10/b/a*D-1/5/b/a*D*Poisson-1/2*b/a^3*D-1/2/b^3*D*a,1/10/a*D-1/10/a*D*Poisson-1/2/b^2*D*a,1/2*D*b/a^2-1/10/b*D+1/10/b*D*Poisson, -7/10/b/a*D+1/5/b/a*D*Poisson+1/2*b/a^3*D-1/b^3*D*a,...
-1/10/a*D+1/10/a*D*Poisson-1/b^2*D*a, 1/2*D*b/a^2-2/5/b*D*Poisson-1/10/b*D,7/10/b/a*D-1/5/b/a*D*Poisson+b/a^3*D+1/b^3*D*a,zeros(1,5)];
K1=[K1;-1/10/a*D+1/10/a*D*Poisson+1/2/b^2*D*a , 1/15*b/a*D-1/15*b/a*D*Poisson+1/3/b*a*D, 0, 1/10/a*D-1/10/a*D*Poisson+1/b^2*D*a , -1/15*b/a*D+1/15*b/a*D*Poisson+2/3/b*a*D , 0, ...
-1/10/a*D-2/5/a*D*Poisson-1/b^2*D*a, 4/15*b/a*D-4/15*b/a*D*Poisson+4/3/b*a*D, zeros(1,4)];
K1=[K1;-1/2*D*b/a^2+1/10/b*D-1/10/b*D*Poisson , 0, 1/3*b/a*D+1/15/b*a*D-1/15/b*a*D*Poisson, 1/2*D*b/a^2-2/5/b*D*Poisson-1/10/b*D, 0, 2/3*b/a*D-4/15/b*a*D+4/15/b*a*D*Poisson,...
2/5/b*D*Poisson+D*b/a^2+1/10/b*D, -D*Poisson, 4/3*b/a*D+4/15/b*a*D-4/15/b*a*D*Poisson,zeros(1,3)];
K1=[K1; -7/10/b/a*D+1/5/b/a*D*Poisson+1/2*b/a^3*D-1/b^3*D*a, -1/10/a*D+1/10/a*D*Poisson-1/b^2*D*a, 2/5/b*D*Poisson-1/2*D*b/a^2+1/10/b*D, 7/10/b/a*D-1/5/b/a*D*Poisson-1/2*b/a^3*D-1/2/b^3*D*a,...
1/10/a*D-1/10/a*D*Poisson-1/2/b^2*D*a, -1/2*D*b/a^2+1/10/b*D-1/10/b*D*Poisson, -7/10/b/a*D+1/5/b/a*D*Poisson-b/a^3*D+1/2/b^3*D*a, 1/10/a*D+2/5/a*D*Poisson-1/2/b^2*D*a,...
-D*b/a^2-1/10/b*D+1/10/b*D*Poisson, 7/10/b/a*D-1/5/b/a*D*Poisson+b/a^3*D+1/b^3*D*a,zeros(1,2)];
K1=[K1;1/10/a*D-1/10/a*D*Poisson+1/b^2*D*a, -1/15*b/a*D+1/15*b/a*D*Poisson+2/3/b*a*D, 0, -1/10/a*D+1/10/a*D*Poisson+1/2/b^2*D*a, 1/15*b/a*D-1/15*b/a*D*Poisson+1/3/b*a*D, 0, ...
1/10/a*D+2/5/a*D*Poisson-1/2/b^2*D*a, -4/15*b/a*D+4/15*b/a*D*Poisson+2/3/b*a*D, 0,-1/10/a*D-2/5/a*D*Poisson-1/b^2*D*a, 4/15*b/a*D-4/15*b/a*D*Poisson+4/3/b*a*D,0];
K1=[K1;2/5/b*D*Poisson-1/2*D*b/a^2+1/10/b*D, 0, 2/3*b/a*D-4/15/b*a*D+4/15/b*a*D*Poisson, 1/2*D*b/a^2-1/10/b*D+1/10/b*D*Poisson, 0, 1/3*b/a*D+1/15/b*a*D-1/15/b*a*D*Poisson,...
D*b/a^2+1/10/b*D-1/10/b*D*Poisson, 0, 2/3*b/a*D-1/15/b*a*D+1/15/b*a*D*Poisson, -D*b/a^2-2/5/b*D*Poisson-1/10/b*D, D*Poisson, 4/3*b/a*D+4/15/b*a*D-4/15/b*a*D*Poisson];
stiffness_e = K1 + K1.' - diag(diag(K1));
M1=[3454,zeros(1,11)];
M1=[M1;922*b,320*b*b,zeros(1,10)];
M1=[M1;-922*a,-252*a*b,320*a*a,zeros(1,9)];
M1=[M1;1226,398*b,-548*a,3454,zeros(1,8)];
M1=[M1;398*b,160*b*b,-168*a*b,922*b,320*b*b,zeros(1,7)];
M1=[M1;548*a,168*a*b,-240*a*a,922*a,252*a*b,320*a*a,zeros(1,6)];
M1=[M1;394,232*b,-232*a,1226,548*b,398*a,3454,zeros(1,5)];
M1=[M1;-232*b,-120*b*b,112*a*b,-548*b,-240*b*b,-168*a*b,-922*b,320*b*b,zeros(1,4)];
M1=[M1;232*a,112*a*b,-120*a*a,398*a,168*a*b,160*a*a,922*a,-252*a*b,320*a*a,zeros(1,3)];
M1=[M1;1226,548*b,-398*a,394,232*b,232*a,1226,-398*b,548*a,3454,zeros(1,2)];
M1=[M1;-548*b,-240*b*b,168*a*b,-232*b,-120*b*b,-112*a*b,-398*b,160*b*b,-168*a*b,-922*b,320*b*b,0];
M1=[M1;-398*a,-168*a*b,160*a*a,-232*a,-112*a*b,-120*a*a,-548*a,168*a*b,-240*a*a,-922*a,252*a*b,320*a*a];
mass_e = (M1 + M1.' - diag(diag(M1)))*a*b*Thickness*rho/6300;
stiffness=zeros(GDof,GDof);
mass=zeros(GDof,GDof);
for i=1:numberElements
indice=elementNodes(i,:);
elementDof=[3*indice(1)-2 3*indice(1)-1 3*indice(1) 3*indice(2)-2 3*indice(2)-1 3*indice(2)...
3*indice(3)-2 3*indice(3)-1 3*indice(3) 3*indice(4)-2 3*indice(4)-1 3*indice(4)];
stiffness(elementDof,elementDof)=stiffness(elementDof,elementDof)+stiffness_e;
mass(elementDof,elementDof)=mass(elementDof,elementDof)+mass_e;
end