function [hb,lmark] = decomposition(node,elem,level,N0)
%--------------------------------------------------------------------------
% DECOMPOSITION gives the hierarchical structure of nodes
%   hb(i,2) ---- hb(i,1) --- hb(i,3)
%   lmark(k): the index of k-th level end at lmark(k)
%   thus lmark(k-1)+1:lmark(k) is the node index in the k-th level

% Long Chen
% 11-28-2007
%--------------------------------------------------------------------------
N=size(node,1);
hb=zeros(N,3);
lmark(level+1)=N;
for k=level+1:-1:2
    [node,elem,newhb] = coarsening(node,elem,N0);
    Nnew=size(newhb,1);
    lmark(k-1)=lmark(k)-Nnew;
    hb(lmark(k-1)+1:lmark(k),:)=newhb;
end

%--------------------------------------------------------------------------
function [node,elem,newhb] = coarsening(node,elem,N0)
%--------------------------------------------------------------------------
% Construct data structure
%--------------------------------------------------------------------------
N = size(node,1); NT = size(elem,1);
dualEdge = sparse(elem(:,[1,2,3]),elem(:,[2,3,1]),[1:NT,1:NT,1:NT]);
valence = accumarray(elem(:),ones(3*NT,1),[N 1]);

newestNode = unique(elem((elem(:,1)>N0),1));
isGood = (valence(newestNode) == 2) | (valence(newestNode) == 4);
goodNode = newestNode(isGood);
marker(goodNode)=1;

newhb=zeros(N,3);
%--------------------------------------------------------------------------
% Coarsen the mesh by removing marked nodes
%--------------------------------------------------------------------------
for t = 1:NT
    p = elem(t,1);
    if (p > 0) % otherwise t is already removed
        if (marker(p) == 1) % p is a marked good point
            brother = dualEdge(elem(t,2),p);
            elem(t,1) = elem(t,2);
            elem(t,2) = elem(t,3);
            elem(t,3) = elem(brother,2);
            elem(brother,1) = 0; % brother is going to be deleted
            if (newhb(p,2)==0) 
                newhb(p,2) = elem(t,2);
                newhb(p,3) = elem(t,3);
            end
        end % endif for all markd good nodes
    end
end % end of for loop on all elements
newhb = [goodNode, newhb(goodNode,2:3)];
%--------------------------------------------------------------------------
% Clear element and eta arrarys
%--------------------------------------------------------------------------
ix = (elem(:,1) == 0); elem(ix,:) = [];