function [node,elem,bdEdge,HB] = uniformrefine(node,elem,bdEdge,HB)
% UNIFORMREFINE divide each traingle into 4 similar triangles. 
%
% [node,elem,bdEdge] = uniformbisect(node,elem,bdEdge)
% [node,elem] = uniformbisect(node,elem)
%
% bdEdge in the input and output is optional and may be omitted.
%--------------------------------------------------------------------------
% Copyright (C) 2008 Long Chen. See COPYRIGHT.txt for details.
%--------------------------------------------------------------------------

%------------------- Construct data structure -----------------------------
totalEdge = sort([elem(:,[2,3]); elem(:,[3,1]); elem(:,[1,2])],2);
[edge, i2, j] = unique(totalEdge,'rows');
N = size(node,1); NT = size(elem,1); NE = size(edge,1);
elem2edge = reshape(j,NT,3);
%------------------- Add new nodes --------------------------
node(N+1:N+NE,:) = (node(edge(:,1),:)+node(edge(:,2),:))/2; 
HB(:,1) = (N+1:N+NE)'; HB(:,2) = edge(:,1); HB(:,3) = edge(:,2);
edge2newNode = int32((N+1:N+NE)'); %edge2newNode = int32(edge2newNode);
%--------------------------------------------------------------------------
% refine each triangle into four triangles by regular refinement
%     3
%    / \
%   5 - 4
%  / \ / \
% 1 - 6 - 2
%--------------------------------------------------------------------------
t=1:NT;
p(t,1:3) = elem(t,1:3);
p(t,4:6) = edge2newNode(elem2edge(t,1:3));
elem(t,:) = [p(t,1), p(t,6), p(t,5)];
elem(NT+1:2*NT,:) = [p(t,6), p(t,2), p(t,4)];
elem(2*NT+1:3*NT,:) = [p(t,5), p(t,4), p(t,3)];
elem(3*NT+1:4*NT,:) = [p(t,4), p(t,5), p(t,6)];
%-------------------- Update boundary edges-------------------------------
if (nargin==3) && (~isempty(bdEdge))
    bdEdge(NT+1:2*NT,[1 3]) = bdEdge(t,[1 3]); 
    bdEdge(2*NT+1:3*NT,[1 2]) = bdEdge(t,[1 2]); 
    bdEdge(3*NT+1:4*NT,1) = 0;
    bdEdge(t,1) = 0;
else
    bdEdge = [];
end