% Zero stiffness matrix

K = zeros(Nv, Nv);
b = zeros(Nv, 1);


% Loop over elements and calculate residual and stiffness matrix

for ii = 1:Nt,
  
  kn(1) = tri2nod(1,ii);
  kn(2) = tri2nod(2,ii);
  kn(3) = tri2nod(3,ii);
    
  xe(1) = xy(1,kn(1));
  xe(2) = xy(1,kn(2));
  xe(3) = xy(1,kn(3));

  ye(1) = xy(2,kn(1));
  ye(2) = xy(2,kn(2));
  ye(3) = xy(2,kn(3));

  
  % Derivatives of node 1's interpolant 
  dNdxi(1,1) = -1.0; % with respect to xi1
  dNdxi(1,2) = -1.0; % with respect to xi2
  
  % Derivatives of node 2's interpolant
  dNdxi(2,1) =  1.0; % with respect to xi1
  dNdxi(2,2) =  0.0; % with respect to xi2

  % Derivatives of node 3's interpolant
  dNdxi(3,1) =  0.0; % with respect to xi1
  dNdxi(3,2) =  1.0; % with respect to xi2
  
  % Sum these to find dxdxi (note: these are constant within an element)
  dxdxi = zeros(2,2);
  for nn = 1:3,
    dxdxi(1,:) = dxdxi(1,:) + xe(nn)*dNdxi(nn,:);
    dxdxi(2,:) = dxdxi(2,:) + ye(nn)*dNdxi(nn,:);
  end
  
  % Calculate determinant for area weighting
  J = dxdxi(1,1)*dxdxi(2,2) - dxdxi(1,2)*dxdxi(2,1);
  A = 0.5*abs(J); % Area is have the Jacobian
  
  % Invert dxdxi to find dxidx using inversion rule for a 2x2 matrix
  dxidx = [ dxdxi(2,2)/J, -dxdxi(1,2)/J; ...
	   -dxdxi(2,1)/J,  dxdxi(1,1)/J];
  
  % Calculate dNdx 
  dNdx = dNdxi*dxidx;

  % Add contributions to stiffness matrix for node 1 weighted residual
  K(kn(1), kn(1)) = K(kn(1), kn(1)) + (dNdx(1,1)*dNdx(1,1) + dNdx(1,2)*dNdx(1,2))*A;
  K(kn(1), kn(2)) = K(kn(1), kn(2)) + (dNdx(1,1)*dNdx(2,1) + dNdx(1,2)*dNdx(2,2))*A;
  K(kn(1), kn(3)) = K(kn(1), kn(3)) + (dNdx(1,1)*dNdx(3,1) + dNdx(1,2)*dNdx(3,2))*A;
  
  % Add contributions to stiffness matrix for node 2 weighted residual
  K(kn(2), kn(1)) = K(kn(2), kn(1)) + (dNdx(2,1)*dNdx(1,1) + dNdx(2,2)*dNdx(1,2))*A;
  K(kn(2), kn(2)) = K(kn(2), kn(2)) + (dNdx(2,1)*dNdx(2,1) + dNdx(2,2)*dNdx(2,2))*A;
  K(kn(2), kn(3)) = K(kn(2), kn(3)) + (dNdx(2,1)*dNdx(3,1) + dNdx(2,2)*dNdx(3,2))*A;
  
  % Add contributions to stiffness matrix for node 3 weighted residual
  K(kn(3), kn(1)) = K(kn(3), kn(1)) + (dNdx(3,1)*dNdx(1,1) + dNdx(3,2)*dNdx(1,2))*A;
  K(kn(3), kn(2)) = K(kn(3), kn(2)) + (dNdx(3,1)*dNdx(2,1) + dNdx(3,2)*dNdx(2,2))*A;
  K(kn(3), kn(3)) = K(kn(3), kn(3)) + (dNdx(3,1)*dNdx(3,1) + dNdx(3,2)*dNdx(3,2))*A;
  
end


% Loop over boundary edges and account for bc's

for ii = 1:Nbc,
  
  bnum = bedge(3,ii) + 1;
    
  kn(1) = bedge(1,ii);
  kn(2) = bedge(2,ii);

  xe(1) = xy(1,kn(1));
  xe(2) = xy(1,kn(2));

  ye(1) = xy(2,kn(1));
  ye(2) = xy(2,kn(2));

  ds = sqrt((xe(1)-xe(2))^2 + (ye(1)-ye(2))^2);
  
  K(kn(1), kn(1)) = K(kn(1), kn(1)) + hwall(bnum)*ds*(1/3);
  K(kn(1), kn(2)) = K(kn(1), kn(2)) + hwall(bnum)*ds*(1/6);
  b(kn(1))        = b(kn(1))        + hwall(bnum)*ds*0.5*Tcool(bnum);
  
  K(kn(2), kn(1)) = K(kn(2), kn(1)) + hwall(bnum)*ds*(1/6);
  K(kn(2), kn(2)) = K(kn(2), kn(2)) + hwall(bnum)*ds*(1/3);
  b(kn(2))        = b(kn(2))        + hwall(bnum)*ds*0.5*Tcool(bnum);
  
end


% Solve for temperature

Tsol = K\b;