#! /usr/bin/env python3
#
def 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.
#
  import numpy as np

  element_num, element_order = element_node.shape
#
#  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 = np.zeros ( [ segment_num, 2 ] )
  s = 0
  for e in range ( 0, element_num ):
    j = element_order - 1
    for jp1 in range ( 0, element_order ):
      segment[s,0] = element_node[e,j]
      segment[s,1] = element_node[e,jp1]
      s = s + 1
      j = jp1
#
#  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
  boundary_segments = np.zeros ( [ boundary_segment_num, 2 ] )
  for s in range ( 0, segment_num ):
    s0 = segment[s,0]
    s1 = segment[s,1]
    segment[s,0] = s1
    segment[s,1] = s0
    segment2 = np.unique ( segment, axis = 0 )
    if ( segment2.shape[0] == segment_num ):
      boundary_segment_num = boundary_segment_num + 1
      boundary_segments = np.append ( boundary_segments, [ [ s0, s1 ] ], axis = 0 )
    segment[s,0] = s0
    segment[s,1] = s1
#
#  Reorder the segments to form a loop.
#
  for b1 in range ( 0, boundary_segment_num - 1 ):
    for b2 in range ( b1 + 1, boundary_segment_num ):
      if ( boundary_segments[b2,0] == boundary_segments[b1,1] ):
        t0 = boundary_segments[b1+1,0]
        t1 = boundary_segments[b1+1,1]
        boundary_segments[b1+1,0] = boundary_segments[b2,0]
        boundary_segments[b1+1,1] = boundary_segments[b2,1]
        boundary_segments[b2,0] = t0
        boundary_segments[b2,1] = t1
        continue

  return boundary_segments

def mesh_boundary_test ( ):

#*****************************************************************************80
#
## mesh_boundary_test() tests mesh_boundary().
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    19 December 2022
#
#  Author:
#
#    John Burkardt
#
  import platform

  print ( '' )
  print ( 'mesh_boundary_test():' )
  print ( '  Python version: ' + platform.python_version ( ) )
  print ( '  Test mesh_boundary()' )

  mesh_boundary_test_square ( )
  mesh_boundary_test_beam ( )
  mesh_boundary_test_ell ( )
  mesh_boundary_test_serendipity ( )
#
#  Terminate.
#
  print ( '' )
  print ( 'mesh_boundary_test():' )
  print ( '  Normal end of execution.' )

  return

def mesh_boundary_test_square ( ):

#*****************************************************************************80
#
## mesh_boundary_test_square() tests a simple quadrilateral mesh.
#
#  Discussion:
#
#    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 ] ]
#
#    The boundary segments array should be computed as: 
#
#      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
#
  import numpy as np

  print ( '' )
  print ( 'mesh_boundary_test_square():' )
  print ( '  Boundary of a 4-quadrilateral mesh.' )
  element_node = np.array ( [ \
    [ 4, 5, 2, 1 ], \
    [ 5, 6, 3, 2 ], \
    [ 7, 8, 5, 4 ], \
    [ 8, 9, 6, 5 ] ] )

  print ( '' )
  print ( '  element_node array:' )
  print ( element_node )

  boundary_segments = mesh_boundary ( element_node )

  print ( '' )
  print ( '  boundary segments returned by mesh_boundary():' )
  print ( boundary_segments )

  return

def mesh_boundary_test_beam ( ):

