#!python
#!python

# edit the list of surfaces and the start id
surface_list = [1,5,6]
start_sideset_id = 1

import cubit
from sets import Set

# finds a quad/face with the given 4 nodes 
def find_face(face_list, opp_face_nodes):
  if len(face_list) == 1:
    opposite_face = face_list[0]
  else:
    for f in face_list:
      face_conn = cubit.get_connectivity("Quad", f)
      if len(opp_face_nodes & Set(face_conn)) == 4:
        opposite_face = f
        break
  return opposite_face

# finds a chord of hexes extending from the given nodes
def get_hex_chord(nodes):
  hex_sheet1 = Set(cubit.get_hex_sheet(nodes[0], nodes[1]))
  hex_sheet2 = Set(cubit.get_hex_sheet(nodes[1], nodes[2]))
  hex_chord = list(hex_sheet1 & hex_sheet2)
  return hex_chord

# finds the nodes at the opposite end of a hex chord
def find_opposing_face_nodes(hex_chord, nodes):
  hex_map = {}
  #
  # place the connectivity in a map using the hex id as the key
  for hex in hex_chord:
    hex_map[hex] = (tuple(cubit.get_connectivity("Hex", hex)))
  #
  # this is a pretty intense n^2 search.  
  #There may be a more efficient way.  
  # Algorithm description:
  # 1) given the nodes associated with a face in a hex
  # 2) find the hex with those nodes
  # 3) intersect the list of nodes in the face with the 
  #    list of nodes in the face.  This list of nodes corresponds
  #    with the opposite face.
  # 4) continue this process until all hexes have been traversed.
  opp_face_nodes = Set(nodes)
  zero_count = 0
  while zero_count < len(hex_map):
    for hex in hex_map:
      if len(hex_map[hex]) > 0:
        test_face = opp_face_nodes ^ Set(hex_map[hex])
        if len(test_face) == 4:
          opp_face_nodes = test_face
          last_hex = hex
          hex_map[hex] = () # done with this hex; connectivity null
          zero_count = 0    # have to start over again
          break
      else:
        zero_count += 1
  #
  return last_hex, opp_face_nodes

# start main routine.  
# This routine finds the opposite face of a hex chord and
# creates a sideset with the faces of both ends of the chords.  
# There are guarantees if you run this on an irregular volume where
# hex chords may wrap back onto the same geometric surface.

for surface_id in surface_list:
  surface_quads = cubit.get_surface_quads(surface_id)
  #
  for quad in surface_quads:
    #
    # get a chord of hexes starting from the given quad
    start_quad_nodes = cubit.get_connectivity("Quad", quad)
    hex_chord = get_hex_chord(start_quad_nodes)
    # 
    # draw the hex chords
    command = "draw hex " + str(list(hex_chord))
    command = command.replace("[","") 
    command = command.replace("]","") 
    cubit.cmd(command)
    #
    # find the last hex and quad nodes the chord
    last_hex, opp_quad_nodes =  \
             find_opposing_face_nodes(hex_chord, start_quad_nodes)
    #
    # get all the dimension 2 (quad) elements off the hex
    quad_list = cubit.get_sub_elements("hex", last_hex, 2)
    #
    # find the quad associated with nodes on 
    # the end of the chord
    opposite_quad = find_face(quad_list, opp_quad_nodes)
    #
    #print "opposite_quad is: ", opposite_quad
    #
    command = "sideset " + str(start_sideset_id) + " face " \
                             + str(quad) + " " + str(opposite_quad)
    cubit.cmd(command)
    start_sideset_id += 1
#