clear; close all
%
%  Problem 13.34
%
thick = 1;      %% unit thick out of plane
width = 2;      %% m
depth = 1;      %% m
Tleft = 100;
Ttop = 500;
Tbot = 100;     %% C
K = 10;         %% W/m*C
h = 20;         %% W/m^2*C
Tinf = 20;      %% C
snodes = [0      0   %% m
          width  0
          width  depth
          0      depth];
edges = [1 2
         2 3
         3 4
         4 1];
hdata.hmax = 
[nodes,elems] = mesh2d(snodes,edges,hdata);
nnodes = size(nodes,1);
nelems = size(elems,1);
ndof = nnodes;
%
%  Find nodes on edges
%
bot   = find(findedge(nodes,snodes,[1 2]));
right = find(findedge(
top   = find(findedge(
left  = find(findedge(
%
%  Constrain temps
%
bc = [left Tleft*ones(size(left))
      top  Ttop*ones(size(top))
      bot  Tbot*ones(size(bot))];
%
%  Average at corners
%
topleft = intersect(top,left);
ind = find(bc(:,1)==topleft);  bc(ind,2) = (Ttop+Tleft)/2;
botleft = 
ind = 
%
%  Solve
%
bigt = zeros(ndof,1);
bcdof = bc(:,1);
bigt(bcdof) = bc(:,2);
bigk = zeros(ndof,ndof);
bigf = zeros(ndof,1);
for i=1:nelems
    n1 = elems(i,1);  no1 = find(n1==right);
    n2 = elems(i,2);  no2 = find(n2==right);
    n3 = elems(i,3);  no3 = find(n3==right);
    nout = [no1 no2 no3];
    dof = [n1 n2 n3];
    if length(nout) == 2
        conv = [right(nout)' h Tinf];
        [kc,fc] = triahk(nodes,n1,n2,n3,K,thick,conv);
        bigf(dof) = bigf(dof) + fc;
    else
        kc = triahk(nodes,n1,n2,n3,K,thick);
    end
    bigk(dof,dof) = bigk(dof,dof) + kc;
end
freedof = 1:ndof;
freedof(bcdof) = [];
RHS = bigf(freedof) - bigk(freedof,bcdof)*bigt(bcdof);
bigt(freedof) = bigk(freedof,freedof)\RHS;
%
%  Plot
%
for i=1:nelems
    nn = elems(i,[1 2 3 1]);
    xx = nodes(nn,1);
    yy = nodes(nn,2);
    patch(xx,yy,bigt(nn));
    hold on
end
axis equal
%
%  Recover fluxes, draw quivers
%
qx = zeros(nelems,1);
qy = qx;
xc = zeros(nelems,1);
yc = xc;
for i=1:nelems
    n1 = elems(i,1);
    n2 = elems(i,2);
    n3 = elems(i,3);
    nn = [n1 n2 n3];
    DB = triahf(nodes,n1,n2,n3,K,thick);
    tt = bigt(nn);
    q = -DB*tt;
    qx(i) = q(1);
    qy(i) = q(2);
    xc(i) = sum(nodes(nn,1))/3;
    yc(i) = sum(nodes(nn,2))/3;
end
quiver(xc,yc,qx,qy,'k');
%
%  Add up total convection losses
%
Q = h*depth*(mean(bigt(right)) - Tinf);  %% roughly
title(sprintf('Tinf = %d, loss = %8.0f watts/m loss to convection',Tinf,Q));