#*****************************************************************************80
#
## mesh_boundary_test_beam() tests a beam quadrilateral mesh.
#
#  Discussion:
#
#    This example is available as mfem file beam-quad.mesh
#
#    Consider the following simple quadrilateral mesh:
#
#      9--10--11--12--13--14--15--16--17
#      |   |   |   |   |   |   |   |   |
#      0---1---2---3---4---5---6---7---8
#
#    for which the element array is:
#
#      element_node = [
#        0,  1, 10,  9 ...
#        1,  2, 11, 10 ...
#        2,  3, 12, 11 ...
#        3,  4, 13, 12 ...
#        4,  5, 14, 13 ...
#        5,  6, 15, 14 ...
#        6,  7, 16, 15 ...
#        7,  8, 17, 16 ]
#
#    The boundary segments array should be computed as: 
#
#      boundary_segments = [
#        [  0,  1 ]
#        [  1,  2 ],
#        [  2,  3 ],
#        [  3,  4 ],
#        [  4,  5 ],
#        [  5,  6 ],
#        [  6,  7 ],
#        [  7,  8 ],
#        [  8, 17 ]
#        [ 17, 16 ],
#        [ 16, 15 ],
#        [ 15, 14 ],
#        [ 14, 13 ],
#        [ 13, 12 ],
#        [ 12, 11 ],
#        [ 11, 10 ],
#        [ 10,  9 ],
#        [  9,  0 ] ]
#
#  Licensing:
# 
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    18 December 2022
#
#  Author:
#
#    John Burkardt
#
  import numpy as np

  print ( '' )
  print ( 'mesh_boundary_test_beam():' )
  print ( '  Boundary of an 8-quadrilateral mesh.' )

  element_node = np.array ( [ \
    [ 0,  1, 10,  9 ], \
    [ 1,  2, 11, 10 ], \
    [ 2,  3, 12, 11 ], \
    [ 3,  4, 13, 12 ], \
    [ 4,  5, 14, 13 ], \
    [ 5,  6, 15, 14 ], \
    [ 6,  7, 16, 15 ], \
    [ 7,  8, 17, 16 ] ] )

  print ( '' )
  print ( '  element_node array:' )
  print ( element_node )

  boundary_segments = mesh_boundary ( element_node )

  print ( '' )
  print ( '  boundary segments returned by mesh_boundary():' )
  print ( boundary_segments )

  return

def mesh_boundary_test_ell ( ):

#*****************************************************************************80
#
## mesh_boundary_test_ell() tests a triangular mesh.
#
#  Discussion:
#
#    Consider the following triangular mesh of the L-shaped region:
#
#      5--10--15
#      | \ | \ |
#      4---9--14
#      | \ | \ |
#      3---8--13--18--21--24
#      | \ | \ | \ | \ | \ |
#      2---7--12--17--20--23
#      | \ | \ | \ | \ | \ |
#      1---6--11--16--19--22
#
#    for which the element array is:
#
#      element_node = [
#        1,  6,  2
#        7,  2,  6
#        2,  7,  3
#        8,  3,  7
#        3,  8,  4
#        9,  4,  8
#        4,  9,  5
#       10,  5,  9
#        6, 11,  7
#       12,  7, 11
#        7, 12,  8
#       13,  8, 12
#        8, 13,  9
#       14,  9, 13
#        9, 14, 10
#       15, 10, 14
#       11, 16, 12
#       17, 12, 16
#       12, 17, 13
#       18, 13, 17
#       16, 19, 17
#       20, 17, 19
#       17, 20, 18
#       21, 18, 20
#       19, 22, 20
#       23, 20, 22
#       20, 23, 21
#       24, 21, 23 ]
#
#    The boundary segments array should be computed as: 
#
#      boundary_segments = [
#        [  1,  6 ]
#        [  6, 11 ],
#        [ 11, 16 ],
#        [ 16, 19 ],
#        [ 19, 22 ],
#        [ 22, 23 ],
#        [ 23, 24 ],
#        [ 24, 21 ],
#        [ 21, 18 ]
#        [ 18, 13 ],
#        [ 13, 14 ],
#        [ 14, 15 ],
#        [ 15, 10 ],
#        [ 10,  5 ],
#        [  5,  4 ],
#        [  4,  3 ],
#        [  3,  2 ],
#        [  2,  1 ],
#
#  Licensing:
# 
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    18 December 2022
#
#  Author:
#
#    John Burkardt
#
  import numpy as np

  print ( '' )
  print ( 'mesh_boundary_test_ell():' )
  print ( '  Boundary of an 24-triangle mesh.' )

  element_node = np.array ( [ \
     [ 1,  6,  2 ], \
    [  7,  2,  6 ], \
    [  2,  7,  3 ], \
    [  8,  3,  7 ], \
    [  3,  8,  4 ], \
    [  9,  4,  8 ], \
    [  4,  9,  5 ], \
    [ 10,  5,  9 ], \
    [  6, 11,  7 ], \
    [ 12,  7, 11 ], \
    [  7, 12,  8 ], \
    [ 13,  8, 12 ], \
    [  8, 13,  9 ], \
    [ 14,  9, 13 ], \
    [  9, 14, 10 ], \
    [ 15, 10, 14 ], \
    [ 11, 16, 12 ], \
    [ 17, 12, 16 ], \
    [ 12, 17, 13 ], \
    [ 18, 13, 17 ], \
    [ 16, 19, 17 ], \
    [ 20, 17, 19 ], \
    [ 17, 20, 18 ], \
    [ 21, 18, 20 ], \
    [ 19, 22, 20 ], \
    [ 23, 20, 22 ], \
    [ 20, 23, 21 ], \
    [ 24, 21, 23 ] ] )

  print ( '' )
  print ( '  element_node array:' )
  print ( element_node )

  boundary_segments = mesh_boundary ( element_node )

  print ( '' )
  print ( '  boundary segments returned by mesh_boundary():' )
  print ( boundary_segments )

  return

