function [k,dof,m] = triak(nodes,n1,n2,n3,th,E,nu,rho)
%function [k,dof,m] = triak(nodes,n1,n2,n3,th,E,nu,rho)
%
%  Create elem stiffness matrix etc.
%  for a 2-d triangle plane stress element
%  2 DOF per node, 6 per elem
%
%  nodes = node list
%  n1, n2, n3 = nodes for this elem
%  E,nu = elastic modulus, Poisson's ratio
%  th = thickness
%  rho = mass density (or weight) -- optional
%
%  k = element stiffness matrix
%  DB = D matrix times B matrix (for stress calc)
%  dof = dof list
%  mass matrix
%  mass = mass of element (or weight, depending on your rho)
%
error(nargchk(7,8,nargin))
%
%  Check user input
%
if(length(unique([n1 n2 n3])) < 1)
    error(sprintf('Dup node for this element: %d %d %d',n1,n2,n3));
end;
%
%  Corner node coordinates
%
x = nodes([n1 n2 n3],1);
y = nodes([n1 n2 n3],2);
i=1; j=2; m=3;  %% Text p. 272
A = det([1 x(i) y(i)
         1 x(j) y(j)
	     1 x(m) y(m)])/2;
if abs(A) < eps
  fprintf('Element connecting nodes %d, %d, %d has zero area\n',n1,n2,n3);
  error('Cannot proceed with stiffness calc.');
end;
%
%  Stuff for B matrix
%
alpha(i) = x(j)*y(m) - y(j)*x(m);
alpha(j) = x(m)*y(i) - y(m)*x(i);
alpha(m) = x(i)*y(j) - y(i)*x(j);
beta(i)  = y(j) - y(m);
beta(j)  = y(m) - y(i);
beta(m)  = y(i) - y(j);
gamma(i)  = x(m) - x(j);
gamma(j)  = x(i) - x(m);
gamma(m)  = x(j) - x(i);
%
%  B matrix p. 277
%
Bi = [beta(i)  0
      0        gamma(i)
      gamma(i) beta(i)];
Bj = [beta(j)  0
      0        gamma(j)
      gamma(j) beta(j)];
Bm = [beta(m)  0
      0        gamma(m)
      gamma(m) beta(m)];
B = [Bi Bj Bm]/(2*A);
%
%  Watch out for reverse numbering
%
A = abs(A);
%
%  D matrix p. 269
%
D = (E/(1-nu^2)) * [ 1  nu   0
                     nu 1    0
        		     0  0    (1-nu)/2];
%
%  elem stiffness matrix  eq. 6.2.52
%
k = th*A*B'*D*B;
%
%  Give 'em dof & we're done
%
dof = [2*n1-1 2*n1 2*n2-1 2*n2 2*n3-1 2*n3];
%
%  Almost, mass too
%
if nargin == 8
    %
    % elem mass matrix, eq. 16.7.9
    %
    m = zeros(6,6);
    for i=1:6
        m(i,i:2:6)=1;
    end;
    m = m+m';
    mass = rho * A * th;
    m = m*mass/12;
end;