source: inundation/parallel/mg2ga.py @ 1855

#########################################################
#   
#  
#  Read in a data file stored in the mg_cell data format.
#
#  mg_cell is a code written by Linda Stals for parallel
# mulitgrid solver based on adaptive finite elements. The
# reason this code is used is because it gives a grid
# partition.
#
#  This is only intended as a temporary file, once an
# automatic grid partitioner has been incorporated this
# file will become redundant.
#
#  Authors: Linda Stals and Matthew Hardy, June 2005
#
#
#########################################################

from string import atoi
from string import atof
from string import split


#########################################################
#  Read in a mg_cell file. It is assumed that the file
# has a very specific format. An example file is given in
# test_3l_2c.out.
#
#  *) The mg_cell data-structure is a node-based
# data-structure, consequently the same triangle will
# appear three times (mg_cell outputs all of the triangles
# connected to a given node).
#
#  *) The node labels are given in terms of the global ID,
# which may be thought of as a string (e.g. 12_2), this
# needs to be converted into an integer format for the GA
# data structure.
#
#  *) It is necessary to specify the boundary conditions,
# which has been done in mg_cel by outputing the edges
# with a tag. If the tag is negative the edge is an
# interior edge, if the tag is positive it is a boundary
# edge. Different boundary conditions will be given
# different tags.
#
# -------------------------------------------------------
#
#  *) The information returned by this routine is the
# nodes, triangles, boundary edges and the number of
# triangles to be assigned to each processor. All of the
# node labels are still stored by the global ID.
#
#########################################################
def readmg(f):

    # read the nuber of processors

    f.readline()
    no_proc = atoi(f.readline().split()[0])

    # skip over additional information not needed
    
    for i in range(5):
        f.readline();

    # nodes list

    nodes = []

    # boundary edge list

    edges = []

    # triangle list

    triangles = []

    # know first nodes_per_proc[0] belong to cell 0, next
    # nodes_per_proc[1] belong to cell 1 etc.

    nodes_per_proc = [] 
    triangles_per_proc = []

    # loop over the processors

    for q in range(no_proc):
    
        # read the nodes

        no_nodes = atoi(f.readline().split()[1])
        nodes_per_proc.append(no_nodes)
        f.readline()

        for i in range(no_nodes):
            line = f.readline().split()
            line[1:3] = map(atof, line[1:3])
            nodes.append(line[0:3])

        # read the edges
        
        f.readline()
        no_edges = atoi(f.readline().split()[1])
        f.readline()

        for i in range(no_edges):
            line = f.readline().split()
            if (atoi(line[2]) >= 0):
                edges.append(line)
        
        # read the triangles

        f.readline(), f.readline(), f.readline()
        no_triangles = atoi(f.readline().split()[1])
        f.readline()

        # remove any triangles that have been repeated
        
        for i in range(no_triangles):
            line = f.readline().split()
            line.sort()
            if (len(triangles)==0):
                triangles.append(line)
                continue
            
            for t in range(len(triangles)):
                if (triangles[t]==line): 
                    break
            if (triangles[t]!=line): 
               triangles.append(line)

        # keep a record of how many triangles belong to each processor
        
        no_triangles=len(triangles)-sum(triangles_per_proc)
        triangles_per_proc.append(no_triangles)

        # skip extra space
    
        f.readline()
        f.readline()
        f.readline()
        f.readline()

    # return the nodes, triangles, boundary edges and the number of triangles
    # to be assigned to each processor
    
    return nodes, triangles, edges, triangles_per_proc


#########################################################
# Convert the format of the data to that used by
# pyvolution
#
#
# *) Change the nodes global ID's to an integer value,
#starting from 0.
#
# *) The triangles and boundary edges must also be
# updated accordingly.
#
# -------------------------------------------------------
#
# *) The nodes, triangles and boundary edges defined by
# the new numbering scheme are returned
#
#########################################################

def restructure(nodes, triangles, edges):
    Nnodes =len(nodes)
    Ntriangles = len(triangles)
    
    index={} #a dictionary mapping existing node ids to the new ids 0,1,2,...

    GAnodes = []
    
    for node_idnew in range(Nnodes): #move through the list of nodes, renumbering to 0,1,2, ...
        index[nodes[node_idnew][0]]=node_idnew #make the dictionary entry, for later use by triangles
        GAnodes.append(nodes[node_idnew][1:3])
