source: inundation-numpy-branch/parallel/pmesh_divide.py @ 3031

#########################################################
#
#
#  Read in a data file and subdivide the triangle list
#
#
#  The final routine, pmesh_divide_metis, does automatic
# grid partitioning. Once testing has finished on this
# routine the others should be removed.
#
#  Authors: Linda Stals and Matthew Hardy, June 2005
#  Modified: Linda Stals, Nov 2005
#            Jack Kelly, Nov 2005
#
#
#########################################################

from os import sep
from sys import path
from math import floor

from Numeric import zeros, Float, Int, reshape, argsort


#########################################################
#
# If the triangles list is reordered, the quantities
# assigned to the triangles must also be reorded.
#
# *) quantities contain the quantites in the old ordering
# *) proc_sum[i] contains the number of triangles in
# processor i
# *) tri_index is a map from the old triangle ordering to
# the new ordering, where the new number for triangle
# i is proc_sum[tri_index[i][0]]+tri_index[i][1]
#
# -------------------------------------------------------
#
# *) The quantaties are returned in the new ordering
#
#########################################################

def reorder(quantities, tri_index, proc_sum):

    # Find the number triangles

    N = len(tri_index)

    # Temporary storage area

    index = zeros(N, Int)
    q_reord = {}

    # Find the new ordering of the triangles

    for i in range(N):
        bin = tri_index[i][0]
        bin_off_set = tri_index[i][1]
        index[i] = proc_sum[bin]+bin_off_set

    # Reorder each quantity according to the new ordering

    for k in quantities:
        q_reord[k] = zeros((N, 3), Float)
        for i in range(N):
            q_reord[k][index[i]]=quantities[k].vertex_values[i]
    del index

    return q_reord


#########################################################
#
# Do a simple coordinate based division of the grid
#
# *) n_x is the number of subdivisions in the x direction
# *) n_y is the number of subdivisions in the y direction
#
# -------------------------------------------------------
#
# *)  The nodes, triangles, boundary, and quantities are
# returned. triangles_per_proc defines the subdivision.
# The first triangles_per_proc[0] triangles are assigned
# to processor 0, the next triangles_per_proc[1] are
# assigned to processor 1 etc. The boundary and quantites
# are ordered the same way as the triangles
#
#########################################################

def pmesh_divide(domain, n_x = 1, n_y = 1):

    # Find the bounding box
    
    x_coord_min = domain.xy_extent[0]
    x_coord_max = domain.xy_extent[2]
    y_coord_min = domain.xy_extent[1]
    y_coord_max = domain.xy_extent[3]

    # Find the size of each sub-box

    x_div = (x_coord_max-x_coord_min)/n_x
    y_div = (y_coord_max-y_coord_min)/n_y

    # Initialise the lists
    
    tri_list = []
    triangles_per_proc = []
    proc_sum = []
    for i in range(n_x*n_y):
        tri_list.append([])
        triangles_per_proc.append([])
        proc_sum.append([])
        tri_list[i] = []

    # Subdivide the triangles depending on which sub-box they sit
    # in (a triangle sits in the sub-box if its first vectex sits
    # in that sub-box)

    tri_index = {}
    N = domain.number_of_elements
    for i in range(N):
        t = domain.triangles[i]

        x_coord = domain.centroid_coordinates[i][0]
        bin_x = int(floor((x_coord-x_coord_min)/x_div))
        if (bin_x == n_x):
            bin_x = n_x-1

        y_coord = domain.centroid_coordinates[i][1]
        bin_y = int(floor((y_coord-y_coord_min)/y_div))
        if (bin_y == n_y):
            bin_y = n_y-1

        bin = bin_y*n_x + bin_x
        tri_list[bin].append(t)
        tri_index[i] = ([bin, len(tri_list[bin])-1])

    # Find the number of triangles per processor and order the
    # triangle list so that all of the triangles belonging to
    # processor i are listed before those belonging to processor
    # i+1

    triangles = []
    for i in range(n_x*n_y):
        triangles_per_proc[i] = len(tri_list[i])
        for t in tri_list[i]:
            triangles.append(t)

    # The boundary labels have to changed in accoradance with the
    # new triangle ordering, proc_sum and tri_index help with this

    proc_sum[0] = 0
    for i in range(n_x*n_y-1):
        proc_sum[i+1]=proc_sum[i]+triangles_per_proc[i]

    # Relabel the boundary elements to fit in with the new triangle
    # ordering

    boundary = {}
    for b in domain.boundary:
        t =  tri_index[b[0]]
        boundary[proc_sum[t[0]]+t[1], b[1]] = domain.boundary[b]

    quantities = reorder(domain.quantities, tri_index, proc_sum)

    # Extract the node list
    
    nodes = domain.coordinates.copy()
    ttriangles = zeros((len(triangles), 3), Int)
    for i in range(len(triangles)):
        ttriangles[i] = triangles[i]
        
    return nodes, ttriangles, boundary, triangles_per_proc, quantities


