source: inundation/parallel/pmesh_divide.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 math import floor
from Numeric import zeros, Float, Int

#########################################################
#
# 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 pmesh_divide_linda(f, Domain, n_x = 1, n_y = 1):

    # read in the pmesh

    domain = pmesh_to_domain_instance(f, Domain)

    # 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]

    rect = domain.xy_extent

    # 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 arrays

    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 = {}
    for i in range(len(domain.triangles)):
        t = domain.triangles[i]

        x_coord = domain.coordinates[t[0]][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.coordinates[t[0]][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[0], t[1], t[2]])
        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]

    # extract the node list

    nodes = []
    for n in domain.coordinates:
        nodes.append([n[0], n[1]])

    return nodes, triangles, boundary,  triangles_per_proc, rect

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

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()

    return nodes, triangles, boundary, triangles_per_proc, quantities

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
    from Numeric import argsort
    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])

    #print tri_list

    #print tri_index

    # 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()

    return nodes, triangles, boundary, triangles_per_proc, quantities