function element_neighbor = neighbor_elements_q4_mesh ( element_num, ...
  element_node )

%*****************************************************************************80
%
%% NEIGHBOR_ELEMENTS_Q4_MESH determines element neighbors in a Q4 mesh.
%
%  Discussion:
%
%    A quadrilateral mesh of a set of nodes can be completely described by
%    the coordinates of the nodes, and the list of nodes that make up
%    each element.  However, in some cases, it is necessary to know
%    element adjacency information, that is, which element, if any,
%    is adjacent to a given element on a particular side.
%
%    This routine creates a data structure recording this information.
%
%    The primary amount of work occurs in sorting a list of 4 * ELEMENT_NUM
%    data items.
%
%    Note that COL is a work array allocated dynamically inside this
%    routine.  It is possible, for very large values of ELEMENT_NUM,
%    that the necessary amount of memory will not be accessible, and the
%    routine will fail.  This is a limitation of the implementation of
%    dynamic arrays in FORTRAN90.  One way to get around this would be
%    to require the user to declare ROW in the calling routine
%    as an allocatable array, get the necessary memory explicitly with
%    an ALLOCATE statement, and then pass ROW into this routine.
%
%    Of course, the point of dynamic arrays was to make it easy to
%    hide these sorts of temporary work arrays from the poor user%
%
%    This routine was revised to store the edge data in a column
%    array rather than a row array.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 February 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer ELEMENT_NUM, the number of elements.
%
%    Input, integer ELEMENT_NODE(4,ELEMENT_NUM), the nodes
%    that make up each element.
%
%    Output, integer ELEMENT_NEIGHBOR(4,ELEMENT_NUM), lists the
%    neighboring element on each side of a given element, or -1 if there is
%    no neighbor.
%
  element_order = 4;

  col = zeros(4,4*element_num);
%
%  Step 1.
%  From the list of nodes for element Q, of the form: (I,J,K,L)
%  construct the four neighbor relations:
%
%    (I,J,1,Q) or (J,I,1,Q),
%    (J,K,2,Q) or (K,J,2,Q),
%    (K,L,3,Q) or (L,K,3,Q)
%    (K,I,4,Q) or (I,K,4,Q)
%
%  where we choose (I,J,1,Q) if I < J, or else (J,I,1,Q)
%
  for q = 1 : element_num

    i = element_node(1,q);
    j = element_node(2,q);
    k = element_node(3,q);
    l = element_node(4,q);

    if ( i < j )
      col(1:4,4*(q-1)+1) = [ i, j, 1, q ]';
    else
      col(1:4,4*(q-1)+1) = [ j, i, 1, q ]';
    end

    if ( j < k )
      col(1:4,4*(q-1)+2) = [ j, k, 2, q ]';
    else
      col(1:4,4*(q-1)+2) = [ k, j, 2, q ]';
    end

    if ( k < l ) 
      col(1:4,4*(q-1)+3) = [ k, l, 3, q ]';
    else
      col(1:4,4*(q-1)+3) = [ l, k, 3, q ]';
    end

    if ( l < i )
      col(1:4,4*(q-1)+4) = [ l, i, 4, q ]';
    else
      col(1:4,4*(q-1)+4) = [ i, l, 4, q ]';
    end

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

  icol = 1;

  while ( 1 )

    if ( 4 * element_num <= icol )
      break
    end

    if ( col(1,icol) ~= col(1,icol+1) | col(2,icol) ~= col(2,icol+1) )
      icol = icol + 1;
      continue
    end

    side1 = col(3,icol);
    q1    = col(4,icol);
    side2 = col(3,icol+1);
    q2    = col(4,icol+1);

    element_neighbor(side1,q1) = q2;
    element_neighbor(side2,q2) = q1;

    icol = icol + 2;

  end

  return
end