function [ node_num2, element_num2, node_xy2, element_node2, ...
  element_neighbor2 ] = triangulation_refine_local ( node_num, ...
  element_num, node_xy, element_node, element_neighbor, element )

%*****************************************************************************80
%
%% TRIANGULATION_REFINE_LOCAL locally refines a triangulation.
%
%  Discussion:
%
%    To use and to update the element neighbor array efficiently, I here
%    assume that the definition is as follows:
%
%      ELEMENT_NEIGHBOR(I,E) is the index of the element which is adjacent
%      to element E, and opposite to the I-th node, whose index is
%      ELEMENT_NODE(I,E).
%
%                 N3
%                 |\
%                 | \
%      Neighbor 2 |  \ Neighbor 1
%                 | E \
%                 |    \
%                 N1----N2
%
%                 Neighbor 3
%
%    Unfortunately, I don't think I have always followed this convention,
%    nor do I think this is a standard convention.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    30 August 2011
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer NODE_NUM, the number of nodes.
%
%    integer ELEMENT_NUM, the number of elements.
%
%    real NODE_XY(2,NODE_NUM), the point coordinates.
%
%    integer ELEMENT_NODE(3,ELEMENT_NUM), the nodes that make 
%    up the elements.
%
%    integer ELEMENT_NEIGHBOR(3,ELEMENT_NUM), the element neighbors on 
%    each side.  If there is no element neighbor on a particular side, the 
%    value of should be negative. 
%
%    integer ELEMENT, the index of the element to be refined.
%
%  Output:
%
%    integer NODE_NUM2, the number of nodes.
%
%    integer ELEMENT_NUM2, the number of elements.
%
%    real NODE_XY2(2,NODE_NUM2), the point coordinates.
%
%    integer ELEMENT_NODE2(3,ELEMENT_NUM2), the nodes that make 
%    up the elements.
%
%    integer ELEMENT_NEIGHBOR2(3,ELEMENT_NUM2), the element neighbors on 
%    each side.  If there is no element neighbor on a particular side, the 
%    value should be negative. 
%

%
%  1: Determine the three new nodes.
%
  node_num2 = node_num + 3;
  n1 = element_node(element,1);
  n2 = element_node(element,2);
  n3 = element_node(element,3);
  n12 = node_num + 1;
  n23 = node_num + 2;
  n31 = node_num + 3;

  node_xy2 = node_xy;

  node_xy2(n12,1:2) = 0.5 * ( node_xy(n1,1:2) + node_xy(n2,1:2) );
  node_xy2(n23,1:2) = 0.5 * ( node_xy(n2,1:2) + node_xy(n3,1:2) );
  node_xy2(n31,1:2) = 0.5 * ( node_xy(n3,1:2) + node_xy(n1,1:2) );
%
%  2: Revise the element information.  Where there were 4 elements,
%  there will now be 10 (unless the element is missing neighbors).
%
  ea1 = element_neighbor(element,1);
  eb1 = element_neighbor(element,2);
  ec1 = element_neighbor(element,3);

  element_num2 = element_num;

  element_num2 = element_num2 + 1;
  e1 = element_num2;

  element_num2 = element_num2 + 1;
  e2 = element_num2;

  element_num2 = element_num2 + 1;
  e3 = element_num2;

  if ( 0 < ea1 )
    element_num2 = element_num2 + 1;
    ea2 = element_num2;
  else
    ea2 = -1;
  end

  if ( 0 < eb1 )
    element_num2 = element_num2 + 1;
    eb2 = element_num2;
  else
    eb2 = -1;
  end

  if ( 0 < ec1 )
    element_num2 = element_num2 + 1;
    ec2 = element_num2;
  else
    ec2 = -1;
  end

  element_node2 = zeros(element_num2,3);
 
  element_node2(1:element_num,1:3) = element_node(1:element_num,1:3);

  element_node2(element,1:3) = [ n23, n31, n12 ];
  element_node2(e1,1:3)      = [ n1, n12, n31 ];
  element_node2(e2,1:3)      = [ n2, n23, n12 ];
  element_node2(e3,1:3)      = [ n3, n31, n23 ];

  if ( 0 < ea1 )
    na = outlier ( element_node(element,1:3), element_node(ea1,1:3) );
    element_node2(ea1,1:3) = [ n2, na, n23 ];
    element_node2(ea2,1:3) = [ na, n3, n23 ];
  end

  if ( 0 < eb1 )
    nb = outlier ( element_node(element,1:3), element_node(eb1,1:3) );
    element_node2(eb1,1:3) = [ n3, nb, n31 ];
    element_node2(eb2,1:3) = [ nb, n1, n31 ];
  end

  if ( 0 < ec1 )
    nc = outlier ( element_node(element,1:3), element_node(ec1,1:3) );
    element_node2(ec1,1:3) = [ n1, nc, n12 ];
    element_node2(ec2,1:3) = [ nc, n2, n12 ];
  end
%
%  3: Revise the neighbor information.
%
  element_neighbor2 = zeros(element_num2,3);
  element_neighbor2(1:element_num,1:3) = element_neighbor(1:element_num,1:3);

  element_neighbor2(element,1:3) = [ e1, e2, e3 ];

  element_neighbor2(e1,1:3) = [ element, eb2, ec1 ];
  element_neighbor2(e2,1:3) = [ element, ec2, ea1 ];
  element_neighbor2(e3,1:3) = [ element, ea2, eb1 ];

  if ( 0 < ea1 )
    i = i4vec_find ( 3, element_neighbor(ea1,1:3), element );
    i = i4_wrap ( i + 1, 1, 3 );
    ea3 = element_neighbor(ea1,i);
    i = i4_wrap ( i + 1, 1, 3 );
    ea4 = element_neighbor(ea1,i);
    element_neighbor2(ea1,1:3) = [ ea2, e2, ea3 ];
    element_neighbor2(ea2,1:3) = [ e3, ea1, ea4 ];
    i = i4vec_find ( 3, element_neighbor(ea4,1:3), ea1 );
    element_neighbor2(ea4,i) = ea2;
  end

  if ( 0 < eb1 )
    i = i4vec_find ( 3, element_neighbor(eb1,1:3), element );
    i = i4_wrap ( i + 1, 1, 3 );
    eb3 = element_neighbor(eb1,i);
    i = i4_wrap ( i + 1, 1, 3 );
    eb4 = element_neighbor(eb1,i);
    element_neighbor2(eb1,1:3) = [ eb2, e3, eb3 ];
    element_neighbor2(eb2,1:3) = [ e1, eb1, eb4 ];
    i = i4vec_find ( 3, element_neighbor(eb4,1:3), eb1 );
    element_neighbor2(eb4,i) = eb2;
  end

  if ( 0 < ec1 )
    i = i4vec_find ( 3, element_neighbor(ec1,1:3), element );
    i = i4_wrap ( i + 1, 1, 3 );
    ec3 = element_neighbor(ec1,i);
    i = i4_wrap ( i + 1, 1, 3 );
    ec4 = element_neighbor(ec1,i);
    element_neighbor2(ec1,1:3) = [ ec2, e1, ec3 ];
    element_neighbor2(ec2,1:3) = [ e2, ec1, ec4 ];
    i = i4vec_find ( 3, element_neighbor(ec4,1:3), ec1 );
    element_neighbor2(ec4,i) = ec2;
  end

  return
end
 
