function triangle_neighbor = triangulation_order6_neighbor_triangles ( ...
  triangle_num, triangle_node )

%% TRIANGULATION_ORDER6_NEIGHBOR_TRIANGLES determines triangle neighbors.
%
%  Discussion:
%
%    A triangulation of a set of nodes can be completely described by
%    the coordinates of the nodes, and the list of nodes that make up
%    each triangle.  However, in some cases, it is necessary to know
%    triangle adjacency information, that is, which triangle, if any,
%    is adjacent to a given triangle on a particular side.
%
%    This routine creates a data structure recording this information.
%
%    The primary amount of work occurs in sorting a list of 3 * TRIANGLE_NUM
%    data items.
%
%  Modified:
%
%    01 January 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer TRIANGLE_NUM, the number of triangles.
%
%    Input, integer TRIANGLE_ORDER(6,TRIANGLE_NUM), the nodes that make up 
%    each triangle.
%
%    Output, integer TRIANGLE_NEIGHBOR(3,TRIANGLE_NUM), the three triangles that
%    are direct neighbors of a given triangle.  TRIANGLE_NEIGHBOR(1,I) is 
%    the index of the triangle which touches side 1, defined by nodes 2 and 3,
%    and so on.  TRIANGLE_NEIGHBOR(1,I) is negative if there is no neighbor
%    on that side.  In this case, that side of the triangle lies on the 
%    boundary of the triangulation.
%

%
%  Step 1.
%  From the list of vertices for triangle T, of the form: (I,J,K)
%  construct the three neighbor relations:
%
%    (I,J,1,T) or (J,I,1,T),
%    (J,K,2,T) or (K,J,2,T),
%    (K,I,3,T) or (I,K,3,T)
%
%  where we choose (I,J,1,T) if I < J, or else (J,I,1,T)
%
  for tri = 1 : triangle_num

    i = triangle_node(1,tri);
    j = triangle_node(2,tri);
    k = triangle_node(3,tri);

    if ( i < j )
      row(3*(tri-1)+1,1:4) = [ i, j, 1, tri ];
    else
      row(3*(tri-1)+1,1:4) = [ j, i, 1, tri ];
    end

    if ( j < k )
      row(3*(tri-1)+2,1:4) = [ j, k, 2, tri ];
    else
      row(3*(tri-1)+2,1:4) = [ k, j, 2, tri ];
    end

    if ( k < i )
      row(3*(tri-1)+3,1:4) = [ k, i, 3, tri ];
    else
      row(3*(tri-1)+3,1:4) = [ i, k, 3, tri ];
    end

  end
%
%  Step 2. Perform an ascending dictionary sort on the neighbor relations.
%  We only intend to sort on columns 1 and 2; the routine we call here
%  sorts on columns 1 through 4 but that won't hurt us.
%
%  What we need is to find cases where two triangles share an edge.
%  Say they share an edge defined by the nodes I and J.  Then there are
%  two rows of ROW that start out ( I, J, ?, ? ).  By sorting ROW,
%  we make sure that these two rows occur consecutively.  That will
%  make it easy to notice that the triangles are neighbors.
%
  row = i4row_sort_a ( 3*triangle_num, 4, row );
%
%  Step 3. Neighboring triangles show up as consecutive rows with
%  identical first two entries.  Whenever you spot this happening,
%  make the appropriate entries in TRIANGLE_NEIGHBOR.
%
  triangle_neighbor(1:3,1:triangle_num) = -1;

  irow = 1;

  while ( 1 )

    if ( 3 * triangle_num <= irow )
      break
    end

    if ( row(irow,1) ~= row(irow+1,1) | row(irow,2) ~= row(irow+1,2) )
      irow = irow + 1;
      continue
    end

    side1 = row(irow,3);
    tri1 = row(irow,4);
    side2 = row(irow+1,3);
    tri2 = row(irow+1,4);

    triangle_neighbor(side1,tri1) = tri2;
    triangle_neighbor(side2,tri2) = tri1;

    irow = irow + 2;

  end