def mesh_boundary_test_serendipity ( ):

#*****************************************************************************80
#
## mesh_boundary_test_serendipity() tests a serendipity mesh.
#
#  Discussion:
#
#    Consider the following mesh:
#
#      1---2---3---4---5
#      |       |       |
#      6       7       8
#      |       |       |
#      9--10--11--12--13
#              |       |
#             14      15
#              |       |
#             16--17--18
#
#    Note that the elements are "quadratic", which in this case means that
#    each edge involves three nodes.
#
#    The element array is:
#
#      element_node = [
#        [  9, 10, 11,  7,  3,  2,  1,  6 ],
#        [ 11, 12, 13,  8,  5,  4,  3,  7 ],
#        [ 16, 17, 18, 15, 13, 12, 11, 14 ] ]
#
#    The boundary segments array should be computed as: 
#
#      boundary_segments = [
#        [  1,  6 ],
#        [  6,  9 ],
#        [  9, 10 ]
#        [ 10, 11 ],
#        [ 11, 14 ],
#        [ 14, 16 ],
#        [ 16, 17 ],
#        [ 17, 18 ],
#        [ 18, 15 ],
#        [ 15, 13 ],
#        [ 13,  8 ],
#        [  8,  5 ],
#        [  5,  4 ],
#        [  4,  3 ],
#        [  3,  2 ],
#        [  2,  1 ] ]
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    19 December 2022
#
#  Author:
#
#    John Burkardt
#
  import numpy as np

  print ( '' )
  print ( 'mesh_boundary_test_serendipity():' )
  print ( '  Boundary of a mesh involving 3 serendipity elements.' )

  element_node = np.array ( [ \
    [  9, 10, 11,  7,  3,  2,  1,  6 ], \
    [ 11, 12, 13,  8,  5,  4,  3,  7 ], \
    [ 16, 17, 18, 15, 13, 12, 11, 14 ] ] )

  print ( '' )
  print ( '  element_node array:' )
  print ( element_node )

  boundary_segments = mesh_boundary ( element_node )

  print ( '' )
  print ( '  boundary segments returned by mesh_boundary():' )
  print ( boundary_segments )

  return

def sortrows ( x ):

#*****************************************************************************80
#
## sortrows() lexically sorts the rows of an array.
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:
#
#    14 October 2022
#
#  Author:
#
#    John Burkardt
#
#  Input:
#
#    array x[]: the array to be sorted.
#
#  Output:
#
#    array x[]: the sorted array.
#
  import numpy as np

  x = x[ np.lexsort ( x.T[::-1] ) ]

  return x

def sortrows_test ( ):

#*****************************************************************************80
#
## sortrows_test() tests sortrows().
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:
#
#    14 October 2022
#
#  Author:
#
#    John Burkardt
#
  from numpy.random import default_rng
  import numpy as np

  rng = default_rng ( )

  print ( '' )
  print ( 'sortrows_test():' )
  print ( '  sortrows() lexically sorts the rows of an array.' )
  print ( '' )
  x = rng.integers ( low = 1, high = 10, size = [ 10, 3 ], endpoint = True )
  print ( '' )
  print ( '  Initial array:' )
  print ( x )

  x = sortrows ( x )

  print ( '' )
  print ( '  Array after sortrows() was applied:' )
  print ( x )

  return

def timestamp ( ):

#*****************************************************************************80
#
## timestamp() prints the date as a timestamp.
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:
#
#    21 August 2019
#
#  Author:
#
#    John Burkardt
#
  import time

  t = time.time ( )
  print ( time.ctime ( t ) )

  return

if ( __name__ == '__main__' ):
  timestamp ( )
  mesh_boundary_test ( )
  timestamp ( )