#        nodes[node_idnew][0]=node_idnew #renumber the node

    GAedges = {}
    for e in edges:
        GAedges[index[e[0]]]=[]
    for e in edges:
        GAedges[index[e[0]]].append([index[e[1]],e[2]])
    
    #Now loop over the triangles, changing the node ids and orientation. 
    for t in range(Ntriangles):
        for i in range(3):
            n=index[triangles[t][i]] #use the dictionary to get the new node id
            triangles[t][i]=n #relabel the node

            
    del(index) #delete the dictionary

    return GAnodes, triangles, GAedges


#########################################################
#
# Reorientate the triangles
#
#
# *) The ordering of the nodes in the triangles may need
# be changed to ensure the correct orientation.
#
# -------------------------------------------------------
#
# *) The modified list of triangles is returned
#
#########################################################
def reorient(nodes, triangles):
    

    Nnodes =len(nodes)
    Ntriangles = len(triangles)
    
    #Now loop over the triangles, changing the node ids and orientation. 
    for t in range(Ntriangles):
        x=[] #x-coordinates of the three nodes (vertices)
        y=[] #y-coordinates of the three nodes (vertices)
        for i in range(3):
            n=triangles[t][i] #use the dictionary to get the new node id
            x.append(nodes[n][0])
            y.append(nodes[n][1])
        if ((x[1]-x[0])*(y[2]-y[0])-(x[2]-x[0])*(y[1]-y[0])<0): #need to change the orientation
            triangles[t][2]=triangles[t][0]#swap the 0th and 2nd nodes
            triangles[t][0]=n
    
    del(x)
    del(y)
    
    return triangles


#########################################################
#
# Put the boundary edge information into a form compatible
# with the GA datastructure
#
# *) The GA datastructure assumes that the boundary edge
# information is stored in a dictionary. The dictionary
# key is a list [ti, ei], where ti is the triangle number
# and ei (0 <= ei < 2) is the edge number for the triangle
# ti. The value assigned to the given key is the boundary
# tag.
#
# -------------------------------------------------------
#
# *) The boundary edge information is returned.
#########################################################
def mesh_boundary(triangles, edges):

    # intialise the boundary dictionary
    
    boundary = {}

    # loop over all of the triangles
    
    for i in range(len(triangles)):

        # given a particular triangle
        
        t = triangles[i]

        # is the first vertex a boundary node?
        
        if (edges.has_key(t[0])):

            # check to see if the triangles has a boundary
            # edge 
            
            for b in edges[t[0]]:
                if (b[0] == t[1]):
                   boundary[i, 2] = b[1]
                elif (b[0] == t[2]):
                    boundary[i, 1] = b[1]

        # is the second vertex a boundary node?
        
        if (edges.has_key(t[1])):

            # check to see if the triangles has a boundary
            # edge
            
            for b in edges[t[1]]:
                if (b[0] == t[0]):
                   boundary[i, 2] = b[1]
                elif (b[0] == t[2]):
                    boundary[i, 0] = b[1]

        # is the third vertex a boundary node?
        
        if (edges.has_key(t[2])):

            # check to see if the triangles has a boundary
            # edge
            
            for b in edges[t[2]]:
                if (b[0] == t[0]):
                   boundary[i, 1] = b[1]
                elif (b[0] == t[1]):
                    boundary[i, 0] = b[1]

    # return the boundary information
    
    return boundary


#########################################################
#
# Read the nodes, triangles and boundary information from
# given file
#
# *) The information in the file is stored using the
# mg_cell format.
#
# *) See the documentation of the previous functions.
#
# -------------------------------------------------------
#
# *The information returned by this routine includes the
# nodes, triangles, boundary edges and the number of
# triangles to be assigned to each processor.
#
#########################################################

def mg2ga(f):

    # read in the information from the file
    
    [nodes, triangles, edges, triangles_per_proc] = readmg(f)

    # change the node numbering scheme to be compatible
    # with the GA datastructure
    
    [nodes, triangles, edges] = restructure(nodes, triangles, edges)

    # makesure the triangles are orientated in the correct
    # way
    triangles = reorient(nodes, triangles)

    # collect the boundary edge information
    
    boundary = mesh_boundary(triangles, edges)

    # return the information
    
    return nodes, triangles, boundary, triangles_per_proc