function boundary_segments = mesh_boundary ( element_node )

%*****************************************************************************80
%
%% mesh_boundary() determines the segments forming a mesh boundary.
%
%  Discussion:
%
%    This function is only appropriate for 2D meshes constructed from
%    a number of polygonal elements.  All the elements should be the
%    same type, such as triangles or quadrilaterals.  The elements should
%    be 'linear', that is, a triangle should be described by three
%    vertices, with no additional nodes along the edges or in the interior.
%    The mesh must be 'conformal', that is, there should not be any
%    'hanging nodes'.  The vertices of each element should be listed
%    in counter clockwise order.
%
%    The mesh should not include holes.
%
%  Example:
%
%    Consider the following simple quadrilateral mesh:
%
%      1---2---3
%      |   |   |
%      4---5---6
%      |   |   |
%      7---8---9
%
%    for which the element array is:
%
%      element_node = [
%        [ 4, 5, 2, 1 ],
%        [ 5, 6, 3, 2 ],
%        [ 7, 8, 5, 4 ],
%        [ 8, 9, 6, 5 ] ]
%
%    This function will return the boundary segments array: 
%
%      boundary_segments = [
%        [ 1, 4 ]
%        [ 4, 7 ],
%        [ 7, 8 ],
%        [ 8, 9 ],
%        [ 9, 6 ],
%        [ 6, 3 ],
%        [ 3, 2 ],
%        [ 2, 1 ] ]
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    18 December 2022
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer element_node[element_num,element_order]: element_node[i,j]
%    lists the index or label of the j-th node of the i-th element.  
%    If we think of the values as indexes, they can run from 0 to node_num-1,
%    or from 1 to node_num.  But we can also think of them as simply
%    labels, in which case they can be any set of node_num distinct values.
%
%  Output:
%
%    integer boundary_segments[boundary_segment_num,2]: a list of pairs of
%    nodes which define the line segments forming the boundary.  The line
%    segments are given in sequence, so that the boundary is traced in a
%    counter clockwise direction.
%
  [ element_num, element_order ] = size ( element_node );
%
%  Determine the number of segments.
%
  segment_num = element_num * element_order;
%
%  Break each element into segments.  
%  Here is where a dictionary might be better.
%
  segment = zeros ( segment_num, 2 );
  s = 0;
  for e = 1 : element_num
    j = element_order;
    for jp1 = 1 : element_order
      s = s + 1;
      segment(s,1) = element_node(e,j);
      segment(s,2) = element_node(e,jp1);
      j = jp1;
    end
  end
%
%  Lexically sort the rows of the array.
%
  segment = sortrows ( segment );
%
%  From the list of segments, 
%  find rows that are not matched by a reversed copy.
%  These are the boundary segments.
%
  boundary_segment_num = 0;
  for s = 1 : segment_num
    segment(s,1:2) = segment(s,2:-1:1);
    segment2 = unique ( segment, 'rows' );
    if ( size ( segment2, 1 ) == segment_num )
      boundary_segment_num = boundary_segment_num + 1;
      boundary_segments(boundary_segment_num,1:2) = segment(s,2:-1:1);
    end
    segment(s,1:2) = segment(s,2:-1:1);
  end
%
%  Reorder the segments to form a loop.
%
  for b1 = 1 : boundary_segment_num - 1
    for b2 = b1 + 1 : boundary_segment_num
      if ( boundary_segments(b2,1) == boundary_segments(b1,2) )
        temp = boundary_segments(b1+1,1:2);
        boundary_segments(b1+1,1:2) = boundary_segments(b2,1:2);
        boundary_segments(b2,1:2) = temp;
        continue
      end
    end
  end

  return
end