#########################################################
#
# Do a coordinate based division of the grid using the
# centroid coordinate
#
# *) n_x is the number of subdivisions in the x direction
# *) n_y is the number of subdivisions in the y direction
#
# -------------------------------------------------------
#
# *)  The nodes, triangles, boundary, and quantities are
# returned. triangles_per_proc defines the subdivision.
# The first triangles_per_proc[0] triangles are assigned
# to processor 0, the next triangles_per_proc[1] are
# assigned to processor 1 etc. The boundary and quantites
# are ordered the same way as the triangles
#
#########################################################
def pmesh_divide_steve(domain, n_x = 1, n_y = 1):
    
    # Find the bounding box
    x_coord_min = domain.xy_extent[0]
    x_coord_max = domain.xy_extent[2]
    y_coord_min = domain.xy_extent[1]
    y_coord_max = domain.xy_extent[3]


    # Find the size of each sub-box

    x_div = (x_coord_max-x_coord_min)/n_x
    y_div = (y_coord_max-y_coord_min)/n_y


    # Initialise the lists
    
    tri_list = []
    triangles_per_proc = []
    proc_sum = []
    for i in range(n_x*n_y):
        tri_list.append([])
        triangles_per_proc.append([])
        proc_sum.append([])
        tri_list[i] = []

    # Subdivide the triangles depending on which sub-box they sit
    # in (a triangle sits in the sub-box if its first vectex sits
    # in that sub-box)

    tri_index = {}
    N = domain.number_of_elements
 
    # Sort by x coordinate of centroid

    sort_order = argsort(argsort(domain.centroid_coordinates[:,0]))

    x_div = float(N)/n_x

    for i in range(N):
        t = domain.triangles[i]
        
        bin = int(floor(sort_order[i]/x_div))
        if (bin == n_x):
            bin = n_x-1

        tri_list[bin].append(t)
        tri_index[i] = ([bin, len(tri_list[bin])-1])

    # Find the number of triangles per processor and order the
    # triangle list so that all of the triangles belonging to
    # processor i are listed before those belonging to processor
    # i+1

    triangles = []
    for i in range(n_x*n_y):
        triangles_per_proc[i] = len(tri_list[i])
        for t in tri_list[i]:
           triangles.append(t)
        

    # The boundary labels have to changed in accoradance with the
    # new triangle ordering, proc_sum and tri_index help with this

    proc_sum[0] = 0
    for i in range(n_x*n_y-1):
        proc_sum[i+1]=proc_sum[i]+triangles_per_proc[i]

    # Relabel the boundary elements to fit in with the new triangle
    # ordering

    boundary = {}
    for b in domain.boundary:
        t =  tri_index[b[0]]
        boundary[proc_sum[t[0]]+t[1], b[1]] = domain.boundary[b]

    quantities = reorder(domain.quantities, tri_index, proc_sum)

    # Extract the node list
    
    nodes = domain.coordinates.copy()

    # Convert the triangle datastructure to be an array type,
    # this helps with the communication

    ttriangles = zeros((len(triangles), 3), Int)
    for i in range(len(triangles)):
        ttriangles[i] = triangles[i]
        
    return nodes, ttriangles, boundary, triangles_per_proc, quantities




###
#
# Divide the mesh using a call to metis, through pymetis.


path.append('..' + sep + 'pymetis')

from metis import partMeshNodal

def pmesh_divide_metis(domain, n_procs):
    
    # Initialise the lists
    # List, indexed by processor of # triangles.
    
    triangles_per_proc = []
    
    # List of lists, indexed by processor of vertex numbers
    
    tri_list = []
    
    # List indexed by processor of cumulative total of triangles allocated.
    
    proc_sum = []
    for i in range(n_procs):
        tri_list.append([])
        triangles_per_proc.append(0)
        proc_sum.append([])

    # Prepare variables for the metis call
    
    n_tri = len(domain.triangles)
    if n_procs != 1: #Because metis chokes on it...
        n_vert = len(domain.coordinates)
        t_list = domain.triangles.copy()
        t_list = reshape(t_list, (-1,))
    
        # The 1 here is for triangular mesh elements.
        
        edgecut, epart, npart = partMeshNodal(n_tri, n_vert, t_list, 1, n_procs)
        # print edgecut
        # print npart
        # print epart
        del edgecut
        del npart
        # Assign triangles to processes, according to what metis told us.
        
        # tri_index maps triangle number -> processor, new triangle number
        # (local to the processor)
        
        tri_index = {}
        triangles = []        
        for i in range(n_tri):
            triangles_per_proc[epart[i]] = triangles_per_proc[epart[i]] + 1
            tri_list[epart[i]].append(domain.triangles[i])
            tri_index[i] = ([epart[i], len(tri_list[epart[i]]) - 1])

        # print tri_list
        # print triangles_per_proc
        
        # Order the triangle list so that all of the triangles belonging
        # to processor i are listed before those belonging to processor
        # i+1

        for i in range(n_procs):
            for t in tri_list[i]:
                triangles.append(t)
            
        # The boundary labels have to changed in accoradance with the
        # new triangle ordering, proc_sum and tri_index help with this

        proc_sum[0] = 0
        for i in range(n_procs - 1):
            proc_sum[i+1]=proc_sum[i]+triangles_per_proc[i]

        # Relabel the boundary elements to fit in with the new triangle
        # ordering

        boundary = {}
        for b in domain.boundary:
            t =  tri_index[b[0]]
            boundary[proc_sum[t[0]]+t[1], b[1]] = domain.boundary[b]

        quantities = reorder(domain.quantities, tri_index, proc_sum)
    else:
        boundary = domain.boundary.copy()
        triangles_per_proc[0] = n_tri
        triangles = domain.triangles.copy()
        
        # This is essentially the same as a chunk of code from reorder.
        
        quantities = {}
        for k in domain.quantities:
            quantities[k] = zeros((n_tri, 3), Float)
            for i in range(n_tri):
                quantities[k][i] = domain.quantities[k].vertex_values[i]
        
    # Extract the node list
    
    nodes = domain.coordinates.copy()
    
    # Convert the triangle datastructure to be an array type,
    # this helps with the communication

    ttriangles = zeros((len(triangles), 3), Int)
    for i in range(len(triangles)):
        ttriangles[i] = triangles[i]
        
    return nodes, ttriangles, boundary, triangles_per_proc, quantities