function mesh_display ( prefix, node_show, element_show )

%*****************************************************************************80
%
%% mesh_display() displays a mesh.
%
%  Discussion:
%
%    mesh_display() displays the nodes and elements of a mesh.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    03 January 2023
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    character prefix, the common prefix for the node and element files:
%    * prefix_nodes.txt,     the node coordinates.
%    * prefix_elements.txt,  the nodes that make up each element.
%
%    integer node_show, indicates the node visibility:
%    0: do not show the nodes;
%    1:        show the nodes;
%    2:        show the nodes, and label them.
%
%    integer element_show, indicates the element visibility:
%    0: do not show the element;
%    1:        show the element;
%    2:        show the element, and color them.
%    3:        show the element, and color them, and label them.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'mesh_display():\n' );
  fprintf ( 1, '  MATLAB/Octave version "%s"\n', version ( ) );
  fprintf ( 1, '  Read a NODE_NUMx2 node coordinate dataset.\n' );
  fprintf ( 1, '  Read an ELEMENT_NUMxELEMENT_ORDER element dataset.\n' );
  fprintf ( 1, '  Display an image of the nodes and elements.\n' );
%
%  Read the node data.
%  Temporarily transpose the data.
%
  node_filename = strcat ( prefix, '_nodes.txt' );
  node_xy = load ( node_filename );
  [ node_num, dim_num ] = size ( node_xy );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Read node coordinate data in "%s".\n', node_filename );
  fprintf ( 1, '  Spatial dimension DIM_NUM = %d\n', dim_num );
  fprintf ( 1, '  Number of nodes NODE_NUM  = %d\n', node_num );
%
%  Read the element data.
%
  element_filename = strcat ( prefix, '_elements.txt' );
  element_node = load ( element_filename );
  [ element_num, element_order ] = size ( element_node );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Read element node data in "%s".\n', element_filename );
  fprintf ( 1, '  Element order      = %d\n', element_order );
  fprintf ( 1, '  Number of elements = %d\n', element_num );
  fprintf ( 1, '\n' );
%
%  Detect and correct zero-based node indexing.
%
  element_node = mesh_base_one ( node_num, element_node );
%
%  Display a plot.
%
  mesh_display_matlab ( prefix, node_num, node_xy, element_order, ...
    element_num, element_node, node_show, element_show );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'mesh_display():\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  return
end
function element_node = mesh_base_one ( node_num, element_node )

%*****************************************************************************80
%
%% mesh_base_one() ensures that the element definition is one-based.
%
%  Discussion:
%
%    The ELEMENT_NODE array contains nodes indices that form elements.
%    The convention for node indexing might start at 0 or at 1.
%
%    Since a MATLAB program will naturally assume a 1-based indexing, it is
%    necessary to check a given element definition and, if it is actually
%    0-based, to convert it.
%
%    This function attempts to detect 0-based node indexing and correct it.
%
%    Thanks to Mark Asch for pointing out an error, 01 December 2009!
%
%  Licensing:
%
%    This code is distributed under the MIT license. 
%
%  Modified:
%
%    22 October 2022
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer NODE_NUM, the number of nodes.
%
%    integer ELEMENT_NODE(ELEMENT_NUM,ELEMENT_ORDER), the element definitions.
%
%  Output:
%
%    integer ELEMENT_NODE(ELEMENT_NUM,ELEMENT_ORDER), the updated definitions.
%
  node_min = min ( min ( element_node ) );
  node_max = max ( max ( element_node ) );

  if ( node_min == 0 && node_max == node_num - 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'mesh_base_one():\n' );
    fprintf ( 1, '  The element indexing appears to be 0-based!\n' );
    fprintf ( 1, '  This will be converted to 1-based.\n' );
    element_node = element_node + 1;
  elseif ( node_min == 1 && node_max == node_num )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'mesh_base_one():\n' );
    fprintf ( 1, '  The element indexing appears to be 1-based!\n' );
    fprintf ( 1, '  No conversion is necessary.\n' );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'mesh_base_one(): Warning!\n' );
    fprintf ( 1, '  The element indexing is not of a recognized type.\n' );
  end

  return
end
function mesh_display_matlab ( prefix, node_num, node_xy, element_order, ...
  element_num, element_node, node_show, element_show )

%*****************************************************************************80
%
%% mesh_display_matlab() displays a mesh.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    03 January 2023
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    string PREFIX, the file name prefix that identifies the data.
%
%    integer NODE_NUM, the number of points.
%
%    real NODE_XY(NODE_NUM,2), the nodes.
%
%    integer ELEMENT_ORDER, the order of the elements.
%
%    integer ELEMENT_NUM, the number of elements.
%
%    integer ELEMENT_NODE(ELEMENT_NUM,ELEMENT_ORDER), lists, 
%    for each element, the indices of the points that form the vertices 
%    of the element.
%
%    logical NODE_SHOW:
%    0, do not show the nodes.
%    1, show the nodes.
%    2, show the nodes, and label them.
%
%    logical ELEMENT_SHOW, 
%    0, do not show the elements.
%    1, show the elements.
%    2, show the elements, and color them.
%    3, show the elements, and color them, and label them.
%

