function K = rect2d(E,nu,t,P1,P2,P3)

%   RECT2D
%
%   This command generates the 8 x 8 stiffness matrix for a two
%   dimensional rectangular plate element, subjected to in-plane
%   forces only, in global coordinates. We emphatically remind the
%   user that this is not a general quadrilateral plate element, 
%   but simply a rectangular one. The local y coordinate is assumed
%   to run along the line connecting points P1 and P2, and the local
%   x direction runs along the line connecting points P1 and P4.
%   It does not matter whether the user defines the corners of the
%   rectangle going clockwise, or counterclockwise, around the
%   element, as long as things are done consistently. The
%   formulation being used here is described in detail on
%   Przemieniecki, J.S., "Theory of Matrix Structural Analysis",
%   McGraw-Hill, 1968 (Section 5.8, pp. 89-102). It should also be
%   pointed out the element is generated using the 
%   linear-edge-displacement assumption. The syntax is:
%              K = rect2d(E,nu,t,P1,P2,P3)
%   where: E is the Young's modulus;
%          nu is the Poisson's ratio;
%          t is the thickness of the plate;
%          P1, P2 and P3 are vectors of the {x,y} coordinates of
%              three vertices of the rectangle (in the global
%              coordinate system), taken sequentially either
%              clockwise or counterclockwise around the element.

a = norm(P3-P2);                 % Find direction cosines.
lqr = (P3(1)-P2(1))/a;
mqr = (P3(2)-P2(2))/a;
b = norm(P2-P1);
lpq = (P2(1)-P1(1))/b;
mpq = (P2(2)-P1(2))/b;

lam = [lqr mqr; lpq mpq];        % Build up lambda matrix.
zz = [0 0; 0 0];
lambda = [lam zz zz zz; zz lam zz zz; zz zz lam zz; zz zz zz lam];

beta = b/a;                            % Create local matrix.
kl(1,1) = 4*beta + 2*(1-nu)/beta;
kl(2,1) = 1.5*(1+nu);
kl(2,2) = 4/beta + 2*(1-nu)*beta;
kl(3,1) = 2*beta - 2*(1-nu)/beta;
kl(3,2) = -1.5*(1-3*nu);
kl(3,3) = kl(1,1);
kl(4,1) = -kl(3,2);
kl(4,2) = -4/beta + (1-nu)*beta;
kl(4,3) = -kl(2,1);
kl(4,4) = kl(2,2);
kl(5,1) = -2*beta - (1-nu)/beta;
kl(5,2) = -kl(2,1);
kl(5,3) = -4*beta + (1-nu)/beta;
kl(5,4) = kl(3,2);
kl(5,5) = kl(1,1);
kl(6,1) = -kl(2,1);
kl(6,2) = -2/beta - (1-nu)*beta;
kl(6,3) = -kl(3,2);
kl(6,4) = 2/beta - 2*(1-nu)*beta;
kl(6,5) = kl(2,1);
kl(6,6) = kl(2,2);
kl(7,1) = kl(5,3);
kl(7,2) = -kl(3,2);
kl(7,3) = kl(5,1);
kl(7,4) = kl(2,1);
kl(7,5) = kl(3,1);
kl(7,6) = kl(3,2);
kl(7,7) = kl(1,1);
kl(8,1) = kl(3,2);
kl(8,2) = kl(6,4);
kl(8,3) = kl(2,1);
kl(8,4) = kl(6,2);
kl(8,5) = -kl(3,2);
kl(8,6) = kl(4,2);
kl(8,7) = -kl(2,1);
kl(8,8) = kl(2,2);
for i = 1:7
   for j = i+1:8
      kl(i,j) = kl(j,i);
   end
end
klocal = E*t/(12*(1-nu^2))*kl;

K = lambda'*klocal*lambda;            % Create matrix in global coordinates.


