function DB = triahf(nodes,n1,n2,n3,K,t)
%
%  function [k,fconv] = triahf(nodes,n1,n2,n3,K,t,conv);
%
%  Triangle element for thermal conductivity
%
%  Inputs:
%
%  n1, n2, n3     node numbers
%  nodes          list of all nodes (x, y)
%  K              conductivity (isotropic)
%  conv           optional vector [m n h Tinf];
%                 m,n = node numbers defining side with conv
%                 h = convection coefficient
%                 Tinf = free-stream temperature
%
%  Outputs:
%
%  k              element conductivity matrix, global coord
%  fconv          load vector for convection, if conv given
%
error(nargchk(6,6,nargin));
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;
Asign = sign(A);
A = abs(A);
%
%  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. 556
%
B = [beta(i)  beta(j)  beta(m)
     gamma(i) gamma(j) gamma(m)]/(2*A);
B = Asign*B;
%
%  D matrix p. 480
%
D = [K 0; 0 K];
DB = D*B;