function [k,fconv] = triahk(nodes,n1,n2,n3,K,t,conv);
%
%  function [k,fconv] = triahk(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,7,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. 480
%
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];
%
%  elem stiffness matrix  eq. 6.2.52
%
k = t*A*B'*D*B;
%
%  Convection matrix and vector if needed
%
if nargin > 6
    m1 = conv(1);
    m2 = conv(2);
    nn = [n1 n2 n3];
    im(1) = find(m1==nn);
    im(2) = find(m2==nn);
    if isempty(im(1)) || isempty(im(2))
        error('Conv nodes not found');
    end
    x1 = nodes(m1,1);  y1 = nodes(n1,2);
    x2 = nodes(m2,1);  y2 = nodes(n2,2);
    len = sqrt( (x2-x1)^2 + (y2-y1)^2 );
    h = conv(3);
    Tinf = conv(4);
    k(im,im) = k(im,im) + (h*t*len/6) * [2 1; 1 2];
    fconv = zeros(3,1);
    fconv(im) = (h*t*Tinf*len/2) * [1;1];
end;