function sig = trias(nodes,n1,n2,n3,th,E,nu,d)
%function sig = trias(nodes,n1,n2,n3,th,E,nu,d)
%
%  Inputs
%     nodes -- node list
%     n1, n2, n3 -- nodes for this tria
%     th -- thickness
%     E, nu -- elastic properties
%     d -- six displacement values, global coord.
%
%   Output
%     sig -- vector of stresses sig_x sig_y tau_xy
%
error(nargchk(8,8,nargin));
%
%  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;
%
%  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];
sig = D*B*d;   %% eq. 6.5.28;