function [edge2tri, bedge2tri] = SetupEdgeList(tri2nod, bedge)

% This routine determines the edge data structures from a
% triangular mesh connectivity pattern (stored in tri2nod)
% and a list of boundary edge nodes (stored in bedge).  The
% resulting data structures are edge2tri and bedge2tri which store
% the following info:
%
% edge2tri(1, i) = index of first  node of edge i
% edge2tri(2, i) = index of second node of edge i
% edge2tri(3, i) = index of triangle to left of edge i
% edge2tri(4, i) = index of triangle to right of edge i
%
% The edge2tri nodes are oriented so that the triangle in 
% edge2tri(3,i) is to the left when moving from node 1 to node 2
%
% bedge2tri(1, i) = index of first node for bedge i
% bedge2tri(2, i) = index of second node for bedge i
% bedge2tri(3, i) = index of the triangle for bedge i
% bedge2tri(4, i) = boundary ID for bedge i
% 
% The bedge2tri nodes are oriented so that the triangle is to the left
% when moving from node 1 to node 2
%
% Note the edges in bedge2tri are stored in the same order as they
% appear in the input array bedge2tri
%

% Get sizes
[Ntemp, Nt] = size(tri2nod);
[Ntemp, Nbe] = size(bedge);

% Determine number of internal edges
Ne = (3*Nt - Nbe)/2;

% Determine maximum node number
Nv = max(max(tri2nod));

% Set up arrays to hold edge lists
edge2tri  = zeros(4,Ne);
bedge2tri = zeros(4,Nbe);

% Copy over contents of bedge to bedge2tri
bedge2tri(1:2,:) = bedge(1:2,:);
bedge2tri(4,:)   = bedge(3,:);

% Set up work space 
nod2edge = sparse(Nv,Nv);

% Set interior triangle index in bedge2tri list
for i = 1:Nbe,

  % Get node numbers
  k1 = bedge2tri(1,i);
  k2 = bedge2tri(2,i);
  
  % Find min and max node numbers
  kmin = min(k1, k2);
  kmax = max(k1, k2);
  
  % Check if edge does not exist in nod2edge list
  ktri = nod2edge(kmin,kmax);
  if (ktri ~= 0), % Edge exists: Error
    fprintf('Multiple boundary edge. ERROR: ktri ~= 0.  ktri = %i\n',ktri);
    return;
  else, % Edge does not exist, insert it
    nod2edge(kmin,kmax) = -i; % Add edge to list noting triangle
  end
end

% Loop over triangles and find all edges
ie = 0;
for i = 1:Nt,
  for j = 1:3,
    
    % Get node numbers
    k1 = tri2nod(         j, i);
    k2 = tri2nod(mod(j,3)+1, i);
    
    % Find min and max node numbers
    kmin = min(k1, k2);
    kmax = max(k1, k2);
    
    % Check if edge does not exist in nod2edge list
    ktri = nod2edge(kmin,kmax);
    if (ktri == 0), % Edge does not exist
      nod2edge(kmin,kmax) = i; % Add edge to list noting triangle
    elseif (ktri > 0), % Edge exists, insert into edge2tri list
      ie = ie + 1;
      edge2tri(1,ie) = kmin;
      edge2tri(2,ie) = kmax;
      if (kmin == k1),
        edge2tri(3,ie) = i;
        edge2tri(4,ie) = ktri;
      else
        edge2tri(3,ie) = ktri;
        edge2tri(4,ie) = i;
      end
      nod2edge(kmin,kmax) = 0; 
    else, % Boundary edge exists, insert into bedge2tri list
      bedge2tri(1,-ktri) = k1;
      bedge2tri(2,-ktri) = k2;
      bedge2tri(3,-ktri) = i;
      nod2edge(kmin,kmax) = 0;
    end
    
  end
end

% Check that no non-zero entries exist in nod2edge
nz = nnz(nod2edge);
if (nz > 0),
  fprintf('ERROR: non-zero entries remain in nod2edge\n');
end