%
%  Clear the graphics screen.
%
  clf ( );
%
%  Get the scale.
%
  xy_min(1) = min ( node_xy(:,1) );
  xy_max(1) = max ( node_xy(:,1) );

  xy_min(2) = min ( node_xy(:,2) );
  xy_max(2) = max ( node_xy(:,2) );

  xy_range(1:2) = xy_max(1:2) - xy_min(1:2);

  margin = 0.025 * max ( xy_range(1), xy_range(2) );

  x_min = xy_min(1) - margin;
  x_max = xy_max(1) + margin;
  y_min = xy_min(2) - margin;
  y_max = xy_max(2) + margin;
%
%  Hold the graphics screen.
%
  hold ( 'on' );
%
%  Draw the element faces.
%
  if ( 2 <= element_show )

    face_color = 'g';

    for element = 1 : element_num
      p = element_node(element,1:element_order);
      patch ( node_xy(p,1), node_xy(p,2), face_color );
    end

  end
%
%  Draw the element lines.
%
  if ( 1 <= element_show )

    line_color = 'r';

    for element = 1 : element_num
      p = [ element_node(element,1:element_order), element_node(element,1) ];
      line ( 'XData', node_xy(p,1), 'YData', node_xy(p,2), 'Color', line_color );
    end

  end
%
%  Label the elements.
%
  if ( 3 <= element_show )

    for element = 1 : element_num

      ave_x = 0.0;
      ave_y = 0.0;

      for i = 1 : element_order
        node = element_node(element,i);
        ave_x = ave_x + node_xy(node,1);
        ave_y = ave_y + node_xy(node,2);
      end

      ave_x = ave_x / element_order;
      ave_y = ave_y / element_order;

      text ( ave_x, ave_y, num2str ( element ) );

    end

  end
%
%  Draw the boundary in red.
%
  boundary_segments = mesh_boundary ( element_node );
  boundary_num = size ( boundary_segments, 1 );
  for i = 1 : boundary_num 
    n1 = boundary_segments(i,1);
    n2 = boundary_segments(i,2);
    bx = [ node_xy(n1,1), node_xy(n2,1) ];
    by = [ node_xy(n1,2), node_xy(n2,2) ];
    plot ( bx, by, 'r-', 'linewidth', 3 );
  end
%
%  Draw the nodes.
%
  if ( 1 <= node_show )
    plot ( node_xy(:,1), node_xy(:,2), 'o', 'markersize', 4, ...
      'markerfacecolor', 'black' );
  end
%
%  Label the nodes.
%
  if ( 2 <= node_show )
    abit = margin;
    for j = 1 : node_num
      text ( node_xy(j,1) + abit, node_xy(j,2) + abit, num2str(j) );
    end
  end

  xlabel ( '<-- X -->', 'Fontsize', 16 );
  ylabel ( '<-- Y -->', 'Fontsize', 16 );
%
%  The TITLE function will interpret underscores in the title.
%  We need to unescape such escape sequences!
%
  title_string = s_escape_tex ( prefix );
  title ( title_string, 'Fontsize', 16 );

% axis ( [ x_min, x_max, y_min, y_max ] );
  axis ( 'equal' );
%
%  Release the graphics screen.
%
  hold ( 'off' );

  pause ( 5 );
  filename = sprintf ( '%s.png', prefix );
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );

  return
end
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.
%
  boundary_loop = 0;

  for b1 = 1 : boundary_segment_num
    another = false;
    for b2 = b1 + 1 : boundary_segment_num
      if ( boundary_segments(b2,1) == boundary_segments(b1,2) )
        another = true;
        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
    if ( ~ another )
      boundary_loop = boundary_loop + 1;
    end
  end

  return
end
function s2 = s_escape_tex ( s1 )

%*****************************************************************************80
%
%% s_escape_tex() de-escapes TeX escape sequences.
%
%  Discussion:
%
%    In particular, every occurrence of the characters '\', '_',
%    '^', '{' and '}' will be replaced by '\\', '\_', '\^',
%    '\{' and '\}'.  A TeX interpreter, on seeing these character
%    strings, is then likely to return the original characters.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    19 January 2007
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    string S1, the string to be de-escaped.
%
%  Output:
%
%    string S2, a copy of the string, modified to avoid TeX escapes.
%
  s1_length = length ( s1 );

  s1_pos = 0;
  s2_pos = 0;
  s2 = [];

  while ( s1_pos < s1_length )

    s1_pos = s1_pos + 1;

    if ( s1(s1_pos) == '\' | ...
         s1(s1_pos) == '_' | ...
         s1(s1_pos) == '^' | ...
         s1(s1_pos) == '{' | ...
         s1(s1_pos) == '}' )
      s2_pos = s2_pos + 1;
      s2 = strcat ( s2, '\' );
    end

    s2_pos = s2_pos + 1;
    s2 = strcat ( s2, s1(s1_pos) );

  end

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% timestamp() prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end
 
