clear; close all
%
%  Coolant tubes
%
K = 0.1;            %% kW/m*degC
depth = 0.01;       %% 10 mm
radius = 0.002;     %% 2 mm hole radius
spacing = 0.02;     %% 20 mm
thick = 1;          %% unit out-of-plane dim
Thole = 20;         %% 20C coolant temp
top_flux = 500;     %% 500 kW/m^2 (!)
%
%  Create mesh
%
snode1 = [0        0
          spacing  0
          spacing  depth
          0        depth];
edge1 = [1 2
         2 3
         3 4
         4 1];
xc = spacing/2;
yc = depth/2;
[snode2,edge2] = hole([xc yc],radius,20,4);
hdata.hmax = spacing/20;
hdata.hmax = spacing/25;
snodes = [snode1; snode2];
edges = [edge1; edge2];
[nodes,elems] = mesh2d(snodes,edges,hdata);
%
%  Generate global conductance matrix bigk
%
nnodes = size(nodes,1);
nelems = size(elems,1);
ndof = nnodes;
bigk = zeros(ndof,ndof);
for i=1:nelems
    n1 = elems(i,1);
    n2 = elems(i,2);
    n3 = elems(i,3);
    kc = triahk(nodes,n1,n2,n3,K,thick);
    dof = [n1 n2 n3];
    bigk(dof,dof) = bigk(dof,dof) + kc;
end;
%
%  Flux thru top
%
bigf = zeros(ndof,1);
top_nodes = findedge(nodes,snode1,[3 4]);  %% nodes on top edge
top_nodes = find(top_nodes);
ntop = length(top_nodes);
for i=1:ntop-1
    n1 = top_nodes(i);
    n2 = top_nodes(i+1);
    elen = nodes(n1,1) - nodes(n2,1);
    elen = abs(elen);
    nn = [n1 n2];
    bigf(nn) = bigf(nn) + elen*thick*top_flux/2;  %% half edge flux each node
end;
%
%  Prescribed temps on hole
%
bigt = zeros(ndof,1);
hole_nodes = findedge(nodes,snodes,edge2);
hole_nodes = find(hole_nodes);
bcdof = hole_nodes;
bigt(hole_nodes) = Thole;
freedof = 1:ndof;
freedof(bcdof) = [];
kff = bigk(freedof,freedof);
kfb = bigk(freedof,bcdof);
RHS = bigf(freedof) - bigk(freedof,bcdof)*bigt(bcdof);
bigt(freedof) = bigk(freedof,freedof)\RHS;
%
%  Temps & fluxes
%
figure;
hh = zeros(nelems);
for i=1:nelems
    nn = elems(i,:);
    xx = nodes(nn,1);
    yy = nodes(nn,2);
    hh(i) = patch(xx,yy,bigt(nn));
end;
hold on
axis equal
set(gca,'XTick',[],'YTick',[]);
title(sprintf('Max temp %4.1f ^oC',max(bigt)));
%qx = zeros(nelems,1);  qy=qx;
%qy = zeros(nelems,1);
%xc = zeros(nelems,1);  yc=xc;
%for i=1:nelems
%    nn = elems(i,:);
%    xx = nodes(nn,1);
%    yy = nodes(nn,2);
%    DB = triahf(nodes,n1,n2,n3,K,thick);
%    q = -DB*bigt(nn);
%    qx(i) = q(1);
%    qy(i) = q(2);
%    xc(i) = sum(nodes(nn,1))/3;
%    yc(i) = sum(nodes(nn,2))/3;
%end;
%hq = quiver(xc,yc,qx,qy,'k');
colorbar