function etoe = mesh_etoe ( e_order, e_num, etov )

%*****************************************************************************80
%
%% MESH_ETOE determines the element-to-element connectivity.
%
%  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.
%
%    This routine was modified to use columns instead of rows.
%
%  Example:
%
%    The input information from ETOV:
%
%    E           ETOV
%    --------   ---------------
%     1         3      4      1
%     2         3      1      2
%     3         3      2      8
%     4         2      1      5
%     5         8      2     13
%     6         8     13      9
%     7         3      8      9
%     8        13      2      5
%     9         9     13      7
%    10         7     13      5
%    11         6      7      5
%    12         9      7      6
%    13        10      9      6
%    14         6      5     12
%    15        11      6     12
%    16        10      6     11
%
%    The output information in TRIANGLE_NEIGHBOR:
%
%    E          Neighboring E
%    --------  ---------------------
%
%     1        -1     -1      2
%     2         1      4      3
%     3         2      5      7
%     4         2     -1      8
%     5         3      8      6
%     6         5      9      7
%     7         3      6     -1
%     8         5      4     10
%     9         6     10     12
%    10         9      8     11
%    11        12     10     14
%    12         9     11     13
%    13        -1     12     16
%    14        11     -1     15
%    15        16     14     -1
%    16        13     15     -1
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 September 2017
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer E_ORDER, the order of the elements.
%
%    Input, integer E_NUM, the number of triangles.
%
%    Input, integer E_TO_V(E_ORDER,E_NUM), the nodes that 
%    make up each triangle.
%
%    Output, integer E_TO_E(E_NUM,E_ORDER), 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 element E, of the form: (V1,V2,...VO)
%  construct the neighbor relations:
%
%    (VO,V1,1,E)
%    (V1,V2,2,E)
%    (V2,V3,3,E)
%    ...
%    (VO-1,VO,O,E)
%
%  but sort the first two compoents.
%
  r = zeros ( e_order * e_num, 4 );
  
  row = 0;

  for e = 1 : e_num

    v2 = etov(e_order,e);

    for o = 1 : e_order
      row = row + 1;
      v1 = v2;
      v2 = etov(o,e);
      r(row,1) = min ( v1, v2 );
      r(row,2) = max ( v1, v2 );
      r(row,3) = o;
      r(row,4) = e;
    end

  end
%
%  Step 2. Sort the data.
%
  r = sortrows ( r );
%
%  Step 3. Neighboring elements show up as consecutive rows with
%  identical first two entries.
%
  etoe(1:e_num,1:e_order) = -1;

  row = 1;

  while ( true )

    if ( e_num * e_order <= row )
      break
    end

    if ( r(row,1) ~= r(row+1,1) || r(row,2) ~= r(row+1,2) )
      row = row + 1;
    else
      s1 = r(row,3);
      e1 = r(row,4);
      s2 = r(row+1,3);
      e2 = r(row+1,4);
      etoe(e1,s1) = e2;
      etoe(e2,s2) = e1;
      row = row + 2; 
    end

  end

  return
end
